CN106373165B - 断层合成图像重建方法和系统 - Google Patents
断层合成图像重建方法和系统 Download PDFInfo
- Publication number
- CN106373165B CN106373165B CN201610796748.6A CN201610796748A CN106373165B CN 106373165 B CN106373165 B CN 106373165B CN 201610796748 A CN201610796748 A CN 201610796748A CN 106373165 B CN106373165 B CN 106373165B
- Authority
- CN
- China
- Prior art keywords
- mrow
- voxel
- image
- msub
- value
- 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
Links
Classifications
-
- G06T12/30—
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
-
- G06T12/20—
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/421—Filtered back projection [FBP]
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/424—Iterative
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2211/00—Image generation
- G06T2211/40—Computed tomography
- G06T2211/436—Limited angle
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明涉及一种断层合成图像重建方法和系统。上述方法包括：S10、获取第n‑1三维图像的未更新区域，计算第m参考图像与所述第n‑1三维图像的体素差值，对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值S20、获取第n‑1补偿图像；S30、对第n角度的投影数据进行重建，得到第n三维图像；S40、判断n是否等于n0，若否，则将n更新为n+1，返回步骤S10，若是，则进入步骤S50；S50、将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，进入步骤S60，若是，则将第m重建图像结果作为最终的断层合成图像；S60、根据第m重建图像结果和迭代次数m对第m参考图像进行加权处理，获得第m+1参考图像，并将m更新为m+1，返回步骤S10。
Description
技术领域
本发明涉及断层合成图像处理技术领域，特别是涉及一种断层合成图像重建方法和系统。
背景技术
乳腺X射线摄影系统是乳腺疾病诊断和筛查的重要手段，其利用压迫板和支撑板将乳房等被测目标进行按压，并通过软X射线对上述被测目标的组织照射，通过探测器对穿透的X射线进行接收和处理，得到压迫后被测目标结构的二维投影图像。然而，上述乳腺X射线摄影系统中，被测目标的三维结构在二维投影图像上无法避免组织重叠和相互遮掩，而且依然无法获得对乳房等被测目标的纵向深度的结构信息。
为实现上述被测目标三维信息的获取并解决其中的组织重叠等问题，目前已经研发出基于数字断层合成技术(Digital Tomosynthesis，DTS)的乳腺断层合成系统(DigitalBreast Tomosynthesis，DBT)，此类系统在常规的乳腺X射线摄影系统基础上，增设X线球管独立旋转运动结构，在压迫器、被测目标(乳房等)、探测器保持固定条件下，将球管沿着圆弧轨道旋转运动，在有限角度范围内对按压的被测目标中的组织进行多次低剂量曝光，然后读取探测器上若干个探测器单元获取一系列二维投影数据，最终利用图像重建方法生成按压的被测目标三维图像。DBT获得的三维图像能解决二维数字乳腺摄影图像(DigitalMammogram，DM)的组织重叠和假阳性等问题，并可与同等压迫状态下的DM图像实现良好的对比和结合，在一定程度上提高了乳腺疾病的诊断效率。
若直接应用滤波反投影FBP和FDK等解析重建方法进行DBT图像重建，欠采样的有限稀疏角度投影对滤波器设计较为困难，投影信息的缺失难以获得精确满意的断层合成图像。而传统的反投影重建(BP)算法，虽能快速得到基本的三维结构信息，但重建结果缺乏精细的细节，无法直接应用于临床。
并行代数重建技术(Simultaneous Algebraic Reconstruction Technique，SART)是目前DBT重建相对精准有效的方法。然而，在DBT成像系统中，圆弧运动的X线球管和探测器法线形成的夹角越大，获得的纵向深度信息越丰富。然而角度过大会导致物体成像区域(Field of Vision，FOV)的部分区域超出探测器的有效范围，从而造成投影信息截断和缺失，形成投影缺失区域，且该投影缺失区域的大小会随着扫描角度和物体尺寸的变化而变化。在SART重建过程中依次在各个角度下对三维图像体素进行更新操作时，对应角度下投影缺失区域的体素将无法得到反投影更新，形成未更新区域，且相邻角度的投影截断差异将导致相应的断层合成图像中各个未更新区域边缘附近像素值出现多次急剧突变，这种现象最终会造成在重建得到的断层合成图像沿扫描方向侧面区域出现多条带状伪影，从而影响被测目标的断层合成图像重建的质量和效果。
发明内容
基于此，有必要针对传统重建方案影响被测目标断层合成图像质量的技术问题，提供一种新的断层合成图像重建方法和系统。
一种断层合成图像重建方法，包括如下步骤：
S10、获取第n-1三维图像的未更新区域，计算第m参考图像与所述第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，m为当前迭代次数，初始值设定为1；
S20、根据所述补偿基值对第n-1三维图像的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像；
S30、利用并行代数重建技术根据第n-1补偿图像对第n角度的投影数据进行重建，得到第n三维图像；
S40、判断n是否等于n0，若否，则将n更新为n+1，返回步骤S10，若是，则进入步骤S50；
S50、将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，进入步骤S60，若是，则将第m重建图像结果作为最终的断层合成图像；
S60、根据第m重建图像结果和迭代次数m对第m参考图像进行加权处理，获得第m+1参考图像，并将m更新为m+1，返回步骤S10。
一种断层合成图像重建系统，包括：
获取模块，用于获取第n-1三维图像的未更新区域，计算第m参考图像与所述第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，m为当前迭代次数，初始值设定为1；
修正模块，用于根据所述补偿基值对第n-1三维图像的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像；
重建模块，用于利用并行代数重建技术根据第n-1补偿图像对第n角度的投影数据进行重建，得到第n三维图像；
判断模块，用于判断n是否等于n0，若否，则将n更新为n+1，返回执行获取模块，若是，则执行更新模块；
更新模块，用于将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，执行设置模块，若是，则将第m重建图像结果作为最终的断层合成图像；
设置模块，用于将第m重建图像结果和迭代次数m对第m参考图像进行加权处理，得到第m+1参考图像，并将m更新为m+1，返回执行获取模块。
上述断层合成图像重建方法和系统，可以获取第n-1三维图像的未更新区域，计算第m参考图像与第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值，根据上述补偿基值对第n-1三维图像的体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像，以对第n角度的投影数据进行重建，得到第n三维图像，并循环上述过程，直至获取第n0三维图像，将第n0三维图像作为第m重建图像结果，并以此结果加权处理获得第m+1参考图像，继续对各个角度的投影数据进行迭代重建直至迭代次数达到预设次数后，将第m重建图像结果作为最终的断层合成图像；其可以消除断层合成图像沿扫描方向侧面区域的多条带状伪影，提高被测目标断层合成图像的重建效果。
附图说明
图1为一个实施例的断层合成图像重建方法流程图；
图2为一个实施例的DBT系统示意图；
图3为一个实施例的X射线射向被测目标示意图；
图4为一个实施例的数字乳房模型示意图；
图5为一个实施例的第一参考图像示意图；
图6为一个实施例的本发明所重建的DBT断层合成图像示意图；
图7为一个实施例的传统方案所重建的DBT断层合成图像示意图；
图8为一个实施例的断层合成图像重建系统结构示意图。
具体实施方式
下面结合附图对本发明的断层合成图像重建方法和系统的具体实施方式作详细描述。
参考图1，图1所示为一个实施例的断层合成图像重建方法流程图，包括如下步骤：
S10、获取第n-1三维图像的未更新区域，计算第m参考图像与所述第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，即DBT系统扫描的总投影数；n的初始值为2，m为当前迭代次数，初始值设定为1；
DBT系统(乳腺断层合成系统)可以如图2所示，其中压迫板、被测目标和探测器保持相对固定，X射线球管沿着弧形轨道依次在预先设定的n0个角度向被测目标发射X射线，以对相应的被测目标进行扫描，采集若干张低剂量曝光图像进行断层合成图像重建。如图2所示，X射线与探测器在某些夹角(如较大角度夹角)曝光时，被测目标会被投射至探测器外部，导致被测目标该区域的投影信息截断和缺失，形成相应的投影缺失区域，造成相应重建图像的信息缺失。利用并行代数重建技术对在进行第n角度下的投影更新修正得到的第n三维图像，该图像存在相应的缺失区域。上述缺失区域在该角度下投影数据的修正过程中得不到投影数据的更新修正，形成相应的未更新区域，由于第n三维图像是根据第n-1补偿图像和第n角度下的投影数据进行更新修正所得，因此第n三维图像的未更新区域与第n角度下的投影数据的截断缺失区域相对应，即第n三维图像的未更新区域即为第n角度下的投影数据的截断缺失区域。若对探测器接收的投影数据直接进行SART重建，依次对各个角度下对三维图像体素进行更新操作时，对应缺失区域的体素将无法得到更新，且相邻角度的缺失区域差异将导致缺失区域边缘像素值出现多次急剧突变，造成在重建图像沿扫描方向侧面区域出现多条带状伪影。上述被测目标可以包括DBT系统扫描的被测器官。
图3所示为第n-1角度下的X射线射向被测目标示意图，如图3所示，若X射线球管从探测器的第一端向第二端运动，获取最靠近探测器第二端的X射线101与探测器103所在平面所成的角度104，被测目标102与探测器103的位置关系(包括被测目标102的第一端距探测器103第一端的距离，包括被测目标102的第二端距探测器103第二端的距离)，还可以获取被测目标102的尺寸(如宽度、高度等)，根据光线投射关系和相关目标几何参数便可计算第n-1角度下的截断缺失区域的空间信息，包括图中的缺失区域105(未投影到探测器的区域)对应的体素在被测目标中的层范围、行范围和列范围，从而可确定在第n-1角度重建所得的第n-1三维图像的未更新区域。
上述步骤可以使用反投影(BP)重建方法等能够获取被测目标基本三维结构的算法对投影数据进行重建，以得到被测目标的初始参考图像，即第一参考图像。利用反投影(BP)重建方法能很好地利用投影数据的全局信息进行计算，得到较好的全局三维信息且不会出现带状伪影，利用上述反投影(BP)重建方法对各个角度下的投影数据进行重建，可以得到被扫描物体(被测目标)的基本三维结构，以进行SART重建时的未更新区域的修正补偿。
S20、根据所述补偿基值对第n-1三维图像的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像；
上述步骤中可以将补偿基值叠加至第n-1三维图像未更新区域的体素值上，以实现第n-1三维图像未更新区域体素的修正补偿。
S30、利用并行代数重建技术根据第n-1补偿图像对第n角度的投影数据进行重建，得到第n三维图像；
上述步骤中，可以先对第n-1补偿图像按角度n进行前向投影，得到第n估计投影数据，再对上述第n估计投影数据和第n角度下的投影数据进行差值加权等处理，以对前向投影和差值加权得到的结果进行反投影后，对第n-1补偿图像各体素进行加权修正，得到更新的第n三维图像。上述第n角度下的投影数据包括第n个角度的X射线照射被测目标时，探测器接收的投影数据。
S40、判断n是否等于n0，若否，则将n更新为n+1，返回步骤S10，若是，则进入步骤S50；
S50、将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，进入步骤S60，若是，则将第m重建图像结果作为最终的断层合成图像；
S60、根据第m重建图像结果和迭代次数m对第m参考图像进行加权处理，获得第m+1参考图像，并将m更新为m+1，返回步骤S10。
本实施例中，可以从X射线起始曝光的投影角度至X射线终止曝光的投影角度，依次对各个角度下被测目标对应的投影数据进行如步骤S10至S30所述的重建操作，当未达到迭代终止条件(一般为预设的迭代次数，亦可采用如图像梯度范数等其他终止条件)重复上述重建过程，以使最终得到的断层合成图像，与相应被测目标的三维图像信息保持较高的一致性。上述迭代次数可以根据断层合成图像的重建精度要求进行设置，比如设置为5或者10等值。
在上述迭代过程中，在第一次迭代过程时对于第一角度的投影数据，可以直接对上述第一角度的投影数据进行重建，以获取第一三维图像，利用光线投射关系和扫描几何参数确定第一三维图像的未更新区域，再根据上述第一三维图像和相应的未更新区域计算获得第一补偿图像，并根据该补偿图像和第二角度下的投影数据进行第二三维图像的重建。若在后续各次迭代过程(除第一次外)中，可以以上一次迭代过程所得到的第m重建图像结果作为初始三维图像，并以第m重建图像结果和第m参考图像按迭代次数进行加权计算，得到第m+1参考图像，依据上述初始三维图像依次各个角度的投影数据进行如步骤S10至S30所述的重建；即以上一次迭代过程所得到的第m重建图像结果作为第一三维图像，根据上述第一三维图像和获得的第m+1参考图像进行对第二三维图像的重建。
本实施例提出的断层合成图像重建方法，可以计算获取第n-1三维图像的未更新区域，计算第m参考图像与第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值，根据上述补偿基值对第n-1三维图像未更新区域各体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像，根据第n角度的投影数据进行重建，得到第n三维图像，并循环上述过程，直至获取第n0三维图像，将第n0三维图像作为第m重建图像结果，并此次结果加权处理获得第m+1参考图像，继续对各个角度的投影数据进行迭代重建直至迭代达到预设的迭代次数后，将第m重建图像结果作为最终的断层合成图像；其可以消除断层合成图像沿扫描方向侧面区域的多条带状伪影，提高被测目标断层合成图像的图像质量和重建效果。
在一个实施例中，上述步骤S10之前还可以包括：
对探测器接收的多张投影数据进行反投影重建，将反投影重建得到的三维图像设置为第一参考图像；
利用并行代数重建技术对第一角度的投影数据进行重建，得到第一三维图像；
获取第一三维图像的未更新区域，计算第一参考图像与所述第一三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；
将所述补偿基值与第一参考图像对第一三维图像的体素值进行修正补偿，根据修正补偿得到的体素值获取第一补偿图像；
利用并行代数重建技术根据第一补偿图像对第二角度的投影数据进行重建，得到第二三维图像。
上述第一三维图像的未更新区域为第一角度下探测器未接收到投影数据的区域，即该角度下投影截断缺失区域。第一参考图像与所述第一三维图像的体素差值即为上述第一参考图像中各个体素的体素值分别与第一三维图像中相应体素的体素值之差。
本实施例根据第一三维图像和第一参考图像进行区域补偿，并对第二角度的投影进行第二三维图像的重建，可以提高重建的第二三维图像的准确性，并消除第一角度下未更新区域的体素突变。
在一个实施例中，上述对探测器接收的多张投影数据进行反投影重建，将反投影重建得到的三维图像设置为第一参考图像的过程可以包括：
对探测器接收的多张不同角度的投影数据根据探测器位置、射线源位置，物体尺寸等三维空间关系和光线投射原理，对目标三维图像各个体素进行反投影重建，将反投影重建得到的三维图像设置为第一参考图像。
上述DBT向所述被测目标发射X射线时，探测器可以接收不同角度下X射线穿透物体的多张投影数据。
作为一个实施例，上述对探测器接收所述X射线的不同角度下投影数据进行反投影重建，将反投影重建得到的三维图像设置为第一参考图像的过程可以包括：
对探测器接收所述X射线的投影数据进行反投影(BP)重建方法处理，根据在不同角度下的X射线源、扫描物体(被测目标)、探测器等空间信息，确定各个体素在该角度下X射线穿透路径，并获得对应路径下的投影图像的像素坐标，依次叠加所有角度对应像素坐标所在的像素值并加权处理后，即可获取被测目标的基本三维结构，将所述基本三维结构设置为第一参考图像。
上述反投影(BP)重建方法能很好地利用投影数据的全局信息进行计算，得到较好的全局三维信息，不会出现带状伪影，利用上述反投影(BP)重建方法对探测器接收的投影数据进行重建，可以得到被扫描物体(被测目标)的基本三维结构，作为第一参考图像，继而迭代次数不断加权更新，以进行SART重建的未更新区域的修正补偿。
在一个实施例中，上述步骤S30可以包括：
对第n-1补偿图像按角度n进行前向投影，得到第n估计投影数据；
根据第n估计投影数据和第n角度下的投影数据进行像素的差值加权；
将上述前向投影和差值加权得到的结果反投影处理后，对第n-1补偿图像各体素进行加权修正，得到更新的第n三维图像。
本实施例可以对第n-1补偿图像按角度n进行前向投影，模拟当前角度n下球管发射X射线并穿透扫描物体投影至探测器上，计算投影图像各个像素(探测器元)对应X射线的穿透路径，并依次将在该穿透路径上对应n-1补偿图像三维空间上各个体素按对应路径权重进行加权累加，得到第n估计投影数据；根据第n估计投影数据和第n角度下的投影数据进行差值加权，将估计投影和采集投影对应各个像素的像素值进行差值处理，并按该像素点(探测器元)所在的X射线穿透扫描物体的总长度对该像素差进行加权处理；对前向投影和差值加权得到的结果进行反投影加权修正，按角度n下对应球管、探测器、扫描物体的空间位置，获取投影差值结果各坐标对应的X射线，按该X射线路径获得被穿透扫描物体对应的空间三维坐标，将上述投影差值结果按所在路径上各个体素的路径权重，加权叠加至第n-1补偿图像对应三维坐标像素上的像素值，得到更新修正的第n三维图像。
本实施例在探测器有限尺寸和大角度扫描的DBT重建条件下，可通过全局投影补偿修正的图像进行SART重建，克服缺失投影数据造成的重建图像扫描方向侧面区域带状伪影，保证了重建的断层合成图像的纵向深度信息，同时提高图像重建的质量。
在一个实施例中，上述步骤S60可以包括：
对所述第m重建图像结果的体素值和第m参考图像的体素值代入加权更新公式，计算第m+1参考图像的体素值；其中，所述加权更新公式为：
Refm+1(i,j,k)＝(Refm(i,j,k)+(m-1)pResm(i,j,k))/(1+(m-1)p)，
其中，Refm+1(i,j,k)表示第m+1参考图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Resm(i,j,k)表示第m重建图像结果中第i行第j列第k层的体素值，m表示当前迭代次数，p为更新权重系数。上述更新权重系数p可以设置为大于0的值。
本实施例在迭代过程中，根据第m重建图像结果对第m参考图像按迭代次数进行加权更新，保证初始参考图像的基本结构的同时，根据迭代次数逐步提高参考图像的质量，进一步提高SART重建过程未更新区域补偿修正的精确性。在一个实施例中，计算第m参考图像与所述第n-1三维图像的体素差值的过程包括：
将所述第m参考图像的体素值和第n-1三维图像的体素值分别代入体素差公式计算第n-1三维图像的体素差值；其中，所述体素差公式为：
Diffn-1(i,j,k)＝Imgn-1(i,j,k)-Refm(i,j,k)，
其中，Imgn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Diffn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层的体素差值，m表示当前迭代次数。
本实施例计算第m参考图像与第n-1三维图像的体素差值，以对未更新区域进行插值，可以消除第m参考图像与相应投影数据重建的第n三维图像的局部体素突变带来的误差。
作为一个实施例，上述根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值的过程包括：
将所述体素差值代入邻域均值公式计算未更新区域的邻域修正值；
将所述未更新区域的邻域修正值代入插值拟合公式计算未更新区域的补偿基值；
所述邻域均值公式为：
所述插值拟合公式为：
其中，Ndiffn-1(i，k)表示第n-1三维图像中第i行第k层的邻域修正值，J表示第n-1三维图像中体素的边界列数，JBn-1(i，k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，N表示预设邻域的体素列数，Comn-1(i，j，k)表示第n-1三维图像中第i行第j列第k层的补偿基值。
上述邻域可以根据被测目标的体素列数以及相应未更新区域的体素列数进行预设，通常情况下，在对第n-1三维图像进行第n角度的投影数据更新时，相应邻域可以设置为第n-1三维图像中，从上一角度投影未更新区域(即第n-1角度的投影数据无法更新的对应区域)的边界线开始，向非未更新区域侧的若干列体素；若上述邻域设置为从未更新区域边界线开始，向非未更新区域侧的两列体素，则N＝2。
本实施例根据邻域列数N、邻域修正值Ndiff(i,k)，可逐层按扫描轨迹方向拟合当前未更新区域各体素的补偿曲线。上述拟合方式可包括线性拟合，也可使用曲线拟合等方式计算未更新区域各体素的补偿值。
在一个实施例中，上述根据所述补偿基值对第n-1三维图像的体素值进行修正补偿的过程包括：
将所述补偿基值与第n-1三维图像的体素值分别代入修正补偿公式进行运算；所述修正补偿公式为：
其中，Fixn-1(i，j，k)表示第n-1补偿图像第i行第j列第k层的体素值，JBn-1(i，k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，Imgn-1(i，j，k)表示第n-1三维图像中第i行第j列第k层体素的体素值，Comn-1(i，j，k)表示第i行第j列第k层的补偿基值，Refm(i，j，k)表示第m参考图像第i行第j列第k层的体素值，m表示当前迭代次数。
本实施例采用叠加修正的方式，可保证叠加后所得到的未更新区域体素值得到全局投影信息的补偿，避免了传统方案下投影数据截断缺失造成的体素值突变。
作为一个实施例，上述将所述补偿基值与第n-1角度对应的叠加三维图像的体素值进行叠加的过程还可以采用均值滤波叠加方式等其他叠加方式。上述均值滤波叠加可为：
其中，Fixn-1(i，j，k)为叠加所得第n-1补偿图像第i行第j列第k层的体素值，Comn-1(i，j，k)表示第n-1三维图像第i行第j列第k层的补偿基值，Refm(i，o，k)表示第m参考图像第i行第o列第k层的体素值，L均值滤波的长度，m表示当前迭代次数。
上述局部插值方法，结合第一参考图像利用邻域拟合和局部插值的方法对未更新区域进行插值补偿，利用投影数据的全局信息获得信息较为完全的补偿三维图像，上述补偿三维图像可以用于SART重建下一次的投影更新操作，能很好地解决缺失邻域边界像素值突变导致的灰度值剧烈变化的问题，避免带状伪影的产生。
本发明提供的断层合成图像重建方法，通过预设角度范围扫描待重建物体，分别获取各角度下对应的投影数据，对获取的投影数据进行反投影重建方法处理，获取被扫描物体的基本三维结构，并作为初始参考图像(第一参考图像)；根据机械扫描几何信息和光线投射关系，分别对每个角度下投影数据反向投影处理，精确计算对应的未更新区域的分布；对投影数据进行SART重建，在当前角度下结合第一参考图像和对上次角度重建图像未更新区域进行局部插值，对未更新区域的体素进行数据修正；根据修正后的三维图像进行前向投影差值加权和反投影加权修正，完成当前投影下三维图像的更新，继而完成所有角度下的投影数据更新，得到当次迭代结果；根据迭代次数和迭代结果对参考图像更新，继续迭代重建得到最终的DBT断层合成图像。其在保证探测器尺寸有限和大角度扫描的前提下，可克服缺失投影数据造成的重建图像沿扫描方向侧面区域带状伪影，保证了DBT图像纵向深度信息，同时提高图像重建的质量。
下面结合一具体应用场景，对本发明上述断层合成图像重建方法进行说明。
选择如图4所示数字乳房模型作为被测目标(模型矩阵大小为400*600*50，体素大小为0.5*0.5*1mm)。DBT重建(断层合成图像重建)步骤可以包括：
1)对乳房模型进行DBT扫描，X射线球管运动范围由-30°到+30°，每隔3°进行曝光和投影数据采集，共采集21张投影数据Projn，n∈[1，21]，假设探测器矩阵大小为600*800，像素大小为0.4*0.4mm，对应得到的投影数据矩阵为600*800*21；
2)利用反投影(BP)重建方法同时对采集的21张投影数据Projn进行反投影重建，得如图5所示参考图像，反投影(BP)重建方法重建得到的得到基本的三维结构，用以作上述被测目标的第一参考图像Ref1；
3)预设初始三维图像Img0为0，根据第1角度的投影数据Proj1进行SART的前向投影和差值加权后，反投影至重建域加权修正后得到第1重建图像Img1；
4)从第2角度开始，即n＞1时，根据设定的投影几何和光线投射原理，计算得第n-1三维图像Imgn-1的未更新区域，计算第m参考图像Refm-1与所述第n-1三维图像Imgn-1的体素差值，根据体素差值对未更新区域Unn-1进行插值拟合，得到所述未更新区域各体素的补偿基值Comn-1；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，此处n0＝21，m为当前迭代次数，初始值设定为1；
所述邻域均值公式为：
所述插值拟合公式为：
其中，Ndiffn-1(i，k)表示第n-1三维图像第i行第k层的邻域修正值，J表示第n-1三维图像中体素的边界列数，JBn-1(i，k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，N表示预设邻域的体素列数，此处N＝3，Comn-1(i，j，k)表示第i行第j列第k层的补偿基值，补偿基值的拟合采用线性方式进行拟合。
5)利用所述修正补偿公式和补偿基值Comn-1对第n-1三维图像Imgn-1的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像Fixn-1；所述修正补偿公式为：
其中，Fixn-1(i，j，k)表示第n-1补偿图像第i行第j列第k层的体素值，JBn-1(i，k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，Imgn-1(i，j，k)表示第n-1三维图像中第i行第j列第k层体素的体素值，Comn-1(i，j，k)表示第n-1三维图像第i行第j列第k层的补偿基值，Refm(i，j，k)表示第m参考图像第i行第j列第k层的体素值，m表示当前迭代次数。
6)利用并行代数重建SART技术根据第n-1补偿图像Fixn-1对第n角度的投影数据Projn进行重建，即得到第n三维图像Imgn；
对第n-1补偿图像Fixn-1按角度n进行前向投影，模拟当前角度n下球管发射X射线并穿透扫描物体投影至探测器上，计算投影图像各个像素(探测器元)对应X射线的穿透路径，并依次将在该穿透路径上对应n-1补偿图像Fixn-1三维空间上各个体素按对应路径权重进行加权累加，得到第n估计投影数据Aprojn；
根据第n估计投影数据AProjn和第n角度下的投影数据Projn进行差值加权，将估计投影和采集投影对应各个像素的像素值进行差值处理，并按该像素点(探测器元)所在的X射线穿透扫描物体的总长度对该像素差进行加权处理；
对前向投影和差值加权得到的结果DProjn进行反投影加权修正，按角度n下对应球管、探测器、扫描物体的空间位置，获取投影差值结果DProjn各坐标对应的X射线，按该X射线路径获得被穿透扫描物体对应的空间三维坐标，将上述投影差值结果按所在路径上各个体素的路径权重，加权叠加至第n-1补偿图像Fixn-1对应三维坐标像素上的像素值，得到更新修正的第n三维图像Imgn。
7)判断n是否等于21，若否，则将n更新为n+1，返回步骤4)依次进行各角度下的未更新区域计算、三维图像补偿、投影数据重建，若是，则完成所有投影数据的更新重建，进入步骤下一步骤8)；
8)将第n三维图像Imgn作为第m重建图像结果Resm，判断m是否等于5(预设的迭代次数)，若是，则将第m重建图像结果作为最终的断层合成图像，若否，则对对第m重建图像结果Resm的体素值和第m参考图像Refm的体素值代入加权更新公式，计算第m+1参考图像Refm+1；其中，所述加权更新公式为：
Refm+1(i,j,k)＝(Refm(i,j,k)+(m-1)pResm(i,j,k))/(1+(m-1)p)，
其中，Refm+1(i,j,k)表示第m+1参考图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Resm(i,j,k)表示第m重建图像结果中第i行第j列第k层的体素值，m表示当前迭代次数，p为大于0的更新权重系数，此应用场景设置为1。
图6所示为应用本发明方案所重建的DBT断层合成图像，图7所示为应用传统方案所重建的DBT断层合成图像，如图7所示，采用传统的SART重建方法，由于探测器尺寸受限和角度过大，导致某些角度下投影数据缺失截断，导致该角度反投影的成像FOV出现缺失区域，在各投影更新迭代操作时，将产生缺失区域边界体素突变，形成多带状伪影，从而影响相应断层合成图像重建的效果。将应用本发明提供的断层合成图像重建方法所重建得到的断层合成图像(如图6)与图7所示的DBT断层合成图像进行比较，可以表明，利用本发明的断层合成图像重建方法，被扫描物体(被测目标)的重建图像与理想模体的结构信息基本一致；相对于边界出现多带状伪影的断层合成图像(如图7)，通过本发明上述,断层合成图像重建方法所得到的断层合成图像，显然更为清晰准确，有效解决上述大角度扫描和探测器尺寸受限导致的带状伪影问题，而且无需增加外围部件实现，执行效率高，稳定性强。
参考图8，图8所示为一个实施例的断层合成图像重建系统结构示意图，包括：
获取模块10，用于获取第n-1三维图像的未更新区域，计算第m参考图像与所述第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，m为当前迭代次数，初始值设定为1；
修正模块20，用于根据所述补偿基值对第n-1三维图像的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像；
重建模块30，用于利用并行代数重建技术根据第n-1补偿图像对第n角度的投影数据进行重建，得到第n三维图像；
判断模块40，用于判断n是否等于n0，若否，则将n更新为n+1，返回执行获取模块，若是，则执行更新模块；
更新模块50，用于将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，执行设置模块，若是，则将第m重建图像结果作为最终的断层合成图像；
设置模块60，用于将第m重建图像结果和迭代次数m对第m参考图像进行加权处理，得到第m+1参考图像，并将m更新为m+1，返回执行获取模块。
在一个实施例中，上述重建模块进一步用于：
对第n-1补偿图像按角度n进行前向投影，得到第n估计投影数据；
根据第n估计投影数据和第n角度下的投影数据进行像素的差值加权；
将上述前向投影和差值加权得到的结果反投影处理后，对第n-1补偿图像各体素进行加权修正，得到更新的第n三维图像。
在一个实施例中，上述设置模块进一步用于：
对所述第m重建图像结果的体素值和第m参考图像的体素值代入加权更新公式，计算第m+1参考图像的体素值；其中，所述加权更新公式为：
Refm+1(i,j,k)＝(Refm(i,j,k)+(m-1)pResm(i,j,k))/(1+(m-1)p)，
其中，Refm+1(i,j,k)表示第m+1参考图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Resm(i,j,k)表示第m重建图像结果中第i行第j列第k层的体素值，m表示当前迭代次数。
本发明提供的断层合成图像重建系统与本发明提供的断层合成图像重建方法一一对应，在所述断层合成图像重建方法的实施例阐述的技术特征及其有益效果均适用于断层合成图像重建系统的实施例中，特此声明。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
Claims (6)
1.一种断层合成图像重建方法，其特征在于，包括如下步骤：
S10、获取第n-1三维图像的未更新区域，计算第m参考图像与所述第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，m为当前迭代次数，初始值设定为1；对被测目标的投影数据进行重建，得到被测目标的初始参考图像，即第一参考图像；在第一次迭代过程时，对第一角度的投影数据进行重建，获取第一三维图像，利用光线投射关系和扫描几何参数确定第一三维图像的未更新区域，根据所述第一三维图像和相应的未更新区域计算获得第一补偿图像，并根据该补偿图像和第二角度下的投影数据进行第二三维图像的重建；
所述计算第m参考图像与所述第n-1三维图像的体素差值的过程包括：
将所述第m参考图像的体素值和第n-1三维图像的体素值分别代入体素差公式计算第n-1三维图像的体素差值；其中，所述体素差公式为：
Diffn-1(i,j,k)＝Imgn-1(i,j,k)-Refm(i,j,k)；
其中，Imgn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Diffn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层的体素差值，m表示当前迭代次数；
所述根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值的过程包括：
将所述体素差值代入邻域均值公式计算未更新区域的邻域修正值；
将所述未更新区域的邻域修正值代入插值拟合公式计算未更新区域的补偿基值；
所述邻域均值公式为：
<mrow>
<msub>
<mi>Ndiff</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>N</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mrow>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</munderover>
<msub>
<mi>Diff</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mi>N</mi>
</mrow>
所述插值拟合公式为：
<mrow>
<msub>
<mi>Com</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>Ndiff</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mfrac>
<mrow>
<mi>j</mi>
<mo>-</mo>
<mi>J</mi>
</mrow>
<mrow>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>J</mi>
</mrow>
</mfrac>
<mo>,</mo>
<mi>j</mi>
<mo>&Element;</mo>
<mo>(</mo>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>(</mo>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
</mrow>
<mo>)</mo>
<mo>,</mo>
<mi>J</mi>
<mo>&rsqb;</mo>
<mo>,</mo>
</mrow>
其中，Ndiffn-1(i,k)表示第n-1三维图像中第i行第k层的邻域修正值，J表示第n-1三维图像中体素的边界列数，JBn-1(i,k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，N表示预设邻域的体素列数，Comn-1(i,j,k)表示第n-1三维图像第i行第j列第k层的补偿基值；
S20、根据所述补偿基值对第n-1三维图像的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像；
S30、利用并行代数重建技术根据第n-1补偿图像对第n角度的投影数据进行重建，得到第n三维图像；
S40、判断n是否等于n0，若否，则将n更新为n+1，返回步骤S10，若是，则进入步骤S50；
S50、将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，进入步骤S60，若是，则将第m重建图像结果作为最终的断层合成图像；
S60、根据第m重建图像结果和迭代次数m对第m参考图像进行加权处理，获得第m+1参考图像，并将m更新为m+1，返回步骤S10；所述步骤S60包括：
对所述第m重建图像结果的体素值和第m参考图像的体素值代入加权更新公式，计算第m+1参考图像的体素值；其中，所述加权更新公式为：
Refm+1(i,j,k)＝(Refm(i,j,k)+(m-1)pResm(i,j,k))/(1+(m-1)p)，
其中，Refm+1(i,j,k)表示第m+1参考图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Resm(i,j,k)表示第m重建图像结果中第i行第j列第k层的体素值，m表示当前迭代次数，p为更新权重系数。
2.根据权利要求1所述的断层合成图像重建方法，其特征在于，所述步骤S10之前还包括：
对探测器接收的多张投影数据进行反投影重建，将反投影重建得到的三维图像设置为第一参考图像；
利用并行代数重建技术对第一角度的投影数据进行重建，得到第一三维图像；
获取第一三维图像的未更新区域，计算第一参考图像与所述第一三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；
将所述补偿基值与第一参考图像对第一三维图像的体素值进行修正补偿，根据修正补偿得到的体素值获取第一补偿图像；
利用并行代数重建技术根据第一补偿图像对第二角度的投影数据进行重建，得到第二三维图像。
3.根据权利要求1所述的断层合成图像重建方法，其特征在于，所述步骤S30包括：
对第n-1补偿图像按角度n进行前向投影，得到第n估计投影数据；
根据第n估计投影数据和第n角度下的投影数据进行像素的差值加权；
将上述前向投影和差值加权得到的结果反投影处理后，对第n-1补偿图像各体素进行加权修正，得到更新的第n三维图像。
4.根据权利要求1所述的断层合成图像重建方法，其特征在于，所述根据补偿基值对第n-1三维图像的体素值进行修正补偿的过程包括：
将所述补偿基值与第n-1三维图像的体素值分别代入修正补偿公式进行运算；所述修正补偿公式为：
<mrow>
<msub>
<mi>Fix</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Img</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>j</mi>
<mo>&Element;</mo>
<mo>&lsqb;</mo>
<mn>0</mn>
<mo>,</mo>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Ref</mi>
<mi>m</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>Com</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>j</mi>
<mo>&Element;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mi>J</mi>
<mo>&rsqb;</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
其中，Fixn-1(i,j,k)表示第n-1补偿图像第i行第j列第k层的体素值，JBn-1(i,k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，Imgn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层体素的体素值，Comn-1(i,j,k)表示第n-1三维图像第i行第j列第k层的补偿基值，Refm(i,j,k)表示第m参考图像第i行第j列第k层的体素值。
5.一种断层合成图像重建系统，其特征在于，包括：
获取模块，用于获取第n-1三维图像的未更新区域，计算第m参考图像与所述第n-1三维图像的体素差值，根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值；其中，n为大于1且小于或等于n0的整数，n0为发射X射线的总角度数，m为当前迭代次数，初始值设定为1；对被测目标的投影数据进行重建，得到被测目标的初始参考图像，即第一参考图像；在第一次迭代过程时，对第一角度的投影数据进行重建，获取第一三维图像，利用光线投射关系和扫描几何参数确定第一三维图像的未更新区域，根据所述第一三维图像和相应的未更新区域计算获得第一补偿图像，并根据该补偿图像和第二角度下的投影数据进行第二三维图像的重建；
所述计算第m参考图像与所述第n-1三维图像的体素差值的过程包括：
将所述第m参考图像的体素值和第n-1三维图像的体素值分别代入体素差公式计算第n-1三维图像的体素差值；其中，所述体素差公式为：
Diffn-1(i,j,k)＝Imgn-1(i,j,k)-Refm(i,j,k)；
其中，Imgn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Diffn-1(i,j,k)表示第n-1三维图像中第i行第j列第k层的体素差值，m表示当前迭代次数；
所述根据体素差值对未更新区域进行插值拟合，得到所述未更新区域各体素的补偿基值的过程包括：
将所述体素差值代入邻域均值公式计算未更新区域的邻域修正值；
将所述未更新区域的邻域修正值代入插值拟合公式计算未更新区域的补偿基值；
所述邻域均值公式为：
<mrow>
<msub>
<mi>Ndiff</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>N</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mrow>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</munderover>
<msub>
<mi>Diff</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mi>N</mi>
</mrow>
所述插值拟合公式为：
<mrow>
<msub>
<mi>Com</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>Ndiff</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mfrac>
<mrow>
<mi>j</mi>
<mo>-</mo>
<mi>J</mi>
</mrow>
<mrow>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>J</mi>
</mrow>
</mfrac>
<mo>,</mo>
<mi>j</mi>
<mo>&Element;</mo>
<mo>(</mo>
<msub>
<mi>JB</mi>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>(</mo>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>k</mi>
</mrow>
<mo>)</mo>
<mo>,</mo>
<mi>J</mi>
<mo>&rsqb;</mo>
<mo>,</mo>
</mrow>
其中，Ndiffn-1(i,k)表示第n-1三维图像中第i行第k层的邻域修正值，J表示第n-1三维图像中体素的边界列数，JBn-1(i,k)表示未更新区域的边界在第n-1三维图像第i行第k层所处的列数，N表示预设邻域的体素列数，Comn-1(i,j,k)表示第n-1三维图像第i行第j列第k层的补偿基值；
修正模块，用于根据所述补偿基值对第n-1三维图像的未更新区域体素值进行修正补偿，根据修正补偿得到的体素值获取第n-1补偿图像；
重建模块，用于利用并行代数重建技术根据第n-1补偿图像对第n角度的投影数据进行重建，得到第n三维图像；
判断模块，用于判断n是否等于n0，若否，则将n更新为n+1，返回执行获取模块，若是，则执行更新模块；
更新模块，用于将第n三维图像作为第m重建图像结果，判断m是否等于预设的迭代次数，若否，执行设置模块，若是，则将第m重建图像结果作为最终的断层合成图像；
设置模块，用于将第m重建图像结果和迭代次数m对第m参考图像进行加权处理，得到第m+1参考图像，并将m更新为m+1，返回执行获取模块；所述设置模块进一步用于：
对所述第m重建图像结果的体素值和第m参考图像的体素值代入加权更新公式，计算第m+1参考图像的体素值；其中，所述加权更新公式为：
Refm+1(i,j,k)＝(Refm(i,j,k)+(m-1)pResm(i,j,k))/1+(m-1)p，
其中，Refm+1(i,j,k)表示第m+1参考图像中第i行第j列第k层的体素值，Refm(i,j,k)表示第m参考图像中第i行第j列第k层的体素值，Resm(i,j,k)表示第m重建图像结果中第i行第j列第k层的体素值，m表示当前迭代次数。
6.根据权利要求5所述的断层合成图像重建系统，其特征在于，所述重建模块进一步用于：
对第n-1补偿图像按角度n进行前向投影，得到第n估计投影数据；
根据第n估计投影数据和第n角度下的投影数据进行像素的差值加权；
将上述前向投影和差值加权得到的结果反投影处理后，对第n-1补偿图像各体素进行加权修正，得到更新的第n三维图像。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610796748.6A CN106373165B (zh) | 2016-08-31 | 2016-08-31 | 断层合成图像重建方法和系统 |
| PCT/CN2016/098738 WO2018040126A1 (zh) | 2016-08-31 | 2016-09-12 | 断层合成图像重建方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610796748.6A CN106373165B (zh) | 2016-08-31 | 2016-08-31 | 断层合成图像重建方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106373165A CN106373165A (zh) | 2017-02-01 |
| CN106373165B true CN106373165B (zh) | 2017-11-28 |
Family
ID=57899881
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610796748.6A Expired - Fee Related CN106373165B (zh) | 2016-08-31 | 2016-08-31 | 断层合成图像重建方法和系统 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN106373165B (zh) |
| WO (1) | WO2018040126A1 (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108010096A (zh) | 2017-11-30 | 2018-05-08 | 上海联影医疗科技有限公司 | Cbct图像重建方法、装置和cbct设备 |
| DE102017223603B4 (de) * | 2017-12-21 | 2020-01-09 | Siemens Healthcare Gmbh | Verfahren zur Rekonstruktion eines dreidimensionalen, mit einer Biplan-Röntgeneinrichtung aufgenommenen Bilddatensatzes, Biplan-Röntgeneinrichtung, Computerprogramm und elektronisch lesbarer Datenträger |
| WO2019047545A1 (zh) * | 2018-05-04 | 2019-03-14 | 西安大医集团有限公司 | 低剂量成像方法及装置 |
| CN109147035B (zh) * | 2018-08-03 | 2023-08-08 | 上海电气集团股份有限公司 | 三维模型的显示方法及系统 |
| EP3660790A1 (en) | 2018-11-28 | 2020-06-03 | Koninklijke Philips N.V. | System for reconstructing an image of an object |
| CN109658465B (zh) * | 2018-12-07 | 2023-07-04 | 广州华端科技有限公司 | 图像重建过程中的数据处理、图像重建方法和装置 |
| CN110335204B (zh) * | 2019-05-07 | 2021-06-01 | 中国人民解放军陆军工程大学 | 一种热成像图像增强方法 |
| CN110706338B (zh) * | 2019-09-30 | 2023-05-02 | 东软医疗系统股份有限公司 | 图像重建方法、装置、ct设备及ct系统 |
| CN112258627B (zh) * | 2020-09-18 | 2023-09-15 | 中国科学院计算技术研究所 | 一种局部断层三维重构系统 |
| CN115147378B (zh) * | 2022-07-05 | 2023-07-25 | 哈尔滨医科大学 | 一种ct图像分析提取方法 |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5133690B2 (ja) * | 2004-10-08 | 2013-01-30 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | ボクセルに依存する補間を用いる画像再構成 |
| CN102419866B (zh) * | 2010-09-27 | 2013-08-21 | 北京农业智能装备技术研究中心 | 用于ct图像断层重建的滤波函数建立方法及断层重建方法 |
| CN102103757B (zh) * | 2010-12-27 | 2012-09-19 | 中国科学院深圳先进技术研究院 | 锥束图像重建方法及装置 |
| CN103971349B (zh) * | 2013-01-30 | 2017-08-08 | 上海西门子医疗器械有限公司 | 计算机断层扫描图像重建方法和计算机断层扫描设备 |
| CN103455989B (zh) * | 2013-09-24 | 2016-06-15 | 南京大学 | 一种结合超声图像提高有限角度ct成像质量的方法 |
| US9600924B2 (en) * | 2014-02-05 | 2017-03-21 | Siemens Aktiengesellschaft | Iterative reconstruction of image data in CT |
| CN105488823B (zh) * | 2014-09-16 | 2019-10-18 | 株式会社日立制作所 | Ct图像重建方法、ct图像重建装置及ct系统 |
| CN105184835B (zh) * | 2015-09-15 | 2018-11-06 | 上海联影医疗科技有限公司 | 乳腺断层图像重建方法和装置 |
-
2016
- 2016-08-31 CN CN201610796748.6A patent/CN106373165B/zh not_active Expired - Fee Related
- 2016-09-12 WO PCT/CN2016/098738 patent/WO2018040126A1/zh not_active Ceased
Also Published As
| Publication number | Publication date |
|---|---|
| WO2018040126A1 (zh) | 2018-03-08 |
| CN106373165A (zh) | 2017-02-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106373165B (zh) | 断层合成图像重建方法和系统 | |
| US8189735B2 (en) | System and method for reconstruction of X-ray images | |
| CN109671128B (zh) | 图像重建过程中的数据处理、图像重建方法和装置 | |
| Yu et al. | Data consistency based translational motion artifact reduction in fan-beam CT | |
| CN103180879B (zh) | 用于从投影数据对对象进行混合重建的设备和方法 | |
| US10722178B2 (en) | Method and apparatus for motion correction in CT imaging | |
| CN106683144A (zh) | 一种图像迭代重建方法及装置 | |
| CN103455989A (zh) | 一种结合超声图像提高有限角度ct成像质量的方法 | |
| CN102783962A (zh) | 提供受抑制的超出测量场伪影的图像数据组的方法及设备 | |
| CN107764846A (zh) | 一种正交直线扫描的cl成像系统及分析方法 | |
| CN105580053B (zh) | 运动补偿的迭代重建 | |
| CN101496064A (zh) | 用于重构图像的方法和用于重构图像的重构系统 | |
| WO2014041889A1 (ja) | Ｘ線ｃｔ装置およびｘ線ｃｔ画像の処理方法 | |
| CN107796834B (zh) | 一种正交电子直线扫描cl成像系统及方法 | |
| CN109920020A (zh) | 一种锥束ct病态投影重建伪影抑制方法 | |
| CN105488823A (zh) | Ct图像重建方法、ct图像重建装置及ct系统 | |
| US12315145B2 (en) | Image processing device, learning device, radiography system, image processing method, learning method, image processing program, and learning program | |
| JP2015231528A (ja) | Ｘ線コンピュータ断層撮像装置及び医用画像処理装置 | |
| KR102348139B1 (ko) | 이중 해상도의 관심 영역 내외 투영 데이터를 이용한 체내 단층 촬영 방법 및 시스템 | |
| US9965875B2 (en) | Virtual projection image method | |
| US10966670B2 (en) | Imaging system and method for dual-energy and computed tomography | |
| US10089758B1 (en) | Volume image reconstruction using projection decomposition | |
| US9237873B2 (en) | Methods and systems for CT projection domain extrapolation | |
| CN1781450A (zh) | 三维图像处理装置 | |
| Dennerlein et al. | Geometric jitter compensation in cone-beam CT through registration of directly and indirectly filtered projections |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | 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 | ||
| 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: 20171128 |