一种基于压缩X射线断层合成的联合优化方法与流程

本发明属于计算成像技术领域，尤其涉及一种基于压缩x射线断层合成的联合优化方法。

背景技术：

目前，x射线计算机断层扫描技术(x-raycomputedtomography，简称ct)已广泛应用于临床诊断、安全检查、工业检测等领域，但是传统的ct会围绕对象进行完全一周的测量，并利用计算机处理的算法，如滤波反投影算法(fileredback-projectionalgorithm，简称fbp)重建对象，这样就大大增加了辐射剂量，采集时间和测量的成本及几何复杂度。x射线断层合成技术(x-raytomosynthesis)作为ct的一种替代技术，从对照射角度的限制出发，进行一组不完全投影，并基于迭代算法重建对象，这种方法有很多优点，例如可大大降低辐射剂量，减小测量数并且有着单一的传感几何。断层合成的不完全测量有利于减少辐射，但是不可避免的会导致不适定逆重构问题。压缩x射线断层合成技术(compressivex-raytomosynthesis)作为一种新兴的技术，减缓了在小测量数和不适定逆重构问题间的权衡，它可以由一组x射线光源产生的二维投影数据中重构出三维的物体，这其中编码孔径用来调制照明结构以降低辐射剂量。通过编码孔径的优化，可以得到最佳的重构图样。

但是，当前的编码孔径优化方法都是只优化编码孔径图样，因此限制了优化的自由度，而自由度的限制会导致无法获得一个最优的重构结果，在断层合成技术中突破这种限制显得至关重要。

技术实现要素：

为解决上述问题，本发明提供一种基于压缩x射线断层合成的联合优化方法，能够显著提高优化自由度，有效改善三维物体的重建效果，图像质量明显提高。

一种基于压缩x射线断层合成的联合优化方法，包括以下步骤：

s1：根据压缩感知理论，获取目标函数ω：

其中，i为各光源依次照射三维物体后，在探测器上得到相应投影的测量值的索引序号，i＝(u-1)×m+v，u＝1,2,...,p，p为光源的数量，u为光源索引，v＝1,2,...,m为探测器的栅格个数，m＝mx×my，mx和my分别为探测器栅格化后沿着x轴和y轴的维度；αu为光源亮暗表征值的中间变量，为光源亮暗表征值的中间变量的优化值；θi为第i个索引序号对应的稀疏表示系数，为第i个索引序号对应的稀疏表示系数的优化值；xu为第u个光源在编码孔径上形成的照射圆域中心横坐标，为第u个光源在编码孔径上形成的照射圆域中心横坐标的优化值；yu为第u个光源在编码孔径上形成的照射圆域中心纵坐标，为第u个光源在编码孔径上形成的照射圆域中心纵坐标的优化值；为将预设的感知矩阵α每一列归一化后得到的感知矩阵，*为转置，i为单位矩阵，f为弗罗贝尼乌斯范数；

s2：采用最陡梯度下降法优化目标函数ω，得到最优的光源索引、光源入射角度以及编码孔径；

s3：采用最优的光源索引、光源入射角度以及编码孔径搭建的断层合成系统照射三维物体，得到三维物体在探测器上的投影；

s4：采用梯度投影重构算法对所述投影进行图像恢复，得到三维物体的重构图样。

进一步地，所述目标函数ω的获取方法具体为：

s101：将三维物体栅格化为nx×ny×nλ的立体图形，其中nx，ny和nλ分别为三维物体沿着x，y和z轴的维度；

s102：获取第u个光源对应的编码孔径的等效传输函数

其中，⊙是点乘运算，πu为第u个光源在编码孔径上投影的光源图形的传输函数，λu为第u个光源入射角度的传输函数，tu为第u个光源对应的编码孔径的传输函数；

s103：获取各光源在探测器上的成像模型y，其中光源的数量用p表示，具体的：

其中，yu为第u个光源的成像模型，cu为等效传输函数对角化后的对角矩阵，hu为预设的第u个光源的断层合成系统矩阵，c为所有光源对应编码孔径的对角矩阵cu集合，h为预设的所有光源对应的断层合成系统矩阵，f为采用光栅扫描的方式扫描三维物体后，三维物体的每个栅格对应的衰减系数向量；

s104：获取成像模型y在完备字典下的稀疏表示：

y＝chψθ

其中，ψ为完备字典，θ为稀疏表示系数；

s105：根据压缩感知理论，获取目标函数ω1：

其中，为稀疏表示系数优化值，a＝cw，w＝hψ，γ为预设的正则项的权重因子；

s106：根据目标函数ω1得到的稀疏表示系数优化值恢复三维物体图像，需要满足以下有限等距约束条件：

其中，δs为限制等距常数；

s107：以限制等距常数δs最小为优化目标，获取目标函数ω2：

其中，为所有光源对应编码孔径的对角矩阵cu集合c经过优化后的值；

s108：根据集合c与等效传输函数的关系，将目标函数ω2变换为目标函数ω3：

s109：对目标函数ω3进行连续化，得到目标函数ω。

进一步地，步骤s109中所述对目标函数ω3进行连续化，得到目标函数ω具体为：

定义如下参数变换式：

其中，πu为光源亮暗表征值，πu＝1表示第u个光源状态为亮，否则状态为暗，tanh为双曲正切函数，ti为编码孔径的二值参数，ti＝1表示第u个光源下编码孔径第v个栅格透光，ti＝0表示第u个光源下编码孔径第v个栅格不透光，编码孔径与探测器的栅格一一对应，λi为光源入射角度的二值参数，λi＝1表示第u个光源下编码孔径第v个栅格在光源照射的圆域内，λi＝0表示第u个光源下编码孔径第v个栅格在光源照射的圆域外，γ为sigmoid函数，(xi,yi)为编码孔径上栅格的坐标；

将上述函数代入目标函数ω3，得到目标函数ω。

进一步地，如果第u个光源状态为亮，则步骤s102中所述第u个光源的传输函数πu为全1矩阵，否则为全0矩阵；

第u个光源入射角度的传输函数λu在编码孔径上表征为一个圆形窗口，且所述圆形窗口表示编码孔径上的曝光区域，在圆形区域以内的像素为1，否则为0。

有益效果：

本发明提供一种基于压缩x射线断层合成的联合优化方法，由于不同的光源索引对应不同的光源位置，则本发明基于压缩断层合成技术，协同优化光源位置、入射角度以及编码孔径；基于压缩感知理论，通过联合优化光源位置、入射角度以及编码孔径，使得断层合成系统中的感知矩阵的相干性最小；结合基于梯度投影重构算法，采用正则化方法来减少收敛误差，并实现对三维物体的均匀感知；与单独优化光源位置、入射角度、编码孔径相比较，本发明提供的方法能够显著提高优化自由度，有效改善三维物体的重建效果，图像质量明显提高。

附图说明

图1为本发明提供的一种基于压缩x射线断层合成的联合优化方法的流程图；

图2为本发明提供的sco的物理模型；

图3为本发明提供的三维物体的四层切片；

图4为本发明提供的采用随机编码孔径对三维物体重构的结果示意图；

图5为本发明提供的采用单独优化光源位置的方法得到的三维物体的重构结果示意图；

图6为本发明提供的采用单独优化入射角度的方法得到的三维物体的重构结果示意图；

图7为本发明提供的采用单独优化编码孔径的方法得到的三维物体的重构结果示意图；

图8为本发明提供的采用联合优化的方法得到的三维物体的重构结果示意图。

具体实施方式

为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。

当前的编码孔径优化方法通常只考虑编码孔径，而x射线光源的位置和入射角度都会对传感矩阵有着重要的影响，因此应该考虑光源位置、入射角度和编码孔径之间的协同作用，进一步提高断层合成系统的重建性能。为此本发明引入光源、角度、编码孔径协同优化(sourceandcodedaperturejointoptimization，简称sco)。

参见图1，该图为本实施例提供的一种基于压缩x射线断层合成的联合优化方法的流程图。一种基于压缩x射线断层合成的联合优化方法，具体包括以下步骤：

s1：根据压缩感知理论，获取目标函数ω：

其中，i为各光源依次照射三维物体后，在探测器上得到相应投影的测量值的索引序号，i＝(u-1)×m+v，u＝1,2,...,p，p为光源的数量，u为光源索引，v＝1,2,...,m为探测器的栅格个数，m＝mx×my，mx和my分别为探测器栅格化后沿着x轴和y轴的维度；αu为光源亮暗表征值的中间变量，为光源亮暗表征值的中间变量的优化值；θi为第i个索引序号对应的稀疏表示系数，为第i个索引序号对应的稀疏表示系数的优化值；xu为第u个光源在编码孔径上形成的照射圆域中心横坐标，为第u个光源在编码孔径上形成的照射圆域中心横坐标的优化值；yu为第u个光源在编码孔径上形成的照射圆域中心纵坐标，为第u个光源在编码孔径上形成的照射圆域中心纵坐标的优化值；为将预设的感知矩阵α每一列归一化后得到的感知矩阵，*为转置，i为单位矩阵，f为弗罗贝尼乌斯范数；

s2：采用最陡梯度下降法优化目标函数ω，得到最优的光源索引、光源入射角度以及编码孔径；

s3：采用最优的光源索引、光源入射角度以及编码孔径搭建的断层合成系统照射三维物体，得到三维物体在探测器上的投影；

s4：采用梯度投影重构算法对所述投影进行图像恢复，得到三维物体的重构图样。

进一步地，下面给出目标函数ω的一种获取方法。

参见图2，该图为本实施例提供的sco的物理模型。目标函数ω的获取方法，具体包括以下步骤：

s101：将三维物体栅格化为nx×ny×nλ的立体图形，其中nx，ny和nλ分别为三维物体沿着x，y和z轴的维度。

s102：获取第u个光源对应的编码孔径的等效传输函数

其中，⊙是点乘运算，πu为第u个光源在编码孔径上投影的光源图形的传输函数，λu为第u个光源入射角度的传输函数，tu为第u个光源对应的编码孔径的传输函数。

需要说明的是，第u个光源的传输函数πu为全1矩阵，否则为全0矩阵；

第u个光源入射角度的传输函数λu在编码孔径上表征为一个圆形窗口，且所述圆形窗口表示编码孔径上的曝光区域，在圆形区域以内的像素为1，否则为0。tu是一个二值矩阵，用来表示编码孔径。

由此可见，本实施例考虑多个优化变量，将光源数量、入射角度、编码孔径整合到编码孔径，与单独优化光源位置、角度、编码孔径相比较，所提出的方法能显著提高优化自由度，从而实现重建效果的改善。

s103：获取各光源在探测器上的成像模型y，其中光源的数量用p表示，具体的：

其中，yu为第u个光源的成像模型，cu为等效传输函数对角化后的对角矩阵，hu为预设的第u个光源的断层合成系统矩阵，c为所有光源对应编码孔径的对角矩阵cu集合，h为预设的所有光源对应的断层合成系统矩阵，为采用光栅扫描的方式扫描三维物体后，三维物体的每个栅格对应的衰减系数向量。

需要说明的是，成像模型的离散形式可以表述为：

y＝hf

其中，h为预设的断层合成系统矩阵。对于每一个光源来说，第u个光源的成像模型表示为：

yu＝cuhuf

则第u个光源的对角矩阵cu选出了第u个光源的断层合成系统矩阵hu中对应于透光部分的行。

s104：获取成像模型y在完备字典下的稀疏表示：

y＝chψθ

其中，ψ为完备字典，θ为稀疏表示系数。

s105：根据压缩感知理论，获取目标函数ω1：

其中，为稀疏表示系数优化值，a＝cw，w＝hψ，γ为预设的正则项的权重因子。

s106：根据目标函数ω1得到的稀疏表示系数优化值恢复三维物体图像需要满足以下有限等距约束条件：

其中，δs为限制等距常数。

s107：以限制等距常数δs最小为优化目标，获取目标函数ω2：

其中，为所有光源对应编码孔径的对角矩阵cu集合c经过优化后的值。

s108：根据集合c与等效传输函数的关系，将目标函数ω2变换为目标函数ω3：

s109：对目标函数ω3进行连续化，得到目标函数ω。

进一步地，对目标函数ω3进行连续化，得到目标函数ω具体为：

定义如下参数变换式：

其中，πu为光源亮暗表征值，πu＝1表示第u个光源状态为亮，否则状态为暗，tanh为双曲正切函数，ti为编码孔径的二值参数，ti＝1表示第u个光源下编码孔径第v个栅格透光，ti＝0表示第u个光源下编码孔径第v个栅格不透光，编码孔径与探测器的栅格一一对应，λi为光源入射角度的二值参数，λi＝1表示第u个光源下编码孔径第v个栅格在光源照射的圆域内，λi＝0表示第u个光源下编码孔径第v个栅格在光源照射的圆域外，γ为sigmoid函数，(xi,yi)为编码孔径上栅格的坐标；

将上述函数代入目标函数ω3，得到目标函数ω。

需要说明的是，圆域中心(xu,yu)之所以能够表征光源的入射角度，是因为光源的入射角度不同，在编码孔径上的投影也不同，则圆域中心(xu,yu)也对应不同，即圆域中心与光源的入射角度一一对应。

需要说明的是，编码孔径有的栅格透光，有的栅格不透光，从而能够对光源进行调制，更大程度的降低x射线的辐射量，辐射量的降低在医学领域显得极为重要；同时，编码孔径挡住了部分光源，压缩探测器采集的信息量，有利于提高运算效率。

图3为三维物体的四层切片。图4为用随机编码孔径，不经过任何优化后对三维物体重构的结果。图5为用单独优化光源位置后，得到的三维物体的重构结果。图6为用单独优化入射角度后，得到的三维物体的重构结果。图7为用单独优化编码孔径后，得到的三维物体的重构结果。图8为经过联合优化方法得到的三维物体重构结果。其中，图3-图8均包括4幅小图，4幅小图分别为以三维物体的x轴和y轴为水平面，沿z轴方向获得的切片。以上五种方法的重构平均psnr值分别为24.65、26.72、26.93、29.47、29.67，由此可见，采用本实施例的联合优化方法重构出的图像的psnr值有明显的提高，也就是说，本发明能够有效提高三维物体的重构质量。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当然可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。