matlab中自带的CT投影与反投影函数(扇束、平行束及其之间的转化)

2020-1-27
作者：老李

matlab中自带的CT投影与反投影函数(扇束、平行束及其之间的转化)

我参考了这篇文章，以及mathwork文档，这些给我帮助很大
链接: linkhttps://blog.csdn.net/nanhuaibeian/article/details/95234813

matlab中的函数的特点

matlab里面自带的反投影函数中，可以去选择滤好几种波器，这是非常棒的

而且有着投影模式转换的函数（完成了数据重排的工作）

扇束投影

函数 fanbeam 可以生成扇形射线束投影：

%% function test fanbeam
%1-26
% B = fanbeam(g,D,param1,val1,parm2,val2,...)
% g：是包含被投影的物体的图像
% 
% D：是从扇形射线束的顶点到旋转中心的距离（单位为像素）
%    规定D大于g的直径的一半，D = K*sqrt(size(g,1)^2+size(g,2)^2)/2
%    其中，K 是大于1的常数(例如，K=1.5到2是合理的值)
% B 的每一列包含扇形射线束传感器每旋转一个角度得到的样本。
% B 中的列数由扇形旋转增量决定。默认情况下B 有 360 列，B中的行数由传感器的数目决定。
% 函数 fanbeam 通过计算对于任意旋转角度覆盖全部图像所需要的射线条数，来求出传感器的数量
% 这个数字强烈的依赖于指定的几何形状(直线或圆弧)
close all;clear;clc;

对于扇束投影方式，我们可以分别生成等弧度的探测器和等长度的探测器。

对于等长度的探测器，他的光子接收器的间距长度是由像素长度进行度量的

我们生成图像并用两种模式来进行投影，而投影的方式是由’FanSensorGeometry’这一个参数来决定的（‘line’,‘arc’），剩下的参数可以在上图中对应找到

%% fanbeam test
% 生成 Shepp-Logan 头部幻影
g2 = phantom('Modified Shepp-Logan',600);
% D 为从扇形射线束的顶点到旋转中心的距离
D = 1.5*hypot(size(g2,1),size(g2,2))/2;
% 生成矩形的扇形射线束投影
% 传感器间距为1个单位(默认间距)，角度增量为0.5度
B1_line = fanbeam(g2,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',1);
% 先转置数组，然后翻转数组
B1_line0 = flipud(B1_line');
B2_arc = fanbeam(g2,D,'FanSensorGeometry','arc','FanSensorSpacing',0.1,'FanRotationIncrement',1);
B2_arc0 = flipud(B2_arc');
figure(2);imshow(B1_line0,[]);colorbar
figure(3);imshow(B2_arc0,[]);colorbar

投影效果:
原图：

投影效果

扇束反投影

反投影时候的参数要把投影时候的参数传入函数之中，比较厉害的一点是，matlab自带的各种滤波器可以运用于其中

%% ifanbeam test
f1 = ifanbeam(B1_line,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',1);
f2 = ifanbeam(B2_arc,D,'FanSensorGeometry','arc','FanSensorSpacing',0.1,'FanRotationIncrement',1);
f3 = ifanbeam(B1_line,D,'FanSensorGeometry','line','FanSensorSpacing',1,'FanRotationIncrement',1,'Filter','Hamming');
f4 = ifanbeam(B2_arc,D,'FanSensorGeometry','arc','FanSensorSpacing',0.1,'FanRotationIncrement',1,'Filter','Hamming');
figure(4);imshow(f1,[]);colorbar
figure(5);imshow(f2,[]);colorbar
figure(6);imshow(f3,[]);colorbar
figure(7);imshow(f4,[]);colorbar

图四和图五是运用了R-L滤波器的效果，R-L滤波器对噪声比较敏感
图六和图七是运用了汉明滤波器的效果

平行束投影与反投影

注意角度这个参数是必须的

%% radon and iradon test
r1 = radon(g2,0:179);
figure;imshow(r1,[]);colorbar
r0 = iradon(r1,0:179);
figure;imshow(r0,[]);colorbar
r01 = iradon(r1,0:179,'Hamming');
figure;imshow(r01,[]);colorbar

以下是效果，第三个图还是用了汉明滤波器

扇束投影到平行束投影的转换

需要把投影时的参数全部传进去
'ParallelCoverage’这个参数有两种‘cycle’和‘halfcycle’指的是一整圈（360度）和半圈（180度）

%% fan2para test
P1_line = fan2para(B1_line,D,'FanRotationIncrement',1,...
    'FanSensorGeometry','line',...
    'FanSensorSpacing',1,...
    'ParallelCoverage','cycle');
P2_arc = fan2para(B1_line,D,'FanRotationIncrement',1,...
    'FanSensorGeometry','line',...
    'FanSensorSpacing',0.1,...
    'ParallelCoverage','halfcycle');
figure(8);imshow(P1_line,[]);colorbar
figure(9);imshow(P2_arc,[]);colorbar

与radon的效果基本上是一致的