基于MATLAB的CT滤波反投影算法的实验研究

已于 2023-12-18 16:12:51 修改
阅读量1.4k 收藏 26
点赞数 1
文章标签： matlab 算法 机器学习
于 2023-04-30 14:37:07 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/fcxgfdjy/article/details/130448400
版权
该文详细展示了如何使用MATLAB中的iradon函数以及自定义FBP算法，结合S-L和R-L滤波器对Shepp-Loganphantom图像进行滤波反投影重建的过程。通过比较不同滤波器的效果，探讨了滤波在图像重建质量中的影响。
摘要由CSDN通过智能技术生成

实验四 滤波反投影算法的实验研究
1、利用”iradon”函数采用S-L与R-L滤波器进行滤波实现滤波反投影算法。

I=phantom(256);
theta=0:1:179;
P=radon(I,theta);
rec=iradon(P,theta,'linear','None');
rec_RL=iradon(P,theta,'Ram-Lak');
rec_SL=iradon(P,theta,'linear','Shepp-Logan');
figure;
subplot(2,2,1);imshow(I,[]),title('原始图像');
subplot(2,2,2);imshow(rec,[]),title('直接反投影图像');
subplot(2,2,3);imshow(rec_RL,[]),title('RL滤波反投影图像');
subplot(2,2,4);imshow(rec_SL,[]),title('SL滤波反投影图像');

在这里插入图片描述

2、利用FBP算法重建图像
1）对椭圆平行束投影数据进行滤波反投影重建，分别采用S-L与R-L滤波器进行滤波

clear;
N=256;
I = phantom('Shepp-Logan',N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
d=1;
%%产生投影数据
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
%%
fh_RL=medfuncRlfilterfunction(N,d);
fh_SL=medfuncSlfilterfunction(N,d);
%%
rec=zeros(N);
for m=1:theta_num
    pm=P1(:,m);
    Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
            n=floor(Xrm);
            t=Xrm-floor(Xrm);
            n=max(1,n);
            n=min(n,N-1);
            p=(1-t)*pm(n)+t*pm(n+1);
            rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
        end
    end
end
%%
rec_RL=zeros(N)
for m1=1:theta_num
    pm1=P1(:,m1);
    pm_RL=conv(fh_RL,pm1,'same');
    Cm1=(N/2)*(1-cos((m1-1)*delta)-sin((m1-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm1=Cm1+(k2-1)*cos((m1-1)*delta)+(k1-1)*sin((m1-1)*delta);
            n1=floor(Xrm1);
            t1=Xrm1-floor(Xrm1);
            n1=max(1,n1);
            n1=min(n1,N-1);
            p1=(1-t1)*pm_RL(n1)+t1*pm_RL(n1+1);
            rec_RL(N+1-k1,k2)=rec_RL(N+1-k1,k2)+p1;
        end
    end
end
%%
rec_SL=zeros(N)
for m2=1:theta_num
    pm2=P1(:,m2);
    pm_SL=conv(fh_SL,pm2,'same');
    Cm2=(N/2)*(1-cos((m2-1)*delta)-sin((m2-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm2=Cm2+(k2-1)*cos((m2-1)*delta)+(k1-1)*sin((m2-1)*delta);
            n2=floor(Xrm2);
            t2=Xrm2-floor(Xrm2);
            n2=max(1,n2);
            n2=min(n2,N-1);
            p2=(1-t2)*pm_SL(n2)+t2*pm_SL(n2+1);
            rec_SL(N+1-k1,k2)=rec_SL(N+1-k1,k2)+p2;
        end
    end
end
%%
figure;
subplot(2,2,1);imshow(I,[]),title('原始图像');
subplot(2,2,2);imshow(rec,[]),title('直接反投影图像');
subplot(2,2,3);imshow(rec_RL,[]),title('RL滤波反投影图像');
subplot(2,2,4);imshow(rec_SL,[]),title('SL滤波反投影图像');

在这里插入图片描述
2）对Shepp-Logan图平行束投影数据进行滤波反投影重建，别采用S-L与R-L滤波器进行滤波

clear;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
d=1;
%%产生投影数据
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
%%
fh_RL=medfuncRlfilterfunction(N,d);
fh_SL=medfuncSlfilterfunction(N,d);
%%
rec=zeros(N);
for m=1:theta_num
    pm=P1(:,m);
    Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
            n=floor(Xrm);
            t=Xrm-floor(Xrm);
            n=max(1,n);
            n=min(n,N-1);
            p=(1-t)*pm(n)+t*pm(n+1);
            rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
        end
    end
end
%%
rec_RL=zeros(N)
for m1=1:theta_num
    pm1=P1(:,m1);
    pm_RL=conv(fh_RL,pm1,'same');
    Cm1=(N/2)*(1-cos((m1-1)*delta)-sin((m1-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm1=Cm1+(k2-1)*cos((m1-1)*delta)+(k1-1)*sin((m1-1)*delta);
            n1=floor(Xrm1);
            t1=Xrm1-floor(Xrm1);
            n1=max(1,n1);
            n1=min(n1,N-1);
            p1=(1-t1)*pm_RL(n1)+t1*pm_RL(n1+1);
            rec_RL(N+1-k1,k2)=rec_RL(N+1-k1,k2)+p1;
        end
    end
end
%%
rec_SL=zeros(N)
for m2=1:theta_num
    pm2=P1(:,m2);
    pm_SL=conv(fh_SL,pm2,'same');
    Cm2=(N/2)*(1-cos((m2-1)*delta)-sin((m2-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm2=Cm2+(k2-1)*cos((m2-1)*delta)+(k1-1)*sin((m2-1)*delta);
            n2=floor(Xrm2);
            t2=Xrm2-floor(Xrm2);
            n2=max(1,n2);
            n2=min(n2,N-1);
            p2=(1-t2)*pm_SL(n2)+t2*pm_SL(n2+1);
            rec_SL(N+1-k1,k2)=rec_SL(N+1-k1,k2)+p2;
        end
    end
end
%%
figure;
subplot(2,2,1);imshow(I,[]),title('原始图像');
subplot(2,2,2);imshow(rec,[]),title('直接反投影图像');
subplot(2,2,3);imshow(rec_RL,[]),title('RL滤波反投影图像');
subplot(2,2,4);imshow(rec_SL,[]),title('SL滤波反投影图像');

在这里插入图片描述

function P = medfuncParallelBeamForwardProjection(theta,N,N_d)
    shep=[   1   .69     .92     0       0       0
        -.8 .6624   08740   0       -.0184  0
        -.2 .1100   .3100   .22     0       -10
        -.2 .1600   .4100   -.22    0       15
        .1  .2100   .2500   0       .35     5
        .1  .0460   .0560   .06       .1      15
        .1  .0260   .0460   .13       -.1     20
        .1  .0460   .0130   -.08    -.605   -5
        .1  .0430   .0230   .03       -.606   3
        .1  .0230   .0460   .06     -.605   0];
    theta_num = length(theta);
    P = zeros(N_d,theta_num); % 存放投影数据
    rho = shep(:,1).';  % 圆对应密度
    be =0.5*N*shep(:,2).'; %椭圆短半轴
    ae = 0.5*N*shep(:,3).'; % 椭圆长半轴
    xe =0.5*N*shep(:,4).'; %圆中心x坐标
    ye = 0.5*N*shep(:,5).'; %圆中心y坐标
    alpha = shep(:,6).'; % 圆旋转角度
    alpha = alpha * pi/180;
    theta = theta * pi/180; % 角度换算成弧度
    TT =-(N_d-1)/2:(N_d-1)/2;% 探测器坐标
    for k1 = 1 : theta_num
        P_theta = zeros(1,N_d);% 存放每一角度投影数据
        for k2 = 1: max(size(xe))
            a =(ae(k2) * cos(theta(k1)-alpha(k2)))^2+(be(k2) * sin(theta(k1)-alpha(k2)))^2;
            temp = a- (TT-xe(k2) * cos(theta(k1))-ye(k2) * sin(theta(k1))).^2;
            ind = temp>0;% 根号内值需为非负
            P_theta(ind)= P_theta(ind) +rho(k2) * (2 * ae(k2) * be(k2) * sqrt(temp(ind)))./a;
        end
        P(:,k1) = P_theta.';
    end
end

function fh_RL=medfuncRlfilterfunction(N,d)
fh_RL=zeros(1,N);
for k1=1:N
    fh_RL(k1)=-1/(pi*pi*((k1-N/2-1)*d)^2);
    if mod(k1-N/2-1,2)==0
        fh_RL(k1)=0;
    end
end
fh_RL(N/2+1)=1/(4*d^2);

function fh_SL=medfuncSlfilterfunction(N,d)
fh_PL=zeros(1,N);
for k1=1:N
    fh_SL(k1)=-2/(pi^2*d^2*(4*(k1-N/2-1)^2-1));
end
看星河的兔子
关注
基于滤波反投影的图像重建算法matlab仿真,R-L滤波和S-L滤波
FPGA/MATLAB学习教程/源码/项目合作开发
10-29 2437
用S-L滤波函数重建的图像中振荡相应较小，对含噪声的数据重建出来的图像质量也较R-L滤波函数重建的图像质量要好。但是，S-L滤波函数重建的图像在高频响应方面不如R-L滤波函数好，这是因为S-L滤波函数在高频段偏离了理想的滤波函数。R-L滤波函数的特点是函数形式简单，重建的图像轮廓清晰。与R-L滤波函数不同的是，S-L滤波函数在频域中不用矩形函数去截断滤波函数，而是用一些比较平滑的窗函数来约束滤波函数。R-L滤波函数的基本出发点是认为实际的二维图像函数总有一个频率上限，所以频域中的滤波函数可表示为。
滤波反投影matlab实现（matlab源程序）
04-13
滤波反投影
平行光滤波反投影matlab程序
09-16
平行光滤波反投影matlab程序，程序不错，可以看看。刚刚数学建模用的程序，感觉适合初学者..................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
Matlab之图像变换技术（十二）
qq_44790423的博客
03-13 2737
图像变换是将图像从空间域变换到变换域。图像变换的目的是根据图像在变换域的某些性质对其进行处理。通常，这些性质在空间城内很难获取。在变换城内处理结束后，将处理结果进行反变换到空间城。 这里将详细介紹图像变换技术，主要包括Radon变换和反变换，傅立叶变换和反变换，离散余弦变换和反变换等。 目录图像RADON变换RADON正变换Radon反变换图像傅里叶变换图像离散余弦变换图像Hadamard变换图像...
matlab iradon函数详解,如何处理iradon函数所得图像
weixin_29534399的博客
03-19 2105
matlab iradon函数小编是个既深情又薄情的人，深情在小编和谁在一起的时候就算出现更好的也没有见异思迁的念头。matlab iradon函数对一组多角度 x光信息作radon逆变换，得到的正方形图歌声形成的空间，任凭年华来去自由，所以依然保护着的人的容颜不曾改和场庞大而没有落幕的恨。M文件大致可以理解为由一系列的语句组成的相对独立的一个运行体。分为M脚本文件与M函数文件。M脚本文件没有参数...
MATLAB滤波反投影算法
10-31
利用MATLAB实现滤波反投影技术，编写的一个程序。。。。
滤波反投影fbp算法matlab
05-16
用matlab实现ct滤波反投影fbp算法，其中投影可以直接用radon变化实现，但代码中手动编写投影部分的代码，使用的滤波器为R-L滤波器
运用滤波反投影的方法对图像进行重建matlab仿真
fpga/matlab/simulink算法仿真工程
11-23 1293
这将意味着，每一幅反投影图像都将有正像素和负像素值，当所有的反投影图像相加形成最终的重建图像时，一些负像素值位置可能变成正像素，而平均值可能不为零，但是，典型地，最终的图像将还是有负像素值。这表明，将对应的投影 g ( ρ , θ ) g(\rho, \theta)g(ρ,θ) 与斜坡滤波器传递函数 s ( ρ ) s(\rho)s(ρ) 的傅里叶反变换进行卷积，可以得到角度 θ \thetaθ 的各个反投影，整个反投影图像可以通过对所有反投影图像积分得到。
【图像重建】基于滤波反投影实现图像重建附matlab代码
matlab_dingdang的博客
09-22 774
图像重建是CT技术的一个研究热点。重建算法的现实是对算法研究的一个重要环节。由Math Works公司推出的MATLAB工具软件具有强大的数学计算和图像处理功能,并为算法提供了一个方便有效的研究和实现的平台。本文在图像重建分析的基础上,运用MATLAB实现了基于扇束的滤波反投影重建算法的计算机模拟。智能优化算法。
滤波反投影算法的matlab仿真
FPGA/MATLAB学习教程/源码/项目合作开发
08-14 2130
本文将详细介绍滤波反投影算法的原理、实现过程，以及其变种——直接滤波反投影和解析法滤波反投影，并探讨其应用领域。：投影数据在经过Radon变换后，进行滤波和反投影操作，以恢复对象的像素值。解析法滤波反投影是滤波反投影的另一种变种，通过解析法得到投影数据的反投影，从而避免了离散化和数值计算带来的误差。直接滤波反投影是滤波反投影算法的一种变体，通过在投影域内直接滤波，而不是在图像域内滤波。将滤波后的投影数据代入解析反投影表达式，得到图像的解析反投影。对投影数据进行滤波，得到滤波后的投影数据。
滤波反投影matlab仿真程序
08-14
本程序用matlab实现滤波反投影，有助于我们理解滤波反投影算法的具体原理
基于matlab的平行束滤波反投影重建算法
07-22
基于matlab的平行束滤波反投影重建算法。基于matlab的平行束滤波反投影重建算法
基于扇形滤波的反投影重建算法附matlab代码.zip
06-18
1.版本：matlab2019a，不会运行可私信 2.领域：【图像重建】 3.内容：基于扇形滤波的反投影重建算法附matlab代码.zip 4.适合人群：本科，硕士等教研学习使用
投影算法：投影算法-matlab开发
05-30
投影算法是一种参数估计方法，用于估计传递函数的参数参考：Astrom 的自适应控制
CT 滤波反投影重建 直接反投影重建
11-29
投影算法调用matlab系统函数 [R,xp] = radon(I,theta); % Compute (forward ) projections 其余为自己编写，通过两种不同插值的方法，实现直接反投影和滤波反投影，最终script实现3种不同投影个数的重建效果。 适合CT重建算法初学者调试学习，直观了解不同算法，不同投影个数引起的不同结果。 phantom 为一个圆形。 是我在CMU的一次课程作业，含源码和文档。
基于matlab的CT图像滤波反投影算法仿真+仿真操作录像
10-15
1.版本：matlab2021a，我录制了仿真操作录像，可以跟着操作出仿真结果 2.领域：滤波反投影 3.内容：基于matlab的CT图像滤波反投影算法仿真，对比滤波反投影前后效果。 4.适合人群：本，硕等教研学习使用
整理好的基于matlab的滤波反投影重建算法源程序.zip
05-22
压缩包中的唯一文件“ksrjacjn.m”很可能是一个MATLAB脚本或函数，实现了滤波反投影算法的一部分或全部。MATLAB文件通常以.m为扩展名，其中.m脚本是可执行的代码，而.m函数则是可以被其他MATLAB代码调用的自定义函数...
FBP滤波反投影算法
12-02
FBP滤波反投影（Filtered Back-Projection，简称FBP）算法是计算机断层扫描（Computed Tomography，CT）图像重建领域的一种经典方法。该算法由Radon变换理论为基础，通过数学模型将从不同角度获取的投影数据转换为横...
滤波反投影图像重建算法
热门推荐
shuangyue的专栏
08-14 3万+
关于CT重建的算法有很多，在这里给大家介绍的是滤波反投影算法，其原理如下： 设f(x,y)表示需要重建的图像，用p(t,θ)表示在角度获取的f(x,y)的一个平行投影，t表示投影射线到对称中心（即旋转中心）的距离，设s为于角度θ下的投影X射线平行的坐标轴，它与t所在的坐标轴垂直，则：         （1） 对p(t,θ)进行一维傅里叶变换     （2） 将（1）代入（2）式得到
shepp-logan 模型
最新发布
11-30
Shepp-Logan 模型是一种在医学图像重建中常用的数学模型。它由两位科学家 Laurence Shepp 和 Benjamin F. Logan 在1974年提出。 Shepp-Logan 模型被广泛应用于计算机断层成像（CT）的重建过程中。CT 是一种通过扫描和测量人体内部组织的 X 射线吸收来生成图像的技术。而重建过程则是逆过程，通过从收集到的数据中恢复出原始的图像。 Shepp-Logan 模型假设人体内的组织结构可以用具有不同密度和形状的简单几何形状来表示。这些几何形状包括椭圆、圆柱等等，它们代表了不同组织的吸收特性和位置。 在重建过程中，使用 Shepp-Logan 模型将扫描数据转换成像素点的数值。这个模型能够估计每个像素点在扫描过程中受到的 X 射线的吸收量，从而可以计算出每个像素点的灰度值。通过这种方式，就能够还原出一张表示人体内部组织结构的图像。 Shepp-Logan 模型的优点是简单且易于实现。它不仅在医学领域中被广泛使用，也被应用于其他领域，如材料科学和工业检测等。然而，这个模型也有一些限制，例如它只能处理二维图像，并且在复杂的情况下可能存在一定的误差。 总的来说，Shepp-Logan 模型在医学图像重建中起到了重要的作用，具有一定的准确性和实用性。随着科学技术的发展，人们对更精确的图像重建方法的研究也在不断进行着。
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

看星河的兔子

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值