2017-03-14 10:55:59 qq278672818 阅读数 3063
使用方法: 1、R = radon(I, theta) 返回亮度图像在角度theta下的Radon变换R。 Radon变换是一幅图像在一个特定的角度下的径向线方向的投影。如果theta是一个标量,R则是一个包含在theta的列向量。如果theta是一个向量,R则是一个矩阵,据真的每一列是对应其中一个theta的Radon变换。如果忽略掉theta,则其默认为0:179. 2、[R,xp] = radon(...) 对应于R中的每一行,返回一个包含径向坐标的向量xp。xp中的径向坐标是沿着X’轴的数值,其为在theta下,X’轴逆时针方向映射来的。两个坐标系的原点为图像的中心点,且为floor((size(I)+1)/2),例如在一个20-by-30的图像中,其中心点为(10,15)。
2016-06-29 13:15:38 lpsl1882 阅读数 11470

  从一个角度,用光源照射对象物体,屏幕上会形成对象物体的影子;如果物体是半透明的,那么影子便有灰度而不是纯黑的,这说明屏幕上的像可以反映物体内部对可见光的衰减作用。我们从落于[0~π]的一系列连续角度照射物体,形成一系列的像,这些像包含物体结构特征信息,基本上可以通过这些像还原物体的形状特征,如果物体是半透明的,那么物体内部的结构也可以还原出来。物体原始形状变换生成这些投影像,称为radon变换;从这些像还原物体形态,称为逆radon变换。人体对可见光是不透明的,但对X光是半透明的,因此CT可以发射X光照射人体,生成人体内部结构的图像信息。
  radon变换的公式是:

xcosθ+ysinθ=ρg(ρ,θ)=f(x,y)δ(xcosθ+ysinθρ)dxdy

本质上就是沿着xcosθ+ysinθ=ρ确定的多条平行射线,建立法线方向为θ的线段上的投影。一个比较简单的实现是,让图像均匀旋转θ角度,然后计算x轴上的投影。假设角度区间[0~π]被分割为n份,投影线段长m,则最终我们获取m×n大小的二维矩阵,称为投影矩阵。

  radon反变换的公式是:

f(x,y)=π0g(xcosθ+ysinθ,θ)dθ
该反变换操作比较简单,但是计算量大,且输出图像模糊有光晕。

  为了得到清晰的图像,我们需要进行频域滤波。首先我们知道二维傅里叶变换对为:

F(u,v)=f(x,y)ej2π(ux+vy)dxdyf(x,y)=F(u,v)ej2π(ux+vy)dudv

这里引入傅里叶切片定理,其中ω是频率分量:
G(ω,θ)=g(ρ,θ)ej2πωρdρ=f(x,y)δ(xcosθ+ysinθρ)ej2πωρdxdydρ=f(x,y)[δ(xcosθ+ysinθρ)ej2πωρdρ]dxdy=f(x,y)ej2πω(xcosθ+ysinθ)dxdy=F(ωcosθ,ωsinθ)

这说明一个投影的一维傅里叶变换,是二维投影矩阵的二维傅里叶变换的一个切片。上述最后一步执行换元操作。
频域逆变换为:

f(x,y)=F(ωcosθ,ωsinθ)dωcosθdωsinθdωcosθdωsinθ=ωdωdθf(x,y)=F(ωcosθ,ωsinθ)ωdωdθ=2π00G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθcos(θ+π=cosθ,sin(θ+π=sinθG(ω,θ+π)=f(x,y)ej2πω(xcosθysinθ)dxdy=f(x,y)ej2π(ω)(xcosθ+ysinθ)dxdy=G(ω,θ)f(x,y)=π00G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+π00G(ω,θ)ej2π(ω)(xcosθ+ysinθ)ωdωdθlet(ωt)f(x,y)=π00G(ω,θ)ej2πω(xcosθ+ysinθ)ωdωdθ+π00G(t,θ)ej2πt(xcosθ+ysinθ)tdtdθ=π0G(ω,θ)ej2πt(xcosθ+ysinθ)|ω|dωdθ

我们可以看出上述式子推导过程中利用了[0~2π]的数据,虽然最后只需要[0~π]的投影数据,但是最终反滤波结果多了一项|ω|。这里面有什么更深层次的数学原理,我确实是不知道的。相比直接空间积分的结果,傅里叶切片重建效果更好的原因,可能有:

  • 我们之前进行radon投影时,看到[0~π]和[π~ 2π]的θ角度下得到的投影是互为镜像,粗看认为信息是重复的,于是贸然舍弃一半,只计算了[0~π]角度的投影数据,这可能隐性地造成了数据不全。
  • 直接进行空间积分变换,可能造成不同角度投影之间发生混淆,导致产生模糊和光晕

在最终式子中,积分计算|ω|项是不可能的,该项不可积。我们可以通过引入窗函数,截断|ω|项或者使用其他近似的窗函数,计算积分计算式并滤波,从而得到一个相对较好的结果。

滤波反投影重建逆radon变化步骤如下:

  • 计算每个投影的一维傅里叶变换
  • 使用截断的|ω|项或者类似的窗口函数进行滤波,得到新的一维傅里叶变换数据
  • 计算傅里叶反变换,获得原图

原图
简单滤波
简单滤波的结果
matlab
Matlab重建结果,注意有条纹

Matlab的radon逆变换效果好,Matlab的牛逼只有研究过的人才懂,自己实现不了的痛也只有研究过的人才懂。。

2018-06-03 21:44:38 corilei 阅读数 5044

radon校正

Radon(拉东)算法是一种通过定方向投影叠加,找到最大投影值时角度,从而确定图像倾斜角度的算法。具体过程如图所示 

                    

拉东变换

若函数F表示一个未知的密度,对F做radon变换,相当于得到F投影后的讯号,举例来说:F相当于人体组织,断层扫描的输出讯号相当于经过radon变换的F。 因此,可以用radon反变换从投影后的密度函数,重建原始的密度函数,它也是重建断层扫描的数学理论基础,另一个被广为人知名词的是三维重建 
radon变换后的讯号称作“正弦图”,因为一个偏离中心的点的radon变换是一个正弦曲线。所以对一些小点的radon变换,会看起来像很多不同振福、相位的正弦函数重叠在一起。


    matlab实现

clear all
clc
close all

[inputfilename,dirname] = uigetfile('*.*');
inputfilename = [dirname, inputfilename];
im = imread(inputfilename); % For example: 'input.jpg'

grayImage = rgb2gray(im);
%%%%%

%%%%% Edge detection and edge linking....
binaryImage = edge(grayImage,'canny'); % 'Canny' edge detector
binaryImage = bwmorph(binaryImage,'thicken'); % A morphological operation for edge linking
%%%%%

%%%%% Radon transform projections along 180 degrees, from -90 to +89....
% R: Radon transform of the intensity image 'grayImage' for -90:89 degrees.
% In fact, each column of R shows the image profile along corresponding angle. 
% xp: a vector containing the radial coordinates corresponding to each row of 'R'.
% Negative angles correspond to clockwise directions, while positive angles
% correspond to counterclockwise directions around the center point (up-left corner).
% R1: A 1x180 vector in which, each element is equal the maximum value of Radon transform along each angle.
% This value reflects the maximum number of pixels along each direction. 
% r_max: A 1x180 vector, which includes corresponding radii of 'R1'.
theta = -90:89;
[R,xp] = radon(binaryImage,theta);
imagesc(theta,xp, R); colormap(jet);
xlabel('theta (degrees)');ylabel('x''');
title('theta方向对应的Radon变换R随着x''的变化图');
colorbar
%%%%%

[R1,r_max] = max(R);
theta_max = 90;
while(theta_max > 50 || theta_max<-50)
    [R2,theta_max] = max(R1); % R2: Maximum Radon transform value over all angles. 
                              % theta_max: Corresponding angle 
    R1(theta_max) = 0; % Remove element 'R2' from vector 'R1', so that other maximum values can be found.
    theta_max = theta_max - 91;
end

correctedImage = imrotate(im,-theta_max); % Rotation correction
correctedImage(correctedImage == 0) = 255; % Converts black resgions to white regions

subplot(1,2,1), subimage(im)
subplot(1,2,2), subimage(correctedImage)

function [bw,qingxiejiao]=radontran(bwbone,bw)%radon倾斜校正
 
 theta=1:90;    
 [R,xp]=radon( bwbone,theta);
 [I0,J]=find(R>=max(max(R)));%找倾斜角
qingxiejiao=90-J;
 bw=imrotate(bw,qingxiejiao,'bilinear','crop');
clc;
clear all;
close all;
[fn pn fi]=uigetfile('*.*','choose a picture');
Img=imread([pn fn]);
imshow(Img);title('Original image');
I = rgb2gray(Img);
I=improve_hist(I);
bw=edge(I,'canny');
theta=1:179;
[R,xp]=radon(bw,theta);
[I0,J]=find(R>=max(R(:)));%J记录了倾斜角
qingxiejiao=90-J;
I1=imrotate(Img,qingxiejiao,'bilinear','crop');
subplot(1,2,1),imshow(Img);title('Original image');
subplot(1,2,2),imshow(I1);title('correct image');

其中

B = imrotate(A,angle)
将图像A(图像的数据矩阵)绕图像的中心点旋转angle度, 正数表示逆时针旋转, 负数表示顺时针旋转。返回旋转后的图像矩阵。
B = imrotate(A,angle,method)
使用method参数可以改变插值算法,method参数可以为下面这三个值:
'nearest':最邻近线性插值(Nearest-neighbor interpolation)
'bilinear': 双线性插值(Bilinear interpolation)
'bicubic': 双三次插值(或叫做双立方插值)(Bicubic interpolation)
B = imrotate(A,angle,method,bbox)
bbox参数用于指定输出图像属性:
'crop': 通过对旋转后的图像B进行裁剪, 保持旋转后输出图像B的尺寸和输入图像A的尺寸一样。
'loose': 使输出图像足够大, 以保证源图像旋转后超出图像尺寸范围的像素值没有丢失。 一般这种格式产生的图像的尺寸都要大于源图像的尺寸。

2019-06-24 23:07:04 nanhuaibeian 阅读数 384

函数 radon 用来对给定的二维矩形数组生成一组平行射线投影。

从而实现前面所说的由投影来重建图像

一、概念

函数 radon 用来对给定的二维矩形数组生成一组平行射线投影。

该函数的基本语法:R = radon(I,theta),其中,I 是一个二维数组,theta 是角度值的一维数组。投影包含在 R 的列中,生成的投影数等于数组 theta 中的角度数。生成的投影长到足以在射线约束旋转时跨越观察的宽度。当射线垂直于矩形数组的主对角线时会出现这种视图。换句话说,对一个大小为 MxN 的输入数组,投影的最小长度为 [M^2 + N^2]^(1/2)
当然,其他角度的投影事实上要短得多,且它们要用零来填充,以便所有投影的长度都相同(如要求的那样,R 应是一个矩形数组)。由函数 radon 返回的实际长度,要稍大于每个像素的单位面积计算的主对角线的长度。
函数 radon 一个更一般的语法:[R,xp] = radon(I,theta),其中 xp 包含沿着 x’ 轴的坐标值,xp 中的值用于标注图轴。

在 CT 算法模拟中,用来生成一幅著名图像(称为 Sheep-Logan头部幻影)的一个有用函数的语法为:P = phantom(def,n)
其中,def 是一个指定所生成头部幻影类型的字符串,n 是行数和列数(默认值为256)

def 的有效值:
‘Sheep-Logan’:计算机断层研究 人员广泛使用的测试图像。这幅图像的对比度很低。
‘Modified Shepp-Logan’:Sheep-Logan 幻影的变体,其对比度得到了改进,因此有更好的视觉效果

二、使用函数 radon

% 创建一个全零的数组,600 x 600
>> g1 = zeros(600,600);
% 将 g1 第100-500行,250-350 列赋值为1
>> g1(100:500,250:350) = 1;
% 生成 Shepp-Logan 头部幻影
>> g2 = phantom('Modified Shepp-Logan',600);
>> imshow(g1),figure,imshow(g2)

在这里插入图片描述
使用半度增量的雷登变换由下面的语句得到:

% 使用 0.5 的角度增量,从 0 度到175.5度为增量
>> theta = 0:0.5:179.5;
% 对给定的二维矩形数组生成一组平行射线投影
>> [R1,xp1] = radon(g1,theta);
>> [R2,xp2] = radon(g2,theta);

此时,R1 的第一列是 θ = 0度,第二列是 θ = 0.5度,等等。
第一列的第一个元素对应于 ρ 的最小的负值,而最后一个元素对应于 ρ 的最大的正值。其他列类似
如果我们想要显示 R1 以便投影从左到右进行,并且第一个投影出现在图像的底部,必须转置并反转数组

% 先转置数组,然后翻转数组
>> R1 = flipud(R1');
>> R2 = flipud(R2');
% 设置显示图像的参数
>> figure,imshow(R1,[],'XData',xp1([1 end]),'YData',[179.5 0])
% 将坐标轴系统的原点由其默认的左上角移到右下角
>> axis xy
% 显示 X 轴和 Y 轴
>> axis on
% 设置 X轴为 θ轴, Y 轴为 ρ 轴
>> xlabel('\rho'),ylabel('\theta')

R1 的雷登变换图
在这里插入图片描述

>> figure,imshow(R2,[],'XData',xp2([1 end]),'YData',[179.5 0])
>> axis xy
>> axis on
>> xlabel('\rho'),ylabel('\theta')

R2 的雷登变换图
在这里插入图片描述
最后的这两幅图显示了结果,在这两幅图像中,每行表示了固定值 θ 的一个完整投影。

2013-12-19 11:32:47 utimes 阅读数 8282

1917年澳大利亚数学家Radon首先论证了通过物体某一平面的投影重建物体该平面两维空间分布的公式。他的公式要求获得沿该平面所有可能的直线的全部投影(无限集合)。所获得的投影集称为Radon变换。由Radon变换进行重建图像的操作则称为逆Radon变换。Radon变换和逆Radon变换对CT成像的意义在于,它从数学原理上证实了通过物体某一断层层面“沿直线衰减分布的投影”重建该层面单位体积,即体素的线性衰减系数两维空间分布的可能性。关于更详细参考[1][2].

clear all;
close all;
clc;

sourcPic = imread('test.jpg');
grayPic =rgb2gray(sourcPic);
ed = edge(grayPic,'canny');
subplot(2,2,1);
imshow(ed),title('sourcPic');

r  = radon(ed,15);
subplot(2,2,2);
plot(r),title('Radon 15');

r1 = radon(ed,30);
subplot(2,2,3);
plot(r1),title('Radon 30');

r2 = radon(ed,45);
subplot(2,2,4);
plot(r2),title('Radon 45');
实验结果

参考资料

[1] Radon Transform From Wikipedia, the free encyclopedia.

[2] Radon Transform From Mathoworld Wolfrom.

=========

关于Image Engineering & Computer Vision的更多讨论与交流,敬请关注本博和新浪微博songzi_tea.


Radon定理与证明

阅读数 448

拉登/Radon变换

阅读数 3159

没有更多推荐了,返回首页