ctsimpson.m

3D--ART重建

%--CTSimpson积分 --good
%--p投影值，K--方向数，S--每方向投影数

function [p,K,S]=CTSimpson(K,S)

%--da间隔角度，xra xrb  x轴起止点，yra yrb  y轴起止点，y轴间隔，h积分步长，p投影向量
da=pi/K;
xra=-sqrt(2)/2;xrb=sqrt(2)/2;
%xra=-0.5;xrb=0.5;
yra=-0.5;yrb=0.5;
d=(yrb-yra)/(S-1);
L=51;
h=(xrb-xra)/(2*(L-1));

p=[];

for k=0:K-1
    for s=0:S-1
        a=k*da;%确定射线方向
        yr=yra+s*d;%确定射线沿着yr轴的位置
        %v1=[cos(a) sin(a);-sin(a) cos(a)];   %旋转向量使积分线以da大小旋转
        %v=inv(v1);
        v=[cos(a) -sin(a);sin(a) cos(a)];
        f1=0;f2=0;
        for l=0:L-1
            xr=xra+(2*l+1)*h;
            xy=v*[xr;yr];
            f1=f1+Agauss(xy(1),xy(2));
        end
        for l=1:L-1
            xr=xra+2*l*h;
            xy=v*[xr;yr];
            f2=f2+Agauss(xy(1),xy(2));
        end
        xy=v*[xra;yr];
        fa=Agauss(xy(1),xy(2));
        xy=v*[xrb;yr];
        fb=Agauss(xy(1),xy(2));
        p(s+1,k+1)=h/3*(fa+4*f1+2*f2+fb);
    end
end