MATLAB 图像处理 头颅模型ct图像重建

主要任务：

1.建立Shepp-Logan模型，使用灰度叠加

2.生成投影数据

3.反投影重建

4.滤波反投影重建

重建过程会产生星状伪迹，因此需要先修正再投影，不使用Radon函数，原理不讲，直接放代码，模型建立根据原函数源码改编，查看函数源码可在命令行窗口输入: >>open 函数

其余过程不是我自己写的，但是时间过久原网址找不到了

%Shepp-Logan模型
n=256;
q=zeros(n);
xax=((0:n-1)-(n-1)/2)/((n-1)/2);  
xg =repmat(xax,n,1);   % x坐标，y坐标旋转90度，形成块矩阵
%            ρ    a     b    x0    y0    α
%        ---------------------------------
ellipse= [    2   .92   .69    0     0     90;
            -.98 .8740 .6624   0  -.0184   90;
            -.02 .3100 .1100  .22    0     72;
            -.02 .4100 .1600 -.22    0     108;
             .01 .2500 .2100   0    .35    90;
             .01 .0460 .0460   0    .1     0;
             .01 .0460 .0460   0   -.1     0;
             .01 .0460 .0230 -.08  -.605   0;
             .01 .0230 .0230   0   -.606   0;
             .01 .0460 .0230  .06  -.605   90   ];
for k = 1:10
    asq = ellipse(k,2)^2;       % a^2
    bsq = ellipse(k,3)^2;       % b^2
    phi = ellipse(k,6)*pi/180;  % α，旋转角度
    x0 = ellipse(k,4);          % x 坐标
    y0 = ellipse(k,5);          % y 坐标
    A = ellipse(k,1);           % 折射指数
    x=xg-x0;                    % 使椭圆居中
    y=rot90(xg)-y0;             % 数组旋转90度
    cosp = cos(phi); 
    sinp = sin(phi);
    idx=find(((x.*cosp + y.*sinp).^2)./asq + ((y.*cosp - x.*sinp).^2)./bsq <= 1); 
    q(idx) = q(idx) + A;
end
[j,k]=size(q);
%归一化
for a=1:j
    for b=1:k
        if q(a,b)<=0.9874
            q(a,b)=0;
        elseif q(a,b)>=1.1644
            q(a,b)=255;
        else
            q(a,b)=(q(a,b)-0.9874)/(1.1644-0.9874)*255;
        end
    end
end
q=uint8(q);


%Sinogram投影数据图
P=q;
[m,n]=size(P);
thm=45;
pre=imrotate(P,thm);
[mmax,nmax]=size(pre);
s=1;
final=zeros(180/s,nmax);%创建图片存储投影后的图。180度，每个角度投影出一条线
t=1;
%对原图旋转一个角度，求线积分
for theta=1:s:179
    protate=imrotate(P,theta);
    pf=sum(protate,1);
    [mreal,nreal]=size(pf);%计算实际尺寸
    %确定起始点
    if ((nmax-nreal)/2-floor(nmax-nreal)/2)==0
        from=floor((nmax-nreal)/2)+1;%总点数为偶数
    else
        from=floor((nmax-nreal)/2)+1;%总点数为奇数
    end
    %确定结束点
    End=floor((nmax-nreal)/2)+nreal;
    %将一个角度的radon变换后的现状图存入结果的某一行
    final(180/s-t,from:End)=pf;%从最底下一行开始存
    t=t+1;
end
final=imrotate(final,90);


%直接反投影
I=pow2(nextpow2(size(final,1))-1);%重构图像大小
p1=zeros(I,I);
for i=1:size(final,2)
    tmp=imrotate(repmat(final(:,i),1,size(final,1)),i-1,'bilinear');
    tmp=tmp(floor(size(tmp,1)/2-I/2)+1:floor(size(tmp,1)/2+I/2),floor(size(tmp,2)/2-I/2)+1:floor(size(tmp,2)/2+I/2));
    p1=p1+tmp;
end
p1=p1/size(final,2);
p1=rot90(p1);


%滤波反投影
N=180;
H=size(final,1);
h=zeros((H*2-1),1);
for i=0:H-1 %RL滤波器，δ=1
    if i==0
        h(H-i)=1/4;
    elseif rem(i,2)==0
        h(H-i)=0;
        h(H+i)=0;
    else
        h(H-i)=-1/(i*pi)^2;
        h(H+i)=-1/(i*pi)^2;
    end
end
x=zeros(H,N);
for i=1:N
    s=final(:,i);
    xx=conv(s',h');
    x(:,i)=xx(H:2*H-1);
end
p2=zeros(I,I);
for i=1:I
    for j=1:I
        for k=1:180
            theta=k/180*pi;
            z=(j-I/2-0.5)*cos(theta)+(I/2+0.5-i)*sin(theta)+(H+1)/2;
            z1=floor(z);
            z2=floor(z+1);
            p2(i,j)=p2(i,j)+(z2-z)*x(z1,k)+(z-z1)*x(z2,k);
        end
    end
end
p2=pi/N*p2;

subplot(2,2,1);imshow(q);title('Shepp-Logan模型');
subplot(2,2,2);imagesc(final);colormap(gray(256));title('投影图');
subplot(2,2,3);imshow(p1,[]);title('直接反投影');
subplot(2,2,4);imshow(p2,[]);title('滤波反投影');