胡拉八扯之Radon变换,断层影像重建

新用户注册 | 用户登陆 | 刷新

CoolingRib


发表文章数: 227
内力值: 285/285
贡献度: 1485
人气: 147

论坛嘉宾

胡拉八扯之Radon变换,断层影像重建 [文章类型: 原创]

看到观星楼有人讨论反问题,提到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-ray scaner一起完成的,显然不干软件什么事。接下来,照理说是要靠计算机
把获得的数据做一个逆Radon变换,就能得到被检测截面的X-ray吸收系数的分布图像了。CT的图像其实就是一个吸收系数的图,类似的B超或者声纳之类的图像是大致是一个弹性模量的图(反射声波)...

但是接下来这里有一个问题就是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.
当L越来越小的时候,频域内周期越来越大,空间分辨率越高。但是当L为有限的时候,分辨率如果用频率来表示的话,从原点(“直流”分量)开始算,由于周期性缘故显然最高到1/2L处。设想一间黑屋子,唯一的光源是一个可调节频率的频闪光源,一台电风扇。
首先,当光源的闪烁频率和风扇的转速(用转/秒来表述)相等的时候,你将看到风扇是停止的,当光源频率高于风扇转速的两倍时,你才能看到风扇正常的转动,如果光源频率介于风扇转速一倍和两倍之间,那你会看到风扇倒着转了。这里的情况被称为频谱混叠。
另外,函数变换本身还带来了坐标平移一类的问题。
实际当中用的最多的是一种叫做滤波反投影的算法来实现断层重建,说穿了关键就针对以上一些问题设计合理的滤波器。

另外值得一提的是,这里用到的数学大概一百年前就有了,但是随着计算机技术的进步,具体实现的时候,出现过不同的版本。譬如说,70年代的商业运行的CT(256X256),带一台长得像电冰箱般的“卷积器(convolver)”,顾名思义,它的滤波器实现大概是用DSP做卷积的(离散的卷积就是一系列的移位连乘连加)。而现代,随着硬件技能的突飞猛进,FFT不成问题了,这个交给CPU在频域内作乘法就
能搞定。退一步说,我甚至怀疑,那个形体巨大的Convolver做卷积的性能恐怕未必能赶上我正在码字的电脑。此刻,它正在运行音乐播放软件foorbar,同时一起运行的还有一堆插件(也可看作卷积器),比如老式电子管音色,教堂混响,耳机模拟现场音效之类的...

以上这些基本上是相关领域的abc,没有深入,基本凭借记忆,说法可能和专门的教材不完全一样,而且很多地方一知半解,肯定会有谬误,大家随便看看不可当真,当然欢迎拍砖

- 错误提示 -

很抱歉,您仍然不具有浏览本相片的权限

发表时间: 2007-08-10, 06:11:26 个人资料

CoolingRib


发表文章数: 227
内力值: 285/285
贡献度: 1485
人气: 147

论坛嘉宾

Re: 胡拉八扯之Radon变换,断层影像重建 [文章类型: 原创]

仔细读了一遍,发现有些地方未写清楚。以下文字可替换:

“另外,检测器只能读出空间上分立的数值,所谓的取样过程就是投影函数乘一个迪拉克函数(假设周期为L)序列。而迪拉克序列变换到频域仍然是一个迪拉克序列,只是周期变成了1/L.
当L越来越小的时候,频域内周期越来越大,空间分辨率越高。但是当L为有限的时候,分辨率如果用频率来表示的话,从原点(“直流”分量)开始算,由于周期性缘故显然最高到1/2L处。设想一间黑屋子,唯一的光源是一个可调节频率的频闪光源,一台电风扇。
首先,当光源的闪烁频率和风扇的转速(用转/秒来表述)相等的时候,你将看到风扇是停止的,当光源频率高于风扇转速的两倍时,你才能看到风扇正常的转动,如果光源频率介于风扇转速一倍和两倍之间,那你会看到风扇倒着转了。这里的情况被称为频谱混叠。

替换为:

另外,检测器只能读出空间上分立的数值,所谓的取样过程就是投影函数乘一个迪拉克函数组成的序列(假设周期为L) 而迪拉克序列变换到频域仍然是一个迪拉克序列,只是周期变成了1/L。投影和取样序列相乘在频域就是卷积,出来的结果就是具有了周期频谱,显然可用的只能是原点(DC)所在的一个周期内的数据。当L越来越小的时候,频谱周期越来越大,空间分辨率越来越高。当L为有限的时候,分辨率如果用频率来表示的话,从原点(“直流”分量)开始算,由于周期性缘故显然最高到1/2L处。

设想一间黑屋子,唯一的光源是一个可调节频率的频闪光源,一台电风扇。假定光源闪烁频率为w,显然理论上能够检测到的风扇转速u将允许加上任意整数个w。比方说,每秒亮一下,你看到了风扇转动了1/4圈,那么你可以认为风扇每秒转动1/4圈,但也可以是5/4圈(多转了一圈,有何不可?),9/4圈...也可以是(反着转)-3/4圈,-7/4圈...原因就是前面说的,用一个脉冲序列(光源频闪)去做取样,必然会得到周期性的频谱。接下来,当光源的闪烁频率和风扇的转速(用转/秒来表述)相等的时候,你将看到风扇是停止的,当光源频率高于风扇转速的两倍时,你才能看到风扇正常的转动,如果光源频率介于风扇转速一倍和两倍之间,那你会看到风扇倒着转了。这里的情况被称为频谱混叠。此类现象生活中常遇到。

- 错误提示 -

很抱歉,您仍然不具有浏览本相片的权限

发表时间: 2007-08-11, 02:51:30 个人资料
您尚未登陆 | 用户登陆