等弧度扇束的投影——matlab实现

时间：2020-1-31
作者：老李

等弧度扇束的投影的投影

工作的步凑分为以下四个次序

% 探测器形式为等弧度探测器

1.找出扇束放射源的位置

2.确定扇束的投射范围

3.表达射线

4.对射线加总起来

代码如下

%% generate fanbeam projection
%% 2020-1-31
% 探测器形式为等弧度探测器
% 1.找出扇束放射源的位置
% 2.确定扇束的投射范围
% 3.表达射线
% 4.对射线加总起来
close all; clear; clc;
%% GenerateModel 
views = 360;%放射源的位置
N = 600;%N为图片的半径
%制图
x1 = 0; y1 = 0; r1 = 600; mu1 = 0.5;
x2 = 200; y2 = 200; r2 = 150; mu2 = 1;
vec = -N:N;
[xx,yy] = meshgrid(vec,-vec);%生成坐标网格
img = zeros(2*N+1,2*N+1);%此时图像大小为（1201，1201）
img((xx - x1).^2 + (yy - y1).^2 <= r1^2) = mu1;
img((xx - x2).^2 + (yy - y2).^2 <= r2^2) = mu2;
figure('name','cirle in circle'); imshow(img, []); title('Original image');
%% Projection by exist function
D = ceil(sqrt(600.^2+600.^2));%半径经过计算为849，but function demand at least 850
D0 = 850;
pic1 = fanbeam(img,D0,'FanSensorGeometry','arc','FanSensorSpacing',1,'FanRotationIncrement',1);
figure;imshow((flipud(pic1')),[]);colorbar;title('fanbeam projection by exist function')
%% 找出放射源的地址
% 使用图像正上方的位置作为起点，逆时针旋转360度
% 每一次投放的角度为180度
V = zeros(views,2);
v0 = [0 D0];%正上方位置（起始位置）
for i =1:views
    theta = pi*i/180;
    V(i,1) = cos(theta).*v0(1) - sin(theta).*v0(2);
    V(i,2) = sin(theta).*v0(1) + cos(theta).*v0(2);
end
%% 确定放射线的表示
% 不需要扩展图片
% 扇束的投射范围是-90到89度
% 以放射源到图像中心的线的直线为第一条线
% 1.表达起始直线
% 2.表达投影
P = zeros(360,180);%投影图的大小
v = 0:179;
beta = 0:359;
kk = -2*D0:2*D0;
lineB_x = zeros(180,4*D0+1);
lineB_y = zeros(180,4*D0+1);
line0x = zeros(1,4*D0+1);
line0y = zeros(1,4*D0+1);
for j = 1:360
    for i = 1:180
        Vtheta = pi*v(i)/180;
        Btheta = pi*beta(j)/180;
        line0x = V(j,1) + kk.*cos(Btheta);
        line0y = V(j,2) + kk.*sin(Btheta);
        lineB_x(i,:) = cos(Vtheta).*(line0x-V(j,1)) - sin(Vtheta).*(line0y-V(j,2)) + V(j,1);
        lineB_y(i,:) = sin(Vtheta).*(line0x-V(j,1)) + cos(Vtheta).*(line0y-V(j,2)) + V(j,2);
    end
    temp = interp2(xx, yy, img, lineB_x, lineB_y, 'linear',0);
    t = temp';
    P(j,:) = sum(t);
end
figure;imshow(P,[]);colorbar;title('fanbeam projection')

效果如下：
原图：

效果：