CT原理与图像重建.doc

Radon变换入门matlab CT原理

简介
图像投影，就是说将图像在某一方向上做线性积分（或理解为累加求和）。如果将图像看成二维函数f(x, y)，则其投影就是在特定方向上的线性积分，比如f(x, y)在垂直方向上的线性积分就是其在x轴上的投影；f(x, y)在水平方向上的线积分就是其在y轴上的投影。通过这些投影，可以获取图像在指定方向上的突出特性，这在图像模式识别等处理中可能会用到。
Radon变换（拉东变换），就是将数字图像矩阵在某一指定角度射线方向上做投影变换。这就是说可以沿着任意角度theta来做Radon变换。
实例
% 实验Radon变换
% By lyqmath
% Dalian University of Technology
% School of Mathematical Sciences
clc; clear all; close all;
I = zeros(256, 256);
[r, c] = size(I);
I(floor(1/5*r:4/5*r), floor(3/5*c:4/5*c)) = 1;
figure;
subplot(2, 2, 1); imshow(I); title('原图像');
[R, xt] = radon(I, [0 45 90]); % 在0、45、90度方向做radon变换
subplot(2, 2, 2); 
plot(xt, R(:, 1));
title('水平方向的radon变换曲线');
subplot(2, 2, 3); 
plot(xt, R(:, 2));
title('45度方向的radon变换曲线');
subplot(2, 2, 4); 
plot(xt, R(:, 3));
title('垂直方向的radon变换曲线');

结果：

总结
由于radon变换将图像变换到按角度投影区域，和有名的hough类似，可以应用与检测直线。个人认为，通过将图像矩阵在多角度做积分投影，再对得到的数据做统计分析，可以确定出图像的一些基本性质。

==================附======================

I = zeros(100,100);
I(25:75, 25:75) = 1;
imshow(I)
[R,xp] = radon(I,[0 45]);
figure; plot(xp,R(:,1)); title('R_{0^o} (x\prime)')
theta = 0:180;
[R,xp] = radon(I,theta);
imagesc(theta,xp,R);
title('R_{\theta} (X\prime)');
xlabel('\theta (degrees)');
ylabel('X\prime');
set(gca,'XTick',0:20:180);
colormap(hot);
colorbar

结果
原图：

变换后：

=====================逆变换=====================
I = zeros(100,100);
I(25:75, 25:75) = 1;
P=I;
imshow(P)
theta1=0:10:170;[R1,xp]=radon(P,theta1);   %存在18个角度投影
theta2=0:5:175;[R2,xp]=radon(P,theta2);    %存在36个角度投影
theta3=0:2:178;[R3,xp]=radon(P,theta3);    %存在90个角度投影
figure,imagesc(theta3,xp,R3);colormap(hot);colorbar;
xlabel('\theta');ylabel('x\prime');  %  定义坐标轴
I1=iradon(R1,10);
I2=iradon(R2,5);
I3=iradon(R3,2);
figure,imshow(I1),title('I1');
figure,imshow(I2),title('I2');
figure,imshow(I3),title('I3');

==========================反变换====================
I = imread('mm.JPG');
I = rgb2gray(I);
P=I;
imshow(P)
theta1=0:10:170;[R1,xp]=radon(P,theta1);   %存在18个角度投影
theta2=0:5:175;[R2,xp]=radon(P,theta2);    %存在36个角度投影
theta3=0:1:179;[R3,xp]=radon(P,theta3);    %存在90个角度投影
figure,imagesc(theta3,xp,R3);colormap(hot);colorbar;
xlabel('\theta');ylabel('x\prime');  %  定义坐标轴
I1=iradon(R1,10); % 开始反变换，还原图像
I2=iradon(R2,5);
I3=iradon(R3,1);
figure,imshow(I1),title('I1');
figure,imshow(I2),title('I2');
figure,imshow(I3),title('I3');
k=input('Your Thresh(0,1):')
thres1=max(I3(:))*k         %maybe 0.5 is perfect
figure,imshow(I3>thres1),title('I4');
iradon函数的使用
iradon:图像处理命令
功能：进行反radon变换
语法：
I=iradon(P,theta)
I=iradon(P,theta, interp,filter,d,n)
[I,h]=iradon(…)
举例：
P=phantom(128);
R=radon(P,0:179);
I=iradon(R,0:179,’nearest’,’Hann’);
imshow(P)
figure,imshow(I)
结果：

相关命令：iradon phantom