[发明专利]一种基于网格匹配的Demons算法的多期CT图像配准方法

摘要

一种基于网格匹配的Demons算法的多期CT图像配准方法，计算机X射线断层扫描设备对人体腹部进行CT增强扫描；获取动脉期和肝静脉期的腹部ＣＴ图像，并对获取的腹部CT图像进行阈值处理，阈值处理后的静脉期图像作为参考图像，动脉期图像作为浮动图像；对多期ＣＴ图像进行基于梯度阈值掩膜的刚性预配准；将刚性预配准的浮动图像作为新的浮动图像，参考图像不变，对其进行网格匹配；将网格匹配后的参考图像作为新的参考图像，刚性预配准的浮动图像作为新的浮动图像，对其进行基于微分同胚Demons方法的形变配准。基于网格匹配的微分同胚Demons形变配准在形变配准前加入网格匹配，使得参考图像和浮动图像的灰度一致，克服灰度差这一难题。

主权项

一种基于网格匹配的Demons算法的多期CT图像配准方法，包括以下步骤：步骤1：计算机X射线断层扫描设备对人体腹部进行CT扫描；步骤2：获取腹部CT图像，每组腹部CT图像包含若干张腹部断层扫描图像；步骤3：阈值处理得到参考图像R，浮动图像M；步骤4：基于梯度阈值掩膜的刚性预配准，得到预配准后的浮动图像；步骤4.1：计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图像序列，初始金字塔层数为0，变换参数默认为0；步骤4.2：取当前层金字塔参考图像和浮动图像，平移变换参数为上一层变换参数的2倍，旋转参数和上一层变换参数一致；步骤4.3：计算参考图像的梯度阈值掩膜图像并在掩膜上随即采样得到一个采样点集合；步骤4.4：对采样点进行刚性变换；步骤4.5：对刚性变换后的采样点进行线性灰度差值，并更新联合概率分布图，然后计算出图像R和M的互信息；步骤4.6：当互信息没有达到最大时，执行步骤4.7；否则执行步骤4.8步骤4.7：梯度优化器根据互信息的梯度值生成下一次刚想变换参数，更新变换矩阵，返回步骤4.4；步骤4.8：如果层数达到最大值，那么循环结束，输出浮动图像Mp；如果层数没有达到最大值，那么让层数的值加1，重新执行步骤4.2，再次进入循环，直至循环到层数达到最大值，输出浮动图像Mp；步骤5：网格匹配，得到网格匹配后的参考图像，即新的参考图像；步骤5.1：以步骤4所得图像Mp和源图像R作为输入，计算网格联合直方图H(i,j)；步骤5.2：对于H(i,j)中的任意i，找出概率最大的配对(i,j)，在参考图像中用j替换i；步骤5.3：如果网格的遍历没有完成，返回步骤5.1；否则输出变换图像R_M；步骤5.4：将源图像R和图像Mp划分网格，在对应的网格上进行直方图匹配，从而对源图像R进行变换，得到网格直方图匹配后的图像R_H；步骤5.4：将源图像R通过网络直方图匹配所得到的变换图像和变换图像R_M进行图像融合，输出变换图像R_I；步骤6：微分同胚Demons形变配准,得到最终的结果图像；步骤6.1：将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算；步骤6.2：输出步骤6.1计算计算结果，得到最终结果图像M_R。所述步骤4.1计算参考图像和浮动图像的两层图像金字塔，用简单的高斯金字塔生成图 像序列，具体步骤如下：步骤4.1.1：设金字塔层数为L，则使用方差为(L/2)²高斯平滑图像；步骤4.1.2：相应的层级对图像进行重采样以缩小到相应大小。所述步骤4.3计算参考图像的梯度阈值掩膜图像并得到采样点，具体步骤如下：步骤4.3.1：首先计算参考图像的梯度图像；步骤4.3.2：将梯度图像进行高斯平滑，即保证梯度阈值掩膜图像中包含部分梯度为0的点；步骤4.3.3：然后设置阈值为Th，该值一般的范围是[1,10]，即阈值处理后将梯度为0的点去除；如果图像中的点的梯度值小于设定阈值Th，则将该点的梯度阈值掩膜设为0；如果图像中的点的梯度值大于等于设定阈值Th，则将该点的梯度阈值掩膜设为1，得到梯度阈值掩膜图像和采样点；步骤4.3.4：在掩膜上随机采样若干数量的点。所述步骤4.7梯度优化器根据互信息的梯度值生成下一次刚想变换参数具体步骤如下：步骤4.7.1：梯度值生成下一次的变换参数的计算公式如下，对于旋转参数，下一次的变换参数X⁽ⁿ⁺¹⁾的计算公式表示为：X⁽ⁿ⁺¹⁾＝γ_n▽F(X⁽ⁿ⁾)对于平移参数，计算公式表示为：X⁽ⁿ⁺¹⁾＝X⁽ⁿ⁾‑γ_n▽F(X⁽ⁿ⁾)互信息的梯度▽F(X⁽ⁿ⁾)表示为其中μ表示X⁽ⁿ⁾，p(l,k；μ)表示联合概率函数，k为参考图像的灰度级，l为浮动图像的灰度级。此处γ_n＞0为松弛因子，具体设置为0.9；梯度下降优化器搜索参数时的步长设置为第一层：(0.001,3)，第二层：(0.0001,1.5)，迭代200次。步骤4.7.2：更新变换矩阵，返回步骤4.4；所述步骤6.1将参考图像R_I和浮动图像Mp作为微分同胚的Demons算法的输入，进行计算，微分同胚的Demons算法的具体步骤如下：步骤6.1.1：预先设定迭代次数N和初始形变t为零向量(t必须属于李群空间，因此初始化为0)；步骤6.1.2：根据下面的公式可得到稠密度形变场u，u属于李代数空间，r为参考图像RI，g为浮动图像Mp；步骤6.1.3：对u做一次高斯卷积：u←K_fluid*u，产生流体映射模型的规则化效果，规则化后仍然属于李群空间；步骤6.1.4：c←tοexp(u)，则c仍然属于李群空间；步骤6.1.5：对c做高斯卷积并更新t，c←tοexp(u)；步骤6.1.6：判断迭代次数是否已经为N，如小于N，则返回(2)；否则输出最优变换t^opt。所述步骤6.1.4计算c的复合运算具体步骤：步骤6.1.4.1：令v＝2^‑nu，选择一个n值，使得max||2^‑nu(p)||≤ε，ε＝0.5；步骤6.1.4.2：对于每一个像素位置，计算v，然后令v＝Id+v；步骤6.1.4.3：对v做n次复合运算，则exp(u)＝vοv...οv。