CN113284205B - Ct迭代重建方法及装置 - Google Patents

Ct迭代重建方法及装置 Download PDF

Info

Publication number
CN113284205B
CN113284205B CN202110443876.3A CN202110443876A CN113284205B CN 113284205 B CN113284205 B CN 113284205B CN 202110443876 A CN202110443876 A CN 202110443876A CN 113284205 B CN113284205 B CN 113284205B
Authority
CN
China
Prior art keywords
integral
value
image
weight matrix
image sequence
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.)
Active
Application number
CN202110443876.3A
Other languages
English (en)
Other versions
CN113284205A (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.)
Neusoft Medical Systems Co Ltd
Original Assignee
Shenyang Advanced Medical Equipment Technology Incubation Center Co 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 Shenyang Advanced Medical Equipment Technology Incubation Center Co ltd filed Critical Shenyang Advanced Medical Equipment Technology Incubation Center Co ltd
Priority to CN202110443876.3A priority Critical patent/CN113284205B/zh
Publication of CN113284205A publication Critical patent/CN113284205A/zh
Application granted granted Critical
Publication of CN113284205B publication Critical patent/CN113284205B/zh
Active 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
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • 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/10016—Video; Image sequence
    • 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]
    • 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/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20224—Image subtraction

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明实施例提供一种CT迭代重建方法及装置。本发明实施例通过将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，预先离线确定所述XY平面权重矩阵并存储，在每次迭代重建过程中包括：读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，根据所述投影数据，确定本次迭代后的输出图像序列，若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列，否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代，在迭代过程中用读取XY平面权重矩阵的操作代替该矩阵的生成操作，减少了迭代耗时，提高了重建速度。

Description

CT迭代重建方法及装置
技术领域
本发明涉及医学图像处理技术领域，尤其涉及一种CT迭代重建方法及装置。
背景技术
CT(Computed Tomography，电子计算机断层扫描)成像是目前应用非常广泛的一种影像学技术。CT解析重建算法和CT迭代重建算法是两种常用的CT图像重建方式。
与CT解析重建算法相比，CT迭代重建算法的优势在于可以结合先验信息对CT成像***建模，适用于信息量不足、噪声较大的情况。CT迭代重建算法的劣势是计算耗时，速度慢。
发明内容
为克服相关技术中存在的问题，本发明提供了一种CT迭代重建方法及装置、CT设备及CT***，提高重建速度。
根据本发明实施例的第一方面，提供一种CT迭代重建方法，将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，预先离线确定所述XY平面权重矩阵并存储；
在每次迭代重建过程中包括：
读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
根据所述投影数据，确定本次迭代后的输出图像序列；
若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
根据本发明实施例的第二方面，提供一种CT迭代重建装置，将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，所述装置中存储预先离线确定的所述XY平面权重矩阵，所述装置包括：
投影模块，用于在每次迭代重建过程中读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
确定模块，用于根据所述投影数据，确定本次迭代后的输出图像序列；
输出模块，用于若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
本发明实施例提供的技术方案可以包括以下有益效果：
本发明实施例，通过将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，预先离线确定所述XY平面权重矩阵并存储，在每次迭代重建过程中包括：读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，根据所述投影数据，确定本次迭代后的输出图像序列，若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列，否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代，通过在迭代过程中用读取XY平面权重矩阵的操作代替原来的生成矩阵的操作，减少了迭代耗时，提高了重建速度。
应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本说明书。
附图说明
此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本说明书的实施例，并与说明书一起用于解释本说明书的原理。
图1是本发明实施例提供的CT迭代重建方法的流程示例图。
图2是本发明实施例提供的内存布局重排示意图。
图3是本发明实施例提供的CT迭代重建装置的功能方块图。
图4是本发明实施例提供的电子设备的一个硬件结构图。
具体实施方式
这里将详细地对示例性实施例进行说明，其示例表示在附图中。下面的描述涉及附图时，除非另有表示，不同附图中的相同数字表示相同或相似的要素。以下示例性实施例中所描述的实施方式并不代表与本发明相一致的所有实施方式。相反，它们仅是与如所附权利要求书中所详述的、本发明实施例的一些方面相一致的装置和方法的例子。
在本发明实施例使用的术语是仅仅出于描述特定本发明实施例的目的，而非旨在限制本发明实施例。在本发明实施例和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。还应当理解，本文中使用的术语“和/或”是指并包含一个或多个相关联的列出项目的任何或所有可能组合。
应当理解，尽管在本发明实施例可能采用术语第一、第二、第三等来描述各种信息，但这些信息不应限于这些术语。这些术语仅用来将同一类型的信息彼此区分开。例如，在不脱离本发明实施例范围的情况下，第一信息也可以被称为第二信息，类似地，第二信息也可以被称为第一信息。取决于语境，如在此所使用的词语“如果”可以被解释成为“在……时”或“当……时”或“响应于确定”。
CT迭代重建过程中每次迭代都需要重新计算前向投影，尤其是3D(3维)前向投影，这个过程非常耗时，导致重建速度慢。
数学上可以将计算3D前向投影表示为三维矩阵和向量的乘法，矩阵的元素表示3D前向投影的权重，向量则是3D前向投影需要的体数据。相关技术中，在每次迭代过程中，都要重新计算一遍三维矩阵，然后用该三维矩阵与体数据的向量相乘，获得投影数据。
下面通过实施例对CT迭代重建方法进行详细说明。
本实施例中，在进行迭代重建之前，将三维前向投影所需的三维权重矩阵分解为XY平面权重矩阵和Z向权重向量，预先离线确定所述XY平面权重矩阵并存储
图1是本发明实施例提供的CT迭代重建方法的流程示例图。如图1所示，本实施例的CT迭代重建方法中，在每次迭代重建过程中可以包括：
S101，读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据。
S102，根据所述投影数据，确定本次迭代后的输出图像序列。
S103，若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
其中，首次迭代过程中的输入图像序列可以是根据非迭代重建算法对CT扫描所得的生数据进行重建得到的图像序列。
3D前向投影的权重可以分解为XY平面权重和Z向权重。本发明实施例中，在进行迭代重建之前，预先离线计算XY平面权重矩阵并存储(例如存储至硬盘)，在每次迭代重建过程中，直接读取已存储的XY平面权重矩阵，并可以利用CUDA(Compute Unified DeviceArchitecture，统一计算设备架构)实时计算Z方向的权重，从而完成3D前向投影计算。
本文中，用AM×N表示XY平面权重矩阵。离线计算AM×N不需要过于考虑计算耗时，因此可以选取建模更加准确的前向投影模型计算AM×N。权重矩阵的行数为M＝ViewNumPerRotation×ChannelNum，列数为N＝Volume_X×Volume_Y,因为矩阵AM×N比较稀疏，因此可以以稀疏矩阵格式储存到硬盘，这样可以减少硬盘空间占用。矩阵AM×N中的元素ai,j表示第j个像素对第i根射线的权重，可以选取适当的前向投影模型，如距离驱动，计算得到。
本实施例中，还可以离线计算并存储焦点沿着射线i到像素j的距离SP(i,j)，SP(i,j)仅与i,j相关。离线计算SP(i,j)并保存至硬盘可以进一步减少迭代重建过程的耗时，提高重建速度。
在应用中，可以离线计算一次XY平面权重矩阵，将计算得到的XY平面权重矩阵存储在CT***中。对任一个受检者，在CT迭代重建过程中，都将存储的XY平面权重矩阵作为已知数据，只需要读取操作，而不需要重新计算。这样，就减少了每次迭代过程的耗时，从而减少了整个重建过程的用时，提高了重建速度。
需要说明的是，在一个CT***中，XY平面权重矩阵对所有受检者都是相同的。因此，XY平面权重矩阵只需要离线计算一次，CT迭代重建所用的时间不包含离线计算XY平面权重矩阵的时间。
可见，本实施例通过预先离线计算并存储XY平面权重矩阵，避免了在每次迭代过程中都重新计算3维权重矩阵，减少了迭代过程中三维前向投影的耗时，从而提高了重建速度。
本实施例中，根据所述投影数据，确定本次迭代后的输出图像序列，可以包括：
计算投影数据与对应的CT扫描生数据的差，得到误差投影；
对误差投影进行反投影，得到误差图像序列；
用本次迭代的输入图像序列的图像减去误差图像序列中的相应图像，得到输出图像序列。
在一个示例中，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，可以包括：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
其中，对输入图像序列沿Z向进行图像积分的公式可以如公式(1)所示：
对
其中，InteVolume(x,y,z)表示积分图像的像素值，(x,y,z)表示第z张图像、第y行、第x列。
Volume(x,y,z)是输入给3D前向投影的图像序列，即输入图像序列，(x,y,z)表示第z张图像，第y行，第x列。
在一个示例中，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，可以包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
其中，最大值对应的图像积分值可以通过在第一图像积分值和第二图像积分值之间进行线性插值得到；最小值对应的图像积分值可以通过在第三图像积分值和第四图像积分值之间进行线性插值得到。
对于第s排检测器，第i跟过中心平面的射线，第j个像素(射线i穿过像素j)，假设Z向一维距离驱动前向投影的积分范围表示为[zmin(i,j,s),zmax(i,j,s)]，其中zmin(i,j,s)是该积分范围的最小值，zmax(i,j,s)是该积分范围的最大值，则zmin(i,j,s)和zmax(i,j,s)的计算公式分别如公式(2)和公式(3)所示：
zmin(i,j,s)＝SD/SP(i,j)*(s-SliceNum/2)*SliceThick/Interval (2)
zmax(i,j,s)＝SD/SP(i,j)*(s-SliceNum/2+1)*SliceThick/Interval (3)
其中，SD为焦点到检测器的距离，SliceThick为CT检测器排的厚度，SliceNum为CT检测器排数，Interval为输入给3D前向投影的图像序列的间隔。
通过线性插值可以得到InteVolume(x,y,z)在zmin(i,j,s)和zmax(i,j,s)的值：InteVolume(x,y,zmin(i,j,s))(即最小值对应的图像积分值)和InteVolume(x,y,zmax(i,j,s))(即最大值对应的图像积分值)。则Z向1D距离驱动前向投影积分等于InteVolume(x,y,zmax(i,j,s))-InteVolume(x,y,zmin(i,j,s))。
这里举例说明如何通过线性插值获得最大值对应的图像积分值和最小值对应的图像积分值。
假设积分图像序列共有10幅图像F1～F10，该10幅图像对应的Z坐标值依次为1、2、3、……10；zmin(i,j,s)＝2.5，zmax(i,j,s)＝7.2。对于zmax(i,j,s)＝7.2，图像F1～F10对应的Z坐标值中与zmax(i,j,s)＝7.2相邻的值为z＝7(对应图像F7)和z＝8(对应图像F8)，则找到元素ai,j在图像F7上对应的图像积分值V7和元素ai,j在图像F8上对应的图像积分值V8，则zmax(i,j,s)对应的图像积分值可以通过在V7和V8之间进行线性插值得到。同理，对于zmin(i,j,s)＝2.5，图像F1～F10对应的Z坐标值中与zmin(i,j,s)＝2.5相邻的值为z＝2(对应图像F2)和z＝3(对应图像F3)，则找到元素ai,j在图像F2上对应的图像积分值V2和元素ai,j在图像F3上对应的图像积分值V3，则zmin(i,j,s)对应的图像积分值可以通过在V2和V3之间进行线性插值得到。
假设Pv,c,s表示第v个投影，第c个检测器通道，第s排检测器的3D前向投影值，Pv,c,s对应的投影值表示为公式(4)：
其中，ρ(s)表示锥束射线与锥束中心平面的夹角，ρ(s)与检测器排号s有关；ai,j为XY权重矩阵AM×N中的元素，表示第j个像素对第i根射线的权重。
Pv,c,s对应的投影值累加得到最终投影值，累加公式如下面的公式(5)所示：
其中，下标v和c与下标i的关系为：i＝v*ChannelNum+c。如此遍历完AM×N的所有元素后便完成了一圈3D前向投影计算。
因为Pv,c,s位于全局显存中，用于保存3D前向投影的计算结果，公式(5)需要频繁读写全局显存，所以比较耗时。
在一个示例中，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，可以包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
本实施例通过对Z向1D距离驱动前向投影进行分块计算，提高了缓存命中率，减少了对全局显存的访问，节约了处理时间。
本实施例中对Z向1D距离驱动前向投影采用分块计算的方式，减少对全局显存的访问。本实施例的分块计算过程可以是：
对于Pv,c,s，0≤s≤SliceNum-1，分配TileSize个线程计算，执行核函数的Block的x分量设置为TileSize，则计算分成了TileNum块。这样，就将TileNum个投影值的所有投影分量分别累加至TileNum个寄存器变量PTemp0～PTempTileNum-1，累加完毕后再进一步将PTemp0～PTempTileNum-1的值赋值给Pv,c,s对应的全局显存。
在一个示例中，所述积分图像序列以Z向优先存储的方式存储在内存中。
本实施例是对内存布局的重排，通过内存布局重排，可以提高CUDA合并访问显存的效率。图2是本发明实施例提供的内存布局重排示意图。图2中，左侧图所示为积分图像序列InteVolume(x,y,z)，相关技术中对该图像序列进行行优先存储[Volume_X,Volume_Y,Volume_Z+1]，如图2中右侧的上图所示。本实施例中改为Z向优先存储[Volume_Z+1,Volume_X,Volume_Y]，如图2中右侧的下图所示。
本发明实施例提供的CT迭代重建方法，通过将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，预先离线确定所述XY平面权重矩阵并存储，在每次迭代重建过程中包括：读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，根据所述投影数据，确定本次迭代后的输出图像序列，若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列，否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代，通过在迭代过程中用读取XY平面权重矩阵的操作代替原来的生成矩阵的操作，减少了迭代耗时，提高了重建速度。
基于上述的方法实施例，本发明实施例还提供了相应的装置、设备及存储介质实施例。
图3是本发明实施例提供的CT迭代重建装置的功能方块图。将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，所述装置中存储预先离线确定的所述XY平面权重矩阵。如图3所示，本实施例中，CT迭代重建装置可以包括：
投影模块310，用于在每次迭代重建过程中读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
确定模块320，用于根据所述投影数据，确定本次迭代后的输出图像序列；
输出模块330，用于若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
在一个示例中，所述投影模块310可以具体用于：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
在一个示例中，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
在一个示例中，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
在一个示例中，所述积分图像序列以Z向优先存储的方式存储在内存中。
本发明实施例还提供了一种CT设备。图4是本发明实施例提供的电子设备的一个硬件结构图。该电子设备例如可以是CT***中的控制台设备。将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，所述电子设备中存储有预先离线确定的所述XY平面权重矩阵。如图4所示，电子设备包括：内部总线401，以及通过内部总线连接的存储器402，处理器403和外部接口404，其中：
所述存储器402，用于存储CT迭代重建逻辑对应的机器可读指令；
所述处理器403，用于读取存储器402上的机器可读指令，并执行所述指令以实现如下操作：
在每次迭代重建过程中包括：
读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
根据所述投影数据，确定本次迭代后的输出图像序列；
若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
在一个示例中，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，包括：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
在一个示例中，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
在一个示例中，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
在一个示例中，所述积分图像序列以Z向优先存储的方式存储在内存中。
本发明实施例还提供一种CT***，包括CT设备和电子设备，其中：
所述CT设备，用于对受检者进行CT扫描，获得生数据：
将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，所述电子设备中存储有预先离线确定的所述XY平面权重矩阵，所述电子设备用于：
在每次迭代重建过程中包括：
读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
根据所述投影数据，确定本次迭代后的输出图像序列；
若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
在一个示例中，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，包括：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
在一个示例中，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
在一个示例中，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
在一个示例中，所述积分图像序列以Z向优先存储的方式存储在内存中。
本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，所述计算机可读存储介质中还存储有预先离线确定的所述XY平面权重矩阵，其中，所述程序被处理器执行时实现如下操作：
在每次迭代重建过程中包括：
读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
根据所述投影数据，确定本次迭代后的输出图像序列；
若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代。
在一个示例中，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，包括：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
在一个示例中，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
在一个示例中，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
在一个示例中，所述积分图像序列以Z向优先存储的方式存储在内存中。
对于装置和设备实施例而言，由于其基本对应于方法实施例，所以相关之处参见方法实施例的部分说明即可。以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的模块可以是或者也可以不是物理上分开的，作为模块显示的部件可以是或者也可以不是物理模块，即可以位于一个地方，或者也可以分布到多个网络模块上。可以根据实际的需要选择其中的部分或者全部模块来实现本说明书方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。
上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。
本领域技术人员在考虑说明书及实践这里申请的发明后，将容易想到本说明书的其它实施方案。本说明书旨在涵盖本说明书的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本说明书的一般性原理并包括本说明书未申请的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本说明书的真正范围和精神由下面的权利要求指出。
应当理解的是，本说明书并不局限于上面已经描述并在附图中示出的精确结构，并且可以在不脱离其范围进行各种修改和改变。本说明书的范围仅由所附的权利要求来限制。
以上所述仅为本说明书的较佳实施例而已，并不用以限制本说明书，凡在本说明书的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本说明书保护的范围之内。

Claims (8)

1.一种CT迭代重建方法，其特征在于，将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，预先离线确定所述XY平面权重矩阵并存储；
在每次迭代重建过程中包括：
读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
根据所述投影数据，确定本次迭代后的输出图像序列；
若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代；
所述利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据，包括：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
2.根据权利要求1所述的方法，其特征在于，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
3.根据权利要求1所述的方法，其特征在于，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
4.根据权利要求1所述的方法，其特征在于，所述积分图像序列以Z向优先存储的方式存储在内存中。
5.一种CT迭代重建装置，其特征在于，将三维前向投影所需的三维权重矩阵视为XY平面权重矩阵和Z向权重向量的乘积，所述装置中存储预先离线确定的所述XY平面权重矩阵，所述装置包括：
投影模块，用于在每次迭代重建过程中读取存储的XY平面权重矩阵，利用所述XY平面权重矩阵对输入图像序列进行三维前向投影，得到投影数据；
确定模块，用于根据所述投影数据，确定本次迭代后的输出图像序列；
输出模块，用于若满足预设的迭代停止条件，将本次迭代后的输出图像序列作为最终的CT重建图像序列；否则，将本次迭代后的输出图像序列作为下一次迭代的输入图像序列，执行下一次迭代；
所述投影模块具体用于：
对输入图像序列沿Z向进行图像积分，得到积分图像序列；
针对所述XY平面权重矩阵的每一元素，遍历每一检测器排，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分；
根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据。
6.根据权利要求5所述的装置，其特征在于，根据所述积分图像序列获取所述元素对应的Z向一维距离驱动前向投影积分，包括：
针对每一检测器排s，读取穿过像素j的射线i到像素j的距离SP(i，j)；其中，距离SP(i，j)是预先离线计算得到并存储的；
根据预设的检测器参数和所述距离SP(i，j)，确定Z向一维距离驱动前向投影的积分范围的最大值和最小值；
从所述积分图像序列对应的Z坐标值中确定与所述最大值相邻的第一值和第二值，以及确定与所述最小值相邻的第三值和第四值；
根据所述第一值对应的积分图像上所述元素对应的第一图像积分值和所述第二值对应的积分图像上所述元素对应的第二图像积分值，获得所述最大值对应的图像积分值；根据所述第三值对应的积分图像上所述元素对应的第三图像积分值和所述第四值对应的积分图像上所述元素对应的第四图像积分值，获得所述最小值对应的图像积分值；
根据所述最大值对应的图像积分值和所述最小值对应的图像积分值，确定所述元素对应的Z向一维距离驱动前向投影积分。
7.根据权利要求5所述的装置，其特征在于，根据所述XY平面权重矩阵中所有元素对应的Z向一维距离驱动前向投影积分以及所述XY平面权重矩阵中的元素，获得投影数据，包括：
将所述投影数据的计算分配至TileSize个线程，每个线程处理TileNum个检测器排对应的投影数据计算；每个线程对应的TileNum个检测器属于TileNum个组；
将TileSize个线程的计算结果按照对应组分别累加到TileNum个寄存器变量中；
累加完毕后，将TileNum个寄存器变量的累加结果写入投影数据p(v,c,s)对应的全局显存。
8.根据权利要求5所述的装置，其特征在于，所述积分图像序列以Z向优先存储的方式存储在内存中。
CN202110443876.3A 2021-04-23 2021-04-23 Ct迭代重建方法及装置 Active CN113284205B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110443876.3A CN113284205B (zh) 2021-04-23 2021-04-23 Ct迭代重建方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110443876.3A CN113284205B (zh) 2021-04-23 2021-04-23 Ct迭代重建方法及装置

Publications (2)

Publication Number Publication Date
CN113284205A CN113284205A (zh) 2021-08-20
CN113284205B true CN113284205B (zh) 2024-01-02

Family

ID=77277262

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110443876.3A Active CN113284205B (zh) 2021-04-23 2021-04-23 Ct迭代重建方法及装置

Country Status (1)

Country Link
CN (1) CN113284205B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105931280A (zh) * 2016-03-29 2016-09-07 中北大学 基于gpu的快速三维ct迭代重建方法
CN111110260A (zh) * 2019-12-24 2020-05-08 沈阳先进医疗设备技术孵化中心有限公司 一种图像重建方法、装置及终端设备
CN111192228A (zh) * 2020-01-02 2020-05-22 沈阳先进医疗设备技术孵化中心有限公司 图像处理方法、装置、ct设备及ct***

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7254209B2 (en) * 2003-11-17 2007-08-07 General Electric Company Iterative CT reconstruction method using multi-modal edge information
US9508163B2 (en) * 2013-06-14 2016-11-29 General Electric Company Accelerated iterative reconstruction
JP6368779B2 (ja) * 2013-06-28 2018-08-01 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． トモシンセシスデータからエッジ保存合成マンモグラムを生成するための方法
JP6215449B2 (ja) * 2014-03-14 2017-10-18 株式会社日立製作所 Ｘ線ｃｔ装置、及び処理装置
US10255697B2 (en) * 2014-11-20 2019-04-09 Koninklijke Philips N.V. Method for generation of synthetic mammograms from tomosynthesis data
US9911208B2 (en) * 2016-04-11 2018-03-06 Toshiba Medical Systems Corporation Apparatus and method of iterative image reconstruction using regularization-parameter control

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105931280A (zh) * 2016-03-29 2016-09-07 中北大学 基于gpu的快速三维ct迭代重建方法
CN111110260A (zh) * 2019-12-24 2020-05-08 沈阳先进医疗设备技术孵化中心有限公司 一种图像重建方法、装置及终端设备
CN111192228A (zh) * 2020-01-02 2020-05-22 沈阳先进医疗设备技术孵化中心有限公司 图像处理方法、装置、ct设备及ct***

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
J. J. Scheins et al.Fully-3D PET Image Reconstruction Using Scanner-Independent, Adaptive Projection Data and Highly Rotation-Symmetric Voxel Assemblies.《IEEE Transactions on Medical Imaging 》.2011,第30卷(第3期),第879 - 892页. *
赵星 ; 胡晶晶 ; 张慧滔 ; 张朋 ; .基于图形处理器的锥束CT快速迭代重建算法.北京理工大学学报.2010,第30卷(第11期),第1310-1314页. *

Also Published As

Publication number Publication date
CN113284205A (zh) 2021-08-20

Similar Documents

Publication Publication Date Title
CN108896943A (zh) 一种磁共振定量成像方法和装置
US8983162B2 (en) Method and apparatus for estimating monte-carlo simulation gamma-ray scattering in positron emission tomography using graphics processing unit
CN109377500B (zh) 基于神经网络的图像分割方法及终端设备
CN110163338B (zh) 具有运算阵列的芯片运算方法、装置、终端及芯片
US20120063658A1 (en) Image reconstruction based on accelerated method using polar symmetries
US9189870B2 (en) Method, computer readable medium and system for tomographic reconstruction
Liu et al. GPU-based branchless distance-driven projection and backprojection
CN110811667B (zh) 一种基于gpu加速的高精度pet重建方法及装置
Cuomo et al. A GPU parallel implementation of the local principal component analysis overcomplete method for DW image denoising
KR101283266B1 (ko) Ｇｐｕ를 이용한 양전자 방출 단층 촬영 영상에서의 몬테카를로 시뮬레이션 감마선 산란 추정 방법 및 장치
Corda-D’Incan et al. Memory-efficient training for fully unrolled deep learned PET image reconstruction with iteration-dependent targets
CN113284205B (zh) Ct迭代重建方法及装置
DE102016105973A1 (de) Bildverarbeitungsvorrichtung, Bildverarbeitungsverfahren und Programm
US8416239B2 (en) Intermediate image generation method, apparatus, and program
CN107730464A (zh) 基于块匹配的图像降噪并行算法
CN106910157A (zh) 一种多级并行的图像重建方法及装置
CN116434303A (zh) 基于多尺度特征融合的人脸表情捕捉方法、装置及介质
DE102022120819A1 (de) Quantisiertes neuronales- netzwerk-training und - inferenzieren
CN115375583A (zh) 一种pet参数图像的增强方法、装置、设备及存储介质
US6992704B2 (en) Image processing apparatus and method, recording medium, and program
CN106657966B (zh) 一种基于cpu多线程的集成成像3d片源快速生成方法
CN116739933A (zh) 一种融合***响应矩阵和神经网络的pet图像生成方法
CN112991482B (zh) 基于gpu的快速重建成像方法、设备及可读存储介质
CN116091219A (zh) 单障碍欧式期权Vega数值计算方法和装置
WO2023228910A1 (ja) 画像処理装置および画像処理方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20240204

Address after: 110167 No. 177-1 Innovation Road, Hunnan District, Shenyang City, Liaoning Province

Patentee after: Shenyang Neusoft Medical Systems Co.,Ltd.

Country or region after: China

Address before: Room 336, 177-1, Chuangxin Road, Hunnan New District, Shenyang City, Liaoning Province

Patentee before: Shenyang advanced medical equipment Technology Incubation Center Co.,Ltd.

Country or region before: China