医学图像重建作业matlab源码

这是一个医学图像重建作业matlab源码，通过滤波反投影实现的，可运行，是医学图像计算课程作业

function  [a]=dps（P）
tic;
P=phantom（256）;
%P=imread（‘C:\0刘金豆\医学图像处理实验重建\lenna.bmp‘）;
%P=rgb2gray（P）;
%P= imresize（P[256 256]）;

[NN]=size（P）;%获取图像P的size，存入N
subplot（231）;imshow（P）;title（[int2str（N）‘*‘int2str（N）‘原始图像‘]）;

%先进行自定义radon变换------------------------------------------------------------
thm=45;                %45度时会出现最大尺寸
pre = imrotate（Pthm）;   %将图像P（图像的数据矩阵）绕图像的中心点旋转thm度，
                         %正数表示逆时针旋转， 负数表示顺时针旋转。返回旋转后的图像矩阵pre。
[mmaxnmax] = size（pre）;
s=1;    %定义步长为1

%创建一个180*nmax的空白图片，用以存储投影后的线状图片
Final = zeros（180/snmax）;%这里180代表180角度，每个角度投影成为一条线
t = 1;
for theta = 1:s:179
%对原图旋转一个角度，求和（线积分）
Protate = imrotate（Ptheta）; %对图像P旋转后的图Protate   
Pf = sum（Protate1）;%对Protate矩阵列求和
[mrealnreal]=size（Pf）;           %计算实际尺寸
%确定起始点
if （nmax - nreal）/2-floor（（nmax - nreal）/2） == 0         %y = floor（x） 函数将x中元素向下取整
    From = floor（（nmax - nreal）/2 + 1）;%总点数为偶数时
   else
    From = floor（（nmax - nreal）/2） + 1;%总点数为奇数时
    end
%确定结束点
End = floor（（nmax-nreal）/2） + nreal;
%将一个角度Radon变换后的线状图存入结果图像的某一行
Final（180/s-tFrom:End） = Pf; %从最底下一行开始存起
%上移一行
t = t + 1;
end
%再逆时针旋转
R=imrotate（Final90）;
subplot（232）;imshow（R[]）;title（‘自定义投影后图像‘）;

z=2*ceil（norm（size（P）-floor（（size（P）-1）/2）-1））+3;% radon变换默认平移点数/角度
e=floor（（z-N）/2）+2;   
R1=R（e:（N+e-1）:）;
[mmnn]=size（R1）;
subplot（233）;
imagesc（R1）;
title（[int2str（mm）‘*‘int2str（nn）‘正弦图‘]）;
colormap（gray）;
colorbar;

%滤波函数构造------------------------------------------------------------
g=1-N:N-1; 
%  d=1;
%  R-L滤波函数
%  for i=1:2*N-1
%         if g（i）==0
%         h（i）=1/（4*（d^2））;
%     else if mod（g（i）2）==0
%             h（i）=0;
%         else 
%             h（i）=（-1）/（pi^2*d^2*（g（i）^2））;
%         end
%     end
%  end

%S-L滤波函数
% d=1;
%  for i=1:2*N-1
%        h（i）=2/（pi^2*d^2*（1-4*g（i）^2））;
%  end

%Parzen滤波函数
 for i=1:2*N-1
     h（i）=（24*pi*g（i）*cos（pi*g（i））-96*sin（pi*g（i））-48*pi*g（i）*cos（pi*g（i）/2）+384*sin（pi*g（i）/2）-2*（pi^3）*（g（i）^3）-72*pi*g（i））/（4*（pi^5）*（g（i）^5））;
 end
 h（N）=0.0435;
 
%将投影与滤波函数卷积----------------------------------------------------
G=zeros（N180）;
for m=1:180      
    for i=1:N
        b=0;
        for n=1:N
            b=b+h（N+n-i）*R1（nm）;
            G（im）=b;
        end
    end
end

%投影滤波后线性内插进行图像重建----------------------------------------------
a=zeros（N）;       %重建图像初始化，每个像素点像素值为0
ns=（N+1）/2;
for m=1:180         %遍历每个投影角度
    r=pi/180;       %将角度换为弧度
    for i=1:N      
        for j=1:N    %遍历原图的每一个像素点
             Xrm=（i-ns）*cos（m*r）+（j-ns）*sin（m*r）+ns-1;  %坐标转换
            if Xrm<1
                n=1;                           %内插运算整数值
                t=abs（Xrm）-floor（abs（Xrm））;     %内插运算小数值
            else
                n=floor（Xrm）;
                t=Xrm-floor（Xrm）;
            end
            if n>N-1
                n=N-1;
            end
            c=（1-t）*G（nm）+t*G（n+1m）;       %内插后的滤波投影值
            a（N+1-ji）=a（N+1-ji）+c;         %像

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        4305  2018-12-05 11:18  医学图像重建作业\FBP_GUI.m
     文件        3786  2019-01-05 10:17  医学图像重建作业\dps.m
     目录           0  2019-01-17 14:14  医学图像重建作业\