CN1284122C - 使用一个或多个微处理器的快速并行锥形线束重建系统和方法 - Google Patents

使用一个或多个微处理器的快速并行锥形线束重建系统和方法 Download PDF

Info

Publication number
CN1284122C
CN1284122C CNB028050894A CN02805089A CN1284122C CN 1284122 C CN1284122 C CN 1284122C CN B028050894 A CNB028050894 A CN B028050894A CN 02805089 A CN02805089 A CN 02805089A CN 1284122 C CN1284122 C CN 1284122C
Authority
CN
China
Prior art keywords
processing unit
point processing
computing machine
data
cone beam
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Lifetime
Application number
CNB028050894A
Other languages
English (en)
Other versions
CN1491404A (zh
Inventor
宁若罗
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Rochester
Original Assignee
University of Rochester
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Rochester filed Critical University of Rochester
Publication of CN1491404A publication Critical patent/CN1491404A/zh
Application granted granted Critical
Publication of CN1284122C publication Critical patent/CN1284122C/zh
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/40—Arrangements for generating radiation specially adapted for radiation diagnosis
    • A61B6/4064—Arrangements for generating radiation specially adapted for radiation diagnosis specially adapted for producing a particular type of beam
    • A61B6/4085—Cone-beams
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/58—Testing, adjusting or calibrating thereof
    • A61B6/582—Calibration
    • A61B6/583—Calibration using calibration phantoms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/40—Imaging
    • G01N2223/419—Imaging computed tomograph
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/428—Real-time
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y10—TECHNICAL SUBJECTS COVERED BY FORMER USPC
    • Y10S—TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y10S378/00—X-ray or gamma ray systems or devices
    • Y10S378/901—Computer tomography program or processor

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Radiology & Medical Imaging (AREA)
  • General Physics & Mathematics (AREA)
  • Pathology (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Biomedical Technology (AREA)
  • Optics & Photonics (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • Animal Behavior & Ethology (AREA)
  • Surgery (AREA)
  • Pulmonology (AREA)
  • Molecular Biology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Mathematical Optimization (AREA)
  • Immunology (AREA)
  • Biochemistry (AREA)
  • Analytical Chemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Pure & Applied Mathematics (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

在具有一个或多个微处理器的计算机上于实践中可接受的时间内执行锥形线束重建。涉及重建的计算被分成将在每一个微处理器的MMX、ALU和SSE单元中执行的计算。对于纯粹的浮点数据，优选情况下，使用MMX单元来调整数据地址和图数据，使用SSE单元执行反投影。数据被z行分区，以便在反投影的每一阶段要被处理的数据容纳在L1高速缓存内。

Description

使用一个或多个微处理器的快速 并行锥形线束重建系统和方法
有关政府利益的声明
本发明在某种程度上受NIH Grants 2R01HL48603-05和1R41HL59703的支持。政府对本发明享有某些权利。
技术领域
本发明涉及医学成像等等中的锥形线束重建的系统和方法，具体来说，涉及在一个或多个微处理器上实现的这样的系统和方法。本发明还可用于无损探伤、单光子发射体层扫描和基于CT的爆炸物检测、微CT或微锥形线束体积CT等等。
背景技术
锥形线束重建吸引了医学成像界的广泛注意。锥形线束重建的示例可以在共同授予的美国专利5,999,587和6,075,836和美国专利申请No.09/589,115和09/640,713中找到，它们的内容被完全引入本说明书中作为参考。
CT(计算机X射线断层造影术)图像重建算法可以分为两个主要类别：滤波反投影(FBP)和迭代重建(IR)。人们对过滤反投影的讨论比较多，因为它精确，并能够快速实现。过滤反投影可以作为准确的重建方法或者作为近似的重建方法来实现，两者都基于拉冬(Radon)变换和/或傅里叶变换。
锥形线束重建过程需要花费大量的时间，并需要许多计算操作。当前，锥形线束重建过程的时间非常长，几乎难以用于临床及其他实际应用场合。假设有一组投影大小N＝512的数据，由于FBP的时间和计算是O(N4)，因此重建需要GELOPS(每秒十亿次浮点运算)级的计算能力。通常，使用改善的算法和更快的计算引擎可以实现快速的锥形线束重建。
现有的用于重建的快速算法是基于傅里叶切片定理或反投影的多分辨率重新采样的。基于傅里叶切片定理的算法使用内插法将傅里叶投影数据从极坐标转换为笛卡儿坐标，从中可以通过反向FFT获得重建。人们为缩短FBP时间做了大量的工作，大多数工作都集中于扇形线束数据上。其中包括莱诺图(linogram)方法和“链接”方法以及用于重新投影的相关的快速方法。有人基于正弦图(sinogram)和“链接”提出了近似法；这样的方法对于二维FBP起作用，并会造成O(N2LogN)复杂性。“链接”方法被扩展到了三维锥形线束FBP；在对每一行中的投影数据进行数据重分箱(rebinning)之后，可以对经数据重分箱的数据应用与二维中的相同方法，对于锥形线束重建，数据处理时间可以缩短到O(N3logN)复杂性。还有人提出了另一种快速算法，对于二维FBP使用快速层次反投影(FHBP)算法，其克服了现有的快速算法的一些缺点。FHBP算法基于拉冬变换的层次分解，并且对于重建需要O(N2logN)计算复杂性。令人遗憾的是，实验证明，对于合理的图像大小，N≈103，相对于比较简单的FBP实现的性能比潜在的N/logN加速少得多。和费尔德坎普(Feldkamp)算法相比，重建质量也会降低。在现实的实施方式中，总的重建时间不仅取决于计算复杂性，而且还取决于循环单位时间。如上所述的使用链接方法的三维锥形线束FBP需要更多的存储器空间来存储“链接”区域。包含内插系数和用于访问“链接”数据的地址信息的链接重建表区域需要O(N3)附加的存储器，并由于存储器存取时间而降低性能。加速小于N/logN。
具有各种并行性和流水线的定制反投影硬件引擎可以将执行速度限制到非常小的程度。硬件可以是基于FPGA的模块或者ASIC模块，定制掩码可编程门阵列，基于单元的IC和场可编程逻辑器件或者带有高速RISC或DSP处理器的插件板。这些板通常使用高速多端口缓冲存储器或者DMA控制器来提高各板之间的数据交换速度。还有一些技术，如矢量计算和预内插投影数据，都和定制引擎一起使用，以便减少重建操作。大多数定制的硬件是为二维FBP重建应用构建的。市场上还不能买到专门为快速锥形线束重建设计的基于单个或多个微处理器的重建引擎。
可以使用多处理器计算机或者多计算机系统来加速锥形线束重建算法。许多大型的并行计算机具有通过高速数据路径互连的紧密耦合的处理器。该多处理器计算机可以是共享存储器计算机或分布式存储器计算机。对大型的和非常昂贵的并行计算机进行了大量的工作。大多数这样的工作都使用基于三维拉冬变换的算法。作为示例，在诸如Cray-3D、Paragon和SPI之类的大型计算机上实现了费尔德坎普算法和两个迭代算法：3D ART和SIRT。在这样的实施方式中，局部数据分区用于费尔德坎普算法和SIRT算法，而全局数据分区用于ART算法。该实施方式是像素驱动的。处理器之间的通信速度对于重建时间非常重要，费尔德坎普实施方式可以在多指令多数据(MIMD)计算机中获得最佳性能。在Intel Paragon和CM5计算机上实现了并行二维FBP。使用定制的加速硬件或大型的并行计算机不是一个经济合算的快速重建解决方案，对于研究工作来说它不便于修改或添加新算法。
在分布式计算环境中，许多计算机可以连接在一起以作为多计算机系统工作。计算任务被分布到每一台计算机。通常，在多计算机系统上运行的并行程序使用诸如MPI(消息传递接口)或PVM(并行虚拟机)之类的一些标准库。在一组与以太网相连接的SunSparc2计算机上测试了并行重建，该实施方式是基于PVM库进行的。费尔德坎普算法是在基于MPI库的非均匀工作站组群上实现的。该实施方式在六个计算机组群上运行，结果显示，负载平衡中的实施方式导致81.8％的处理器利用率，使用异步通信将处理器利用率提高到91.9％。多计算机组群的最大缺点是，通信速度降低了重建速度。由于锥形线束重建涉及大量的数据存储，所以数据通常分布到每一台计算机。计算机需要在反投影阶段交换数据。存储器通信是影响重建速度的比较大的因素。另一个缺点是，不能用多计算机群集获得小型的重建引擎。还有人试图采用诸如COBRA(公共对象请求代理体系结构和规范)之类的分布式计算技术实现锥形线束重建。通常，分布式计算库比直接使用MPI库花费更多的通信时间，因此导致重建速度降低。
除各处理器之间的并行性之外，单处理器用某些微型体系结构技术也可以获得数据和操作并行性。指令级并行性(ILP)是一个处理器和编译器设计技术的家族，可以通过使单机操作并行地执行来提高执行速度。现代的处理器可以将指令执行划分为几个阶段；诸如流水线和分支预测之类的某些技术允许多个指令同时执行。为实现数据处理并行性，一些处理器添加了单指令多数据(SIMD)指令，这可以在一个指令中处理多个数据。这样的处理器包括Intel的带有MMXTM和SSE/SSE2的IA 32体系结构，Motorola的带有AltVecTM的PowerPCTM和带有3DnowTM的AMD Athlon。然而，到目前为止，这样的并行性没有在锥形线束重建中得到利用。
发明内容
综上所述，显而易见，需要一种以在实践中可以接受的速度执行锥形线束重建而不需要定制硬件或大型计算机的技术。因此，本发明的目的是提供一种可以用便宜的广泛存在的设备用于锥形线束重建的系统和方法。
本发明提供一种用于生成表示一个对象的内部的三维图像的系统，该系统包括：射线扫描器，该扫描器通过将射线穿过该对象传递到探测器上来生成投影信号；以及接收投影信号的计算机，用于通过对投影信号执行多种计算来生成三维图像，该计算机包括至少一个定点处理单元和至少一个浮点处理单元，该至少一个定点处理单元与该至少一个浮点处理单元并行操作，该计算机将该多种计算分为将在该至少一个定点处理单元中执行的第一组计算和将在该至少一个浮点处理单元中执行的第二组计算。
本发明提供一种用于生成表示一个对象的内部的三维图像的方法，该方法包括：(a)将线束穿过该对象传递到检测器上来生成投影信号；以及(b)接收投影信号并通过对锥形线束投影信号执行多种计算来生成三维图像；其中，步骤(b)是在计算机上执行的，该计算机包括至少一个定点处理单元和至少一个浮点处理单元，该至少一个定点处理单元与该至少一个浮点处理单元并行操作，该计算机将该多种计算分为将在该至少一个定点处理单元中执行的第一组计算和将在该至少一个浮点处理单元中执行的第二组计算。
为实现上述及其他目的，本发明提供在基于混合型计算(HC)的市场上可买到的PC上的高速CBR的实际实施方式。用多级加速实现费尔德坎普CBR，利用单指令多数据(SIMD)执行HC，并使处理器中的执行单元(EU)有效地工作。可以利用操作系统中的多线程和纤维支持(fiber support)，它们能在多处理器环境中自动实现重建并行性，还能使至硬盘的数据I/O更有效。存储器和缓存访问通过适当的数据划分而得到优化。经过在Intel Pentium III 500Mhz计算机上进行测试和与传统的实施方式进行比较，本发明对于288投影(每一个都具有5122个数据点)可以缩短过滤时间75％以上，对于5123立方体节省60％以上的重建时间，同时还能保持良好的精度，平均误差小于0.08％。所产生的系统是经济合算的，而且速度很快。有效的重建引擎可以用市场上可买到的对称式多处理器(SMP)计算机构建，这种计算机可以随着更新的PC处理器和较高存取速度的存储器的问世很方便地升级，而且费用也比较便宜。
在本发明中，费尔德坎普算法锥形线束重建(FACBR)可以实现高速度，高精度。测试环境是Intel Pentium III 500Mhz，640MB100Mhz存储器。结果显示，对于带有288投影的5123立方体的重建可以在少于20分钟的时间内完成，并保持良好的精度，而老的实施方式需要100分钟以上。使用了多个模拟幻影测试HCFACBR的精度。将重建的图像与模拟的幻像和用传统方法重建的图像进行比较，结果显示，用本发明的方法重建的图像与传统方法图像比较，平均误差少于0.04％，并且比计算机模拟的幻影有较好的精度。可以快速而准确地重建三维对象的线性衰减系数分布。
较高速度的启用SSE-2的Pentium IV和2或4处理器PC预计将来能在几分钟内完成5123立方体FACBR。FACBR是用多级加速和利用SIMD和ILP技术的混合型计算来实现的。存储器和超高速缓存访问通过适当的数据分区而得到优化。与采用大型计算机和计算机群集的实施方式相比，本发明经济合算，并且速度快。市场上可以买到的SMP计算机提供了有效的重建引擎，随着更新的PC处理器升级起来也很方便和便宜。相比之下，定制的硬件代价大，也很难升级。
以下将说明PC上的FACER的高速实施方式，还将说明用于混合型执行(HE)和混合型数据(HD)的技术。利用这些混合型计算功能，可以在通用PC上实现优良的存储器组织与指令优化、高速的费尔德坎普实施方式，并且可获得较高的性价比。HD和HE也可以应用于其他硬件平台上的实施方式以便提高FACBR性能。利用较高时钟频率的处理器和便宜的市场上可以买到的SMP PC，可以获得用昂贵而不方便的定制硬件获得的良好性能。由于使用市场上可以买到的PC来实现高性能，因此，为锥形线束重建设计新的算法和新的系统是方便的，将图像抓取系统和三维呈现系统集成到一个容易配置和升级的单一系统中也是有用处的。
由于Intel x86CPU频率已经提高到GHz级，所以构建基于x86的多处理器高速锥形线束重建计算引擎是切合实际的，也是经济可行的。虽然费尔德坎普算法是一种近似锥形线束重建算法，但是它是一个切实可行的并且有效的三维重建算法，是包括本发明的几种精确的锥形线束重建算法中的基本组成部分。
本发明在单微处理器或多处理器上实现并行处理。使用混合型计算(定点和浮点计算)加速了锥形线束重建，而不会降低重建的精度，也不会提高图像噪声。这些特征对于重建软组织(例如癌检测)特别重要。
附图说明
下面将参考附图详细阐述本发明的优选实施例，其中：
图1显示了在优选实施例中的重建中使用的锥形线束坐标系；
图2显示了Intel x86处理器的体系结构；
图3显示了巳知的FAGBR实施方式的UML图表；
图4显示了根据优选实施例的混合型执行的UML图表；
图5显示了在优选实施例中使用的数据分区方案；以及
图6显示了可以在其上实现本发明的优选实施例的系统的方框图。
具体实施方式
现在将参考附图详细阐述本发明的优选实施例。
将参考图1所示的坐标系说明优选实施例。O-XYZ是世界坐标系。X-Y-Z轴给出了用于重建的立体像素(voxel)的物理坐标。Z轴是旋转轴。t-s轴是旋转的台架X-Y坐标系。S轴始终穿过x射线源并垂直于探测器平面。
多个研究组已经调查了锥形线束几何形状的重建。使用中的最有效的算法是由Feldkamp，L.A.，Davis，L.C.，and Kress，J.W.，″Practical Cone-Beam Algorithm，″J Opt.Soc.Am.A 14：(6)612-619(1984)，and Kak，A.C.，and Slaney，M.，Principles of ComputerizedTomographic Imaging，IEEE Press，1988开发的一种。在此算法中，投影数据被反投影到图像缓冲区上，每一个被检测的射线都被在其方向上反投影。图像缓冲区中的像素增量是投影像素的量。投影必须在反投影之前过滤。
本讨论中使用的坐标如图1所示，并且是相对于x射线源102、探测器104和台架106定义的。假设对象P(θ)的投影以角度θ被探测器坐标u和v索引。重建的立体像素值被物理坐标x、y和z索引。旋转中心是z轴。从x射线焦点到旋转轴的距离是dso。通过缩放投影像素大小，探测器的垂直轴可以移动到z轴。通过随后缩放几何形状，z轴上的像素大小可以成为一体。这些缩放简化了重建算法的计算。费尔德坎普算法属于过滤反投影算法。费尔德坎普算法的实施方式包含下列步骤：
a)向投影数据应用加权和倾斜过滤(ramp filter)；这是通过向每一个P(θ)值应用加权并向带有过滤数据的行或列中的数据应用卷积来进行的。
P θ filter ( i , j ) = ( d so d so 2 + i 2 + j 2 P θ ( i , j ) ) * h ( i ) - - - ( 1 )
b)将数据Pθ filter(u，v)反投影到重建的立体像素f(x，y，z)：
f ( x , y , z ) = ∫ 0 2 π u 2 P θ filter ( u · t , u · z ) dθ
u = d so t d so - s - - - ( 2 )
t＝x·cosθ+y·sinθ
s＝y·cosθ-x·sinθ
如上所述，(t，s)是台架系统中的坐标，这是X-Y轴进行角度为θ的旋转转换。
假设在x、y和z方向上的重建体积是Nx×Ny×Nz立体像素。对于M投影，假设dθ是连续的角度位置之间的角度差，Pθ(i，i)是Ni ×Nj像素。问题的复杂性是大致与立体像素的总数成线性关系并与投影视角的数量成线性关系。立体像素的数量加倍，也会使处理时间大致加倍。每一个方向上的尺寸加倍会导致所需要的处理八倍地增大。视图数量加倍也会使处理加倍。随着立体像素数量的增大，角度视图的数量也必须增加，以保持相同的外围分辨率。这在重建较大的对象或者以较好的分辨率重建时是一个重要因素。对于图像尺寸是N×N的情况，投影的数量M通常应该与N处于相同级别，因此，问题的复杂性是O(N4)。大多数计算集中在反投影中。
图2显示了Intel x86体系结构图表。Intel处理器200具有多个执行单元(EU)以同时执行整型和浮动操作。在Intel PentiumIII中，有两个整型单元(ALU0202和ALU1204)，一个浮动单元(FPU)206，一个MMX单元208，用于并行处理8位和16位整型，一个流式SIMD执行单元(SSE)210，用于并行处理四个32位单精度浮点数据。还提供了地址生成单元212、存储器负载214、存储地址计算单元216、存储器存储218和保留站(reservationstation)220。SSE中的SIMD指令允许在一个指令中进行四个浮动整型操作。Pentium III处理器具有五个流水线222以在指令执行中利用并行性。在每一个时钟周期，处理器核心可以在一个端口上向五个流水线中的任何一个流水线发送零个或一个μop，以便达到每个周期五个μop的最大发行(issue)带宽。每一个流水线都连接到不同的EU。对于用于Intel处理器的C/C++编译器的所有不同版本，普通C/C++代码只能形成利用整型单元ALU和浮动单元FPU的代码。为使用MMX单元和SSE单元的硬件资源，使用特殊的内部或人工编写的汇编程序代码。Pentium III处理器中有两级高速缓存。第1级(L1)高速缓存是芯片内(on-chip)高速缓存子系统，包括两个16 Kbyte路集相关联高速缓存(set associativecache)，对于指令和数据，带有32字节的缓存线长度。数据高速缓存具有八个交织在四个字节边界上的存储体。第2级2(L2)高速缓存是片外的，但在同一个处理器包内。其大小通常介于128K字节到1M字节之间。对于数据访问，L2通常具有4到10个周期的等待时间。当处理器需要获取指令或数据时，L1比L2快得多，L2比访问主存储器快。
可以用多个处理器构建多处理器计算机；现在市场上销售的大多数这样的计算机是带有两个或四个处理器的SMP。通常，在计算机上运行的操作系统具有诸如多线程和多进程之类的技术，以利用多处理器的硬件资源。例如，Microsoft Windows具有Win32线程，而Unix/Linux具有pthread支持。可以预期，优选实施例将最常用Windows NT/2000来实现，这些操作系统具有多线程和纤维支持，这会自动在多处理器环境中重建并行性。
可以使用一些技巧来最大限度地减少反投影循环中操作的实际数量，例如，通过更改x、y和z增大的顺序；向重建立体像素应用特殊边界；预内插投影数据，以允许实际内循环中的尽可能最简单的内插法。单处理器计算机系统的基本性能可以用下列公式来表示
T＝n*CPI*t    (3)
其中T是执行的总时间，n是执行的指令的数量，t是每个指令的时间，CPI是每个指令的周期的数量。缩短时钟时间t是工程设计的关键所在。一般来说，较小的、较快的电路会产生更好的时钟速度。
减少其他两个因素涉及并行性的某些版本。有多个并行性的级别：首先，单个程序可以分为几个组成部分，不同的处理器计算每一个部分，这叫做“程序级并行性”。其次，诸如流水线之类的技术可以通过重叠指令的执行获得更大的吞吐量，这叫做“指令级并行性”。最后，低级别并行性主要对算术逻辑单元的设计人员有意义，对用户相对不可见，这叫做“算术和位级的并行性”。优选实施例主要依赖程序级并行性和指令级并行性。程序级并行性在程序的独立的部分或者在循环的单次迭代中显示。这样的并行性可以通过使用多处理器来加以利用。指令级并行性具有两个基本类型：单个的指令在处理器中重叠(在相同时间执行)，一个给定指令被分解成子操作，然后子操作被重叠。
如费尔德坎普算法所描述的，使用一组M个投影，每一个投影都具有大小为N×N像素，以重建N3立方体。每一个投影都需要N3次循环计算才能执行反投影。M个投影需要M*N3次循环计算。通常，M应该与N处于同一级别才能获得较好的结果。总的实际重建时间可以表示为：
Trecon＝k*tuint*O(N4)                     (4)
如果算法不变，费尔德坎普算法的O(N4)计算复杂性不能降低。然而，由于总时间还取决于因数k和循环单位时间tunit，所以较小的因数k和较短的反投影循环单位时间tunit将缩短重建时间。
可以用如图3所示的统一建模语言(UML)顺序图讲述普通FACBR实施方式。当FACBR进程开始时，M个投影数据一个接一个加载，每一个都根据公式(1)过滤；然后过滤数据用于反投影到每一个立体像素(x，y，z)。
在对所有数据完成反投影之后，重建完成，数据被保存或以三维显示呈现。通常，过滤时间大约是反投影时间的1/15到1/30或更小。在根据图3的采用C/C++代码的实施方式中，在FACBR进程期间，只使用FPU和较少部分ALU，最强大的EU最大功效的被浪费。
巳知在现有技术中，费尔德坎普算法实施方式中有四个可能的并行性形式：像素并行性、投影并行性、射线并行性和操作并行性。在重建进程中，所有立体像素和投影都是彼此独立的，射线可以独立地反投影。每一个投影的过滤和反投影操作都是独立的；甚至可以独立地划分低级别乘法和加法操作。要实现快速锥形线束重建，在FACBR实施方式中应用了下列方法：
a)将FDK反投影过程分为两个阶段：投影图生成以计算(u，t，s)和数据反投影以计算f(x，y，z)。公式(2)显示，(u，t，s)只取决于(x，y)，这样，投影图只需要O(N3)计算时间；因此k和tunit都会减小。
b)使用一些先验(priori)知识生成如球面或柱面的一些边界；对于在边界外面的并且不能重建的一些立体像素，可以跳过计算，从而可以提供较小的k。如果重建的立体像素被可视化为长度为N的立方体，那么，立体像素的完整的数量大约为N3，但利用柱面边界，立体像素的缩小的数量可以将k缩小π/4，而球面边界将使k缩小π/6。
c)FACBR进程要求大容量存储器；对于N＝512，一个投影或一个切片的数据将是1兆字节。这甚至比L2的大小还大。tunit实际考虑了(x，y，z)处的每一个立体像素的存储器存取时间和计算时间。在重建进程期间高速缓存丢失将会降低性能，因此数据应该在存储器中这样排列，以便大多数数据访问都靠近处理器，即，处理器从L1获取的数据比从L2获取更多，从L2获取的数据比从主存储器获取更多。
d)要获得较短的tunit，反投影循环核心需要优化。可以使用人工编写的汇编语言控制EU并行地工作，这是混合型执行(HE)以降低k和tuniT。由于浮点数据始终比定点数据占用更长的计算时间，所以部分中间结果可以以定点数据处理，以便混合型数据(HD)用于降低tunit。
e)当SMP可用时，重建数据和投影数据处理可以分成多个部分并在不同的处理器上运行。对于n处理器SMP，当不能被转换到并发工作的任务的部分为f，则根据Amdahl定律，k值按照n/(1+(n-1)f)的理论加速而减少。多处理器计算机按照多线程实施方式工作，并在多个处理器之间谨慎地分配任务。如上所述，巳知有能够以这样的方式控制多处理器计算机的操作系统。对于单处理器，上下文环境切换将牺牲CPU时间，因此可能会实际降低性能，因此，希望多线程方法只有在SMP可用的情况下才使用。
在具有下面的表1提出的规格的普通PC上实现并测试了上文描述的方法：
表1
 机器   CPU频率   物理内存   SIMD/ILP支持
 Pentium III   500MHz   640MB   MMX，SSE
Microsoft Visual C++6.0和Intel C++4.5用作开发工具。在传统的实施方式中，执行纯粹的浮点(PF)计算，因此，由于编译器，只有FPU被完全使用。固有的和微调的汇编程序代码提供混合型计算方法。即，在使用HD的重建处理期间，使用了浮点和定点计算两者，以便在Pentium III处理器中完全利用EU。
现在描述并行性注意事项。首先是使用混合型执行(HE)和混合型数据(HD)。
HD用于缩小并使ALU单元并行地工作。由于SSE单元与FPU单元和MMX单元无关，故SSE单元可以与ALU单元、MMX单元，甚至还有FPU单元同时工作，由此可以对PF数据或HD实现混合型执行模式。有不同的方法来使用EU以加速FACBR过程。图数据和一些中间结果可以由ALU以定点数据格式处理，重建数据和最后的输出结果可以以浮点格式处理。不同的数据和阶段的该混合型数据格式可以提高EU的效率。对于PF数据，最好的方法是使用MMX单元来调整数据地址和图数据，并使用SSE执行反投影计算。MMX可以为两个或更多的点处理数据地址和图数据，而ALU只能处理一个点。图4中作为UML活动图显示了用于PF的HE方法。
由于Pentium III处理器中的MMX单元只能处理8位和16位整数乘法，因此，为HD数据执行HE不如为PF数据执行HE那么有效。然而，对于诸如Pentium IV处理器中的SSE2之类的新处理器技术，利用HD的混合型执行将获得更大的速度改善。
第二个并行性注意事项是数据分区模式。对于不同轴的重建，重建数据被分区为不同的子单元。图5显示了数据分区方案。数据作为一维数组存储在存储器中，对于z，增大每一个数据点的索引，然后对于x，最后对于y增大索引。数据在z行中进行处理，因为相同的投影数据u值可以用于一个z行。事实上，用于对一个z行中的立体像素执行反投影的投影数据实际上位于两个相邻的u行中，因为四个相邻的点(每一个相邻的u行中有两个)用于为z行中的一个立体像素内插数据。因此，容易将投影行数据和重建z行预取到高速缓存，因为一行只占用4N字节。对于N＝512，整个数据行只需要6K字节，这适合于L1和L2高速缓存。在所有立体像素都被重建之后，执行特殊的现场三维转置操作，以将重建数据改变为一维数组，该数组的索引从x，然后是y，最后是z增大。转置操作是现场执行的，因为整个N3立方体占用4N3字节的存储器。沿着Z轴的对称也可以在重建一个Z行中使用，以节省计算图数据的时间。
为确保改进锥形线束重建速度时的精度，首先使用计算机模拟幻影确定实施方式的重建精度。其次，使用纯粹的浮点实施方式和HD计算实施方式量化重建误差噪声电平和重建图像的一致性，将两个实施方式的重建结果与模拟幻影和实验幻影数据进行比较。在判断HC实施方式不会产生人为的和无法接受的重建误差之后，与普通纯浮点计算重建相比评估MC实施方式的加速。在实际工作中，实验幻影数据也用于评估该实施方式的有效性。
下面的表2中显示的两个模拟幻影用于评估精度。SheppLogan幻影用作一般精度误差比较参考。柱面幻影用于比较不同z位置的精度误差。通常，费尔德坎普算法在中心切片具有最好的结果，而对于两端的切片，精度误差增大。柱面幻影用于检查HD和PF精度误差是否随着到中心的z距离而变化。
表2
                                         带有11个椭圆的Shepp Logan幻影
  a   b   c   x0   y0   z0   α   β   μ
  0.6900   0.920   0.900   0.000   0.000   0.000   0.0   0.00   200
  0.6624   0.874   0.880   0.000   0.000   0.000   0.0   0.00   -0.98
  0.4100   0.160   0.210   -0.22   0.000   -0.25   -72.0   0.00   -0.02
  0.3100   0.110   0.220   0.220   0.000   -0.25   72.0   0.00   -0.02
  0.2100   0.250   0.500   0.000   0.350   -0.25   0.00   0.00   0.01
  0.0460   0.046   0.046   0.000   0.100   -0.25   0.00   0.00   0.02
  0.0460   0.023   0.020   -0.08   -0.605   -0.25   0.00   0.00   0.01
  0.0460   0.023   0.020   0.060   -0.605   -0.25   90.0   0.00   0.01
  0.0560   0.040   0.100   0.060   -0.105   0.625   90.0   0.00   0.02
  0.0560   0.056   0.100   0.000   -0.100   0.625   0.00   0.00   -0.02
  0.0230   0.023   0.023   0.000   -0.605   -0.25   0.00   0.00   0.01
                                              带有11个柱面的柱面幻影
  a   b   c   x0   y0   z0   α   β   μ
  0.6900   0.6900   0.900   0.000   0.000   0.000   0.00   0.00   200
  0.6624   0.6624   0.880   0.000   0.000   0.000   0.00   0.00   -1.00
  0.0800   0.0800   0.800   0.000   0.000   0.000   0.00   0.00   0.05
  0.0800   0.0800   0.800   0.250   0.000   0.000   0.00   0.00   0.20
  0.0800   0.0800   0.800   0.500   0.000   0.000   0.00   0.00   -0.20
  0.0800   0.0800   0.800   -0.25   0.000   0.000   0.00   0.00   0.20
  0.0800   0.0800   0.800   -0.50   0.000   0.000   0.00   0.00   -0.20
  0.0800   0.0800   0.800   0.000   0.25   0.000   0.00   0.00   0.20
  0.0800   0.0800   0.800   0.000   0.50   0.000   0.00   0.00   -0.20
  0.0800   0.0800   0.800   0.000   -0.25   0.000   0.00   0.00   0.20
  0.0800   0.0800   0.800   0.000   0.50   0.000   0.00   0.00   -0.20
总的FACBR时间包含过滤时间和反投影时间；过滤时间只占总时间的一小部分。在用SSE以PF数据格式优化过滤过程之后，下面的表3显示的2885122投影的时间结果被视为是令人满意的，因此不需要为过滤过程考虑HE或HD优化。
表3
  过滤模式   4N ext   SSE 4N   2N ext   SSE 2N
  时间(秒)   146.891   58.163   78.753   35.310
  加速                      2.5255                      2.2303
反投影过程是FACBR的最耗时部分。测试了五个实施方式的加速度；结果在下面的表4中示出。测试使用了2885122投影来重建5123数据。所有重建都用柱面边界运行。程序在Windows NT4.0中运行，并占用95％到98％的处理器时间。第一列是带有边界的传统的PF方法，第二列是带有数据分区的PF，第三列是HD方法，第四列是带有HE的HD数据，第五列是带有SSE加速的PF，最后一列是带有HE的PF。结果显示，HD比传统的实施方式提供了3到3.5的加速，HE-HD提供了4到5的加速，几乎与带有SSE的PF相同。此结果说明了两点：首先，HE-HD不涉及SSE单元，这样，诸如Celeron之类的比较便宜的处理器可用于获取与SSE-PF几乎相同的性能；其次，可通过使用以与SSE处理浮点数据的同样的方式处理定点数据的功能单元获得较高的加速；这样的功能已经出现在Pentium IV处理器中。HE-PF是Pentium III处理器中最有效的方法。使用球面边界的重建，连同结合了MMX和SSE单元的功能的混合型计算，为5123FACBR提供了5到6的加速和15.03分钟的重建时间。对于较高时钟频率的处理器和SMP计算机，时间结果将更好。
表4
  切片和参数   普通FD   PF   HD   HE-HD   SSEPF   HE-PF
  64   时间(秒)   870   618   263   166   169   148
  加速   1.0   1.408   3.308   5.241   5.148   5.878
  28   时间(秒)   I799   1256   529   348   349   308
  加速   1.0   1.432   3.401   5.170   5.155   5.841
  256   时间(秒)   3288   2424   1040   681   744   671
  加速   1.0   1.356   3.149   4.828   4.419   4.900
  512   时间(秒)   6562   4714   2047   1516   1634   1101
  加速   1.0   1.392   3.206   4.328   4.016   5.960
下面的表5用于显示了HE-PD和HE-PF方法的不同切片的有效tunit。由于程序在多处理器操作系统中运行，处理器时间资源随着时间而变化，因此，有效的tunit也随着时间而变化。基本上，tunit随着切片数量增大而变得稳定。当切片数量少于4或者数据不是16字节对准时，处理器不能使用SSE，且tunit比SSE可用时大。因此，单切片的时间可以比其他切片大。
表5
  方法                              切片数量
  1   32   64   128   256   512
  HE-PF   时间(秒)   4.126   68   148   308   671   1101
  tunit(纳秒)   54.6   28.1   30.6   31.9   34.7   28.4
  HE-HD   时间(秒)   5.033   76   166   348   681   1516
  tunit(纳秒)   66.6   31.6   34.4   36.0   35.2   39.2
现在来比较精度。两个图像P和Q之间的相对误差按如下公式来计算
E = 1 N pixel Σ | P i , j - Q i , j | 1 N pixel Σ | Q i , j | × 100 % - - - ( 5 )
首先，为每一个像素计算误差比。然后，计算平均误差比，以便进行感兴趣区域(ROI)的完全比较。
由于HE-PF可以与浮点数据一起使用，因此，与传统的PF方法相比，不会造成额外的精度误差。最大的问题是HD计算是否会产生比较大的精度误差。如果PF重建的图像和幻像之间的相对精度误差是EPF，HD重建的图像和幻像之间的相对误差是EHD，HD重建的图像和PF重建的图像之间的相对误差是EHP，混合型计算误差与整个精度误差的比被定义为：
R HD = | E HD - E PF | E HD × 100 % - - - ( 6 )
确定HD重建的图像相对于模拟的幻像和PF重建的图像的精度误差。对于Shepp Logan幻影，HD图像和PF图像之间的精度误差小于0.03％；对于柱面幻影，EHP小于0.02％。因此，HD图像与PF图像相比保持了良好的精度。EHP对模拟幻像产生总误差百分比的5％以下。这意味着，该算法导致了总误差的95％以上。HD图像具有足够的精度并可与PF图像相比。
图6显示了可以实现本发明的设备，该图是从上文引用的美国专利No.5,999,587的图9复制而来。在标准CT中，通过堆积一系列切片获得三维重建。在卷CT中，可以获得对象的直接重建。现在请参看图6，该图显示了本发明的锥形线束层析成象系统600如何用于获得对象的直接的三维重建。应该理解，以简化方框图的形式说明了锥形线束卷CT扫描设备600。本发明最好可以与这样的锥形线束卷CT扫描设备结合使用以生成对象的三维重建矩阵。基于三维重建矩阵，可以获得所期望的三维显示。
锥形线束卷CT扫描设备使用遍历身体中一组路径的锥形的辐线束604检查身体P。如图6所示，x射线源610和二维检测器611安装在围绕待检查的身体P旋转的龙门架602上。x射线源的工作电压是从常规的高压发生器608以这样的方式获取的，即当向x射线源610施加高电压时其产生所期望的锥形线束。高压发生器608用电源618通过开关616通电。需要时，还可以使用造影溶液(contrast solution)喷射器640。
电源618还为第一电机612提供电源，以便它在围绕身体的轨道上，例如，在靠近框架的箭头显示的顺时针方向驱动龙门架602。电源618由开关620或其他常规控制设备打开，以便启动测量序列。速度控制电路614用于控制龙门架602的转速，并提供表示何时电机712的速度处于进行测量所期望的级别的输出控制信号。还可以利用来自转动控制614的输出来操作开关616，以便高压发生器608只有在龙门架602以进行测量所需要的速度被驱动的情况下才打开。
为了获取如前面所讨论的弧形测量，利用倾斜控制615通过龙门架倾斜电机613使龙门架602倾斜相对小的角度±15°到±30°。该倾斜允许获取有关垂直弧形的弧形投影数据。这样的几何形状导致对应于放大率为1.5的探测器611处的37-60cm场的直径为25-40cm的对象的整套数据。虽然一般来说在标准CT台架中都可以进行台架602的倾斜，但是，要获得弧形投影，必须要对标准CT台架进行最小的修改，以便台架的倾斜、x射线曝光时间和投影获取由系统控制计算机624进行同步，如图6所示。
除上文描述的获取圆形和弧形投影的方法之外，还可以以下列两种方式之一实现圆形加弧形的几何形状。在三种方法中的第一个并且也是优选方法中，台架602倾斜一个小的角度(±15°到±30°)，然后，在台架602倾斜的同时旋转x射线管610和二维探测器611。只有在x射线管610和二维探测器611处于旋转角为0°的情况下才能获取半组弧形投影。当倾斜角变成零时，在预置的旋转角位置将获取圆形投影。当圆形投影获取完成时，台架602将倾斜-15°到-30°。只有在x射线管610和二维探测器611处于旋转角为0°的情况下才能获取另外半组弧形投影。
第二个备选方法是在机械结构方面修改标准CT台架，以便向台架添加两个短弧形轨道，x射线管610和二维探测器611可以在该弧形上移动以获得弧形投影，并可在圆形上移动以获得圆形投影。一个弧构成了x射线管610的轨道，另一个弧是二维探测器611的轨道。两个弧形轨道彼此相隔180°。x射线管610和二维探测器611同步地在弧形轨道上移动，以获得弧形投影。然后，x射线管610和二维探测器611在台架上旋转，以获得圆形投影。
在龙门架602上x射线源610的对面安装了二维探测器611，该探测器具有等于或大于1000∶1的动态范围和小于10％的影像延迟，例如，硒薄膜晶体管(STFT)阵列或硅STFT阵列，以便提供对应于x射线衰减信号模式的二维投影。x射线源610和二维探测器611以这样的方式安装在龙门架602上，以便它们两个可以同步移动。
由x射线源610生成的辐射604的锥形线束通过身体或者接受检查的对象被投射。二维探测器611测量跨锥面沿着线束路径集传输的辐射。
或者，连续几个二维探测器(未显示)可固定地安装在龙门架602的附近，x射线源610安装到龙门架上，以便在龙门架旋转时，锥形辐线束604穿过接受检查的身体P投射，并连续地由系列探测器中的每一个接收。
二维投影获取控制和A/D转换单元626，在从系统控制计算机624(包括时钟622)连续获取的扫描脉冲的控制之下，接收对应于二维探测器611的不同行的输出序列。二维探测器的每一行都包括许多探测单元(至少＞100)。每一探测器单元的输出表示可沿着其中一个相应的线束路径测量的衰减值的线积分。锥形线束604形成一个锥角，足以包括身体的感兴趣的整个区域。因此，对于对象的完整的扫描可以通过只将支撑x射线源610和二维探测器611的龙门架602围绕身体转动以在不同的角位置获得二维投影信号来进行。
模数转换单元626用来将投影信号数字化并将它们保存在三维图像重建阵列处理器628和存储设备630中。三维图像重建阵列处理器628所使用的方法是本申请中描述的本发明的算法和方法。三维图像重建阵列处理器628用来将数字化的投影信号转换为x射线衰减数据矢量。x射线衰减数据矩阵对应于被检查的身体躯干内的空间网格位置处的x射线衰减。矩阵的每一个数据元素都代表x射线衰减值，且元素的位置对应于身体内的相应的三维网格位置。
根据前面讨论的本发明的原理，显示处理器632获取作为三维x射线衰减信号模式存储在存储器630中的数据，并按前面描述的方法处理该数据，然后，所期望的三维图像显示在三维显示设备634上。例如三维图像重建阵列处理器632可以是如上所述的带有一个或多个Intel或Intel兼容的x86级的微处理器的计算机。然而，也可以使用任何能够执行相同或基本上相同的并行操作的处理器。
虽然上文阐述了本发明的优选实施例，但是，阅读了本说明书的那些本技术的普通技术人员能够很容易理解，在本发明的范围内可以实现其他实施例。例如，具体的数值都是说明性的，而不是限制性的，所提及的具体的商用产品也是如此。本发明不是特别针对费尔德坎普算法的，而是可用于有效而最佳地实现任何过滤反投影锥形线束算法。本发明也不特别针对x86处理器；相反，本发明可以与能够实现上文描述的算法的任何处理器一起使用，并能够与以下的任何处理器有特定的用途：该处理器具有可以在一个指令集内处理一个以上单精度32位数据的浮点处理单元，以及可以在一个指令集内处理一个以上的16位或32位数据的定点处理单元。因此，本发明应该理解为只受所附的权利要求的限制。

Claims (33)

1.一种用于生成表示一个对象的内部的三维图像的系统，该系统包括：
射线扫描器，该扫描器通过将射线穿过该对象传递到探测器上来生成投影信号；以及
接收投影信号的计算机，用于通过对投影信号执行多种计算来生成三维图像，该计算机包括至少一个定点处理单元和至少一个浮点处理单元，该至少一个定点处理单元与该至少一个浮点处理单元并行操作，该计算机将该多种计算分为将在该至少一个定点处理单元中执行的第一组计算和将在该至少一个浮点处理单元中执行的第二组计算。
2.根据权利要求1所述的系统，其中：
所述第一组计算包括生成投影图；以及
所述第二组计算包括根据投影图对锥形线束投影信号进行反投影，以便产生图像。
3.根据权利要求2所述的系统，其中：
生成投影图的过程包括将世界坐标系(x，y，z)映射到探测器的坐标系(u，t，s)，其中u独立于z；
所述计算机将锥形线束投影信号组织为z行，其中z可以变化，但x和y保持恒定；以及
该计算机为z行的每一个执行反投影。
4.根据权利要求3所述的系统，其中，所述计算机进一步包括一个足够大而可以保存z行之一的高速缓存。
5.根据权利要求3所述的系统，其中，在所述计算机为所有z行执行反投影以作为多个立体像素形成图像之后，该计算机对图像中的立体像素执行三维转置操作，以将该立体像素组织为x行，其中x可以变化，但y和z保持恒定。
6.根据权利要求1所述的系统，其中：
对象的一个边界已知是一个先验；以及
三维图像只在该边界内生成。
7.根据权利要求1所述的系统，其中，投影信号作为纯粹的浮点数据来处理。
8.根据权利要求1所述的系统，其中，投影信号作为浮点和定点数据的混合来处理。
9.根据权利要求1所述的系统，其中，所述计算机包括微处理器，在该微处理器上实现该至少一个定点处理单元和该至少一个浮点处理单元，该定点处理单元可在一个指令集内处理一个以上16位或32位整型数据，该浮点处理单元可在一个指令集内处理一个以上单精度32位浮点数据。
10.根据权利要求9所述的系统，其中，所述计算机包括多个所述的微处理器，每一个微处理器包括至少一个所述的定点处理单元和至少一个所述的浮点处理单元，该定点处理单元可在一个指令集内处理一个以上16位或32位整型数据，该浮点处理单元可在一个指令集内处理一个以上单精度32位浮点数据。
11.根据权利要求1所述的系统，其中，图像是所述对象的内部的线性衰减系数分布。
12.根据权利要求1所述的系统，其中，辐射扫描器是一个辐射锥形线束扫描器，该辐线束是辐射锥形线束，且所述投影信号是锥形线束投影信号。
13.一种用于生成表示一个对象的内部的三维图像的方法，该方法包括：
(a)将线束穿过该对象传递到检测器上来生成投影信号；以及
(b)接收投影信号并通过对锥形线束投影信号执行多种计算来生成三维图像；
其中，步骤(b)是在计算机上执行的，该计算机包括至少一个定点处理单元和至少一个浮点处理单元，该至少一个定点处理单元与该至少一个浮点处理单元并行操作，该计算机将该多种计算分为将在该至少一个定点处理单元中执行的第一组计算和将在该至少一个浮点处理单元中执行的第二组计算。
14.根据权利要求13所述的方法，其中：
对象的一个边界已知是一个先验；以及
三维图像只在该边界内生成。
15.根据权利要求13所述的方法，其中，投影信号作为纯粹的浮点数据来处理。
16.根据权利要求13所述的方法，其中，投影信号作为浮点和定点数据的混合来处理。
17.根据权利要求13所述的方法，其中，所述计算机包括微处理器，在该微处理器上实现所述至少一个定点处理单元和所述至少一个浮点处理单元，该定点处理单元可在一个指令集内处理一个以上16位或32位整型数据，该浮点处理单元可在一个指令集内处理一个以上单精度32位浮点数据。
18.根据权利要求17所述的方法，其中，所述计算机包括多个微处理器，每一个微处理器包括至少一个所述的定点处理单元和至少一个所述的浮点处理单元，该定点处理单元可在一个指令集内处理一个以上16位或32位整型数据，该浮点处理单元可在一个指令集内处理一个以上单精度32位浮点数据。
19.根据权利要求13所述的方法，其中，所述图像是所述对象的内部的线性衰减系数分布。
20.根据权利要求13所述的方法，其中，所述对象包括软组织。
21.根据权利要求20所述的方法，其中，所述图像用于探测软组织中的癌。
22.根据权利要求13所述的方法，其中，步骤(b)是使用人工编写的汇编语言执行的。
23.根据权利要求13所述的方法，其中，步骤(b)是通过在单微处理器或多处理器上使用混合型计算进行并行处理来执行的，以便加速用于重建软组织的锥形线束重建。
24.根据权利要求13所述的方法，其中，步骤(b)只有在多个处理器可用的情况下才通过多线程执行。
25.根据权利要求13所述的方法，其中，所述线束是锥形线束，所述投影信号是锥形线束投影信号。
26.根据权利要求25所述的方法，其中，步骤(b)是使用过滤的反投影锥形线束重建算法执行的。
27.根据权利要求26所述的方法，其中，过滤的反投影锥形线束重建算法是使用采用了单指令多数据技术的混合型计算来执行的。
28.根据权利要求27所述的方法，其中，过滤的反投影锥形线束重建算法是使用多线程在多个处理器上执行的。
29.根据权利要求25所述的方法，其中，步骤(b)是通过费尔德坎普锥形线束重建算法执行的。
30.根据权利要求25所述的方法，其中：
所述第一组计算包括生成投影图；以及
所述第二组计算包括根据投影图对锥形线束投影信号进行反投影，以便产生重建的图像。
31.根据权利要求30所述的方法，其中：
生成投影图的过程包括将世界坐标系(x，y，z)映射到探测器的坐标系(u，t，s)，其中u独立于z；
所述计算机将锥形线束投影信号组织为z行，其中z可以变化，但x和y保持恒定；以及
该计算机为z行的每一个执行反投影。
32.根据权利要求31所述的方法，其中，所述计算机进一步包括一个足够大而可以保存z行之一的高速缓存。
33.根据权利要求31所述的方法，其中，在所述计算机为所有z行执行反投影以作为多个立体像素形成图像之后，该计算机对该图像中的立体像素执行三维转置操作，以将该立体像素组织为x行，其中x可以变化，但y和z保持恒定。
CNB028050894A 2001-02-16 2002-02-13 使用一个或多个微处理器的快速并行锥形线束重建系统和方法 Expired - Lifetime CN1284122C (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US09/784,331 US6477221B1 (en) 2001-02-16 2001-02-16 System and method for fast parallel cone-beam reconstruction using one or more microprocessors
US09/784,331 2001-02-16

Publications (2)

Publication Number Publication Date
CN1491404A CN1491404A (zh) 2004-04-21
CN1284122C true CN1284122C (zh) 2006-11-08

Family

ID=25132101

Family Applications (1)

Application Number Title Priority Date Filing Date
CNB028050894A Expired - Lifetime CN1284122C (zh) 2001-02-16 2002-02-13 使用一个或多个微处理器的快速并行锥形线束重建系统和方法

Country Status (6)

Country Link
US (1) US6477221B1 (zh)
EP (1) EP1366469B1 (zh)
CN (1) CN1284122C (zh)
AU (1) AU2002251922B2 (zh)
CA (1) CA2438387A1 (zh)
WO (1) WO2002067223A2 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101664583B (zh) * 2009-09-09 2012-05-09 深圳市海博科技有限公司 基于cuda的剂量计算优化方法和系统
CN102877828A (zh) * 2012-09-09 2013-01-16 山西山地物探技术有限公司 一种三维多井联合井地ct成像方法

Families Citing this family (50)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6987831B2 (en) 1999-11-18 2006-01-17 University Of Rochester Apparatus and method for cone beam volume computed tomography breast imaging
US7430271B2 (en) * 2000-11-13 2008-09-30 Digitome Corporation Ray tracing kernel
US20020169680A1 (en) * 2001-05-10 2002-11-14 International Business Machines Corporation Method and apparatus for building commercial distributed computing networks via computer cost subsidization
US6741730B2 (en) * 2001-08-10 2004-05-25 Visiongate, Inc. Method and apparatus for three-dimensional imaging in the fourier domain
US6771733B2 (en) * 2001-08-16 2004-08-03 University Of Central Florida Method of reconstructing images for spiral and non-spiral computer tomography
US6638226B2 (en) * 2001-09-28 2003-10-28 Teratech Corporation Ultrasound imaging system
JP3870105B2 (ja) * 2002-02-22 2007-01-17 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 逆投影方法およびｘ線ｃｔ装置
US6904117B2 (en) * 2002-10-30 2005-06-07 Toshiba Corporation Tilted gantry helical cone-beam Feldkamp reconstruction for multislice CT
DE10304662A1 (de) * 2003-02-05 2004-08-19 Siemens Ag Verfahren zur Erzeugung von Bildern in der Computertomographie mit einem 3D-Bildrekonstruktionsverfahren
AU2004244308B2 (en) * 2003-05-27 2008-04-24 Clean Earth Technologies, Llc Method for fast image reconstruction with compact radiation source and detector arrangement using computerized tomography
US7134036B1 (en) * 2003-12-12 2006-11-07 Sun Microsystems, Inc. Processor core clock generation circuits
US7362843B2 (en) * 2004-09-23 2008-04-22 General Electric Company System and method for reconstruction of cone beam tomographic projections with missing data
WO2006116316A2 (en) 2005-04-22 2006-11-02 University Of Chicago Open source trajectory method and apparatus for interior imaging
US7991242B2 (en) 2005-05-11 2011-08-02 Optosecurity Inc. Apparatus, method and system for screening receptacles and persons, having image distortion correction functionality
EP1886257A1 (en) 2005-05-11 2008-02-13 Optosecurity Inc. Method and system for screening luggage items, cargo containers or persons
US7492858B2 (en) * 2005-05-20 2009-02-17 Varian Medical Systems, Inc. System and method for imaging and treatment of tumorous tissue in breasts using computed tomography and radiotherapy
US7646842B2 (en) * 2005-09-23 2010-01-12 General Electric Company Methods and apparatus for reconstructing thick image slices
US20070132754A1 (en) * 2005-12-12 2007-06-14 Intel Corporation Method and apparatus for binary image classification and segmentation
CN101297325B (zh) * 2005-12-29 2013-04-24 英特尔公司 用于射线跟踪的方法和装置
US20070274435A1 (en) * 2006-02-27 2007-11-29 Ruola Ning Phase contrast cone-beam CT imaging
EP1988826A4 (en) * 2006-02-27 2010-05-19 Univ Rochester METHOD AND DEVICE FOR COMPUTERIZED TOMOGRAPHY DYNAMIC IMAGING WITH CONICAL BEAM
US20070237287A1 (en) * 2006-03-28 2007-10-11 Predrag Sukovic Ct scanner with automatic determination of volume of interest
WO2007124338A1 (en) * 2006-04-19 2007-11-01 Xoran Technologies, Inc. Ct scanner with untracked markers
DE102007020879A1 (de) 2006-05-10 2009-04-02 Gachon University Of Medicine & Science Industry-Academic Cooperation Foundation Verfahren und Vorrichtung für die äußerst schnelle Symmetrie- und SIMD- gestützte Projektion/Rückprojektion für die 3D-PET-Bildrekonstruktion
US7899232B2 (en) 2006-05-11 2011-03-01 Optosecurity Inc. Method and apparatus for providing threat image projection (TIP) in a luggage screening system, and luggage screening system implementing same
CN100386779C (zh) * 2006-06-02 2008-05-07 清华大学 基于通用图形显示卡的被测体正投影与反投影方法
US8494210B2 (en) 2007-03-30 2013-07-23 Optosecurity Inc. User interface for use in security screening providing image enhancement capabilities and apparatus for implementing same
DE102006036327A1 (de) * 2006-08-03 2008-02-14 Siemens Ag Verfahren zum Bereitstellen von 3D-Bilddaten und System zum Aufnehmen von Röntgenbildern
US8217937B2 (en) * 2007-03-28 2012-07-10 The Aerospace Corporation Isosurfacial three-dimensional imaging system and method
EP2219525B1 (en) * 2007-08-23 2017-01-04 Bearf, Llc Improved computed tomography breast imaging and biopsy system
US8023767B1 (en) 2008-03-10 2011-09-20 University Of Rochester Method and apparatus for 3D metal and high-density artifact correction for cone-beam and fan-beam CT imaging
US7940891B2 (en) 2008-10-22 2011-05-10 Varian Medical Systems, Inc. Methods and systems for treating breast cancer using external beam radiation
US8327345B2 (en) * 2008-12-16 2012-12-04 International Business Machines Corporation Computation table for block computation
US8407680B2 (en) * 2008-12-16 2013-03-26 International Business Machines Corporation Operand data structure for block computation
US8285971B2 (en) * 2008-12-16 2012-10-09 International Business Machines Corporation Block driven computation with an address generation accelerator
US8281106B2 (en) * 2008-12-16 2012-10-02 International Business Machines Corporation Specifying an addressing relationship in an operand data structure
US8458439B2 (en) * 2008-12-16 2013-06-04 International Business Machines Corporation Block driven computation using a caching policy specified in an operand data structure
WO2010093357A1 (en) * 2009-02-11 2010-08-19 Tomotherapy Incorporated Target pedestal assembly and method of preserving the target
US7949095B2 (en) 2009-03-02 2011-05-24 University Of Rochester Methods and apparatus for differential phase-contrast fan beam CT, cone-beam CT and hybrid cone-beam CT
KR102067367B1 (ko) 2011-09-07 2020-02-11 라피스캔 시스템스, 인코포레이티드 적하목록 데이터를 이미징/검출 프로세싱에 통합시킨 x-선 검사 방법
EP2822468B1 (en) 2012-03-05 2017-11-01 University Of Rochester Methods and apparatus for differential phase-contrast cone-beam ct and hybrid cone-beam ct
US9364191B2 (en) 2013-02-11 2016-06-14 University Of Rochester Method and apparatus of spectral differential phase-contrast cone-beam CT and hybrid cone-beam CT
WO2014133849A2 (en) 2013-02-26 2014-09-04 Accuray Incorporated Electromagnetically actuated multi-leaf collimator
JP6195980B2 (ja) * 2013-06-05 2017-09-13 エーファウ・グループ・エー・タルナー・ゲーエムベーハー 圧力マップを測定する測定装置及び方法
CN105806858B (zh) * 2014-12-31 2019-05-17 北京固鸿科技有限公司 Ct检测方法和ct设备
EP3772702A3 (en) 2016-02-22 2021-05-19 Rapiscan Systems, Inc. Methods for processing radiographic images
US10482632B2 (en) * 2017-04-28 2019-11-19 Uih America, Inc. System and method for image reconstruction
CN109157215B (zh) * 2018-08-29 2021-09-28 中国医学科学院生物医学工程研究所 一种基于系统矩阵的磁感应磁声电导率图像重建方法
EP4141427A4 (en) * 2020-11-18 2024-06-05 Jed Co., Ltd X-RAY INSPECTION DEVICE
CN116543088B (zh) * 2023-07-07 2023-09-19 有方（合肥）医疗科技有限公司 Cbct图像重建方法及装置

Family Cites Families (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5253171A (en) * 1990-09-21 1993-10-12 General Electric Company Parallel processing method and apparatus based on the algebra reconstruction technique for reconstructing a three-dimensional computerized tomography (CT) image from cone beam projection data
US5257183A (en) 1990-12-21 1993-10-26 General Electric Company Method and apparatus for converting cone beam X-ray projection data to planar integral and reconstructing a three-dimensional computerized tomography (CT) image of an object
US5170439A (en) 1991-06-11 1992-12-08 Picker International, Inc. Cone beam reconstruction using combined circle and line orbits
US5365560A (en) 1991-07-29 1994-11-15 General Electric Company Method and apparatus for acquiring a uniform distribution of radon data sufficiently dense to constitute a complete set for exact image reconstruction of an object irradiated by a cone beam source
US5333164A (en) * 1991-12-11 1994-07-26 General Electric Company Method and apparatus for acquiring and processing only a necessary volume of radon data consistent with the overall shape of the object for efficient three dimensional image reconstruction
US5390226A (en) 1992-07-02 1995-02-14 General Electric Company Method and apparatus for pre-processing cone beam projection data for exact three dimensional computer tomographic image reconstruction of a portion of an object
GB2271261A (en) * 1992-10-02 1994-04-06 Canon Res Ct Europe Ltd Processing image data
US5517602A (en) 1992-12-03 1996-05-14 Hewlett-Packard Company Method and apparatus for generating a topologically consistent visual representation of a three dimensional surface
US5278884A (en) 1992-12-18 1994-01-11 General Electric Company Complete 3D CT data acquisition using practical scanning paths on the surface of a sphere
US5461650A (en) 1993-10-18 1995-10-24 General Electric Company Method and system for pre-processing cone beam data for reconstructing free of interpolation-induced artifacts a three dimensional computerized tomography image
US5400255A (en) 1994-02-14 1995-03-21 General Electric Company Reconstruction of images from cone beam data
US6002738A (en) * 1995-07-07 1999-12-14 Silicon Graphics, Inc. System and method of performing tomographic reconstruction and volume rendering using texture mapping
US5671265A (en) 1995-07-14 1997-09-23 Siemens Corporate Research, Inc. Evidential reconstruction of vessel trees from X-ray angiograms with a dynamic contrast bolus
JPH09149902A (ja) 1995-12-01 1997-06-10 Hitachi Medical Corp Ｘ線断層撮影方法および装置
CA2227531C (en) * 1997-01-20 2003-03-18 Hitachi, Ltd. Graphics processing unit and graphics processing system
US6078638A (en) * 1998-09-30 2000-06-20 Siemens Corporate Research, Inc. Pixel grouping for filtering cone beam detector data during 3D image reconstruction
US6343108B1 (en) * 1999-06-18 2002-01-29 Philips Medical Systems (Cleveland), Inc. Cone beam scanner using oblique surface reconstructions
US6201849B1 (en) * 1999-08-16 2001-03-13 Analogic Corporation Apparatus and method for reconstruction of volumetric images in a helical scanning cone-beam computed tomography system

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101664583B (zh) * 2009-09-09 2012-05-09 深圳市海博科技有限公司 基于cuda的剂量计算优化方法和系统
CN102877828A (zh) * 2012-09-09 2013-01-16 山西山地物探技术有限公司 一种三维多井联合井地ct成像方法

Also Published As

Publication number Publication date
AU2002251922B2 (en) 2008-01-03
US20020154727A1 (en) 2002-10-24
EP1366469B1 (en) 2017-10-04
EP1366469A2 (en) 2003-12-03
CN1491404A (zh) 2004-04-21
WO2002067223A3 (en) 2002-10-24
US6477221B1 (en) 2002-11-05
CA2438387A1 (en) 2002-08-29
WO2002067223A2 (en) 2002-08-29

Similar Documents

Publication Publication Date Title
CN1284122C (zh) 使用一个或多个微处理器的快速并行锥形线束重建系统和方法
AU2002251922A1 (en) System and method for fast parallel cone-beam reconstruction using one or more microprocessors
Noël et al. GPU-based cone beam computed tomography
Chou et al. A fast forward projection using multithreads for multirays on GPUs in medical image reconstruction
CN101596113B (zh) 一种ct并行重建系统及成像方法
US20060109952A1 (en) Fan-beam and cone-beam image reconstruction using filtered backprojection of differentiated projection data
CN100464707C (zh) 三维锥束ct图像重建的处理系统
CN101283913A (zh) Ct图像重建的gpu加速方法
US7209535B2 (en) Fourier space tomographic image reconstruction method
CN1063171A (zh) 再现三维计算机断层扫描图象的方法和装置
Scherl et al. Evaluation of state-of-the-art hardware architectures for fast cone-beam CT reconstruction
CN105118039B (zh) 实现锥束ct图像重建的方法及系统
Pedemonte et al. GPU accelerated rotation-based emission tomography reconstruction
Liu et al. GPU-based branchless distance-driven projection and backprojection
Chen et al. A hybrid architecture for compressive sensing 3-D CT reconstruction
Lu et al. Cache-aware GPU optimization for out-of-core cone beam CT reconstruction of high-resolution volumes
Li Design of an FPGA-based computing platform for real-time three-dimensional medical imaging
Ni et al. Review of parallel computing techniques for computed tomography image reconstruction
Scherl et al. Implementation of the FDK algorithm for cone-beam CT on the cell broadband engine architecture
Yu et al. High-speed cone-beam reconstruction on PC
CN101268950B (zh) 基于cell宽频引擎的螺旋ct精确重建系统
CN2860349Y (zh) 三维锥束ct图像重建的处理系统
Xu et al. Mapping iterative medical imaging algorithm on cell accelerator
Boukhamla et al. Parallelization of filtered back-projection algorithm for computed tomography
Shih et al. Fast algorithm for X-ray cone-beam microtomography

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract

Assignee: Koning Corp.

Assignor: Univ. of Rochester

Contract record no.: 2011990000984

Denomination of invention: System and method for fast parallel cone-beam reconstruction using one or more microprocessor

Granted publication date: 20061108

License type: Exclusive License

Open date: 20040421

Record date: 20111019

Assignee: Corning (Tianjin) medical equipment Co., Ltd.

Assignor: Koning Corp.

Contract record no.: 2011990000985

Denomination of invention: System and method for fast parallel cone-beam reconstruction using one or more microprocessor

Granted publication date: 20061108

License type: Exclusive License

Open date: 20040421

Record date: 20111019

CX01 Expiry of patent term
CX01 Expiry of patent term

Granted publication date: 20061108