在医学图像处理中,往往通过对某切面作多个X射线投影,来获得切面的结构图形,这就是所谓的图像重建。图像重建有好多种方法,但实际上,当人们在处理二维或三维投影数据时,真正有效的重要算法都是以Radon变换或Radon逆变换作为数学基础的。
在图像变换中,Radon变换的实质就是计算图像矩阵在某个方向上的投影。例如图像矩阵如果是二维的,在某个方向上的投影,就是在其垂直方向上的线积分。如下图:
Radon变换的几何示意图——投影
事实上,Radon变换最常用的是其逆变换。因为一个物体的投影可以通过光线照射得到,如果我们能得到该物体不同方向上的投影,就可以对这些投影图像进行Radon逆变换,来重建其二维结构或三维结构。
在MATLAB中,有专门的radon()函数和iradon()函数来实现Radon变换的算法,这也是图像重建过程中必不可少的算法。
(1)在MATLAB图像处理工具箱中提供了radon()函数用于计算指定方向上图像的投影,该函数的调用格式如下:
R = radon(I, theta):计算图像I在矢量theta指定方向上的radon变换;
[R,xp] = radon(...):R的各行返回theta中各方向上的radon变换,xp表示theta方向对应的坐标值。图像I的中心,即floor((size(I)+1)/2),为新坐标轴的中心。
(2)在MATLAB图像处理工具箱中提供了iradon()函数用于实现Radon逆变换,完成图像重建,该函数的调用格式如下:
I = iradon(R, theta):R是投影矩阵;theta可以是一个包含所有扫描角度的向量,且每两个相邻角度等间隔;使用的投影越多,所获得的图像越接近原图像,失真越小;
I = iradon(P, theta, interp, filter, frequency_scaling, output_size):interp是插值函数;filter是滤波函数,通过加窗消去投影过程中产生的高频噪声;frequency_scaling是一个标量值,取值范围[0,1],通过缩放滤波函数的频率修改滤波函数;output_size是一个标量,用来规定重建图像的行数和列数;
[I,H] = iradon(...):H的返回值为滤波器的频率响应。
interp是插值函数,有以下几种差值方式可以选择:
filter是滤波函数,有以下几种滤波器可以选择:
‘Ram-Lak’:频带有限的斜坡函数滤波器,是MATLAB默认的滤波器,对投影中的噪声敏感(R-L滤波函数和S-L滤波函数);
‘Shepp-Logan’:sinc函数*R-L函数。
‘Cosine’:cosine函数*R-L函数 。
‘Hamming’:Hamming函数*R-L函数 。
‘Hann’:hann函数*R-L函数 。
‘None’:没有滤波。
下面通过具体实例来说明Radon变换的含义和操作方法:
代码:
Demo1
Demo2
显示效果:
Demon1.投影——即线积分
Demo1.上述黑白图像的Rando变换
Demo1
可以看出,图像的生成是以(theta,xp)构建一个平面坐标,然后R矩阵中的数据作为该平面各对应点的高度,这样就构成一个三维矩阵,从而显示为上述的Radon变换图像。
Demo2_figure1
Demo2_figure2
Demo2_figure3
可以看到,在角度是0度、30多度、90多度和140多度的位置,有明显的亮点,说明该角度的投影——即垂直于该角度的线积分,值比较大。所以,可以猜测投影方向上存在直线段。例如在90度时,其投影方向为水平方向,有大量的白点,说明原图像中存在大量水平的直线段。我们从原图像可以看到,确实存在较多的水平直线段,其重建图像也表明了这一点。
简单来说,就是Radon变换和我上一篇文章提到的Hough一样,都可以作为检测图像中直线段的方法。
Demo2_figure4
Demo2_figure5
Demo2_figure6
从上面可以看出,确实是投影图像越多,其重建图像失真越小。当投影角度间隔为1度或2度时,重建图像失真较小;而当影角度间隔为5度或10度时,重建图像失真较大,图像信息丢失严重。
「十」MATLAB语言之程序性能优化
「二十九」MATLAB数字图像运算之块操作
「三十二」MATLAB图像变换之Hough变换——寻找图像中的直线段