本发明主要涉及医学图像的重建。更具体地，本发明涉及新的图像重建方法和X射线计算机断层摄影设备，以利用混合型滤波提高重建图像质量和重建效率。

本发明包括各种技术的使用，在以下参考列表记载的文档中参考并描述，在整个说明书中由括号中相应的参考标号表示：参考文献[1]C.R.Crawford和K.E.King，Computerized tomography scanning withsimultaneous patient translation，Med.Phys.，17，967-982，1990年。
[2]P.E.Danielsson，P.Edholm，J.Eriksson和M.Seger，Towards exact3D-econstruction for helical cone-beam scanning of long objects。A new detectorarrangement and a new completeness condition，Proc.on Fully 3D ImageReconstruction in Radiology and Nuclear Med.，141-144，1997年。
[3]L.A.Feldkamp，L.C.Davis和J.W.Kress。Practical cone-beamalgorithm，J.Opt.Soc.Am，vol.1，pp.612-619，1984年。
[4]G.T.Herman和A.Naparstek，Fast image reconstruction based on aRadon inversion formula appropriate for rapidly collected data，SIAM J.Appl.Math.，vol.33，pp.511-533，1977年。
[5]A.Katsevich，Theoretically Exact FBP-Type Inversion Algorithm forSpiral CT，第六次无线电与核医学中的完整三维图像重建国际会议，PacificGrove，CA，USA，2001年10月-11月。
[6]A.Katsevich。An Improved Exact Filtered Backprojection Algorithm forSpiral Computed Tomography，Advances in Applied Mathematics，V.32-4，pp.625-825，2004年5月。
[7]A.Katsevich。Analyis of an Exact Inversion Algorithm for SpiralCone-Beam CT，Physics in Medicine and Biology，2002年，vol.47，pp.2583-2598。
[8]H.Kudo，F.Noo，M.Defrise和R.Clackdoyle，New super-short-scanreconstruction algorithm for fan-beam and cone-beam tomography，IEEENSS-MIC 2002年，M5-3。
[9]H.Kudo，F.Noo，M.Defrise和T.Rodet，New Approximate FilteredBackprojection Algorithm for Helical Cone-Beam CT with Redundant Data，IEEE NSS-MIC 2003年，M14-330。
[10]H.Kudo，T.Rodet，F.Noo和M.Defrise，Exact and approximatealgorithm for helical cone-beam ct with redundant data，Phys.Med.Biol.49(2004年)，2525-2546。
[11]A.C.Kak和M.Slaney，Principles of computerized tomographicimaging，IEEE Press，纽约，1988年。
[12]A.Katsevich和A.Zamyatin，Analysis of a family of exact inversionformulas for cone beam CT，已提交。
[13]A.V.Lakshminarayanan，Reconstruction from divergent ray data，Tech.Rep.92，Dept.of Computer Science，布法罗，SUNY，1975年。
[14]F.Natterer，The mathematics of computerized tomography，纽约：Wiley，1986年。
[15]F.Noo，M.Defrise，R.Clackdoyle，H.Kudo，Image reconstruction fromfan-beam projections on less then a short scan，Phys.Med.Biol.，47，2525-2546，2002年。
[16]D.L.Parker，Optimal short scan convolution reconstruction forfan-beam CT，Med.Phys.9(2)，1982年，254-257。
[17]J.D.Pack，F.Noo和H.Kudo，Investigation of saddle trajectories forcardiac CT imaging in cone-beam geometry，Phys.Med.Biol.，49，No 11，2317-2336，2004年。
[18]M.D.Silver，A method for including redundant data in computedtomography，Med.Phys.，27，773-774，2000年。
[19]M.D.Silver，K.Taguchi和K.S.Han，Field-of-view dependent helicalpitch in multi-slice CT Proc.of SPIE Med.Imag.Conf.4320，839-850，2001年。
[20]K.Taguchi，Temporal resolution and the evaluation of candidatealgorithms for four-dimensional CT，Med.Phys.，30(4)，640-650，2003年。
[21]K.Taguchi，B.S.Chiang和M.D.Silver，A new weighting scheme forcone-beam helical CT to reduce the image noise，Phys.Med.Biol.49，2351-2364，2004年。
[22]K.C.Tam，S.Samarasekera和F.Sauer，Exact cone-beam CT with aspiral scan，Phys.Med.Biol.，43，1015-1024，1998年。
[23]G.Wang，C.R.Crawford，W.A.Kalender，Multi-row-detector andcone-beam spiral/helical CT.IEEE Trans.Med.Imaging 19：817-821，2000年。
[24]G.Wang，T.H.Lin，P.Cheng和D.M.Shinozaki，A general cone-beamreconstruction algorithm，IEEE Trans.Med.Imaging，MI-12，486-96，1993年。
[25]A.Zamyatin，Analysis of Cone Beam Reconstruction in ComputerTomography，PhD Dissertation，中佛罗里达大学，2003年。
由计算机断层摄影(CT)装置创建的重建图像的质量和效率对该CT装置的总体性能来说是很重要的。重建图像时使用的算法影响质量和效率。
图1表示CT装置的一种可能的源轨迹的实例，其中可得到重建函数f(x)，在二维空间中x＝(x1，x2)，在三维空间中x＝(x1，x2，x3)。图1的x3轴用z轴表示。图1所示的螺旋轨迹由方程a(t)＝(R cos t，R sin t，tH/2π)给出，其中H称为螺距(helical pitch)。对于一个图像切片的重建，仅需要有限的投影范围。通常这种投影范围，表示为Λ，集中在由t0＝x3·2π/H限定的图像切片上。可根据s改写该轨迹方程(如a(s)＝R cos s，R sin s，sH/2π)，s∈Λ)其中s＝t-t0。有时投影角β(或λ)代替参数s使用。注意β＝s mod 2π。图1所示的L(s，x)表示从重建点x到焦点a(s)的水平距离。
图1的虚线圆圈表示视野(FOV)。FOV为包含重建对象或病人的圆形区域。通常FOV由其半径rFOV限定。FOV直径的典型值是500mm。感兴趣区(theregion of interest)(ROI)可以是整个FOV或FOV的一部分。
检测器几何形状和数据结构用g(s，Θ)表示在Θ方向上沿着来自源位置a(s)的射线的总衰减(或线积分)。该函数表示发生重建前通过CT扫描获得的数据。如图2所示，γ表示扇形角，α表示圆锥角。在扇形射束的情况下，Θ＝γ。在圆锥射束的情况下，Θ＝(γ，α)。也可参照该检测器表示Θ角。参照该检测器，ΘD＝(u，v)(或在扇形射束的情况下ΘD＝u)，因此ΘD取决于检测器几何形状。实际上，圆锥射束数据是检测器坐标g(s，ΘD)＝g(s，u，v)的函数，但一些理论公式将其作为扇形角和圆锥角g(s，Θ)＝g(s，γ，α)，即无检测器表达式(detector-free notation)的函数。这里使用无检测器表达式。
有许多具有不同几何形状的检测器类型：等角的、等距的、不等距的、平的、圆柱的、球的、倾斜的、旋转的、PI掩蔽的(PI-masked)等等。图2表示等角的圆柱形检测器，图3表示平的等距共线检测器。
卷积为清楚起见，下面使用算子符号。每个算子、卷积或反向投影大体上是独立于其自变量(argument)而定义的。算子用粗斜体表示，以区别于函数。相关算子定义如下。
傅立叶变换表示为G(s，ω)＝FT[g(s，γ)]，其中傅立叶变换应用于第二变量γ。
希耳伯特滤波算子(H)表示为H[g(s，γ)]＝g(s，γ)*h(γ)H[g(s，γ)]＝FT-1[G(s，ω)H(ω)]，其中h(γ)＝-1/(πγ)，和H(ω)＝i sinnω。
斜坡滤波算子(Q)表示为Q[g(s，γ)]＝g(s，γ)*q(γ)Q[g(s，γ)]＝FT-1[G(s，ω)Q(ω)]，其中q(γ)＝FT-1[Q(ω)]，和Q(ω)＝|ω|。
改进的希耳伯特滤波算子(Hm)(“改进的”表示核函数是h(sinγ)而不是h(γ))表示为Hm[g(s，γ)]＝g(s，γ)*h(sinγ)。
改进的斜坡滤波算子(Qm)表示为
Qm[g(s，γ)]＝g(s，γ)*q(sinγ)，其中q(sinγ)=(sinγγ)2q(γ),]]>以及Qm(ω)＝FT[q(sinγ)]。
最后，伴有直流漂移(Qm0)的改进的斜坡滤波算子表示为Qm0[g(s，γ)]＝FT-1[G(s，ω)Qm(ω)+G(s，0)/4π]。
斜坡滤波一般用于滤波的反向投影(FBP)算法。FBP是计算机断层摄影(CT)中主要的重建方法。与希耳伯特滤波相比，具有如下优点：(1)对于同一噪声级，具有更高的分辨率(更清晰的图象)，(2)为了在重建中平衡噪声，对分辨率进行更多的控制。第二个优点对于临床应用是非常重要的。在本发明中，主要贡献是从斜坡滤波的数据重建图像。斜坡滤波保持图像的所有高频分量。希耳伯特滤波为低频校正，使圆锥射束赝像(artifact)减少。
反向投影一般地，反向投影算子作用于数据g(s，Θ)，其中该算子将该数据映射到图像空间。由如下方程表示：BPJ[g(s,Θ)](x)=14π∫Λg(s,Θ)|Θ=Θ(s,x)ds]]>这里，首先在Θ＝Θ(s，x)计算g(s，Θ)，其中Θ(s，x)是来自横跨x的a(s)的射线。然后该方程沿着投影范围Λ上的源轨迹对所有这种射线进行积分。为了找到Θ(s，x)如下方程：γ(s,x)=arcsinx1sins-x2cossL(s,x)]]>和α(s,x)=arctanx3-za(s)L(s,x),]]>其中za(s)是a(s)的z坐标。反向投影通常应用于滤波的数据。
有两种方法加权该反向投影算子。该反向投影算子可由L的倒数或L2的倒数加权。这两种方法在如下两个方程中表示。
BPJL[g(s,Θ)](x)=14π∫Λ1L(s,x)g(s,Θ)|Θ=Θ(s,x)ds]]>BPJL2[g(s,Θ)](x)=14π∫ΛRL2(s,x)g(s,Θ)|Θ=Θ(s,x)ds]]>与1/L相比，1/L2的反向投影加权导致噪声恶化和图像上的点展开函数(BSF)的均匀性。
该反向投影范围由Λ表示。用BPJ算子符号中的上角标表示反向投影范围如下：BPJΛ[g(s,Θ)](x)=12π∫Λg(s,Θ)|Θ=Θ(s,x)ds]]>图4A-4E表示五个不同的重建范围。图4A表示全扫描，图4B表示过扫描，图4C表示短扫描，图4D表示超短扫描，图4E表示一种弹性范围。
冗余补偿加权可使用加权因子w(s，γ)修正数据中的冗余。对扇形射束数据不必在整个区域上扫描。对于扇形射束几何结构，g(s，γ)＝g(s+π+2γ，-γ)。对于扇形射束数据，整个扫描会导致对每个射线计数两次。因此，该数据需要由函数w(s，γ)＝1/2加权。当跨过FOV的任何线与重建段相交至少一次时，满足二维数据充分条件[14]。当重建范围Λ≥π+2arcsin(rFOV/R)时，满足二维数据充条件。对扇形射束数据来说，覆盖Λ＝π+2arcsin(rFOV/R)的投影集称为“最小完整数据集”[16]。如下术语在现有技术中已经使用并在这里使用：短扫描，或半扫描重建：Λ＝π+2Γ，arcsin(rFOV/R)≤Γ＜π/2；超短扫描重建：Λ＜π+2arcsin(rFOV/R)；全扫描重建：Λ＝2π；和过扫描重建Λ＞2π。
短扫描的Parker加权从关系g(s，γ)＝g(s+π+2γ，-γ)，其中-γm≤γ≤γm可以看到，对于精确扇形射束重建，仅π+2γm的重建范围已经足够，且对于区域π-2γ＜s＜π+2γm中的数据，数据g(s，γ)，0＜s＜2γm-2γ是冗余的。这里γm＝arcsin(rFOV/R)为该检测器允许的最大扇形角。参考[16]指出，在Parker加权中，应该以这样的方式加权最小完整数据区的数据，即尽可能对间隔均匀分配。参考[16]提出了如下加权函数：wP(s，γ)+wP(s+π+2γ，-γ)＝1wP(s,γ)=sin2(π4sγm-γ),0≤s≤2γm-2γ1,2γm-2γ≤s≤π-2γcos2(π4π+2γm-sγm+γ),π-2γ≤s≤π+2γm]]>
γm＝arcsin(rFOV/R)为该检测器允许的最大扇形角。如果没有定义，则该加权为零，例如，如果s＜0或s＞π+2γm，则wP(s，γ)＝0。在BPJ算子符号中用上角标P表示反向投影范围π+2γm：BPJP[g(s,Θ)](x)=12π∫-π/2-γmπ/2+γmg(s,Θ)|Θ=Θ(s,x)ds.]]>通用Parker加权(MHS加权)Parker加权是MHS加权的一种特殊情况。在MHS加权中，γm应至少为arcsin(rFOV/R)。通过实质上增大的γm，可获得更大的重建范围，从而获得更佳的噪声特性[18，19]。在wP(β，γ)的方程中，通过用虚拟最大扇形角Γ代替物理最大扇形角γm，得到如下方程：wMHS(β,γ)=sin2(π4βΓ-γ),0≤β≤2Γ-2γ1,2Γ-2γ≤β≤π-2γsin2(π4π+2Γ-βΓ+γ),π-2γ≤β≤π+2Γ]]>在BPJ算子符号中用上角标MHS表示反向投影范围π+2Γ：BPJMHS[g(s,Θ)](x)=12π∫-π/2-Γπ/2+Γg(s,Θ)|Θ=Θ(s,x)ds,]]>其中sin-1γFOVR≤Γ<π2.]]>过扫描加权由于g(s，γ)＝g(s+2π，γ)，过扫描加权用于反向投影范围Λ＝2πn+Δ，其中n＝1，2，......，0＜Δ＜π。加权函数如下。
wOS(β)=12nsin2(π2βΔ),0≤β<Δ1,Δ≤β≤2πncos2(π2β-2πnΔ),2πn<β≤2πn+Δ]]>在过扫描加权中，可使用MHS加权和Parker加权代替三角函数sin2或cos2、多项式3x2-2x3[1]或某些其它光滑函数。在BPJ算子符号中，用上角标OS表示反向投影范围2πn+Δ：BPJOS[g(s,Θ)](x)=12π∫-πn-Δ/2πn+Δ/2g(s,Θ)|Θ=Θ(s,x)ds.]]>Noo的加权Noo的加权的优点在于可使用任意重建范围Λ(s0，s1)，其中s0和s1为该重建段的起点和终点。这种加权具有小于半扫描的重建范围，可用于ROI重建。这种重建范围称为短扫描。该加权函数如下：WN(s,γ)=c(s)Σcompc(scomp,γcomp)]]>例如，在简单的短扫描情况下，Σcompc(scomp,γcomp)=c(s)+c(s+π+2γ).]]>在过扫描情况下：Σcompc(scomp,γcomp)=c(s)+c(s+π+2γ)+c(s+2π)]]>函数c(β)如下：c(β)=cos2π(s-s0-Δs)2Δs,s0≤s≤s0+Δs1,s0+Δs≤s≤s1-Δscos2π(s-s1+Δs)2Δs,s1-Δs≤s≤s1]]>其中Δs为平滑滤波间隔。Noo建议取Δs＝10°[15]。注意具有大Δs(≈50°)的Noo的加权等价于Parker加权。Noo的加权可使用任意反向投影范围Λ，如图4A-4E所示[15]。
准圆锥射束加权扇形射束加权可扩展至圆锥射束数据[21]。一旦计算出扇形射束加权w(β，γ)，作为圆锥角的函数，对其加权，并对其标准化，以得到准圆锥射束加权函数wQ3D(β，γ，α)。必须基于该数据的有效性(精确性)定义这种加权，有效性加权由wVal(α)＝1-α(3x2-2x3)表示，其中且α1＝tan-1[t1D/(2R)]；α2＝tan-1[t2D/(2R)]。
两个圆锥角(α1和α2)定义有效曲线的拐点。为了补偿冗余采样，有效加权wVal(α)与扇形射束加权函数w(β，γ)相结合，然后标准化。有效加权wVal(α)可以是任意值；但是，如果选择参数t1和t2，则wVal(α)是有意义的，从而wVal(α)将全加权分配给有效(测量的)射线和，将较少的加权分配给无效(未测量的)射线和，并对该转换滤波。因此，准圆锥射束加权函数为：
wQ3D(β,γ,α)=w(β,γ)wVal(α)/Σcompw(βcomp,γcomp)wval(αcomp),]]>其中在所有补偿位置求和，以使Σcompw(βcomp,γcomp)wval(αcomp)=1.]]>Tam窗口加权Tam窗口是由来自焦点a(s)的螺旋线的上部和下部投影限定的检测器的一部分。图5表示仅包含非冗余数据的完整的πTam窗口[22，2]。图5表示rH＝1.0的情况。虚线表示当rH＝1.25时物理检测器的边界。这里rH＝H/D，其中D为检测器宽度。πTam窗口的加权函数如下：3πTam窗口的加权函数如下：在BPJ算子符号中，用上角标π仅表示πTam窗口的反向投影，根据螺距，在BPJ算子符号中，用上角标TW表示πTam窗口或3πTam窗口的反向投影：BPJTW[g(s,Θ)](x)=12π∫wT(γ,v)g(s,Θ)|Θ=Θ(s,x)ds.]]>注意Parker和Noo的加权函数最初都是为扇形射束数据引入的和圆锥角的独立。这些加权函数的使用产生圆锥射束的赝像。Tam窗口加权为移位常量或呈投影角独立提供三维加权，即同一加权用于每个投影。
Tam窗口加权的优点是：(1)真实的圆锥射束加权，(2)移位常量，即加权函数wT(γ，v)对所有投影来说是相同的，独立于位置z，(3)执行简单。Tam窗口加权的缺点是：(1)未使用冗余数据(未使用测量的数据的某部分)，(2)Tam窗口是固定的(因此对应于π和3π，仅两个螺距是最佳的)，(3)利用不同的反向投影范围，重建不同的图像象素，产生较少的空间均匀图像。
在图6中，虚线表示当该源沿该螺旋线移动时，两个不同图像点沿该栓测器的路径。粗线表示在反向投影中使用的数据，以重建A点和B点。这样，如从图6所见，A点的重建比B点的重建使用更多的数据。
扩展的Tam窗口加权利用平滑滤波函数，Tam窗口可在z方向上扩展。该扩展的Tam窗口在图7中表示。扩展的Tam窗口加权可失去理论精确度，但具有更高的数据使用率。在BPJ算子符号中，用上角标ETW表示扩展的Tam窗口加权：BPJETW[g(s,Θ)](x)=12π∫wET(γ,v)g(s,Θ)|Θ=Θ(s,x)ds.]]>滤波线的方向图8A-8F表示不同的滤波线。起初，沿检测器行或段(图8A)滤波。对于螺旋轨迹，沿着平行于源位置点的螺旋线的切线的方向(切线)滤波，使圆锥射束赝像减少。由HTan和QTan表示的算子沿切线(图8B)滤波。在旋转滤波中(图8C)，中间滤波线与该螺旋线相切，最顶部和最底部的滤波线是水平的(平的)。它们之间的其它滤波线逐步旋转，以形成平滑滤波线族。由HRot和QRot表示的算子沿旋转线滤波。精确的算法[12]使用Katsevich的滤波线族。可在图8D-8F中看到Katsevich滤波线，沿Katsevich的曲线族滤波表示为HKat和QKat。
在实际操作中，如果检测器行不平行于滤波方向，则滤波需要从检测器点阵变换(rebinning)到滤波器点阵。每个变换包括内插，产生更光滑的图像并可能失去细节。HTan和QTan(及其它)算子的使用使圆锥射束赝像减少，其代价是重建时间增加且分辨率减小。
加权和卷积的顺序在操作中，执行加权和卷积的顺序是非常重要的。冗余加权w(s(x3)，γ)取决于重建切片z位置x3。如果卷积前必须加权，即滤波应用于加权数据子集，则对于每个切片，重建范围中的所有投影需要重加权和重卷积。另一方面，如果加权在卷积之后，则对于所有图像切片，该数据可以仅卷积一次，且每个切片仅需重加权。因此，在后一种情况下，重建所需的卷积数减少了很多。
重建算法斜坡滤波一般用于重建医学图像的FBP算法。扇形射束数据的原始FBP算法具有如下形式：[FBP]f(x)=14π∫-ππ1L2(s,x)Qm[g(s,γ)cos(γ)]γ=γ(s,x)ds,]]>其算子表达式的形式为：
f(x)=BPJL2[Qm[g(s,γ)cosγ]].]]>FBP算法首先在等距共线检测器的情况下发展起来[13]，后来在等角射线[11]的情况下得以扩展[4]。它使用1/L2加权的斜坡滤波。可以看到1/L2的反向投影加权与1/L相比，导致噪声恶化和整个图像上的PSF均匀性。该原始算法起初是为全扫描轨迹设计的，后来扩展至短扫描轨迹[16]：[FBP-P]f(x)=BPJL2P[Qm[wP(s,γ)g(s,γ)cosγ]],]]>再后来[3]是为圆形圆锥射束的几何结构设计的：[FDK]f(x)=BPJL2[Qm[g(s,Θ)cosΘ]],]]>(cosΘ＝cosγcosα)。
卷积前应用短扫描轨迹的加权函数[FBP-P]，即将滤波应用于加权的数据子集。由此，在该重建范围中所有投影的每个切片需要重加权并重卷积。卷积后进行加权的算法是更有效的。这是因为对所有数据图像切片来说，该数据仅卷积一次，然后重加权。因此，如果卷积后加权，重建所需的卷积数会减少很多。因此，卷积后应用加权函数的算法可提供显著的计算优点。
Feldkamp算法扩展至螺旋轨迹[24]，后来扩展至任意扫描段[18]如下：[GFDK]f(x)=BPJL2MHS,OS[Qm[wMHS,OS(s,γ)g(s,Θ)cosΘ]].]]>该算法的流程图如图13A所示。
上述讨论的算法限于圆形或螺旋源轨迹。不限于圆形或螺旋轨迹的算法是优选的。由于它可应用于其它轨迹如鞍形轨迹，因此可提供更通用的算法。
Katsevich[5，7]引入了FBP类型的精确圆锥射束算法。
[Katsevich]f(x)=BPJLπ[HmKat[(∂∂s+∂∂Θ)g(s,Θ)]].]]>它使用该圆锥射束数据的偏微分的改进的希耳伯特变换，而不使用一般的斜坡滤波。该算法的流程图如图13C所示。后来Katsevich将其方程概括为一个，不带/s项，并也发展了3π情况下的算法。在上述算法中，在滤波线的特殊族中必须进行滤波[12]，在卷积前后，需要附加的非一般的(non-trivial)变换步骤。图8D-8F表示Katsevich滤波线族。然而出于实际目的，通用的Feldkamp算法相当适用，且具有相对较少检测器行的现有的扫描仪不需要这种复杂性。但是基于希耳伯特滤波的重建提供了一些Feldkamp算法所不具备的非常好的特性[3]。Katsevich的算法与Feldkamp算法相比，显示出基于希耳伯特滤波的重建提供了一些非常好的特性。
后来引入基于希耳伯特变换重建的扇形射束重建算法[15]：[NDCK-FB]f(x)=BPJLΛ[wN(s,γ)Hm[(∂∂s+∂∂γ)g(s,γ)]].]]>参考[15]指出对于精确重建，仅需覆盖Λ＝π的投影。这使得小于短扫描的重建成为可能。对于短扫描重建，需要整个FOV上的投影集(π+2γm)。另外，该数据充分性条件放松。[NDCK-FB]算法的其它优点是卷积后执行冗余扇形射束数据的加权。由于每个切片的数据不必重卷积，这使得[NDCK-FB]算法与任何Feldkamp类型的算法相比更加有效。
在发明人估算噪声流体圆柱形模型的过程中，发现了Noo的算法的另一个优点。可以证明，与Feldkamp算法相比，噪声变化在整个图像上更加均匀。对于Noo的算法，PSF也具有较少的空间变化。这可以解释为反向投影加权是距离的倒数，而不是距离平方的倒数，且减小了所谓放大作用。
基于希耳伯特变换的算法与利用斜坡滤波[8]的算法相比，由于附加的数值微分步骤，引入了更强的平滑滤波。Kudo[8]提出的扇形和圆锥射束数据的算法包括斜坡和希耳伯特过滤两个滤波：[KNDC]f(x)=BPJL2Λ[wN(s,γ)Qm[g(s,γ)cosγ]+∂wN(s,γ)∂γHm[g(s,γ)cosγ]].]]>在参考[10]中，将[NDCK-FB]和[KNDC-FB]两算法概括为具有平的检测器的圆锥射束的情况如下：[KRND-1]f(x)=BPJLΛ[wN(s,γ)HmTan[(∂∂s+∂∂Θ)g(s,Θ)]]]]>和[KRND-2]f(x)=BPJL2Λ[wN(s,γ)QmTan[g(s,Θ)cosΘ]+∂wN(s,γ)∂γHmTan[g(s,Θ)cosΘ]].]]>Kudo的算法有几个缺点。首先，其平方倒数加权和仅对特殊加权函数起作用的缺点类似于Feldkamp的算法。其次，Kudo的算法包括对加权函数取偏微分，这使得该算法缺乏实用性。除了Katsevich的算法，上述算法对圆形轨迹上的扇形射束扫描是精确的，对圆锥射束扫描是近似的。
公知的重建算法有几个缺点。上述讨论算法中没有一个具有所有下列有益的方面：(1)1/L反向投影加权，(2)卷积后应用冗余补偿加权函数，(3)独立于所使用的加权函数的类型的算法，(4)独立于源轨迹的算法，以及(5)混合型过滤滤波的结合(斜坡型滤波和希耳伯特型滤波)。

因此，为了克服相关技术中重建算法的问题，本发明寻求提供一种方法、系统和计算机程序产品，以在扫描对象的计算机断层摄影重建点确定图像数据值。
因此提供了一种方法、系统和计算机程序产品，以在扫描对象的计算机断层摄影重建点确定图像数据值，包括：(1)获取扫描对象的投影数据；(2)用一维斜坡型滤波器滤波获取的投影数据，以生成斜坡滤波数据；(3)将具有距离倒数加权的反向投影算子作用于该斜坡滤波数据，以在CT图像的重建点生成图像数据值。
另外，按照本发明的一个实施例，上述方法还包括：将投影减法作用于获取的投影数据，以生成减法数据；对该减法数据应用希耳伯特型滤波器，以生成希耳伯特滤波数据；将投影加法作用于希耳伯特滤波数据和斜坡滤波数据，以生成滤波的数据；将冗余加权作用于滤波的数据，以生成加权数据，其中应用该反向投影算子的步骤作用于该加权数据，以在CT图像的重建点生成图像数据值。
按照本发明的另一方面，提供了一种X射线计算机断层摄影(CT)系统，以在重建点确定图像数据值，包括：(1)CT扫描单元，用于生成扫描对象的投影数据，该扫描单元包括生成X射线的X射线源和具有检测器元件的检测器，该检测器用于产生该投影数据；以及(2)处理器，包括：滤波单元，用于对该投影数据应用斜坡型滤波器，以生成斜坡滤波数据；以及反向投影单元，用于对该斜坡滤波数据应用具有距离倒数加权的反向投影算子，以在重建点生成图像数据值。
另外，按照本发明的一个实施例，上述系统还包括：投影减法单元，用于将投影减法作用于该投影数据，以生成减法数据；希耳伯特滤波单元，用于对该减法数据应用希耳伯特型滤波器，以生成希耳伯特滤波数据；投影加法单元，用于将投影加法作用于希耳伯特滤波数据和斜坡滤波数据，以生成滤波的数据；加权单元，用于将冗余加权作用于滤波的数据，以生成加权数据，其中该反向投影单元用于将具有距离倒数加权的反向投影算子作用于由该加权单元产生的加权数据。
通过下列附图和优选实施例的具体描述，本发明的其它方法、系统和计算机程序产品对本领域普通技术人员来说是显而易见的。

容易获得本发明更完整的说明及其许多附带的优点，结合附图参照以下具体描述，更容易理解这些内容，其中：图1表示具有螺旋几何结构的源轨迹；图2表示等角圆柱形检测器；图3表示共线检测器；图4A-4E分别表示全扫描、过扫描、短扫描、超短扫描和弹性重建范围；图5表示πTam窗口；图6表示Tam窗口重建范围；图7表示扩展的Tam窗口；图8A-8F表示滤波线；图9A表示按照本发明，在全扫描的重建点确定图像数据值的方法；图9B是执行图9A所示方法的系统；图10A表示按照本发明，在弹性投影范围、超短投影范围、短扫描或过扫描的重建点确定图像数据值的方法；图10B是执行图10A所示方法的系统；图11是描述可在本发明中使用的不同源轨迹、投影范围、加权函数以及滤波方向的表格；图12表示CT设备；以及图13A-13D表示四种重建算法。

本发明针对扇形射束数据的精确算法和圆锥射束数据的准精确算法。将这些算法用于由CT装置产生的图像重建中。下面给出本发明的实施例。
扇形射束数据的2π投影范围(全扫描)方程：
f(x)=14π∫-ππ1L(s,x)Qm0[g(s,γ)]γ=γ(s,x,)ds]]>f(x)=BPJL2π[Qm0[g(s,γ)]].]]>扇形射束数据的弹性投影范围，超短扫描，短扫描或过扫描方程：f(x)=12π∫Λw(s,γ)L(s,x)[Qm0[g(s,γ)]+Hm[∂∂sg(s,γ)]]γ=γ(s,x)ds]]>f(x)=BPJLΛ[w(s,γ)(Qm0[g(s,γ)]+Hm[∂∂sg(s,γ)])]]]>圆锥射束数据的2π投影范围(全扫描)方程：f(x)=14π∫-ππ1L(s,x)Qm0[g(s,γ,α)]γ=γ(s,x),α=α(s,x)ds]]>f(x)=BPJL2π[Qm0[g(s,Θ)]]]]>圆锥射束数据的弹性投影范围，超短扫描，短扫描或过扫描方程：f(x)=12π∫Λw(s,γ)L(s,x)[Qm0[g(s,γ,α)]+Hm[∂∂sg(s,γ,α)]]γ=γ(s,x),α=α(s,x)ds]]>f(x)=BPJLΛ[w(s,γ)(Qm0[g(s,Θ)]+Hm[∂∂sg(s,Θ)])]]]>本发明的上述实施例针对等角检测器给出。本发明的算法不仅限于等角检测器，还对其它检测器的几何结构起作用如共线检测器、等距检测器、非等距检测器、平面检测器、圆柱形检测器、球形检测器、倾斜的和PI掩蔽检测器。
另外，本发明的算法独立于源轨迹。源轨迹不需限于圆形或螺旋轨迹。例如可使用鞍形轨迹。
图9A表示全扫描的本发明的重建算法的一个实施例。在步骤900，通过公知方法获取CT投影数据。函数g(s，γ)或g(s，Θ)分别为图9A中步骤900表示的扇形射束数据或圆锥射束数据的投影数据。
在步骤901，该数据用斜坡型滤波器滤波，生成滤波的数据。该斜坡型滤波器可以是斜坡滤波器、改进的斜坡滤波器或者具有DC偏置改进的斜坡滤波器。
在步骤902，1/L加权的反向投影应用于该斜坡滤波的数据，以在CT图像中产生图像数据值。
在步骤903，确定在该CT图像中是否有需要重建的其它图像数据值(CT数)。如果有其它图像数据值，重复步骤902，直至在该CT图像中没有需要重建的其它图像数据值。
在步骤904，通过按照重建点设置该图像数据值，从步骤902输出的图像值用于生成CT图像。
图9B表示在CT图像中重建图像数据值的系统。CT扫描单元951产生投影数据。该CT扫描单元可向处理器960或向存储单元952传输该投影数据。处理器960用于直接从CT扫描单元951或从存储单元952存取接收该投影数据。处理器960包括滤波单元961，反向投影单元962和输出单元963。滤波单元961用于对该投影数据应用斜坡型滤波器，以生成斜坡滤波的数据。反向投影单元962用于使用反向投影算子对斜坡滤波的数据应用1/L加权，以在重建点上生成图像数据值。然后输出单元963向显示器970、存储单元971或两者输出图像数据值。斜坡滤波单元961和反向投影单元962均可存储并恢复来自存储单元952的信息。
图10A表示本发明的重建算法的另一个实施例。对弹性范围、超短扫描、短扫描或过扫描的扇形射束和圆锥射束方程使用g(s，γ)或g(s，Θ)分别表示图10A中步骤1000表示的扇形射束数据或圆锥射束数据的投影数据。
在步骤1001，投影减法应用于投影数据，生成减法数据。投影减法是在本发明的上述算法中对偏微分项的应用。
在步骤1002，希耳伯特型滤波器应用于该减法数据，以生成希耳伯特滤波数据。在步骤1002，可使用希耳伯特滤波器或改进的希耳伯特滤波器。
在步骤1003，斜坡型滤波器应用于投影数据，生成斜坡滤波数据。
在步骤1004，希耳伯特滤波数据和斜坡滤波数据相结合，生成滤波数据。
在步骤1005，冗余加权函数应用于该滤波数据。冗余加权函数w(s，γ)在所述重建算法中不是特定的，可在Parker加权、通用Parker加权(MHS)、Noo的加权、过扫描加权、准圆锥射束加权以及Tam窗口加权{wP(s，γ)，wMHS(s，γ)，wN(s，γ)，wOS(s，γ)，w3D(s，γ)，wT(γ，v)}中自由选择。
图11的表格包括使用哪种加权函数，用哪种类型的射束数据。图11表示本发明的算法对于扇形射束数据和圆锥射束数据，如何作用于不同的源轨迹、不同的检测器几何结构，可使用什么投影范围，可使用什么加权函数以及可使用什么滤波方向。投影范围Λ决定加权函数的选择。对扇形射束数据，超短扫描使用Noo加权；短扫描使用Parker加权、MHS加权或Noo加权；过扫描使用OS加权或Noo加权；弹性扫描使用Noo加权。除此之外当该投影范围是短扫描或过扫描时，可使用Q3D加权，圆锥射束数据与扇形射束数据相同。Tam窗口加权也可与圆锥射束数据一起使用。
在步骤1006，该加权数据通过距离的倒数加权作用于反向投影算子，以生成图像数据值。
在步骤1007，确定在该CT图像中是否有需要重建的其它图像数据值。如果有其它图像数据值，重复步骤1005-1007，直至在该CT图像中没有其它图像数据值。
在步骤1008，输出该图像数据值，通过按照重建点设置该图像数据值，生成CT图像。
图10B表示在CT图像中处理数据图像值的系统。CT扫描单元951获取扇形射束数据或圆锥射束数据。该CT扫描单元可向处理器960或向存储单元952传输该射束数据。处理器960直接从CT扫描单元951或从存储单元952存取接收该射束数据。在该实施例中，处理器960包括斜坡滤波装置961，投影减法单元964，希耳伯特滤波单元965，投影附加单元966，加权单元967，反向投影单元962以及输出单元963。投影减法单元964对来自CT扫描单元951或存储单元952的射束数据应用投影减法，以生成减法数据。然后希耳伯特滤波单元965对该减法数据应用希耳伯特型滤波器，以生成希耳伯特滤波数据。斜坡滤波单元961对来自CT扫描单元951或存储单元952的射束数据应用斜坡型滤波器，以生成斜坡滤波数据。投影附加单元966将该斜坡滤波数据与该希耳伯特滤波数据结合，以生成滤波数据。加权单元967对该滤波数据应用加权函数，以生成加权数据。反向投影单元962对该加权数据应用1/L加权的反向投影算子，以按照重建点生成数据图像值。输出单元963向显示器970、存储单元971或两者输出该数据图像值。希耳伯特滤波单元965，投影附加单元966，加权单元967，反向投影单元962以及输出单元963均可存储并恢复来自存储单元952的信息。
可以在扫描模式下应用这里描述的实施例，其中X射线源在圆形、螺旋形或鞍形轨迹上移动。该算法独立于该检测器的几何结构，即可使用等角的、等距的、非等距的、平的、圆柱的、球的、倾斜的、旋转的、PI掩蔽的等等。该方程独立于所使用的滤波线的类型；水平的，切向的、旋转的、Katsevich或任意其它合理的滤波线族。可使用超短扫描、短扫描、过扫描以及满足该充分条件的任意轨迹。
本发明与其它CT图像重建算法的比较重建算法可应用于不同的源轨迹，不同的检测器几何结构使用不同的滤波方向和不同的加权函数。但是，重建流和主要步骤对每个算法是唯一的。为了更清楚地看到该算法的主要特征，可使用算子表达式。下面对具有冗余加权螺旋圆锥射束几何结构的方程进行比较，该加权表示最实用的关注区。对等角检测器几何结构重写所有的方程。表示以下算法的流程图在图13A-D中表示，按照本发明的算法在图10A中表示。
[GFDK]f(x)=BPJL2MHS,OS[Qm[wMHS,OS(s,γ)g(s,Θ)cosΘ]]]]>[NDCK]f(x)=BPJLΛ[wN(s,γ)Hm[(∂∂s+∂∂Θ)g(s,Θ)]]]]>[KRND]f(x)=BPJL2Λ[wN(s,γ)Qm[g(s,Θ)cosΘ]+∂wN(s,γ)∂γHm[g(s,Θ)cosΘ]]]]>[Katsevich]f(x)=BPJLπ[HmKat[(∂∂s+∂∂Θ)g(s,Θ)]]]]>[本发明的一个实施例]f(x)=BPJLΛ[w(s,γ)(Q0m[g(s,Θ)]+Hm[∂∂sg(s,Θ)])]]]>空间均匀性取决于该反向投影加权。1/L反向投影加权的算法具有更佳的空间均匀性。因此，NDCK、Katsevich和本发明的算法产生比GFDK和KRND具有更佳的空间均匀性的图像。另一方面，在Katsevich的算法中，不同的图像象素具有不同的投影范围，导致更差的空间均匀性。图10和图13A-D的方框表示每个算法的反向投影加权。
由于本段的讨论不具有圆锥射束赝像，因此与[Katsevich]无关。所有其它算法仅对具有圆形轨迹的扇形射束几何结构是精确的(即在二维空间中)。通过扩展到三维空间，它们还可为各个相关小圆锥角和螺距提供相当良好的结果。但是，该算法在三维空间中不能同样良好地执行。圆锥射束赝像遮盖在物体的一个侧面上，并在另一个侧面上闪烁。[NDCK]和[KRND]伴有切向滤波，可减少但不会消除该圆锥射束赝像。由于切向滤波(以及旋转滤波)还可应用于[GFDK]和本发明的算法，因此不能视为仅[NDCK]和[KRND]的优点。由于[NDCK]和本发明的算法包含/s项，即相邻投影之间的微分，因此具有较少的圆锥射束赝像。这就补偿了螺旋数据的不一致性。注意圆锥射束赝像在Katsevich的算法中不出现，该算法也具有/s项。
圆锥射束赝像的主要原因之一是对螺旋圆锥射束数据应用扇形射束冗余加权。使用Tam窗口冗余加权函数使圆锥射束赝像显著减少。表1表示不同算法可使用哪些加权。
表1.不同算法的冗余加权注意仅[NDCK]和本发明使用所有加权。
当对CT重建算法相互比较时，不管是否有弹性重建范围和软件实施的简化，因子诸如体积重建的速度是很重要的。
体积重建速度主要由在切片重建环中执行多少操作确定。在图13A-D中，[GFDK]的切片重建环包含滤波，这意味着在每个图像切片中对相同的投影重卷积多次。所有其它算法更加有效：每个投影仅卷积一次(对于[NKDC]和Katsevich)或两次(对于[KRND]和图10A中本发明的一个实施例)；不需要重卷积。但是，[KRND]切片重建环比[NKDC]、[Katsevich]和本发明更复杂。注意反向投影是反向投影中最需要计算的部分，每个切片周期的仅一个反向投影是商用CT重建设备所强烈期望的。
弹性重建范围意味着源轨迹的任意子集可用于精确重建，该轨迹在xy平面上的投影满足二维数据充分条件。在结构上[NKDC]、[KRND]和本发明具有弹性重建范围。[Katsevich]由于使用Tam窗口加权，因此没有弹性重建范围。弹性重建范围还意味着超短扫描的可能性。具有这种可能性的算法为：[NKDC]、[KRND]和本发明。
实施的简化对商用CT重建算法是关键的。一个重要的准则是，每个步骤(滤波、加权、反向投影)可表示为单独的模块。另外，[NKDC]中的角度微分取决于它如何实施，并会变得相当复杂。在本发明中，仅有的微分/s为简单的投影减法。
数字数据微分导致分辨率降低。在[GFDK]、[KRND]和本发明中，由于图像的主要部分从斜坡滤波数据重建，因此没有分辨率损失，且不经过微分步骤。
表2表示在此情况下所述算法的主要特征：
表2性能比较图12表示X射线计算机断层摄影装置，可用于获得由本发明的方法所处理的数据。由台架1组成的投影数据测量系统包括X射线源3，用于生成近似圆锥形的X射线流的圆锥射束，还包括二维阵列型X射线检测器5，包括以二维方式排列的多个检测器元件，即多个以一维多行堆叠方式排列的元件。X射线源3和二维阵列型X射线检测器5安装在旋转环2上，位于床6的滑动板上的主体(subject)的相对面上。二锥阵列型X射线检测器5安装在旋转环2上。每个检测器元件与一个通道相对应。来自X射线源3的X射线通过X射线滤波器4指向主体。通过二维阵列型X射线检测器5，对已经通过该主体的X射线作为电信号进行检测。
X射线控制器8向高电压发生器7提供触发信号。高电压发生器7按照接收触发信号的时序向X射线源3典型地为X射线管提供高电压。这导致X射线从X射线源3发出。台架/床控制器9同步地控制台架1的旋转环2的旋转和床6的滑动板的滑动。系统控制器10组成整个系统的控制中心，控制X射线控制器8和台架/床控制器9，如从该主体所见，使X射线源3执行所谓螺旋扫描，其中它沿着螺旋路线移动。特别地，旋转环2以固定角速度连续旋转，同时该滑动板以固定速度移动，X射线以固定间隔角连续或间歇地从X射线源3发出。
二维阵列型X射线检测器5的输出信号由每个通道的数据收集单元11放大并转化为数字信号，以生成投影数据。从数据收集单元11输出的投影数据提供给重建处理单元12。重建处理单元12利用该投影数据找到反向投影数据，该反向投影数据在每个体素中反射该X射线的吸收。如第一实施例所述，在使用X射线的圆锥射束的螺旋扫描系统中，成像区域(有效视野)呈绕旋转轴中心半径为ω圆柱形。重建处理单元12在该成像区域中定义了多个体素(三维象素)，并为每个体素找到反向投影数据。利用该反向投影数据编码的三维图像数据或断层摄影数据发送到显示装置14，在那里在视觉上显示为三维图像或断层摄影图像。该扫描装置包括单元1-11。
为了这种描述的目的，我们应该定义一种图像，以表示物理场景，其中该图像已经由某种成像技术生成。成像技术的实例可包括电视机或CCD照相机或X射线，声纳或超声波成像装置。记录图像的初始媒介可以是电子固态装置，光学记录胶片或某种其它装置，例如可光激发的荧光。然后记录的图像可通过电子结合(如在CCD信号的情况下)或机械/光手段(如在对光学记录胶片进行数字化或对来自可光激发的荧光的数据进行数字化的情况下)转化为数字形式。
利用一般的通用计算机或按照本发明的教导的可编程的微处理器可方便地实施本发明的所有实施例，这对计算机领域的普通技术人员是显而易见的。合适的软件可以由基于本发明教导的普通编程人员容易地制定，这对软件领域的普通技术人员是显而易见的。特别地，计算机外壳可容纳主板，包含CPU、存储器(如DRAM、ROM、EPROM、EEPROM、SRAM、SDRAM以及FlashRAM)及其它可选的特定目的的逻辑器件(如ASICS)或可配置逻辑器件(如GAL和可再编程FPGA)。该计算机还包括多个输入装置，(如键盘和鼠标)，以及控制显示器的显示卡。另外，该计算机可包括软盘驱动器；其它可移动的媒介装置(如小型盘、磁带和可移动的磁光媒介)；还包括硬盘或其它固定的高密度的媒介驱动器，并利用合适的装置总线连接(如SCSI总线、增强IDE总线或超DMA总线)。该计算机还可包括小型盘阅读器，小型盘读/写单元或小型盘光盘机，它们可连接于同一装置总线或连接于另一个装置总线。例如，典型地，在计算机系统上执行的计算机程序包括获得该主体投影数据的装置，生成斜坡滤波数据的装置，用于通过一维斜坡滤波获得的投影数据，还包括对所述斜坡滤波数据应用反向投影处理的装置，其中用对应于从焦点到重建点的距离的倒数的加权处理对该反向投影进行处理，以在CT图像中的重建点生成图像数据值。
与本发明相关的计算机可读媒介包括小型盘、硬盘、软盘、磁带、磁光盘、PROMs(如EPROM、EEPROM、Flash EPROM)、DRAM、SRAM、SDRAM等。存储在这些计算机可读媒介的任何一个或结合体上，本发明包括控制计算机硬件和使计算机与用户交互的软件。这种软件可包括但不限于装置驱动器、操作系统和用户应用程序，如开发工具。本发明的计算机程序产品包括存储计算机程序指令的任何计算机可读媒介(如计算机代码装置)，当由计算机执行时，引起该计算机执行本发明的方法。本发明的计算机代码装置可以是任意可注释或可执行的代码机构，包括但不限于描述语言、注释、动态链接库、Java类以及完整的可执行程序。此外，为了更好的性能、稳定性和/或更低的成本，可对本发明的处理部分进行分配(如(1)多重CPUs或(2)至少一个CPU和至少一个可配置逻辑器件)。例如，可在第一计算机上选择轮廓或图像，并发送至第二计算机进行远程诊断。
本发明还可通过应用特定集成电路的准备工作或通过相互连接合适的常规元件电路网络得以实施，这对本领域普通技术人员来说是显而易见的。
按照本发明的图像数据源可以是任何合适的图像获取装置，如X射线机器或CT设备。另外，如果获取的数据未转换成数字形式，则可将其数字化。可选地，正在获取和处理的图像数据源可以是存储由图像获取装置产生的数据的存储器，该存储器可以是本地或远程的，其中可使用数据通信网络，如PACS(图片存档计算机系统)，以存取按照本发明处理的图像数据。
根据上述教导，本发明的许多修改和变化都是可能的。因此可以理解，在所附的权利要求书的范围内，可以对本发明在特定描述之外实施。