ctsimpson3d.m

用于computed tomography 的算法测试计算投影数据的simpson积分函数

%--CTSimpson3D积分 ,对三维火焰温度场进行投影的数值模拟 --2005/9/16
%--p1投影值（S×K×S矩阵），K--投影方向数，S--每方向投影数（射线数），此处沿Z轴所取层数也取值S

function [p1,K,S]=CTSimpson3D(K,S)

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

for zs=0:S-1
    zr=dz*zs;
    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);
            f1=0;f2=0;
            for l=0:L-1
                xr=xra+(2*l+1)*h;
                xy=v*[xr;yr];
                f1=f1+ThrDMod(xy(1),xy(2),zr);
            end
            for l=1:L-1
                xr=xra+2*l*h;
                xy=v*[xr;yr];
                f2=f2+ThrDMod(xy(1),xy(2),zr);
            end
            xy=v*[xra;yr];
            fa=ThrDMod(xy(1),xy(2),zr);
            xy=v*[xrb;yr];
            fb=ThrDMod(xy(1),xy(2),zr);
            p1(s+1,k+1,zs+1)=h/3*(fa+4*f1+2*f2+fb);
        end
    end
end