完全三维图像重建的一个解析公式资料下载.pdf
- 文档编号:5967099
- 上传时间:2023-05-05
- 格式:PDF
- 页数:8
- 大小:426.65KB
完全三维图像重建的一个解析公式资料下载.pdf
《完全三维图像重建的一个解析公式资料下载.pdf》由会员分享,可在线阅读,更多相关《完全三维图像重建的一个解析公式资料下载.pdf(8页珍藏版)》请在冰点文库上搜索。
1984年Feldkamp2等人受22D扇束直接卷积算法的启发,利用微分求增量的方法得到近似重建公式,正如作者在文中指出他们的公式只是一种近似方法,理论上并不具有任何严格性。
1985年Smith3在理论上作了进一步的探讨,得到了三维图像重建的充分、必要条件,并推导出后来称之为Smith重建过程的算法。
1986年Natterer在TheMathematicsofComputer2izedTomography中从微分算子的观点仔细研究了重建公式。
最终Grangeat4在1987年得到了锥束投影的精确的重建公式。
我们注意到,在以后的研究过程中有关32D重建主要是沿着Smith2Grangeat的路子进行的,可参看文7、8、9。
本文受Grangeat的启发,给出了一个不同的具有几何直观的推导过程,这在理论和实用计算中都具有很重要的价值。
如图2,假设所探测的物体位于坐标原点附近,并且在探测锥之内。
如果用f(x,y,z)表示该物体的像素值(在SPECT中为放射药物的密度分布)。
这时在探头上观察到的值为:
g(,)=0f(a()+t)dt,S2,+。
这里(,)指沿方向从点a()出发的射线,+为扫描轨道参数变化的范围。
我们的目标就是如何利用g(,)(对充分多的(,)来确定函数f(x,y,z)。
如果记1为探头上的一条直线,令P为点a()与1组成的平面,假如能得到f(x,y,z)在此平面上的积分,记此平面的法方向为S2,距原点的距离为s,并令此积分值为Rf(,s),则Rf(,s)=X=sf(X)(X-s)dX,XR3,S2。
这里X=s表示方向为距原点距离为s的超平面,()表示Dirac函数,这时有如下的重建公式:
f(X)=182s2929s2Rf(,s)d。
这就是32DRadon逆变换公式。
对于锥束投影情形,Rf(,s)一般不能直接计算得到,但可以获得99sRf的值,因此离散化计算可以分成以下几个步骤:
11求g(,)关于特定方向的导数然后进行加权线积分从而计算出99sRf;
21利用差分方法计算929s2Rf(,s);
31球面积分2背投影求和:
s2929s2Rf(,X)d。
从上面的讨论中可以看出,计算步骤基本上与平行束类似,只是二维情形需进行卷积运算,详细的讨论可参看Herman11,而三维情形的局部性质避免了这种奇异积分型的卷积运算。
595第5期尤江生等:
完全三维图像重建的一个解析公式13-DRadon逆变换公式22DRadon逆变换公式是众所周知的,本节将简要推导32DRadon逆变换公式,并指出22D与32D之间的本质区别。
在讨论之前,将对密度函数f(x,y,z)作适当的假设,约定f(x,y,z)=0当x2+y2+z2R02,且f(x,y,z)为二次连续可微函数,当然这些假定仅仅出于某些技术上的原因。
记S2为R3中的单位球面,R3中的点表示为X=(x,y,z),如果S2,定义:
Rf(,s)=X=sf(X)d,其中d=(X-s)dX,()为一维Dirac函数。
在球坐标系中,可以写成:
=(coscos,cossin,sin),其中-2,2,0,2),由此可知:
Rf(-,s)=Rf(,-s)。
对固定的,Rf(,s)可以看成实轴R1上的函数,根据前面的假设Rf(,s)为s的二次连续可微函数,对函数Rf(,s)求关于变量s的一维Fourier变换:
Rf(,)=-Rf(,s)e-2isds=-dsX=sf(X)e-2isd=-dsX=sf(X)e-2iXd=R3f(X)e-2iXdX。
记f(x,y,z)为f(x,y,z)的三维Fourier变换,则f(x,y,z)=R3f(x,y,z)e-2i(x,y,z)(x,y,z)dxdydz。
显然有关系:
Rf(,)=f(),但是由于Rf(-,)=Rf(,-)。
因此总可以假设0,而(coscos,cossin,sin)可以看作(x,y,z)在Fourier空间中的极坐标表示。
由此由Fourier逆变换公式可以看出:
f(x,y,z)=R3f(x,y,z)e2i(x,y,z)(x,y,z)dxdydz=-20d?
2-?
2d0Rf()2cose2i(x,y,z)d=-0d?
2cosd-Rf()2e2ia(x,y,z)d=140d?
2cosd92Rf9s2(,s)s=(x,y,z)。
回顾22DRadon逆变换公式:
695北京大学学报第32卷f(r,)=1220d-dlrcos(-)-l99lP(l,),注意到公式中包含卷积形式的积分,也就是Hilbert变换,因此第二项积分可以认为是99lP(l,)沿直线1的某种形式的平均,它具有整体性质,从Smith3中可以看出,Smith的重建过程关键在于建立这种形式的平均与32DRadon逆变换公式之间的关系。
另一方面我们注意到32DRadon逆变换公式中仅仅出现项92Rf9s2,这可由函数Rf(,s)局部地唯一确定,因此具有局部性质,这就是22D与32D之间的本质区别,其实这是由维数的奇偶性引起的。
2锥束投影的重建公式由节1中知道,重建计算的关键在于面积分的计算,这也是三维重建的困难之处。
对于锥束投影,将探测数据限制在一平面上它可以看作极坐标系中函数沿径向的积分,但节1中所讨论的面积分是在笛卡尔坐标系中表示的。
因而如何从探测数据获得面积分就成为大家所讨论的焦点。
最早的FDK方法是利用微分求增量的办法,只能是一种近似的算法,不可能成为精确的重建公式。
另外,最近Wu在其博士论文中提出的重排计算也只是属于一种近似修正法,而且在实际计算中修正因子在很大程度上依赖于要重建的函数,这样一来修正因子的计算就带有很大的不确定性。
受Grangeat1的启发,将给出适用于任何轨道的重建计算过程。
图3扫描与探测Fig.3Scanneranddetector选取一固定的坐标系(O,x,y,z),假设点源运动的轨道为#,并且#位于球心在坐标原点半径为R0的球BR0以外,在实际设计中,对不同的点源S应对应不同的探测平面P。
不失一般性,假设探测平面P垂直于向量SO,当然也可以有其他形式的探测平面,但对从S出发的射线来看,线积分是相同的,因此只要建立探测平面与点源之间的固定位置关系即可,以后就在此坐标系中来讨论问题。
令D为直线SO与探测平面P的交点,在平面P内建立坐标系(D,u,),其中u,为P内795第5期尤江生等:
完全三维图像重建的一个解析公式两单位正交向量,并令=SO?
SO,于是由u,可以建立过D点的新的直角坐标系。
设AP,在A点探测到的值记为Xf(S,A),则:
Xf(S,A)=0f(S+tk)dt,其中k=SA?
SA。
直观地讲,对过S点的任何平面P1,在P1内建立以S为原点的极坐标系,则获得的数据正好是原密度函数关于径向变量的积分。
以平面极坐标系中的函数g(r,)为例,探测值相当于积分g(r,)dr,要从这样的积分获得g(r,)在整个平面上的积分,当然是不可能的,因为g(r,)平面积分的形式为g(r,)rdrd,积分号内要包含由坐标变换引起的加权因子r。
所以限制在一个平面上,要获得g(r,)rdrd仅仅知道g(r,)dr是不够的。
Wu在3中曾经讨论过一种重排计算方法,但作者并未指出如何重组g(r,)dr来获得g(r,)rdrd,特别对于多聚焦形式的探测,这一过程就更为复杂。
值得庆幸的是重建公式(3)f(x,y,z)=1420d2-2cosd92Rf9s2(,s)s=(x,y,z)告诉我们,虽然Rf(,s)不能直接计算得到,但只要能算出9Rf9s(,s)或92Rf9s2(,s)就可以完全重建出原密度函数f(x,y,z)。
这就是Grangeat最初的想法,以下就按照这种思路来计算9Rf9s(,s)。
首先建立如下的引理:
引理f(r,)为平面极坐标系中的函数,令()=0f(r,)dr,则有:
9()9=099pf(r,)rdr。
这里99pf(r,)指沿方向(-sin,cos)的方向导数。
接着前面的论述,如图4所示,对于探测平面P中的任何直线l,考虑过l与S的平面P1,令C为D在l上的投影,M为点D在直线SC上的投影,直线DA平行直线AC,且A为A在直线AD上的投影,线段QM平行直线AC,线段QB平行MC,B为点B在直线AD上的投影,从图中可以看出MD平行QA,DM为平面P1的法向量,平面SBA垂直于直线l。
当Al,易知Xf(S,A)=0f(S+tk)dt,k=SA?
计算Xf(S,A)沿方向n1=BABA的方向导99pXf(S,A)。
由极限关系易知:
99pXf(S,A)=lim0Xf(S,A+n1)-Xf(S,A)。
895北京大学学报第32卷图4Radon变换一阶导数与锥束投影的几何关系Fig.4RelationsbetweenthefirstderivativeofRadontransformandthecone2beamprojections假设为SA与SB之间的夹角,这时由前面的引理的与图3中各线段之间的关系容易看出99Xf(S,A)=099nf(S+tk)tdt,这里99nf(S+tk)表示f(S+tk)沿平面P1的法方向n的导数。
从理论上来讲,只要Xf(S,A)是已知的,就可以求出99Xf(S,A),但具体的解析表达式不容易推导。
我们的目标就是要建立99Xf(S,A)与99pXf(S,A)之间的关系,确切地说就是求9p9的表达式。
在图4中令SD=d,AC=u0,CD=0,ABB=0。
利用各种几何关系可以计算得到:
SA2=u02+02+d2,SB2=d2+(u0-0tg0)2,AB2=021+u0202(02+d2)2,在三角形SBA中,由余弦定理知:
AB2=SA2+SB2-2SASBcos,这时假如在点A沿方向BA有增量,则:
BA+2=(u0+sin0)2+(0+cos0)2+d2+SB2-2(u0+sin0)2+(0+cos0)2+d2)1?
2SBcos(+),取极限后可以得到:
9p9=SA2SBsinSAAB-(0cos0+u0sin0)(SA-SBcos)。
这样对任意固定的Al,有关系:
099nXf(S+tk)tdt=99pXf(S,A)9p9。
记为SA与SC之间的夹角,易知:
d=SCSA2dl。
995第5期尤江生等:
完全三维图像重建的一个解析公式99nRf(n,s)=-99pXf9p9SCSA2dl=-99pXfW(u0,0,d)dl,(u0,0)为点A在P1内的坐标,W(u0,0,d)可以看作是加权因子,W(u0,0,d)的具体表达式如下:
W(u0,0,d)=SCSBsinSAAB-(0cos0+u0sin0)(SA-SBcos)。
这个公式告诉我们,对所有过S的平面P我们都可以求出积分P1fd关于法方向的一阶导数,由此在理论上二阶导数也是可求的,具体计算的时候可用差分方法来求。
根据重建公式(3)要完全重建原函数,必须对所有的s以及角度知道99sRf(,s)的值,但对每个点源S,实际只能求得(n,OSn)点处的值。
这里顺便提一下,从这件事实不难看出Tuy的充分条件是显然的,同时Smith所论述的充分必要条件也是自然的。
下面谈一下对于不同的轨道所确定的精确重建区域,#为点源的扫描轨道,定义8=S#(n,OSn),进一步我们定义如下的集合:
Radon=XR3使得9B(XO)8这里9B(XO)表示以线段XO为直径的球面。
重建公式(3)告诉我们对所有Radon中的点都可以精确地算出f(X)的值,同时从Radon的定义中易知,对于平面轨道来说,只能精确计算出f(X)在轨道所在平面上的值,这就是说要想完全重建原函数,仅由平面轨道探测是不可能的,因而在以后的发展过程中人们设计螺旋轨道6、双环轨道7以及园+直线轨道10来解决这一问题。
3结论我们给出一个不同于Grangeat的方法来计算Radon变换的一阶导数,它的关键在于文中引理的应用。
与其相应的数值计算步骤可以按照引言中描述的程序进行。
进一步还论述了22D与32D之间的本质区别,以及由此引起的精确重建区域问题。
因此本文的结果可以看作是Grangeat工作的进一步发展。
这里给出的仅是理论公式,在以后的文章中将具体讨论离散化计算格式以及模拟实验。
参考文献1TuyHK.Aninversionformulaforcone2beamreconstruction.SIAMJApplMath,1983,43:
5465522FeldkampLA,DavisLC,KressJW.Practicalcone2beamalgorithm.JOptSocAm,1984,1:
6126193SmithBD.Imagereconstructionfromcone2beamprojections:
necessaryandsufficientconditionsandre2constructionmethods.IEEETransMedImag,1985,4:
14254GrangeatP.Mathematicalframeworkofcone2beam3DreconstructionviathefirstderivativeofRadontransform.LectureNotesinMathematics,1497:
66975NattererF.TheMathematicsofComputerizedTomography.NewYork:
Wiley,1986006北京大学学报第32卷6YanXH,LeahyRM.Derivationandanalysisofafilteredbackprojectionalgorithmforcone2beampro2jectiondata.IEEETransMedImag,1991,10:
4624727KudoH,SaitoT.Feasiblecone2beamscanningmethodforexactreconstructionin3Dtomography.JOptSocAm,1990,7:
216921838DefriseM,ClackR.Acone2beamreconstructionalgorithmusingshiftvariantfilteringandcone2beambackprojection.IEEETransMedImag,1994,13:
1861959KudoH,SaitoT.Derivationandimplementationofacone2beamreconstructionalgorithmfornonplanarorbits.IEEETransMedImag,1994,13:
19621110ZengGL,GullbergGT.Acone2beamtomographyalgorithmfororthogonalcircle2and2lineorbit.PhysMedBiol,1992,37:
56357711ImageReconstructionfromProjections:
Implementationandapplications,editedbyHermanGT.NewYork:
Springer2Verlag,1980AnAnalyticFormulaforFully3DImageReconstructionYOUJiangshengBAOShanglian(InstituteofHeavyIonPhysics,PekingUniversity,Beijing,100871)AbstractUsing3DRadoninversionformulaandalemmaofradialintegralinpolarcoordinatsystem,weobtainanewmethodtocomputethefirstderivativeofRadontransform,whichisdifferentfromthetreatmentofGrangeat.Thenewformulaisconvenienttodealwiththediscreteorbit.Keywordscone2beam;
Radoninversionformula;
Smith2Grangeatreconstructionprogram106第5期尤江生等:
完全三维图像重建的一个解析公式
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 完全 三维 图像 重建 一个 解析 公式
![提示](https://static.bingdoc.com/images/bang_tan.gif)