【图像处理】平行线投影radon变换

原创 2016年06月29日 13:15:38

  从一个角度,用光源照射对象物体,屏幕上会形成对象物体的影子;如果物体是半透明的,那么影子便有灰度而不是纯黑的,这说明屏幕上的像可以反映物体内部对可见光的衰减作用。我们从落于[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的牛逼只有研究过的人才懂,自己实现不了的痛也只有研究过的人才懂。。

版权声明:本文为博主原创文章,未经博主允许不得转载。

相关文章推荐

Refactor反编译C#程序

两篇很不错reflector的教程,很有用~ Reflactor就不介绍了,.net下的免费反编译工具。Reflactor在某些情况下是很必须的,当遇到bug的时候,可以通过Reflactor看其内...

ajax入门 不要畏惧 很简单 进了门一切都好学多了

以前总是听别人说ajax是多么的好,然后自己就去借了本书看,哇塞感觉好难哦,什么介绍JavaScript,html,css,还有很多一些东西。看的那个难啊,然后就是硬着头皮把它给看完了,但是就是缺少了...
  • y75475
  • y75475
  • 2017年04月28日 17:37
  • 1053

【opencv+python】图像处理之二、几何变换(仿射与投影)的原理

该系列文章为 OpenCV+Python Tutorials的学习笔记 代码托管在Github 转载请注明: http://blog.csdn.net/a352611/article/detail...
  • a352611
  • a352611
  • 2016年05月15日 17:55
  • 2292

【opencv+python】图像处理之二、几何变换(仿射与投影)的应用

该系列文章为 OpenCV+Python Tutorials的学习笔记 代码托管在Github 转载请注明: http://blog.csdn.net/a352611/article/detail...
  • a352611
  • a352611
  • 2016年05月15日 20:34
  • 6602

图像处理之(直方图)反向投影

反向投影是一种记录像素点或xiangsud

OpenCV之imgproc 模块. 图像处理(4)直方图均衡化 直方图计算 直方图对比 反向投影 模板匹配

直方图均衡化 目标 在这个教程中你将学到: 什么是图像的直方图和为什么图像的直方图很有用用OpenCV函数 equalizeHist 对图像进行直方图均衡化 原...

Python3与OpenCV3.3 图像处理(十三)--反射投影

一、什么是反射投影 简单的说就是通过给定的直方图信息,在图像找到相应的像素分布区域 二、反射投影的应用 物体跟踪、定位物体等 三、示例代码 import ...

Win8 Metro(C#)数字图像处理--2.39二值图像投影

[函数名称]   二值图像投影         ImageProjection(WriteableBitmap src) [算法说明] [函数代码] /// ///...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:【图像处理】平行线投影radon变换
举报原因:
原因补充:

(最多只允许输入30个字)