专利名称：投影与背投影方法与相应的实现算法的制作方法
技术领域：
本发明一般涉及重投影-背投影过程，更具体地，涉及重投影-背投影技术/算法，该技术/算法包括新的插值与数据存取方案，能够导致比现有技术更快的速度、较少的非正常现象，较低的噪声，以及较高的空间分辨率。
背景技术：
在计算X线断层照相术(computed tomography)中，将N维图像转换为N维组线积分集的操作被称为前向投影(forward projection)或重投影(reprojection)。这种操作的最明显的例子是生成物体的X光图像的物理过程。在对数转换之后，能够很好地将X光图像近似为物体的线性衰减系数的分布的线积分投影。实际上，为了X线断层照相术模拟或当进行递归重构时，需要前向投影仪。
相反的操作被称为背投影(backprojection)。这用于过滤背投影与递归重构之中，形成了当今重构算法的主体。
存在许多重投影与背投影方法。在一种方法中，每条X光线以线表示，并且每条线与每个像素的交叉长度被用作加权因子。另外的技术对X光线交叉的每一行或列，在两个像素之间进行线性插值(参看图1)。后两种方法为光线驱动方法。
在投影的情况下，循环通过所有投影线，并且对每一投影线都通过图像加权并对图像像素值求和，以近似光线积分。背投影被定义为相反操作加权因子保持相同，但检测器值被加权并分配给图像像素。
另一种技术为像素驱动方法，一般用于过滤背投影(参看图2)。循环通过所有图像像素，并对每一图像像素，划一条线连接源与该图像像素。然后确定该线与检测器阵列的交叉。在最靠近该交叉的两个检测器值之间进行线性插值，并将结果分配给该图像像素。重投影被定义为相反操作。左右检测器格(bin)的权重由下式给定
ωl=dr-ddr-dl]]>等式(1)ωr=d-dldr-dl]]>其中d为交叉的位置，dr与dl为交叉的右侧与左侧的第一个检测器格中心。
存在有其他方法，诸如基于球型基本函数(spherical basic function)的方法以及使用最近邻居或不插值的方法。
重投影与背投影操作计算量繁重，除了诸如在CT等等中所使用的模拟与重构技术的必要的部分。大部分现有方法可以细分为光线驱动与像素驱动方法。光线驱动与像素驱动方法两者都有的缺点在于它们引入了非正常信号，第一种(即光线驱动方法)在背投影时，第二种(即像素驱动方法)在重投影时。二者都有的另一缺点在于在每一视图重投影/背投影中所使用的数据的百分比。
例如，在具有比检测器格尺寸小很多的像素的图像的光线驱动投影的情况下，只有一部分像素对该角度的投影有作用。对像素驱动背投影的相反情况也如此。在递归重构时，其中要求重投影与背投影方法，可以考虑光线驱动重投影与像素驱动背投影的结合，以避开前面的问题。但是，即使当这是可能的时，人们常常更喜欢用匹配的重投影器-背投影器对。实际上，选择重投影器-背投影器方法的一个重要标准是速度。
两个限制速度的主要因素是计算复杂度与数据存取时间。对于光线驱动方法，计算相对简单。因此对于小数据量，光线驱动方法大大快于像素驱动方法。然而，在较大数据量上，数据存取时间变得更重要，而在这一阶段像素驱动方法开始受益于它的顺序图像存取时间，而光线驱动方法几乎是随机地存取数据。对于3D锥型光束情况，数据集合变得更大，因此数据存取时间变得更重要。
对于这些技术以及与其关联使用的各种装置的进一步公开，可以参照于1998年12月8日以Kawai等人名义公布的美国专利5848114，2002年2月26日以Besson名义公布的美国专利6351514，2002年1月15日以Besson名义公布的美国专利6339632。此处包括了这些专利作为参考。
发明内容
更具体地，本发明的第一方面在于一种图像处理方法，包括将像素网格中的像素投影到具有多个格的检测器上，或反之亦然；动态调整像素与检测器格之一的方形窗口的尺寸，使相邻窗口形成该检测器的检测器格和该像素之一之上的连续阴影；以及确定每一像素在该检测器每一检测器格上的效果，或反之亦然。
本发明的第二方面在于一种图像处理方法，包括将像素网格的每一像素的边沿，以该像素网格中预定的像素线性序列，投影到通过该网格的预定的线上，该网格由从源到检测器所投影的光线交叉；将检测器的每一格的边沿投影到该预定的线上；根据像素边沿与检测器格边沿在预定线上的投影，确定每一像素对该检测器阵列的格的作用，或反之亦然。
本发明的第三方面在于一种图像处理方法，包括确立包含图像像素的像素网格，这些像素以图像行与列排列；连续映射图像像素与检测器的检测器格之间的各个迁移，该检测器已检测到来自源的射线，包括将检测器格迁移投影到预定线上；将像素迁移映射到该预定线上；以及基于相邻投影之间的距离，将该检测器格与像素之一用该预定线上段长度加权。
本发明的第四方面在于一种计算机可读介质，载有计算机可执行的程序，该程序用来处理图像，所述程序配置为指挥计算机将像素网格中的像素投影到具有多个格的检测器上，或反之亦然；动态调整像素与检测器格之一的方形窗口的尺寸，使相邻窗口形成该检测器的检测器格与该像素之一之上的连续阴影；确定每一像素在该检测器每一检测器格上的效果，或反之亦然。
本发明的第五方面在于一种计算机可读介质，载有计算机可执行的程序，该程序用来处理图像，所述程序配置为指挥计算机将像素网格的每一像素的边沿，以该像素网格中预定的像素线性序列，投影到通过该网格的预定的线上，该网格由从源到检测器所投影的光线交叉；将检测器的每一格的边沿投影到该预定的线上；根据像素边沿与检测器格边沿在预定线上的投影，确定每一像素对该检测器阵列的格的作用，或反之亦然。
本发明的第六方面在于一种计算机可读介质，载有计算机可执行的程序，该程序用来处理图像，所述程序配置为指挥计算机确立包含图像像素的像素网格，这些像素以图像行与列排列；连续映射图像像素与检测器的检测器格之间的各个迁移，该检测器已检测到来自源的射线，包括将检测器格迁移投影到预定线上；将像素迁移映射到该预定线上；以及基于相邻投影之间的距离，将该检测器格与像素之一用该预定线上段长度加权。
图1为具有线性插值的光线驱动重投影-背投影的示意图，其中对于投影线所交叉的每一行或列，在两个相邻像素之间进行线性插值；图2为具有线性插值的像素驱动重投影-背投影的示意图，其中连接源与图像像素的线确定与探测器阵列的交叉，并且其中在两个相邻探测器格之间进行线性插值；图3示出均匀视图的光线驱动背投影，显示的结果中引入了高频非正常信号，因为有些像素比其邻近像素更新得更快；图4表示均匀盘(disk)的像素驱动投影图，其中引入了高频非正常信号，因为有些检测器格比其邻近格更新得更快；图5为像素驱动线性插值方法的示意图，其中由于所投影的方型窗口的不规则重叠，某些检测器格将比其他格有更多作用，从而导致高频震荡；图6示出像素驱动线性插值方法，其中调整方形窗口的宽度，使其总是相邻；图7示出距离驱动重投影器-背投影器，其中检测器格界面与像素界面都映射到x轴，并且其中导致的段长度在投影与背投影中被用作加权因子；图8示出距离驱动重投影器-背投影器，提供了像素界面pi与检测器界面di的交织模式的进一步的视图；图9表示均匀盘的距离驱动投影，其中完全消除了高频非正常信号；图10示出均匀视图的距离驱动背投影，其中完全消除了高频非正常信号；图11示出对SUN E4500计算机的每一背投影的时间对数据大小的图。
具体实施例方式
为更好地理解本发明的实施例，需要更详细地解释一下现有技术。在图1、2、6、7中，网格表示像素图像重构网格，其固定于三维坐标系统之中，根据响应于从源到检测器(两者都被大致示出)所投射的光线而取得的数据，将像素映射到该三维坐标系统之上。这些网格中的每一个方块表示一个像素。
如上所述，光线驱动与像素驱动方法都有的一个缺点是它们都引入了高频非正常信号，一种在背投影中，另一种在重投影中。图3示出一均匀视图的光线驱动背投影。干扰模式的原因在于有些像素比其他像素更新得更快。当像素大小与检测器格大小相比是小的时候，非正常信号问题将更严重，而在像素大小与检测器格大小相比是大的时候，这一问题将消失。
图4示出均匀盘的像素驱动投影的窦腔X线照相(sinogram)线。作为示例，在计算X线断层照相术中，所测量的数据集合(窦腔X线照相)由大量的视图(投影)构成。每一视图相应于整个检测器阵列的测量，所以每一视图包括大量的检测器格(投影线)。一幅典型的窦腔X线照相包括1000检测器格/投影线的1500个视图/投影。
如上所述，干扰模式的原因在于有些检测器格比其他格更新得更频繁。另外，当检测器格大小与像素大小相比是小的时候，非正常信号问题将更明显，而在检测器格大小与像素大小相比是大的时候，这一问题将消失。在本例中，只是作为示例，进行了重投影与背投影，采用平面二维扇形光束，放大率1.76，256×256像素，256个探测器格，360°上256个视图，和126°的任意开始角度。
这两种方法都有的另一缺点在于每一视图投影/背投影中数据的使用情况。为解释方便，假设对具有比探测器格大小大很多的像素的图像的光线驱动投影(参看图5)。只有部分像素对该角度的投影有作用。类似地，在对具有比探测器格大小小很多的像素的图像的像素驱动背投影中，在每一视图中只使用了一部分探测器值。这导致了低劣的噪声性能。在递归重构中，这也可能导致差的收敛特性。
在选择投影器-背投影器方法时，一个非常重要的标准是计算速度。两个限制计算速度的主要因素是计算复杂度与数据存取时间。对于光线驱动方法，计算相对简单。因此对于小数据量，光线驱动方法快于像素驱动。然而，在较大数据量上，数据存取时间变得更重要。在这些条件下，像素驱动方法由于其降低存取时间的固有的顺序图像存取而开始显示较好的处理速度，而光线驱动方法需要高得多的随机存取，因为它跳过大数据块并由此脱离数据所存储的顺序方式。结果是处理延迟。
然而，对于3D锥型光束情况，数据集合变得更大，并且这些效果变得更加重要。
a)像素驱动与光线驱动投影器-背投影器适应图5与6分别示出展示现有像素驱动方法的缺点的特征，以及本发明的实施例，其中适应或改变像素驱动方法以防止高频非正常信号。
更具体地，定位与检测器阵列交叉。在交叉处，假定一具有与像素值相等的区域的迪拉克脉冲。将其与宽度等于检测器格大小矩形窗口进行卷积。通过在两个相邻检测器格上积分结果，取得权重。结果用于权重的下列公式ωl=dm-(d-(dr-dl)/2)dr-dl]]>ωr=(d+(dr-dl)/2)-dmdr-dl]]>等式(2)dm=dl+dr2,]]>其中dm为位于dl与dr之间中心的界面的位置。这与等式1相同，其显示了该表达式的对等式。希望通过投影像素的一个均匀行，来达到在与该行相应的投影范围(除由于交叉的位置变化引起通路长度的稍微变化)上的基本均匀的投影。然而，由于所投影的方型窗口的不规则重叠，某些检测器格将比其他格有更多作用，从而导致高频震荡。
这通过以下解决根据本发明的适应的光线驱动的实施例，通过调整方形窗口或像素阴影的宽度，使其总是相邻，并消除缝隙，使它们实际上变得连续。这由图6中的灰色阴影区域示出，并且可以表示为ωl=max(min(dm,d+W/2)-(d-W/2)W,0)]]>ωr＝1-ωl等式(3)W＝Δp·M·cosαd/Δd，其中W为方形窗口的新的宽度，Δp为像素大小，Δd为检测器格大小，M为放大率，αd为投影线的角度。cosαd如果可用cosαdm近似，则可以预先计算。然而窗口宽度W不能大于检测器格大小dr-dl，因为那样其可能覆盖多于两个检测器格。
当然可以通用算法，以允许覆盖多个检测器格，例如通过使用while-loop。然而，这却带来了这一问题减少非正常信号的优点不能抵消算法复杂度的增加。
在对像素驱动技术的适应中，动态调整应用于像素而不是格。
更具体地，对在光线驱动背投影中所引入的非正常信号有类似理由。这导致下列权重用于校正的算法ωl=max(min(pm,p+W/2)-(p-W/2)W,0)]]>ωr＝1-ωl等式(4)W＝Δd/M/cosαp/Δp，其中p交叉位置，pl与pr为该交叉左右的第一个像素中心。然而在本例中，窗口宽度W不能大于像素大小pr-pl，因为那样其可能覆盖多于两个像素。
假定这些适应的算法的速度与原来算法相当。两种适应的方法都完全消除了原来算法所导致的图3与4所示的非正常信号。
b)距离驱动投影-背投影在此实施例中，本发明基于检测器阵列向图像行或列的连续映射，或反之亦然，并且更具体地基于沿投影线方向的映射。为了快速计算，所有检测器位置与图像位置都投影到任意选取的线上，该线可能是(例如)图像的x或y轴。
有了这一点，图像数据为顺序存取，类似于像素驱动方法计算简单，并且类似于光线驱动方法不引入非正常信号，并且在每一视图中都均匀地使用了所有数据。新算法在硬件与软件中都可以实现，其简单，速度块，完全使用了数据而降低了噪声，并且没有引入非正常信号。
更具体地，图7示出了该技术的实施例，其基于检测器阵列向图像行(或列)的连续映射，或反之亦然，并且更具体地基于沿投影线方向的映射。如上所述，为了快速计算，x轴(y轴)用做像素与检测器格的相对位置的基准。为定义图像像素与检测器格的连续映射，不同于对中心操作，使用了像素之间与检测器格之间的迁移。首先，将所有检测器格迁移投影到x轴(或y轴或任意确定的轴)。接着循环通过所有图像行(或列)，并将像素迁移投影到该轴。从图像中读取值，具有投影之间所定义的适当的段长度权重，并根据情况分配给检测器格或像素。
图8示出像素界面pi、检测器界面di、检测器值dij、以及像素值pij的交织模式的更详细的视图。在本实施例中，当前行对光线总和dij的作用可写为d23＝p12d34＝p12等式(5)d45=(p2-d4)·p12+(d5-p2)·p23d5-d4]]>
而对背投影有p12=((d2-p1)·d12+(d3-d2)·d23+(d4-d3)·d34+(p2-d4)·d34)p2-p1]]>p23=(d5-p2)·d45+(d6-d5)·d56+(p3-d6)·d67p3-p2]]>等式(6)图9表示均匀盘的距离驱动投影，相当于图4中的像素驱动投影的结果。将会看到，使用本技术，如同所修改的像素驱动投影器与光线驱动投影器一样，完全消除了高频震荡。
图10示出距离驱动投影相当于图3中的光线驱动背投影的结果。再一次地，使用该方法，如同像素驱动投影器与所适应的光线驱动投影器一样，完全消除了高频非正常信号。
为了比较性能，因为投影与背投影的计算时间非常类似，所以集中于背投影。图像与窦腔X线照相均选择为n×n像素。
图11示出在对于SUN E4500(10 UltraSPARC-II，400Mhz，8Mb快速存储器，10GB RAM)在使用三种方法时，每一背投影所需的时间对比数据大小的图。对于小数据量，算术过程成为瓶颈，因为所有的数据可以放入高速缓存。此处像素驱动方法明显较差，而距离驱动方法接近于光线驱动方法。对所有三种算法都进行了同样的优化。对于大数据集合，存储器存取时间变得更加重要，因为现在整个图像不再能放入高速缓存。只有光线驱动方法真正受此影响，因为存储器存取不是顺序的。这解释了光线驱动方法曲线的斜率。对于更大的数据集合，像素驱动与距离驱动方法具有更大优势，即它们可以在硬件中实现。光线驱动则不能，因为硬件背投影器一般不能同时访问所有的存储器。
上面公开的距离驱动投影-背投影方法汇总如下。然而，为更好地理解本技术的实质，将首先大致列出未修改的像素驱动与光线驱动技术。
像素驱动技术处理所有的图像像素(*)，并且对每一图像像素，执行以下步骤确定连接源与图像像素中心的线。
找到该线与检测器阵列的交叉。
确定检测器格的中心最接近该交叉的两个检测器格。
对于背投影通过这两个检测器格之间的线性插值，计算在该交叉上的值，并将该值分配给图像像素。
对于(重)投影使用与背投影中相同的权重，将图像像素的值分配给这两个检测器格。
光线驱动技术处理所有投影线(**)(所有视图中)通过连接源与检测器格的中心定义投影线。
对每一投影线，执行以下步骤对于(重)投影重置投影总和。
处理所有图像行(***)，并对每一图像行(***)，执行以下步骤计算该投影线与该图像行(***)(的中心线)的交叉。
确定此行(***)中其中心最接近该交叉的两个图像像素。
对于(重)投影通过这两个图像像素之间的线性插值，计算该交叉上的值，并将该值加到投影合。
对于背投影使用与(重)投影中相同的权重，将检测器格的值加到这两个图像像素上。
对于(重)投影将投影总和分配到检测器格上。
距离驱动技术处理所有视图，并且对每一视图，执行以下步骤对于每一检测器格确定该检测器格的边沿确定通过连接该检测器格边沿与X光源的线。
计算该线与x轴(***)的交叉。
该交叉定义所投影的检测器格边沿。
处理所有图像行，并对每一图像行执行以下步骤处理该行中的所有像素，并对每一像素，执行以下步骤确定该图像像素的左右(***)边沿确定通过连接该像素边沿与x光源的线。
计算此线与x轴(***)的交叉。
该交叉定义所投影的像素边沿。
将所投影的检测器格边沿与所投影的像素边沿排序列表。
从x轴的最左侧开始，并且确定当前像素与当前检测器格。
直到达到最右边沿，执行以下步骤确定哪个是下一个边沿(****)。
更新当前像素或当前检测器格。
计算加权因子，即当前边沿的位置减去前一边沿的位置。
对于(重)投影将当前图像像素的值乘以该加权因子，并将其加到当前当前检测器格。
对于背投影将当前检测器格的值乘以该加权因子，并将其加到当前图像像素。
标记(*)表示/属于“像素驱动”(**)表示/属于“光线驱动”(***)如果投影线的方向是更为水平，则以下转换是必须的′行′<-->′列′′x轴′<-->′y轴′′左′<-->′底′′右′<-->′顶′(****)表示/属于“距离驱动”功能应该注意，对所公开技术的总结是示意性的，不能用来具体地限定本发明的范围，并且虽然上述公开只集中于几种投影与背投影方法，但这些技术的应用并不局限于CT应用。还应该注意到在给定限定性假设之下，修改现有的光线驱动与像素驱动线性插值方法消除了高频非正常信号。然而，距离驱动方法完全消除了高频非正常信号，而没有任何限定性假设，在每一视图中，所有数据均匀地作用于结果投影或背投影，并且该方法还具有有利的计算特性。
另外，虽然讨论了用于二维平面检测器扇形光束CT几何方法，但应该理解这些方法与结论并不局限于此，并且本领域或紧密相关领域的技术人员应该理解，这些概念适用于其他二维与三维(或更多)几何关系，包括(仅作为示例)PET与SPECT几何关系。
权利要求
1.一种图像处理方法，包括将像素网格中的像素投影到具有多个格的检测器上，或反之亦然；动态调整像素与检测器格之一的方形窗口的尺寸，使相邻窗口形成该检测器的检测器格和该像素之一之上的连续阴影；以及确定每一像素在该检测器每一检测器格上的效果或反之亦然。
2.如权利要求1所述的图像处理方法，其中在像素驱动图像形成技术中，像素的方形窗口的宽度W的动态调整用以下等式确定ωl=max(min(dm,d+W/2)-(d-W/2)W,0)]]>ωr＝1-ωlW＝Δp·M·cosαd/Δd，其中W为方形窗口的新的宽度，Δp为像素大小，Δd为检测器格的大小，M为放大率，αd为投影线的角度。
3.如权利要求1所述的图像处理方法，其中用于光线驱动图像形成技术中的像素阴影的动态调整使用以下等式确定ωl=max(min(pm,p+W/2)-(p-W/2)W,0)]]>ωr＝1-ωlW＝Δd/M/cosαp/Δp，其中p为光线在检测器上的交叉位置，pl与pr为该交叉左右的第一个像素中心。
4.一种图像处理方法，包括将像素网格的每一像素的边沿，以该像素网格中预定的像素线性序列，投影到通过该网格的预定的线上，该网格由从源到检测器所投影的光线相交；将检测器的每一格的边沿投影到该预定的线上；根据像素边沿与检测器格边沿在预定线上的投影，确定每一像素对该检测器阵列的格的作用或反之亦然。
5.一种图像处理方法，包括确立包含图像像素的像素网格，这些像素以图像行与列排列；连续映射图像像素与检测器的检测器格之间的各个迁移，该检测器已检测到来自源的射线，包括将检测器格迁移投影到预定线上；将像素迁移投影到该预定线上；以及基于相邻投影之间的距离，将该检测器格与像素之一用该预定线上段长度加权。
6.一种计算机可读介质，载有计算机可执行的程序，该程序用来处理图像，所述程序配置为指挥计算机将像素网格中的像素投影到具有多个格的检测器上，或反之亦然；动态调整像素与检测器格之一的方形窗口的尺寸，使相邻窗口形成该检测器的检测器格与该像素之一之上的连续阴影；确定每一像素在该检测器每一检测器格上的效果，或反之亦然。
7.如权利要求6所述的计算机可读介质，其中在像素驱动图像形成技术中，像素的方形窗口的宽度W的动态调整使用以下等式确定ωl=max(min(dm,d+W/2)-(d-W/2)W,0)]]>ωr＝1-ωlW＝Δp·M·cosαd/Δd，其中W为方形窗口的新的宽度，Δp为像素大小，Δd为检测器格大小，M为放大率，αd为投影线的角度。
8.如权利要求6所述的计算机可读介质，其中对于光线驱动图像形成技术中，像素阴影的动态调整使用以下等式确定ωl=max(min(pm,p+W/2)-(p-W/2)W,0)]]>ωr＝1-ωlW＝Δd/M/cosαp/Δp，其中p为光线在检测器上的交叉位置，pl与pr为该交叉左右的第一个像素中心。
9.一种计算机可读介质，载有计算机可执行的程序，该程序用来处理图像，所述程序配置为指挥计算机将像素网格的每一像素的边沿，以该像素网格中预定的像素线性序列，投影到通过该网格的预定的线上，该网格由从源到检测器所投影的光线相交；将检测器的每一格的边沿投影到该预定的线上；根据像素边沿与检测器格边沿在预定线上的投影，确定每一像素对该检测器阵列的格的作用，或反之亦然。
10.一种计算机可读介质，载有计算机可执行的程序，该程序用来处理图像，所述程序配置为指挥计算机确立包含图像像素的像素网格，这些像素以图像行与列排列；连续映射图像像素与检测器的检测器格之间的各个迁移，该检测器已检测到来自源的射线，包括将检测器格迁移投影到预定线上；将像素迁移映射到该预定线上；以及基于相邻投影之间的距离，将该检测器格与像素之一用该预定线上段长度加权。
全文摘要
公开了相对于像素/检测器格(bin)的投影与背投影光线的方法，以减轻/消除高频非正常信号(artifact)。前两个方法分别为像素驱动与光线驱动线性插值技术的适应。在这些技术中，每一像素/格的窗口或阴影被动态调整，并被投影到检测器格/像素上，以消除阴影之间的缝隙。这使在给定检测器格上的每一像素(或反之亦然)的效果能够被适当的加权。第三种方法为距离驱动技术，其中像素与检测器格的迁移被分别投影到共同的轴上。这允许用较短的计算时间与改善了正常的图像来确定每一像素/检测器格对每一检测器格/像素的作用。
文档编号A61B6/03GK1524248SQ03800609
公开日2004年8月25日 申请日期2003年4月15日 优先权日2002年4月15日
发明者布鲁诺·K·B·德曼, 萨米特·K·巴苏, K 巴苏, 布鲁诺 K B 德曼 申请人:通用电气公司