CN104952043B - 图像滤波方法及ct 系统 - Google Patents

图像滤波方法及ct 系统 Download PDF

Info

Publication number
CN104952043B
CN104952043B CN201410118326.4A CN201410118326A CN104952043B CN 104952043 B CN104952043 B CN 104952043B CN 201410118326 A CN201410118326 A CN 201410118326A CN 104952043 B CN104952043 B CN 104952043B
Authority
CN
China
Prior art keywords
image
filtering
pixel
calculated
dimension
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 - Fee Related
Application number
CN201410118326.4A
Other languages
English (en)
Other versions
CN104952043A (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.)
Hitachi Ltd
Original Assignee
Hitachi Ltd
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 Hitachi Ltd filed Critical Hitachi Ltd
Priority to CN201410118326.4A priority Critical patent/CN104952043B/zh
Priority to US15/129,260 priority patent/US9984433B2/en
Priority to JP2016555324A priority patent/JP6259531B2/ja
Priority to PCT/CN2015/075123 priority patent/WO2015144069A1/zh
Publication of CN104952043A publication Critical patent/CN104952043A/zh
Application granted granted Critical
Publication of CN104952043B publication Critical patent/CN104952043B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/20—Image enhancement or restoration using local operators
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • 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
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04N—PICTORIAL COMMUNICATION, e.g. TELEVISION
    • H04N5/00—Details of television systems
    • H04N5/30—Transforming light or analogous information into electric information
    • H04N5/32—Transforming X-rays
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • Signal Processing (AREA)
  • Image Processing (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明涉及一种图像滤波方法及CT系统，利用多核处理器，并行地对图像进行滤波处理。首先，根据图像的维数及滤波处理的预定的邻域范围，确定滤波处理的多个计算方向。然后，按每个计算方向，由多个线程中的每个线程分别针对图像中的一行像素，以该行的每个像素作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行预定的滤波计算，并将滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果。最后，将针对全部计算方向得到的滤波结果累加，从而获得图像滤波结果。由此，能够减少现有并行计算中的重复计算，大幅提高图像滤波的并行计算速度。

Description

图像滤波方法及CT系统
技术领域
本发明涉及一种图像滤波方法及CT系统，尤其涉及利用并行计算的图像滤波方法及CT系统。
背景技术
图像滤波，即消除图像中的无用噪声，是图像预处理中不可缺少的操作。滤波器是图像处理中最关键的组成部分之一，无论是在图像变换、图像增强和图像恢复中都是相当重要。选择不同的滤波器可以实现不同的处理效果，例如低通滤波器可用于平滑图像，高通滤波器可用于边缘提取等。随着信息技术的发展，我们需要处理的数据量大幅增加，这就对图像滤波的处理速度提出了更高的要求。
除了直接用于图像处理外，滤波器还可以作为一些迭代计算的正则化约束条件。例如，计算机断层成像术（CT）广泛应用于医学影像领域。它是利用计算机技术对被测物体断层扫描图像进行重建获得3维断层图像的扫描方式。该扫描方式是通过单一轴面的射线穿透被测物体，根据被测物体各部分对射线的吸收与透过率不同，由计算机采集透过射线并通过3维重构成像。解析重建和迭代重建是CT图像重建的两种基本方法，其中迭代重建在保证图像质量恒定的前提下能够大幅降低辐射剂量，有利于作为未来的发展方向的低辐射剂量CT。在低辐射剂量CT图像的迭代重建中，每次迭代计算中都要进行滤波正则化计算，由于医疗图像数据量大，滤波的速度将直接影响CT图像重建的速度。因此，对大规模图像滤波处理进行加速十分必要。
另一方面，以图像处理器（GPU）为代表的多核处理器的高性能数值运算能力在近些年发展迅速。NVIDIA公司于2007年正式发布的CUDA（Computer Unified DeviceArchitecture：计算统一设备架构），使用一种类C语言（支持现有C语言基础上，进行了部分扩展），使得开发工作 更加易于掌握。GPU不再局限于图形处理，也可应用于通用数值计算中，特别适用于并行度高数值运算量大的运算。
因此，近年来开始利用以GPU为代表的多核处理器的并行计算来加速大规模图像滤波处理。在传统的并行处理方法中，将图像中的每一个像素点的滤波计算作为基本计算单元交给每一个线程完成。但是，这不但会产生大量的重复计算，造成计算量和计算时间的无谓增加，而且会产生全局存储器的大量重复访问，造成全局存储器的访问次数和访问时间的无谓增加。由此，现有技术在利用并行计算来实现图像滤波处理时存在图像滤波速度低下的技术问题。
发明内容
本发明针对现有技术中的上述技术问题，其目的在于，提供一种基于以图像处理器（GPU）为代表的多核处理器的图像滤波方法及CT系统，在利用并行计算实现快速的大规模图像滤波处理中，通过减少现有并行计算中的重复计算，减少计算量，大幅提高图像滤波的并行计算速度。
为了实现上述目的，本发明提供一种图像滤波方法，利用多核处理器，并行地对图像进行滤波处理，其特征在于，包括以下步骤：计算方向确定步骤，根据所述图像的维数及所述滤波处理的预定的邻域范围，确定所述滤波处理的多个计算方向；各方向滤波计算步骤，按所确定的多个计算方向中的每个计算方向分别进行下述处理：由多个线程中的每个线程分别针对所述图像中的一行像素，以该行的每个像素作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行预定的滤波计算，并将所述滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果，其中所述邻域像素是位于对象像素的所述邻域范围中的像素；以及图像滤波结果获得步骤，将通过所述各方向滤波计算步骤针对所确定的全部计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果。
根据本发明的图像滤波方法，在利用并行计算的图像滤波处理中，能够减少现有并行计算中的重复计算，减少计算量，从而大幅提高图像滤波的并行计算速度，并提高图像滤波处理的速度和实用性。
本发明的图像滤波方法也可以是：在所述各方向滤波计算步骤中，多个线程中的每个线程相互并行地进行滤波计算。
在此，通过多个线程中的每个线程相互并行地进行滤波计算，能够进一步提高图像滤波处理的速度。
本发明的图像滤波方法也可以是：在所述各方向滤波计算步骤中，每个线程分别针对所述图像中的一行像素，从该行的起始像素开始直到该行的最终像素为止依次作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行所述滤波计算，并将所述滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果。
在此，在针对各计算方向的滤波计算中，每个线程将一行像素中的起始像素到最终像素依次作为对象像素。由此，能够在该行的各像素之间充分地相互利用滤波计算结果，大幅提高图像滤波的并行计算速度。
本发明的图像滤波方法也可以是：在所述计算方向确定步骤中确定了相对于所述图像中的像素排列方向倾斜的计算方向的情况下，在所述各方向滤波计算步骤中，按该倾斜的计算方向进行处理时，每个线程所计算的各行像素的像素数量相等。
在此，在按倾斜的计算方向进行滤波计算时，通过合理设定由每个线程所计算的各行像素，使得各行像素的像素数量相等，能够最大程度地在各个线程之间使计算量平均化，而且有利于各个线程并行计算，从而提高了图像滤波的并行计算速度。
本发明的图像滤波方法也可以是：在针对2维图像进行滤波处理、且滤波处理的邻域范围的半径为1的情况下，在所述计算方向确定步骤中，所确定的计算方向的数量为4，在所述图像滤波结果获得步骤中，将通过所述各方向滤波计算步骤针对所确定的4个计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果。
由此，通过多核处理器并行地针对2维图像进行滤波处理，能够减少现有并行计算中的重复计算，减少计算量，从而大幅提高2维图像滤波的并行计算速度，并提高2维图像滤波处理的速度和实用性。
本发明的图像滤波方法也可以是：在针对3维图像进行滤波处理、且滤波处理的邻域范围的半径为1的情况下，在所述计算方向确定步骤中， 所确定的计算方向的数量为13，在所述图像滤波结果获得步骤中，将通过所述各方向滤波计算步骤针对所确定的13个计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果。
由此，通过多核处理器并行地针对3维图像进行滤波处理，能够减少现有并行计算中的重复计算，减少计算量，从而大幅提高3维图像滤波的并行计算速度，并提高3维图像滤波处理的速度和实用性。
本发明的图像滤波方法也可以是：在针对3维图像进行滤波处理、且滤波处理的邻域范围的半径为r的情况下，通过所述计算方向确定步骤确定了多个计算方向之后进行下述处理：（1）从全局存储器向所述多核处理器的共享内存读入构成所述3维图像的多个2维图像中的起始2r+1个2维图像，并对该2r+1个2维图像执行所述各方向滤波计算步骤和所述图像滤波结果获得步骤，（2）将读入的2r+1个2维图像中的前r+1个2维图像的各像素的滤波结果从共享内存写入全局存储器，（3）从全局存储器向共享内存读入构成所述3维图像的多个2维图像中的后续r+1个2维图像，并对由该后续r+1个2维图像和在（2）中未从共享内存写入全局存储器的r个2维图像组成的2r+1个2维图像，执行所述各方向滤波计算步骤和所述图像滤波结果获得步骤，（4）重复处理（2）和处理（3），直到将构成所述3维图像的全部2维图像的图像滤波结果从共享内存写入全局存储器，并将写入全局存储器的所述全部2维图像的图像滤波结果作为所述3维图像的图像滤波结果。
一般而言，全局存储器不具有缓存，读写速度较慢，而多核处理器的共享内存具有缓存这一优势，读写速度快。通过上述处理，充分利用共享内存的缓存机制，减少全局存储器的访问次数，大大减少处理器内存的访问时间，从而提高了图像滤波处理的速度和实用性。
为了实现上述目的，本发明还提供一种CT系统，通过X射线对扫描对象进行扫描，输出所述扫描对象的CT图像，其特征在于，具有：CT扫描仪，通过X射线对扫描对象进行扫描，获得所述扫描对象的投影图像；CT图像重建装置，具有多核处理器，根据所述投影图像重建CT图像，在为了重建CT图像而利用所述多核处理器对图像并行地进行的滤波处理中，根据图像的维数及所述滤波处理的预定的邻域范围，确定所述滤波处理的 多个计算方向，按所确定的多个计算方向中的每个计算方向，由多个线程中的每个线程分别针对所述图像中的一行像素，以该行的每个像素作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行预定的滤波计算，并将所述滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果，将针对所确定的全部计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果，其中所述邻域像素是位于对象像素的所述邻域范围中的像素；CT图像输出装置，输出由所述CT图像重建装置重建的CT图像。
根据本发明的CT系统，在CT图像重建中利用并行计算的图像滤波处理中，能够减少现有并行计算中的重复计算，减少计算量，从而大幅提高图像滤波的并行计算速度，并提高图像滤波处理的速度，使得图像滤波更好地应用在以CT系统为代表的医疗图像处理等大规模数据计算中。
本发明的CT系统也可以是：所述CT图像重建装置所具有的所述多核处理器为GPU即图像处理器。
图像处理器（GPU）的高性能数值运算能力在近些年发展迅速，特别适用于并行度高数值运算量大的运算。通过使用GPU来并行地实现CT图像重建中的滤波处理，能够更大程度地提高图像滤波处理的速度。
本发明不限于上述图像滤波方法及CT系统，也可以通过其他方式实现。例如，通过包含上述图像滤波方法的CT图像重建方法、CT图像生成方法等，也能够实现本发明的目的。另外，通过基于软件模块或硬件架构的图像滤波装置或CT图像重建装置，执行与上述图像滤波方法的各步骤所对应的处理，也能够实现本发明的目的。
附图说明
图1是表示本发明的图像滤波方法的流程图。
图2A、图2B是本发明中2维和3维图像滤波的计算方向示意图。
图3A～图3E是本发明各计算方向上线程的任务分配示意图。
图4是本发明针对3维图像的一个线程块的并行处理流程示意图。
图5A、图5B是本发明3维图像滤波的内存管理示意图。
图6是本发明的CT系统的系统构成图。
具体实施方式
为了便于理解本发明，首先对图像滤波处理进行整体的说明。
图像的滤波方法很多，主要可以分为频率域法和空间域法两大类。其中，空间滤波方法是一类直接的滤波方法，它在处理图像时直接对图像灰度作运算。本发明涉及针对空间滤波的并行计算方法。图像滤波利用计算位置像素（中心像素、对象像素）的像素值与其邻域内各个像素值按照滤波公式进行计算，最后得到该计算位置像素的滤波结果。在此，邻域是本领域的技术用语，又称滤波核，指的是对计算位置像素（中心像素、对象像素）的滤波有贡献（影响）的区域，例如可以具有半径r的邻域范围。其中，将位于计算位置像素（中心像素、对象像素）的邻域范围中的像素称为邻域像素。
在图像滤波处理中，不同滤波器的滤波公式不相同。例如，在CT图像的迭代运算中，正则化计算是将邻域像素与计算位置像素的像素值相减，然后利用差值进行后续运算。在后续的说明中，以邻域像素与计算位置像素的像素值相减这一运算作为基本运算，举例介绍并行处理方法。但本发明同样可以应用于其他滤波核中，不限制于减法。
滤波计算就是计算邻域内各个像素对中心像素影响的总和。邻域内各个像素的影响大小可以由权重值来控制。在一般的图像滤波中，邻域像素对计算位置像素的影响权重，主要取决于各个邻域像素相对于中心像素的位置，或者是两者之间的距离。一般来讲，邻域内像素a对像素b的影响，与像素b对像素a的影响应该是相同的。在传统的并行图像滤波处理中，将图像中的每一个像素点的滤波计算作为基本计算单元交给每一个线程完成。此时，在像素a成为计算位置像素（中心像素）时，计算图像滤波结果内像素b对像素a的影响。之后，在像素b成为计算位置像素（中心像素）时，又计算图像滤波结果内像素a对像素b的影响。如上所述，像素a对像素b的影响与像素b对像素a的影响一般相同，因此在传统的并行图像滤波处理中进行了大量的重复计算，造成了运算量和运算时间的无谓增加。
另一方面，在利用多核处理器的并行计算中，每一个线程处理自己的 任务，各个线程可以并行进行计算。所有的线程，分成一系列线程块，在一个线程块内的线程可以访问它们的共享内存，共享内存内有缓存，读写速度快，但是存储大小有限。因此，一般并行计算中，一个线程块首先从全局存储器（例如显存）中读入一部分输入数据到共享内存中，然后在共享内存中进行计算，最后再将结果从共享内存中拷贝回全局存储器，以利用共享内存读写速度快的优点。在传统的并行图像滤波处理中，将图像中的每一个像素点的滤波计算作为基本计算单元交给每一个线程完成。此时，例如在处理3维图像中的某一片图像（例如第1片图像）中的像素时，需要利用第2片图像中的邻域像素进行滤波计算，因此需要从全局存储器读取该第2片图像。之后，在处理第2片图像中的像素时，也需要从全局存储器读取该第2片图像。之后，在处理第3片图像中的像素时，也需要利用第2片图像中的邻域像素进行滤波计算，因此需要再次从全局存储器读取该第2片图像。也就是说，在邻域半径r为1时，第2片图像需要被读取3次。因此，在传统的并行图像滤波处理中对全局存储器进行了大量的重复访问，造成访问次数和访问时间的无谓增加。
针对以上现状，本发明提供一种图像滤波方法，利用以GPU为代表的多核处理器，并行地对图像进行滤波处理。图1是表示本发明的图像滤波方法的流程图。如图1所示，本发明的图像滤波方法包括计算方向确定步骤S1、各方向滤波计算步骤S2和图像滤波结果获得步骤S3。
在计算方向确定步骤S1中，根据图像的维数及滤波处理的预定的邻域范围，确定滤波处理的多个计算方向。在此，计算方向是将相对于对象像素相互对称的至少两个邻域像素连接而成的方向。在滤波计算中，邻域内的像素都是相对于中心像素（对象像素）中心对称的，例如右上像素（x+1,y+1）与左下像素（x-1，y-1）关于中心像素（x,y）对称，那么右上像素与左上像素都位于通过中心像素的斜45度的直线上，这样的一条直线我们称之为一个计算方向。邻域内的像素根据图像维数及邻域大小的不同可以划分为几个不同计算方向。图2A、图2B是本发明中2维和3维图像滤波的计算方向示意图。如图2A所示，作为2维图像，在邻域半径为1（即邻域范围为3*3像素）的情况下，邻域像素共8个，如箭头所示共有A、B、 C、D这4个计算方向。如图2B所示，作为3维图像，在邻域半径为1（即邻域范围为3*3*3像素）的情况下，邻域像素共26个，如箭头所示共有13个计算方向。
在各方向滤波计算步骤S2中，按所确定的多个计算方向中的每个计算方向分别进行下述处理：由多个线程中的每个线程分别针对图像中的一行像素，以该行的每个像素作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行预定的滤波计算，并将滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果。在此，以滤波计算是基于减法运算为例进行说明，公式例如下式，该滤波计算可用于CT图像迭代重建的正则项计算中。
其中，x,y,z是滤波点（像素）在图像中的坐标值。μ是滤波前像素灰度值。μ'是滤波后像素灰度值。mx、my、mz是滤波核在x、y、z方向上的邻域的半径大小，滤波计算使用的邻域大小为2*m+1，邻域的中心点即为要滤波的点（像素）。当mz=0时，3维滤波变为2维滤波。另外，Nrst是滤波邻域范围内各点的权值，是滤波函数.该函数通常具有一定的对称性。
在图像滤波结果获得步骤S3中，将通过各方向滤波计算步骤S2针对所确定的全部计算方向得到的滤波结果按图像的每个像素累加，从而获得图像的图像滤波结果。
以下具体说明本发明的图像滤波方法的几个具体例。
首先说明针对2维图像进行滤波处理的例子。在针对2维图像进行滤波处理、且滤波处理的邻域范围的半径r为1的情况下，如上所述，在计算方向确定步骤S1中确定的计算方向的数量为4。
然后，读入输入图像（2维图像）到GPU全局存储器中，全局存储器中的输出图像初始化为零。接下来进行任务分配。图3A～图3D是本发明针对2维图像在各计算方向上线程的任务分配示意图。在此，假设2维图像的大小为X*Y（例如512*512），输入图像在XY平面上可以分成几个矩形区域，每个矩形区域对应的2维图像由一个线程块处理。假设一个线程 块中含有32个线程，原图像在XY平面上可以分成256个大小为32*32的矩形区域，每个线程块处理一个矩形区域对应的图像（大小为32*32），矩形区域的大小可以根据GPU中处理器的数量进行调整，以保证并行效率。各个线程块并行执行。在上述4个计算方向中，并行处理其中一个计算方向上的所有滤波计算，内存读写都在共享内存中进行。每个线程计算图像中沿该计算方向的一行，它依次读取这一行中在邻域大小范围内的两个像素点，按照滤波公式进行计算，并将计算结果累加到输出图像相应的两个像素点的结果中。
如图3A所示，在图3A所示的方向（即图2A中的A方向）的计算中，一个线程块中的每个线程如箭头直线所示，并行计算沿y方向的1列，线程1计算第1列，线程2计算第2列，线程32处理第32列，即完成2维图像A方向的计算。其中，每个线程分别针对图像中的一行像素（在此为y方向的一列），例如从该行的起始像素开始直到该行的最终像素为止依次作为对象像素，在对象像素与该对象像素在该计算方向（A方向）上的各邻域像素之间进行滤波计算，并将滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果。例如，可以利用下式计算滤波后像素灰度值μ'：
如图3B所示，在图3B所示的方向（即图2A中的B方向）的计算中，每个线程并行计算沿x方向的1行，并行方法类似于A方向。例如，可以利用下式计算滤波后像素灰度值μ'：
如图3C所示，在图3C所示的方向（即图2A中的C方向）的计算中，每个线程如箭头直线所示，并行计算沿斜率为45度方向的一条直线上的滤波。线程1计算点(32,1)至点(1,32)连线上的各点。线程2计算点（32,2）至点（2,32）连线上的各点，及点（1,2）至点（2,1）连线上的各点。不同于A方向、B方向的是，由于矩形图像的原因，线程2处理的两段直线不连续， 其实可以将第二段直线看做第一段直线在该矩形范围内的延伸。也就是说，在计算方向确定步骤S1中确定了相对于图像中的像素排列方向（在此为A方向、B方向）倾斜的计算方向的情况下，在各方向滤波计算步骤S2中，按该倾斜的计算方向（例如C方向）进行处理时，使得每个线程所计算的各行像素的像素数量相等。例如，可以利用下式计算滤波后像素灰度值μ'：
如图3D所示，图3D所示的方向（即图2A中的D方向）的计算类似于C方向，每个线程并行计算沿斜率为135度方向的一条直线上的滤波。例如，可以利用下式计算滤波后像素灰度值μ'：
在上述A方向至D方向的计算中，线程块中的每个线程优选相互并行地进行滤波计算。例如，在多个线程之间，各线程同时以各自计算的各行像素中从起始像素开始次序相同的像素（例如第n个像素）作为对象像素。这样，在各行像素数量相同的情况下，能够实现各个线程相互并行地进行滤波计算。
经过上述计算，通过各方向滤波计算步骤S2计算上述4个计算方向上的滤波结果。然后，在图像滤波结果获得步骤S3中，将针对上述4个计算方向得到的滤波结果按2维图像的每个像素累加，从而获得2维图像的图像滤波结果。最后，将全部线程块的全部图像滤波结果从共享内存写入全局存储器中，作为2维图像的图像滤波结果。在此，每一个线程块均按上述方法处理2维图像数据中的一个矩形区域，各个线程块并行执行，直到所有图像数据处理完毕。
接下来说明针对3维图像进行滤波处理的例子。在针对3维图像进行滤波处理、且滤波处理的邻域范围的半径r为1的情况下，如上所述，在计算方向确定步骤S1中确定的计算方向的数量为13。在此，继续参照图3A～图3D并追加图3E、图4、图5A、图5B进行说明。在此，假设3维图像 的大小为X*Y*Z（例如512*512*512），可以看成是Z个X*Y的2维图像片。
然后，读入输入图像（3维图像）到GPU全局存储器中，全局存储器中的输出图像初始化为零。接下来进行任务分配。输入图像在XY平面上可以分成几个矩形区域，每个矩形区域对应的3维图像体由一个线程块处理。假设一个线程块中含有32个线程，原图像在XY平面上可以分成256个大小为32*32的矩形区域，每个线程块处理一个矩形区域对应的图像体（大小为32*32*Z）。与2维图像时相同，矩形区域的大小可以根据GPU中处理器的数量进行调整，以保证并行效率。各个线程块并行执行。
图4是本发明针对3维图像的一个线程块的并行处理流程示意图。图5A、图5B是本发明3维图像滤波的内存管理示意图。在图5A、图5B中，下角标U表示构成3维图像的多个2维图像中的上一片图像的贡献，M表示中间图像的贡献，L表示下一片图像的贡献，输出为U+M+L。如图4、图5A、图5B所示，在步骤S31中，从全局存储器向多核处理器的共享内存（共享存储）读入构成3维图像的多个2维图像中的起始2r+1个2维图像，并对该2r+1个2维图像执行各方向滤波计算步骤S2和图像滤波结果获得步骤S3。在此r为1，因此从全局存储中拷贝要处理的图像体的前3片图像（第一片至第三片）到共享内存。此部分图像为该线程块要处理的输入图像，并在共享内存中分配一块相同大小的输出图像，初始值为零。在步骤S32中，在上述13个计算方向中，并行处理其中一个计算方向上的所有滤波计算，内存读写都在共享内存中进行。每个线程计算图像中沿该计算方向的一行，它依次读取这一行中在邻域大小范围内的两个像素点，按照滤波公式进行计算，并将计算结果累加到输出图像相应的两个像素点的结果中。
其中，在上述13个计算方向之中，2维计算方向上的计算与上述2维图像时相同，例如在图3A所示的A方向的计算中，线程1的具体操作为，首先读取第1列的前两个像素，即x坐标为1和2的两个像素，将两点像素值相减，带入上述滤波公式，计算将结果累加至μ'1,1,1，计算将结果累加至μ'1,2,1，然后读入x坐标为3的点，将μ1,2,1与μ1,3,1相减，计算将结果累加至μ'1,2,1，计算
将结果累加至μ'1,3,1。以此类推直至完成整列计算。其他2维方向上的计算在此不做赘述。
对于3维计算方向上的并行方法，可以由2维计算方向上的并行方法类推得到，思路一致，只是计算直线方向由2维变为3维。在此，举出3维方向之中的一个计算方向为例进行说明。图3E是本发明针对3维图像在一个3维计算方向（Z方向）上线程的任务分配示意图。如图3E所示，每个线程计算沿箭头直线所示的一行，例如从该行的起始像素开始直到该行的最终像素为止依次作为对象像素，在对象像素与该对象像素在Z方向上的各邻域像素之间进行滤波计算（例如像素值的相减），并将滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果。
在步骤S33中，判断是否完成了所有计算方向的计算。如果判断结果为否，则返回步骤S32。如果判断结果为是，即上述13个计算方向的计算都完成之后，共享内存中2r+1（在此为3）片图像之间的对彼此的滤波贡献计算完成。其中，前r+1片（在此为第1片和第2片）的滤波计算已全部完成，而剩下的r片（在此为第3片）图像还需计算第4片图像对其的贡献，因此第3片图像的滤波计算还未完成。此时，在步骤S34中，将读入的2r+1个（3片）2维图像中的前r+1个（第1片和第2片）2维图像的各像素的滤波结果（1M+L、2U+M+L）从共享内存写入全局存储器，而第3片图像的中间结果（3U+M）仍然保留在共享内存中以用于后续计算。
接着，在步骤S35中，从全局存储器向共享内存读入构成3维图像的多个2维图像中的后续r+1个（在此为第4片和第5片）2维图像，并将这2片图像的输出图像初始化为0。然后，对由该后续r+1个2维图像和之前未从共享内存写入全局存储器的r个2维图像组成的2r+1个（在此为第3片至第5片）2维图像，执行各方向滤波计算步骤S2和图像滤波结果获得步骤S3。具体与步骤S32～S34相同，在此不再赘述。在步骤S36中，判断针对3维图像体是否全部处理完成。如果判断结果为否，则返回步骤S32，重复上述处理直到z片2维图像处理完成，即将构成3维图像的全部2维图像的图像滤波结果从共享内存写入全局存储器。
在此，每一个线程块均按上述方法处理3维图像数据中的一个3维图像体，各个线程块并行执行，直到所有图像数据处理完毕。最后，结束计 算，将写入全局存储器的全部图像滤波结果作为3维图像的图像滤波结果。
接下来，参照图6说明本发明的CT系统。图6是本发明的CT系统的系统构成图。如图6所示，CT系统通过X射线对扫描对象进行扫描，输出扫描对象的CT图像，具有CT扫描仪、CT图像重建装置和CT图像输出装置。
CT扫描仪通过X射线对扫描对象进行扫描，获得扫描对象的投影图像。具体而言，通过单一轴面的射线穿透扫描对象，根据扫描对象各部分对X射线的吸收与透过率不同，获得扫描对象的投影图像。
CT图像重建装置具有多核处理器，根据投影图像重建CT图像。为了重建CT图像，CT图像重建装置利用多核处理器对图像并行地以如上所述的方式进行滤波处理。其中，CT图像重建装置所具有的多核处理器可以是GPU。
CT图像输出装置输出由CT图像重建装置重建的CT图像。在此，CT图像输出装置例如是CT图像显示装置，由液晶显示器或触摸屏等显示设备实现，显示由CT图像重建装置重建的CT图像。另外，CT图像输出装置也可以是打印设备等，打印由CT图像重建装置重建的CT图像。当然，CT图像输出装置只要能够输出CT图像即可，也可以是通用的I/O接口。
根据本发明的上述实施方式，将每一个像素点的滤波计算拆分成几部分，所有线程每次同时完成其中的一部分，直至所有部分完成，累加便得到最终结果。即，本发明将图像滤波计算按照计算方向进行划分，每次所有的线程共同并行完成图像上所有像素在一个计算方向上的计算，直至所有计算方向都完成。前面已经讲过，邻域内像素a对像素b的影响，与像素b对像素a的影响应该是相同的。利用这一点，按照方向计算，像素a与b的计算结果既可以用于像素a的滤波，也可以用于像素b的滤波，因此就避免了传统并行方法中的重复计算。
而且，根据本发明的上述实施方式，对全局存储器无重复读写访问。计算中的所有读写都在共享内存中进行，充分利用了缓存。本发明方法对全局存储器的访问只有读入一次输入图像数据和写出一次输出图像数据。 而传统的并行方法中，每个线程负责处理一个点的所有滤波计算，处理第1片图像时需要第2片图像，处理第3片图像时也需要第2片图像，第2片图像需要被读取3次。因此，本发明中全局存储器中图像读取次数大幅减少，处理速度大幅提升。
在上述实施方式中，以邻域像素与对象像素的像素值相减这一运算作为基本运算说明了图像滤波计算。但是，本发明能够使用的图像滤波计算不限于此，只要是满足互为邻域像素的像素a对像素b的影响与像素b对像素a的影响相同这一原则的图像滤波计算，本发明的图像滤波方法都适用。本发明主要针对的是并行处理方法，而不是具体的图像滤波计算。
在上述各实施方式中，以GPU处理器作为多核处理器的例子。但不限于此，只要是能够并行多线程处理的多核处理器，都能够用于实现本发明。

Claims (8)

1.一种图像滤波方法，利用多核处理器，并行地对图像进行滤波处理，其特征在于，包括以下步骤：
计算方向确定步骤，根据所述图像的维数及所述滤波处理的预定的邻域范围，确定所述滤波处理的多个计算方向；
各方向滤波计算步骤，按所确定的多个计算方向中的每个计算方向分别进行下述处理：由多个线程中的每个线程分别针对所述图像中的一行像素，以该行的每个像素作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行预定的滤波计算，并将所述滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果，其中所述邻域像素是位于对象像素的所述邻域范围中的像素；以及
图像滤波结果获得步骤，将通过所述各方向滤波计算步骤针对所确定的全部计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果，
在针对3维图像进行滤波处理、且滤波处理的邻域范围的半径为r的情况下，通过所述计算方向确定步骤确定了多个计算方向之后进行下述处理：
(1)从全局存储器向所述多核处理器的共享内存读入构成所述3维图像的多个2维图像中的起始2r+1个2维图像，并对该2r+1个2维图像执行所述各方向滤波计算步骤和所述图像滤波结果获得步骤，
(2)将读入的2r+1个2维图像中的前r+1个2维图像的各像素的滤波结果从共享内存写入全局存储器，
(3)从全局存储器向共享内存读入构成所述3维图像的多个2维图像中的后续r+1个2维图像，并对由该后续r+1个2维图像和在(2)中未从共享内存写入全局存储器的r个2维图像组成的2r+1个2维图像，执行所述各方向滤波计算步骤和所述图像滤波结果获得步骤，
(4)重复处理(2)和处理(3)，直到将构成所述3维图像的全部2维图像的图像滤波结果从共享内存写入全局存储器，并将写入全局存储器的所述全部2维图像的图像滤波结果作为所述3维图像的图像滤波结果。
2.如权利要求1所述的图像滤波方法，其特征在于，
在所述各方向滤波计算步骤中，多个线程中的每个线程相互并行地进行滤波计算。
3.如权利要求1所述的图像滤波方法，其特征在于，
在所述各方向滤波计算步骤中，每个线程分别针对所述图像中的一行像素，从该行的起始像素开始直到该行的最终像素为止依次作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行所述滤波计算，并将所述滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果。
4.如权利要求1所述的图像滤波方法，其特征在于，
在所述计算方向确定步骤中确定了相对于所述图像中的像素排列方向倾斜的计算方向的情况下，在所述各方向滤波计算步骤中，按该倾斜的计算方向进行处理时，每个线程所计算的各行像素的像素数量相等。
5.如权利要求1～4中任一项所述的图像滤波方法，其特征在于，
在针对2维图像进行滤波处理、且滤波处理的邻域范围的半径为1的情况下，
在所述计算方向确定步骤中，所确定的计算方向的数量为4，
在所述图像滤波结果获得步骤中，将通过所述各方向滤波计算步骤针对所确定的4个计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果。
6.如权利要求1～4中任一项所述的图像滤波方法，其特征在于，
在针对3维图像进行滤波处理、且滤波处理的邻域范围的半径为1的情况下，
在所述计算方向确定步骤中，所确定的计算方向的数量为13，
在所述图像滤波结果获得步骤中，将通过所述各方向滤波计算步骤针对所确定的13个计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果。
7.一种CT系统，通过X射线对扫描对象进行扫描，输出所述扫描对象的CT图像，其特征在于，具有：
CT扫描仪，通过X射线对扫描对象进行扫描，获得所述扫描对象的投影图像；
CT图像重建装置，具有多核处理器，根据所述投影图像重建CT图像，在为了重建CT图像而利用所述多核处理器对图像并行地进行的滤波处理中，根据图像的维数及所述滤波处理的预定的邻域范围，确定所述滤波处理的多个计算方向，按所确定的多个计算方向中的每个计算方向，由多个线程中的每个线程分别针对所述图像中的一行像素，以该行的每个像素作为对象像素，在对象像素与该对象像素在该计算方向上的各邻域像素之间进行预定的滤波计算，并将所述滤波计算的结果分别以累加的方式保存为该对象像素和各邻域像素的滤波结果，将针对所确定的全部计算方向得到的滤波结果按所述图像的每个像素累加，从而获得所述图像的图像滤波结果，其中所述邻域像素是位于对象像素的所述邻域范围中的像素；
CT图像输出装置，输出由所述CT图像重建装置重建的CT图像，
在针对3维图像进行滤波处理、且滤波处理的邻域范围的半径为r的情况下，所述CT图像重建装置确定了多个计算方向之后进行下述处理：
(1)从全局存储器向所述多核处理器的共享内存读入构成所述3维图像的多个2维图像中的起始2r+1个2维图像，并对该2r+1个2维图像进行所述滤波计算并获得所述图像滤波结果，
(2)将读入的2r+1个2维图像中的前r+1个2维图像的各像素的滤波结果从共享内存写入全局存储器，
(3)从全局存储器向共享内存读入构成所述3维图像的多个2维图像中的后续r+1个2维图像，并对由该后续r+1个2维图像和在(2)中未从共享内存写入全局存储器的r个2维图像组成的2r+1个2维图像，进行所述滤波计算并获得所述图像滤波结果，
(4)重复处理(2)和处理(3)，直到将构成所述3维图像的全部2维图像的图像滤波结果从共享内存写入全局存储器，并将写入全局存储器的所述全部2维图像的图像滤波结果作为所述3维图像的图像滤波结果。
8.如权利要求7所述的CT系统，其特征在于，
所述CT图像重建装置所具有的所述多核处理器为GPU即图像处理器。
CN201410118326.4A 2014-03-27 2014-03-27 图像滤波方法及ct 系统 Expired - Fee Related CN104952043B (zh)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN201410118326.4A CN104952043B (zh) 2014-03-27 2014-03-27 图像滤波方法及ct 系统
US15/129,260 US9984433B2 (en) 2014-03-27 2015-03-26 Image filtering method and CT system
JP2016555324A JP6259531B2 (ja) 2014-03-27 2015-03-26 画像フィルタリング方法及びｃｔシステム
PCT/CN2015/075123 WO2015144069A1 (zh) 2014-03-27 2015-03-26 图像滤波方法及ct系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410118326.4A CN104952043B (zh) 2014-03-27 2014-03-27 图像滤波方法及ct 系统

Publications (2)

Publication Number Publication Date
CN104952043A CN104952043A (zh) 2015-09-30
CN104952043B true CN104952043B (zh) 2017-10-24

Family

ID=54166677

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410118326.4A Expired - Fee Related CN104952043B (zh) 2014-03-27 2014-03-27 图像滤波方法及ct 系统

Country Status (4)

Country Link
US (1) US9984433B2 (zh)
JP (1) JP6259531B2 (zh)
CN (1) CN104952043B (zh)
WO (1) WO2015144069A1 (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106097410B (zh) * 2016-06-21 2018-12-25 沈阳东软医疗系统有限公司 一种扫描数据处理方法和装置
EP3409642A1 (en) 2017-06-01 2018-12-05 Paqell B.V. A process to convert bisulphide to elemental sulphur
CN109521456B (zh) * 2018-09-25 2022-10-21 中国辐射防护研究院 一种基于正则化最小二乘法的γ辐射源项反演方法及系统
CN109727206B (zh) * 2018-12-05 2023-05-02 安徽紫薇帝星数字科技有限公司 一种二值图像中值滤波的快速计算方法以及其实现方法
CN110428483B (zh) * 2019-07-31 2021-09-03 北京长木谷医疗科技有限公司 一种图像处理方法及计算设备
CN115643403A (zh) * 2022-11-11 2023-01-24 上海哔哩哔哩科技有限公司 Av1的滤波方法及装置
CN117094917B (zh) * 2023-10-20 2024-02-06 高州市人民医院 一种心血管3d打印数据处理方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101596113A (zh) * 2008-06-06 2009-12-09 中国科学院过程工程研究所 一种ct并行重建系统及成像方法
CN101630404A (zh) * 2009-07-30 2010-01-20 上海交通大学 弱小目标二维图片噪声滤除中的匹配滤波方法
CN102005037A (zh) * 2010-11-12 2011-04-06 湖南大学 结合多尺度双边滤波与方向滤波的多模图像融合方法
US8391649B2 (en) * 2003-08-01 2013-03-05 Texas Instruments Incorporated Image filter method
CN103208097A (zh) * 2013-01-29 2013-07-17 南京理工大学 图像多方向形态结构分组的主分量分析协同滤波方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5602934A (en) * 1993-09-08 1997-02-11 The Regents Of The University Of California Adaptive digital image signal filtering
JP3860548B2 (ja) * 2003-02-14 2006-12-20 誠 小川 画像処理装置及び画像処理方法
US7916144B2 (en) * 2005-07-13 2011-03-29 Siemens Medical Solutions Usa, Inc. High speed image reconstruction for k-space trajectory data using graphic processing unit (GPU)
US7970191B2 (en) * 2006-09-28 2011-06-28 Siemens Medical Solutions Usa, Inc. System and method for simultaneously subsampling fluoroscopic images and enhancing guidewire visibility
JP2009042909A (ja) * 2007-08-07 2009-02-26 Sanyo Electric Co Ltd 特徴点検出装置、およびそれを搭載した動画像処理装置
US8660330B2 (en) * 2008-06-27 2014-02-25 Wolfram Jarisch High efficiency computed tomography with optimized recursions
US8731269B2 (en) * 2011-10-19 2014-05-20 Kabushiki Kaisha Toshiba Method and system for substantially reducing artifacts in circular cone beam computer tomography (CT)
US9020230B2 (en) * 2011-11-02 2015-04-28 General Electric Company Method and apparatus for iterative reconstruction
CN102609978B (zh) * 2012-01-13 2014-01-22 中国人民解放军信息工程大学 基于cuda架构的gpu加速锥束ct图像重建的方法
JP6118698B2 (ja) * 2013-09-30 2017-04-19 株式会社Ihi 画像解析装置及びプログラム
CN103593850A (zh) * 2013-11-26 2014-02-19 北京航空航天大学深圳研究院 一种在cuda平台上基于递归高斯滤波的sift并行化系统及方法
JP5758533B1 (ja) * 2014-09-05 2015-08-05 楽天株式会社 画像処理装置、画像処理方法、ならびに、プログラム

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8391649B2 (en) * 2003-08-01 2013-03-05 Texas Instruments Incorporated Image filter method
CN101596113A (zh) * 2008-06-06 2009-12-09 中国科学院过程工程研究所 一种ct并行重建系统及成像方法
CN101630404A (zh) * 2009-07-30 2010-01-20 上海交通大学 弱小目标二维图片噪声滤除中的匹配滤波方法
CN102005037A (zh) * 2010-11-12 2011-04-06 湖南大学 结合多尺度双边滤波与方向滤波的多模图像融合方法
CN103208097A (zh) * 2013-01-29 2013-07-17 南京理工大学 图像多方向形态结构分组的主分量分析协同滤波方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于多线程的三维医学影像的重建;高齐新等;《第十四届全国图象图形学学术会议》;20080501;第542-546页 *

Also Published As

Publication number Publication date
WO2015144069A1 (zh) 2015-10-01
JP6259531B2 (ja) 2018-01-10
US9984433B2 (en) 2018-05-29
US20170109860A1 (en) 2017-04-20
JP2017509073A (ja) 2017-03-30
CN104952043A (zh) 2015-09-30

Similar Documents

Publication Publication Date Title
CN104952043B (zh) 图像滤波方法及ct 系统
Coric et al. Parallel-beam backprojection: an FPGA implementation optimized for medical imaging
Xu et al. High-performance iterative electron tomography reconstruction with long-object compensation using graphics processing units (GPUs)
Sabne et al. Model-based iterative CT image reconstruction on GPUs
US20080043024A1 (en) Method for reconstructing an object subject to a cone beam using a graphic processor unit (gpu)
DE102021121187A1 (de) EINRICHTUNG UND VERFAHREN ZUR EFFIZIENTEN GRAFIKVERARBEITUNG EINSCHLIEßLICH STRAHLVERFOLGUNG
CN102831577B (zh) 基于gpu的二维地震图像的快速缩放方法
DE102020129969A1 (de) Verbesserungen der verarbeitung und des caching von graphikverarbeitungseinheiten
CN115330986B (zh) 一种分块渲染模式图形处理方法及系统
Zhao et al. GPU based iterative cone-beam CT reconstruction using empty space skipping technique
Zhao et al. GPU-based 3D cone-beam CT image reconstruction for large data volume
DE102023105565A1 (de) VERFAHREN UND VORRICHTUNG FÜR EFFIZIENTEN ZUGRIFF AUF MEHRDIMENSIONALE DATENSTRUKTUREN UND/ODER ANDERE GROßE DATENBLÖCKE
Buurlage et al. A geometric partitioning method for distributed tomographic reconstruction
Sui et al. Simultaneous image reconstruction and lesion segmentation in accelerated MRI using multitasking learning
Johnston et al. GPU-based iterative reconstruction with total variation minimization for micro-CT
Hansen et al. Synthetic aperture beamformation using the GPU
DE112020000902T5 (de) Datenvorabruf für die grafikdatenverarbeitung
Philip et al. Hybrid CPU-GPU solver for gradient domain processing of massive images
Neophytou et al. GPU-Accelerated Volume Splatting With Elliptical RBFs.
Malkin et al. CUDA-Optimized real-time rendering of a Foveated Visual System
DE102021132981A1 (de) Dreidimensionale tomographie-rekonstruktionspipeline
DE102021110598A1 (de) Schätzen von produktintegralen unter verwendung einer zusammensetzung von verzerrungen
WO2024078049A1 (en) System and method for near real-time and unsupervised coordinate projection network for computed tomography images reconstruction
Christiansen et al. Sinogram Interpolation Inspired by Single-Image Super Resolution
Shi et al. Algorithm of ray casting volume rendering based on CUDA

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C41 Transfer of patent application or patent right or utility model
TA01 Transfer of patent application right

Effective date of registration: 20160624

Address after: Tokyo, Japan, Japan

Applicant after: Hitachi Ltd.

Address before: Tokyo, Japan, Japan

Applicant before: Hitachi Medical Corporation

GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171024

Termination date: 20180327