专利名称：小尺度下热声成像的反卷积重建法的制作方法
技术领域：
本发明属于热声成像技术领域，具体为一种小尺度下热声成像的反卷积重建法。
背景技术：
热声成像是一种新近迅速发展起来的无损医学成像方法[1～3]。在生物组织的热声成像中，一束短脉冲的电磁波(通常是微波或激光)[3～5]照射到组织上，组织快速吸收电磁波能量并膨胀，从而在空气中产生压力波。该压力波叫做热声波，通常是一种超声波，它携带了组织对电磁波吸收的分布特性。使用超声换能器在组织周围扫描探测热声波，然后基于相应的算法[4～8]就可以重建出组织内部的电磁波吸收系数的分布图像。这个过程就是热声成像。热声成像也可被称为光声成像，习惯上如果电磁波是远红外光或微波则被称为热声成像，如果是可见光或近红外光则被称为光声成像[9]。热声成像可被应用于探测组织内部血管分布[10，11]、对组织内血管特征进行测定[3，12]等。因此热声成像在肿瘤检测方面有很大的应用价值[10]，可以为肿瘤检测提供高对比度的成像[13]。
热声成像技术的关键问题是如何通过探测到的热声波重建出组织内电磁波吸收系数的分布，即成像算法。Kruger等提出了三维和二维下的滤波反投影成像算法[2，4]。Wang等提出了三维和二维下的时域重建法、频域重建法[5-7]。Kstli等提出了二维的频域上的快速近似算法[8]。
本发明针对超声换能器做圆球面扫描探测的情况，提出了一种三维下的基于反卷积的热声成像算法反卷积重建法。对于小尺度下的热声成像问题，反卷积重建法是一种快速有效的算法。
发明内容
本发明的目的在于提出一种快速有效的、小尺度下热声成像的反卷积重建法。
本发明提出的一种小尺度下热声成像的反卷积重建法，具体步骤为根据在圆周面上扫描的超声换能器探测到的声压函数p(r0，t)，构造出一个新函数B(r)；B(r)经三维傅里叶变换得到
基于反卷积方法由
计算出
经三维反傅里叶变换得到A′(r)；A′(r)的关于原点对称的函数就是电磁波吸收系数分布A(r)。
下面对各步骤作进一步具体描述。
假定电磁波脉冲均匀地照射在样品组织上，电磁波脉冲函数为阿I(t)，样品的电磁波吸收系数分布为A(r)，在空间r处产生的声压p(r，t)为[2，4，5] 其中β是等压膨胀系数，c是声速，Cp是比热。热声成像算法是一个典型的逆问题，即如何由p(r，t)求出A(r)。
通常假定电磁波脉冲函数为δ函数，即I(t)＝δ(t)。超声换能器在半径为r0的圆球面扫描探测，球心为原点，则p(r，t)中r＝r0。式(1)可化简为[2，5]
整理式(2)得
其中常数令
B(r)就是我们构造的函数。整理式(4)与式(5)得
对于小尺度的待测样品，认为式(6)中|r′|max＜＜r0。考察式(6)中的积分面 结合式(7)与式(8)，可将积分面近似为|r+r′|＝2r0。则式(6)可写为
令A′(r)＝A(-r)，则式(9)可写为
容易证明，如果将A′(r)与B(r)分别看作一个系统的输入和输出，则该系统是线性时不变系统。假设有
显然有
其中a与b是任意常数。这就证明了系统的线性性。又根据
就证明了系统的时不变性。因此该系统是线性时不变系统。则 A′(r)*h(r)＝B(r)(14) 其中h(r)是系统的冲激响应。当系统输入为δ函数时，系统输出就是h(r)。有
A′(r)、h(r)和B(r)经过三维傅里叶变换后分别得到
和
有 但是，直接使用式(16)计算
会有一些问题。首先，
可能存在零值点。其次，在
幅值很小的点，
的误差会被放大。一个常用的解决方法是采用下式来计算[14-15] 其中λ是一个常数，可根据实际情况合理选取[14、15]，这样，
经过三维反傅里叶变换得到A′(r)，A′(r)的关于原点对称的函数就是电磁波吸收系数分布A(r)。
图1、四种待测样品的原始图像和重建结果。图中四行从上到下对应了四种由小到大的待测样品，四列从左到右对应了样品原始图像、时域重建法结果、滤波反投影法结果和反卷积重建法结果。其中原始图像周围的白色圆圈表示超声换能器的扫描轨迹。
图2、r′/r0为不同值时三种算法重建图像的PSNR。其中r0是超声换能器扫描半径，r′是待测样品离扫描圆心的最远距离，点划线对应时域重建法(TDR)，实线对应滤波反投影法(FBP)，虚线对应反卷积重建法(DR)。
具体实施例方式 本发明提出的反卷积重建法的流程概括为根据式(4)和式(5)由声压函数p(r0，t)构造出B(r)；B(r)经三维傅里叶变换得到根据式(15)和式(17)计算出
经三维反傅里叶变换得到最后得到电磁波吸收系数分布A(r)＝A′(-r)。
反卷积重建法的前提条件是待测样品是小尺度的，即样品离原点的最远距离|r|max应远小于超声换能器的扫描半径r0。而事实上，根据后文的仿真实验，只要满足|r|max＜0.3r0，本算法就会有很高的精度。反卷积重建法的主要计算步骤是傅里叶变换和反傅里叶变换，通过应用FFT(快速傅里叶变换)和IFFT(快速反傅里叶变换)算法，使得反卷积重建法有较快的计算速度。
在计算机上使用Matlab进行二维情况下的仿真实验。虽然前文对反卷积重建法的论证是在三维情况下进行的，但将其应用到二维下显然也是成立的。基于超声换能器做圆周扫描探测的情况，将反卷积重建法与当前流行的时域重建法、滤波反投影法进行比较。时域重建法是在数学上严格成立的热声成像算法，但由于计算公式过于复杂，实际应用中都采用其简化近似后的公式。滤波反投影法出现较早，虽然它在数学上是一种近似算法，但是可以完全正确地反映待测组织的轮廓信息[16]。
首先建立待测样品的模型，根据式(2)仿真出扫描圆周上接收到的声压(有160个扫描位置)，然后分别利用三种算法重建出待测样品的图像。图1显示了四种待测样品的原始图像和重建结果。图1中的四行从上到下对应了四种由小到大的待测样品，四列从左到右对应了样品原始图像、时域重建法结果、滤波反投影法结果和反卷积重建法结果。其中原始图像周围的白色圆圈表示超声换能器的扫描轨迹。
从图1可以看出，待测样品尺寸较小时，反卷积重建法有很高的精度(1d，2d)，基本与时域重建法、滤波反投影法相当。而待测样品尺寸接近超声换能器扫描半径时，反卷积重建法的效果较差(3d，4d)。时域重建法和滤波反投影法对各种尺寸的待测样品都有较高的成像精度。不过经过精确对比可以发现，滤波反投影法的重建图像在灰度上的偏差大于时域重建法。
进一步分析三种算法对各种尺寸的待测样品的成像精度。仿真从小到大的正方形待测样品，以重建图像的PSNR(峰值信噪比)为标准对三种算法进行比较。图2显示了r′/r0为不同值时三种算法重建图像的PSNR，其中r0是超声探头扫描半径，r′是待测样品离扫描圆心的最远距离，点划线对应时域重建法(TDR)，实线对应滤波反投影法(FBP)，虚线对应反卷积重建法(DR)。
由图2可见，r′/r0＜0.3时反卷积重建法的精度与时域重建法相当，略好于滤波反投影法；而r′/r0＞0.3时反卷积重建法的精度有明显下降。因此，反卷积重建法适合于对小尺度物体进行成像。
接下来比较三种算法的计算速度。仿真所用计算机为P4 2.0G，512Mb内存。表1显示了三种算法在不同成像分辨率下的计算耗时。表中数据表明，在仿真实验条件下，反卷积重建法的速度最快，约是时域重建法的4～6倍，滤波反投影法的25～100倍。
表1 三种算法在不同成像分辨率下的计算耗时 由实验结果可见，对小尺度物体进行成像时，反卷积重建法具有很高的精度，与时域重建法相当，略好于滤波反投影法。反卷积重建法的速度快于时域重建法，远快于滤波反投影法。因此，对于小尺度下的热声成像问题，反卷积重建法是一种快速有效的算法。
参考文献Song KH，Stoica G，Wang LV.In vivo three-dimensional photoacoustic tomography of awhole mouse head[J]Optics Letters，2006，31(16)2453-2455.Kruger RA，Reinecke DR，Kruger GA.Thermoacoustic computed tomography-technicalconsiderations[J].Medical Physics，1999，26(9)1832-1837.Niederhauser JJ，Jaeger M，Lemor R，et al.Combined ultrasound and optoacoustic system forreal-time high-contrast vascular imaging in vivo[J].IEEE Transactions on Medical Imaging，2005，24(4)436-440.Kruger RA，Liu P，Fang Y，et al.Photoacoustic ultrasound(PAUS)-reconstruction tomogra-phy[J].Medical Physics，1995，22(10)1605-1609.Xu M，Wang LV.Time-domain reconstruction for thermoacoustic tomography in a sphericalgeometry[J].IEEE Transactions on Medical Imaging，2002，21(7)814-822.Xu Y，Feng D，Wang LV.Exact frequency-domain reconstruction for thermoacoustic tomo-graphy-Iplanar geometry[J].IEEE Transactions on Medical Imaging，2002，21(7)823-828.Xu Y，Xu M，Wang LV.Exact frequency-domain reconstruction for thermoacoustic tomogra-phy-IIcylindrical geometry[J].IEEE Transactions on Medical Imaging，2002，21(7)829-833.Kstli KP，Beard PC.Two-dimensional photoacoustic imaging by use of Fourier-transformimage reconstruction and a detector with an anisotropic response[J].Applied Optics，2003，42(10)1899-1908.Yan Hua-gang，Liu Wu.A new technique of medical imaging-thermoacoustic imaging[J].Information of Medical Equipment，2002，736-40+43.Gu Huai-min，Yang Si-hua，Xiang Liang-zhong.Photoacoustic tomography and applicationsin the medical clinic diagnosis[J].Progress in Biochemistry and Biophysics，2006，33(5)431-437.Wang XD，Pang YJ，Ku G，et al.Three-dimensionallaser-induced photoacoustic tomogra-phy ofmouse brain with the skin and skull intact[J].Optics Letters，2003，28(19)1739-1741.Kolkman RGM，Hondebrink E，Steenbergen W，et al.In vivo photoacoustic imaging ofblood vessels using an extreme-narrow aperture sensor[J].IEEE Journal of Selected Topics inQuantun electronics，2003，9(2)343-346.Karabutov AA，Andreev VG，Bell B，et al.Optoacoustic images of early cancer in forwardand backward modes[C].InBoccara AC，Oraevsky AA.Proceedings of SPIE.2001，443413-27.Riad SM.The deconvolution probleman overview[C].InProceeding of the IEEE.1986，74(1)82-85.Bennia A，Riad SM.An optimization technique for iterative frequency-domain deconvolu-tion[J].IEEE Transactions on Instrumentation and Measurement，1990，39(2)358-362.Xu Y，Wang LV.Reconstructions in limited-view thermoacoustic tomography[J].MedicalPhysics，2004，31(4)724-733.
权利要求
1、一种小尺度下热声成像的反卷积重建算法，其特征在于根据在球面面上扫描的超声换能器探测到的声压函数p(r0，t)，构造出一个新函数B(r)；B(r)经三维傅里叶变换得到
基于反卷积方法由
计算出经三维反傅里叶变换得到A′(r)；A′(r)关于原点对称的函数就是电磁波吸收系数的分布A(r)；其中所述的B(r)为
其中常数β是等压膨胀系数，c是声速，Cp是比热，r0为球面半径；而p(r0，t)与A′(r)有如下关系
由反卷积方法计算
的公式为
其中λ是一个常数，
是h(r)经过三维傅里叶变换后得到的，h(r)为
全文摘要
本发明属于热声成像技术领域，具体为一种小尺度下热声成像的反卷积重建算法。该算法通过探测到的声压函数构造出一个新函数，再基于反卷积的方法重建出待测生物体内部的电磁波吸收系数的分布。与当前流行的热声成像算法--时域重建法和滤波反投影法相比，对于小尺度的待测生物体，反卷积重建算法的精度与时域重建法相当，好于滤波反投影法；重建速度快于时域重建法，远快于滤波反投影法。本发明对于小尺度生物组织热声成像问题是一种快速有效的算法。
文档编号G06F17/00GK101099681SQ20071004170
公开日2008年1月9日 申请日期2007年6月7日 优先权日2007年6月7日
发明者汪源源, 驰 张 申请人:复旦大学