radon变换 两维情况下radon变换大致可以这样理解:一个平面内沿不同的直线(直线与原点的距离为d,方向角为alfa)对f(x,y)做线积分,得到的像F(d,alfa)
就是函数f的Radon变换。也就是说,平面(d,alfa)的每个点的像
函数值对应了原始函数的某个线积分值。一个更直观的理解是,假设你的手指被一个很强的平行光源透射,你迎着光源看到的手指图像就是手指的光衰减系塌尘数的三维Radon变换(小小的推广)在给定方向(两个角坐标)的时候的值, 一个最简单而直接的应用就是拿来检测图像里面含有的直线成分,很显然地,任何直线都会导致Randon像在该直线对应(d,alfa)处的极值。 具体的CT断层团逗禅影像重建算法当中其实没怎么用到Radon变换,或者说Radon变换仅仅只有一点点理论上的意义。原因是:CT机做扫描:球管发出X-ray,经过人体,被吸收一部分,进入检测器队列(球管是旋转的,检测器呈扇形分布,很老的和很新的除外,老式的ct做平行扫描,效率低,很新式的什么多层螺旋扫描,我也不知道咋回事)显然检测器读数就是人体的x-ray吸收系数(是空间的函数)对相应路径的线积分,所以这样转一圈下来再把所有的检测器读数值按照(d,alfa)的方式排列一下就算完成了某个被检测截面的Radon变换了,这个过程是人体和X-rayscaner一起完成的,显然不干软件什么事。接下来,照理说是要靠计算机把获得的数据做一个逆Radon变换,就能得到被检测截面的X-ray吸收系数的分布图像了。CT的图像其实就是一个吸收系数的图,类似的B超或者声纳之类的图像是大致是一个d性模量的图(反射声波)... 但是接下来这里有一个问题就是Radon变换是不是可逆,google了下好像是可逆的,我的理解: 1)有另外一种求逆方法,就是解代数方程,简化地说可以大致设想把整个截面离散网格化,每个格子对应一个吸收系数,把每根扫描积分路径经过的格子按照权重(显然透心凉和擦点皮对吸收的贡献不同)作累加,令他们等于相应Radon值(积分变成了加权累加而已)显然设计好的话,这个方程组肯定是有解的(不过运算量会很庞大,比如一个512X512的网格...) 2)工程师不问这么无聊不切实际的问题,所以以前学的时候就压根没想到。 3)最重要的原因,是下面要说的求逆问题,竟然变成了二维的fourier逆变换。所以忘掉Radon变换吧。 有这样一个事实:把某个角度坐标alfa对应的一“条”Radon值(一系列检测器的读数,也实际上就是原始截面吸收系数在方向为alfa+-Pi/2直线簇上的线积分值)作一个fourier变换,得到的就是整个原始被检测截面(吸收系数)的二维fourier像在某条直线上的值(这条直线经过频域的原点并且方向为alfa)如果把所有角度的Radon值作一维Fourier变换,然后按照合适的角度(alfa)经过原点把这些一维fourier像值放在频域平面上,就得到了整个二维fourier像!!!这个其实直观上很容易想象其合理性,还是以手指头为例(不考虑它指向的方向)对着光源看,从左至右,透光率不同产生明暗的变化,亮暗本身是沿前后方向的积分结果决定,但是相邻的亮暗变化却反应了整个手指截面的从左至右这个方向上的频域信息,看到的细节越多,频域的高频分量越多(与前后方向的细节毫无关系,因为被radon积分掉了)。 以上关于CT其实是过分简化的描述,只能提供一个大致的原理。实际情况会有些不同,首先检测器读数是有限空间的,这就是相当于理想的投影函数乘了一个窗函数(某段区间内为一,其他地方为零的函数),在频域内窗函数会“扩散”所以他们频域做
卷积的结果是频域的扩展。也可以说成是,对于非周期函数(包括周期不为无穷大)的fourier级数在边界的间断处只能是平均收敛,“平均”的结果就是在光滑的地方拟合的很好,在间断点处发生振荡。工程中管这个叫做吉布斯(Gibbs)效应,它告诉我们:用有限项指消级数的和去表示一个函数,随着项数的增加,振荡发生的位置会越来越接近间断点,但是它的摆幅不变(写到这忽然觉得它的名字似乎翻译成“挤不死”更贴切)另外,检测器只能读出空间上分立的数值,所谓的取样过程就是投影函数乘一个迪拉克函数组成的序列(假设周期为L)而迪拉克序列变换到频域仍然是一个迪拉克序列,只是周期变成了1/L。投影和取样序列相乘在频域就是卷积,出来的结果就是具有了周期频谱,显然可用的只能是原点(DC)所在的一个周期内的数据。当L越来越小的时候,频谱周期越来越大,空间分辨率越来越高。当L为有限的时候,分辨率如果用频率来表示的话,从原点(“直流”分量)开始算,由于周期性缘故显然最高到1/2L处。 设想一间黑屋子,唯一的光源是一个可调节频率的频闪光源,一台电风扇。假定光源闪烁频率为w,显然理论上能够检测到的风扇转速u将允许加上任意整数个w。比方说,每秒亮一下,你看到了风扇转动了1/4圈,那么你可以认为风扇每秒转动1/4圈,但也可以是5/4圈(多转了一圈,有何不可?),9/4圈...也可以是(反着转)-3/4圈,-7/4圈...原因就是前面说的,用一个脉冲序列(光源频闪)去做取样,必然会得到周期性的频谱。接下来,当光源的闪烁频率和风扇的转速(用转/秒来表述)相等的时候,你将看到风扇是停止的,当光源频率高于风扇转速的两倍时,你才能看到风扇正常的转动,如果光源频率介于风扇转速一倍和两倍之间,那你会看到风扇倒着转了。这里的情况被称为频谱混叠。此类现象生活中常遇到。另外,函数变换本身还带来了坐标平移一类的问题。实际当中用的最多的是一种叫做滤波反投影的算法来实现断层重建,说穿了关键就针对以上一些问题设计合理的滤波器。 另外值得一提的是,这里用到的数学大概一百年前就有了,但是随着计算机技术的进步,具体实现的时候,出现过不同的版本。譬如说,70年代的商业运行的CT(256X256),带一台长得像电冰箱般的“卷积器(convolver)”,顾名思义,它的滤波器实现大概是用DSP做卷积的(离散的卷积就是一系列的移位连乘连加)。而现代,随着硬件技能的突飞猛进,FFT不成问题了,这个交给CPU在频域内作乘法就能搞定。退一步说,我甚至怀疑,那个形体巨大的Convolver做卷积的性能恐怕未必能赶上我正在码字的电脑。此刻,它正在运行音乐播放软件foorbar,同时一起运行的还有一堆插件(也可看作卷积器),比如老式电子管音色,教堂混响,耳机模拟现场音效之类的... 以上这些基本上是相关领域的abc,没有深入,基本凭借记忆,说法可能和专门的教材不完全一样,而且很多地方一知半解,肯定会有谬误,大家随便看看不可当真,当然欢迎拍砖。theta=0:180theta是radon变换的角度,可以自由设置,这里是从0°到180°一共181个角度,如果只求等分的60个角度的radon变换,可以改成theta=0:3:177;如果只求一个特殊角度的radon变换可以写成theta=x这里x就是希望的角度兆察。
R = radon(I, theta)
R是存储radon变换的值,它是一个矩阵,列数是theta的个数,表示每一个角度生拦猜纳成一列,行数是被处理的矩阵(I)对角线的长度。
R一般用于inverse radon transform,程序如下:
I1=iradon(R,theta)
imshow(I1)
iradon命令有简没很多参数,具体使用方法请参考http://www.mathworks.com/access/helpdesk/help/toolbox/images/index.html?/access/helpdesk/help/toolbox/images/iradon.html&http://www.mathworks.com/cgi-bin/texis/webinator/search?pr=Whole_site&db=MSS&prox=page&rorder=750&rprox=750&rdfreq=500&rwfreq=500&rlead=250&sufs=0&order=r&whole=Whole_site&entire_flag=1&is_summary_on=1&ResultCount=10&query=inverse+radon+transform
:)
Radon 变换是平行束对图像的线积分,根据各个角度得到的一系列投影值逆radon重建得到原始图像。变换角度默认是逆时针。r=radon(im,30)得到的是一维数组,平行束与X轴夹角为30度时,升稿脊距原点不同距离的投影线(平行束)上对图像的线积分吵渗。[R,Xp] = RADON(...) XP对应平行束的位置。敬埋
评论列表(0条)