一种半自动分段的顺序优化配准方法与流程

本发明属于医学图像配准领域，具体涉及一种半自动分段的顺序优化配准方法。

背景技术：

传统的骨科手术中，病灶需要主刀医生凭借临床经验，结合术前ct或x光图像，在大脑中构建骨结构的三维模型来定位。过度依赖医生的临床经验，缺乏的医学影像依据。随着计算机辅助手术系统的发展，出现了将术前的ct数据三维重建为体数据，结合术中x光图像进行手术的导航系统。既为手术提供了客观的影像的依据，又使整个手术流程变得可视化。

手术导航的研究关键就是多模态医学图像的配准问题。2d-3d医学图像配准是医学图像配准中的热点。分为基于特征的配准方法，以及基于灰度的配准方法。基于特征的配准方法仅使用了少量的特征，配准时间较短，但需要侵入式标记特征点；基于灰度的配准方法更加稳定、精确，但耗时长。

传统方法均将脊椎结构看作刚体运动，但椎体骨组织各段都存在着相对运动。x光与ct的拍摄的时间不同，所以各段椎体骨组织间存在着相对运动，会影响配准精度。

技术实现要素：

本发明的目的是提供一种半自动分段的顺序优化配准方法，以现有技术中存在的由于骨组织相对运动产生配准误差，分段后配准数量增加导致耗时长的问题。

一种半自动分段的顺序优化配准方法，包括以下步骤：

步骤1，输入ct切片数据，利用css算法提取椎体骨组织的角点来确定初始分割范围；

步骤2，对在初始分割范围具有椎体骨组织的所有ct切片数据进行角点数据处理，得到最终分割范围；

步骤3，利用最终分割范围分割所有ct切片数据得到骨组织的ct体数据，并生成drr图像；

步骤4，将得到的drr图像与被最终分割范围处理过的x光图像进行配准，依次求得6个配准参数rx,ry,rz,tz,tx,ty的初始解；

步骤5，利用步长加速法对6个配准参数的初始解进行优化，得到配准参数的最终解。

进一步地，步骤1提取椎体骨组织的角点确定初始分割范围包括通过选取序列号为ct切片数量1/2的切片确定椎体骨组织初始分割范围。

进一步地，步骤1利用css算法提取椎体骨组织的角点包括：

将ct切片利用canny算子轮廓线提取后，对轮廓线利用不同尺度的高斯核进行平滑，接着对最大尺度进行曲率计算得到角点，并且在其他尺度上跟踪验证角点，尺度为σ的曲线的曲率为:

进一步地，步骤1将最大尺度曲率绝对值大于200且小于10000的点，或比相邻的曲率极小点大两倍以上的作为角点候选点。

进一步地，步骤1所述的初始分割范围左边缘x＝xmin，右边缘x＝xmax，上边缘y＝ymin，下边缘y＝ymax所围成的区域。

进一步地，步骤2进行角点数据处理包括将不超出范围5mm的角点视为范围内的点，将右边缘上不超出3mm的点视为范围内的点。

进一步地，步骤4依次求得6个配准参数rx，ry，rz，tz，tx，ty的初始解包括采用步长加速算法求解旋转参数rx，ry，采用边缘方向求解旋转参数rz，采用相似三角形法求解平移参数tz，采用相位相关法求解平移参数tx，ty。

进一步地，步骤5利用步长加速法对6个配准参数的初始解进行优化，得到配准参数的最终解包括：

在得到了全部配准参数初始解后，按照参数重新生成drr图像，首先对rx，ry进行步长加速优化算法，初始步长为3开始，减少为0.005时搜索结束，搜索rx，ry最优值的[-5，5]的范围，接着对rz在反方向[-1，0]的范围以0.3为步长开始进行搜索，每轮步长减少1/3，当步长小于0.01时停止，接着重新求解更新tz，tx，ty的值，对tz，tx，ty在反方向[-1，1]的范围以0.3为步长开始进行搜索，每轮搜索步长减少1/3，当步长小于0.01时停止。

与现有技术相比，本发明的有益效果是：

(1)本发明所采用的半自动分段方法，克服了骨组织间发生相对运动而产生实验误差的缺点；

(2)本发明所采取的顺序优化方法，克服分段后配准的数量增加导致耗时增加，从而不满足手术要求的问题，将传统配准的时间复杂度o(n⁶)降为了o(6n)；

(3)本发明采用的drr生成方法，减少了生成drr图像的耗时以及生成drr图像的数量；

(4)本发明所采取求解6个配准参数的方法与现有技术相比，减少了配准参数的求解时间并提高了精确度。

附图说明

图1为本发明分段顺序优化配准的流程图；

图2为准旋转参数rz时用到的边缘方向直方图的性质；

图3为drr生成过程体数据沿z轴平移产生的相似三角形；

图4为css角点检测算法使用的开始和结束切片；

图5为利用css角点检测算法对x光数据图像配准区域提取的结果；

图6为本发明实施例得到的各段配准结果。

具体实施方式

下面结合附图与本发明实施例对本发明的技术方案进行更为清楚，完整的描述。

本实施例中采用北医三院脊椎ct图像以及x光图像，截取了病灶部位ct以及x光图像，总共包括三段骨组织。ct图像体数据大小为512×359×238，层厚为0.7mm。x光图像大小为1024×803。分辨率均为0.28346mm/pixel。

本方法所述的体数据是三维图形可视化中的术语，体数据由体素组成。体素是基本体积元素，也可以理解为三维空间内的具有排列和颜色的点或一小块区域。本实施例所采用的ct图像体数据是将所有切片叠加未进行分割的类似于长方体数据，512×359×238指的是图像横向512个像素点，纵向359个像素点，238个切片，x光图像大小为1024×803，即横向1024个像素点，纵向803个像素点。

一种半自动分段的顺序优化配准方法，包括以下步骤：

步骤1，输入ct切片数据，利用css算法提取椎体骨组织的角点确定初始分割范围；

由于ct切片是按照顺序排列的，所以每一个切片都有对应的序列，选取序列号为ct切片数量1/2的切片，即当ct切片的数量为238片，我们此时选取第119号切片，这是因为此时椎体骨组织的面积最大，可以确定椎体骨组织的初始范围。

本发明采用了css角点提取算法进行椎体骨组织的角点定位，通过四个角点来确定骨组织的范围。角点是椎体骨组织边缘曲线上局部曲率最大的点，角点的具体确定是通过将ct切片数据利用canny算子轮廓线提取后，对轮廓线利用不同尺度的高斯核进行平滑，接着对最大尺度进行曲率计算得到角点，并且在其他尺度上跟踪验证角点。尺度为σ的曲线的曲率为:

其中x(u)，y(u)分别为轮廓线上的横坐标与纵坐标的一维函数，u为一维函数的自变量，σ为尺度空间，x′(u，σ)，y′(u，σ)分别为x(u，σ)与y(u，σ)的一阶导数，x″(u，σ)，y″(u，σ)分别为x(u，σ)与y(u，σ)的二阶导数，g(u，σ)为高斯核函数，为卷积运算。

公式(1)为曲率的求解函数，公式(2)、(3)为对横纵坐标的一维函数进行高斯平滑的过程。

本发明将最大尺度曲率绝对值大于200且小于10000的点，或比相邻的曲率极小点大两倍以上的作为角点候选点。将角点候选点在逐渐减少的尺度中寻找极大点，直到最低的尺度。

设各段四个角点的坐标按照顺时针方向分别为(x1，y1)，(x2，y2)，(x3，y3)，(x4，y4)；

x1，x4中选出较小点保存到xmin中，x2，x3中选出较大点保存到xmax中，在上边缘中选择纵坐标最小点保存到ymin中，在下边缘选择纵坐标最大点保存到ymax中，则椎体骨组织的初始范围即x＝xmin，x＝xmax，y＝ymin，y＝ymax所围成的区域。在本实施例中得出的119号切片的三段骨组织初始分割范围分别为(51.86，91.65，3.26，38.81)，(47.29，85.13，43.38，75.35)，(45.34，83.83，80.56，113.18)。

步骤2，对在初始分割范围具有椎体骨组织的所有ct切片数据进行角点数据处理，得到最终分割范围；

如图4，从刚有椎体骨组织出现的的第94号切片开始，到分割范围内的椎体骨组织刚好消失的168号结束，分别使用css算法提取角点并记录角点的坐标值；在椎体骨组织没有完全出现的切片中，如果存在角点ymin≤y≤ymax，则判定该角点为相应段的角点，如果|ymin-y|≤|ymax-y|则为上边缘的角点，否则为下边缘的角点；

在本方法中，对所有ct切片数据进行角点数据处理具体指的是：ct切片数据中每张切片的椎体骨组织之间偏差很小，算法将不超出范围5mm的角点视为范围内的点，将右边缘上不超出3mm的点视为范围内的点。

即如果y≥ymax或y≤ymin，当y-ymax≤5时，则将y的值作为新的ymax，如果存下边缘点的最大纵坐标t≥y，则将t的值作为新的ymax；当ymin-y≤5时，则将y的值作为新的ymin，如果存上边缘点的最小纵坐标t≤y，则将t的值作为新的ymin；

如果x≥xmax或x≤xmin，当x-xmax≤5时，则将x的值作为新的xmax，如果存在右边缘点的最小纵坐标t-y≤3，则将t的值作为新的ymin；当xmin-x≤5时，则将x的值作为新的xmin；

得到的最终分割范围为使用x＝51.21，x＝93.29，y＝3.26，y＝39.79分割第一段，使用x＝45.99，x＝87.74，y＝43.38，y＝76.65分割第二段，使用x＝43.71.21，x＝85.41，y＝79.59，y＝116.44来分割第三段。

步骤3，利用最终分割范围分割所有ct切片数据得到骨组织的ct体数据，并生成drr图像；

遍历到所有骨组织都消失的序列号即168号切片，分别用最终得到的x＝xmin，x＝xmax，y＝ymin，y＝ymax即步骤2中最终分割范围分割94号到168号切片ct切片数据，并按照序列号叠加即为ct体数据。本发明采用bresenham算法生成drr图像，bresenham算法是计算机图形学中为了“显示器(屏幕或打印机)系由像素构成”的这个特性而设计出来的算法。通过3维直线bresenham生成算法可以快速的得到一个体素序列，这个体素序列与光线投射算法中每条x射线所遍历的真实像素点最为接近。在计算投影点ct值的过程中只需要进行少量的浮点运算，不需要进行插值运算和求整运算，大大加快了运算速度。

体数据生成drr图像具体过程为首先确定虚拟点光源的固定位置，然后从光源向外发射虚拟射线，射线穿过3维ct体数据向垂直于射线中轴的投影面进行投射。射线同平面的交点决定了drr图像中各像素点的位置。当射线以一定的步长向ct体数据投射时，可以得到所经过的每个切片上交点的ct值，若交点不在网格位置处，则需要用插值算法计算出该交点的ct值。随着射线的前进，对得到的ct值进行累加，投射结束后，便可以得到一个累加值(累加过程可以加权)。所有的光线投射完毕后，就得到了整个drr图像上每一点的ct累加值，将这些ct累加值映射成像素灰度值，就得到drr图像。

步骤4，将得到的drr图像与被最终分割范围处理过的x光图像进行配准，依次求得6个配准参数rx，ry，rz，tz，tx，ty的初始解；

采用最终分割范围来对x光图像进行处理，即使用步骤2得到的x＝xmin，x＝xmax，y＝ymin，y＝ymax分割本实施例所采用的北医三院脊椎病灶部位的x光图像得到的配准区域作为参考图像，所谓的配准主要是比对处理后的x光图像与步骤3生成的drr图像的相似性，相似性越高配准效果越好。

本发明所采用的配准方法是顺序优化配准算法，并将步长设为2，即隔一张切片进行一次浮点计算；按照不同方法按照顺序求出6个配准参数rx，ry，rz，tz，tx，ty的值，这6个参数分别对应体数据绕x旋转，绕y旋转，绕z旋转，z轴方向平移，x轴方向平移，y轴方向平移，。先求出绕x轴以及y轴的旋转参数rx，ry，接着求出z轴的旋转参数rz，然后求出z轴方向的平移参数tz，最后求出x轴方向以及y轴方向的平移参数tx，ty。以上的求解顺序必须严格单向，这样才能保证各参数间相互影响最小。在得到所有6个参数的初始解后，再对每个参数按照顺序，以不同的步长在精确的范围内进行优化，最终得到配准参数的最终解。

4.1本发明采用了步长加速算法求解绕x，y轴旋转参数rx，ry，步长从30开始，步长小于0.05时搜索结束。为防止局部最优解，进行80次迭代，将每一次的输出结果作为下一次的输入，直到相似性度量不再发生改变为止。采用相似性度量为归一化互信息，计算公式如下：

其中pa(a)，pb(b)为边缘概率分布，pab(a，b)为联合概率分布。

4.2本发明引入了边缘方向求z轴的旋转参数rz。图像边缘的每个像素都对应一个边缘方向，可以将边缘分解为特定方向的像素点，对这些像素点的方向进行统计得到边缘方向直方图。通过将drr图像与x光图像进行小波分解，得到各尺度图像。再对各尺度图像进行轮廓线提取。通过如下公式来计算梯度方向，其中gh(x，y)为水平方向的梯度，gv(x，y)为垂直方向的梯度，θ(x，y)为梯度方向。f(x，y)为对应点的灰度值，对任意x，y有：

gh(x，y)＝f(x+1，y)-f(x-1，y)(5)

gv(x，y)＝f(x，y+1)-f(x，y-1)(6)

使用边缘方向直方图法可以解决drr图像与x光图像光照不同的问题，排除配准参数rz时配准参数tx，ty的影响。本文采用360个单位的边缘方向直方图，参见图2，分别求出x光图像与drr图像的边缘方向直方图，将直方图中峰值作为主方向，rz为两直方图的偏移量。

4.3本发明利用相似三角形法求解z轴方向的平移参数tz。通过css算法分别求出drr图像与x光影像间的左上角与右下角角点的距离ddrr和dxray，利用公式(8)来计算缩放比例。

drr图像的产生过程正好构成了一个相似三角形。当在z轴上放生平移的时，产生了2个相似三角形。我们将其放入一个二维平面中，参见图3。通过相似三角形原理即可得到配准参数tz。

从三角形的左侧顶点向右侧的边做垂线，这条线为光源点到drr图像影像的距离设为m。将体数据对角线的长度设为x，将光源到体数据的距离设为t。由相似三角形的原理可得：

由(8)可得：

将(9)(10)(11)相结合得到如下公式：

公式中的t是在生成drr时设置的，r是通过公式(8)计算得出。

4.4本发明采用相位相关法求解x轴方向以及y轴方向的平移参数tx，ty。本发明只利用它求解平移参数。因为第一步存在得到的峰值不能可靠地对应正确旋转角度的位置。

得到配准参数rx，ry，rz，tz后，重新生成drr图像。分别对drr图像与x光图像进行傅里叶变化(公式13)，并求出互功率谱(公式14)。接着对互功率谱进行逆傅里叶变化，得到相关值的分布，分布中的的峰值即为tx，ty。

f1(u，v)为图像的傅里叶变化，f2(u，v)为平移图像的傅里叶变化，x0，y0为平移量

为f2(u，v)的共轭。

步骤5，利用步长加速法对6个配准参数的初始解进行优化，得到配准参数的最终解。

本发明得到了全部配准参数初始解后，按照参数重新生成drr图像，但此时生成drr图像得到的配准结果具有较大的误差，所以需要对生成drr图像的6个配准参数进行优化，以获得更小的误差，优化的具体过程为首先对rx，ry进行步长加速优化算法，初始步长为3开始，减少为0.005时搜索结束，搜索rx，ry最优值的[-5，5]的范围，接着对rz在反方向[-1，0]的范围以0.3为步长开始进行搜索，每轮步长减少1/3，当步长小于0.01时停止，接着重新求解更新tz，tx，ty的值，对tz，tx，ty在反方向[-1，1]的范围以0.3为步长开始进行搜索，每轮搜索步长减少1/3，当步长小于0.01时停止。见图6，最终得到的最终的配准结果，以tx，ty，tz，rx，ry，rz顺序分别为(10.74，-6.34，5.02，5.21，6.88，10.05)，(7.88，-5.35，4.89，3.91，4.34，6.98)，(6.79，3.77，5.05，5.16，4.22，-4.04)。

本实施例选取第一段骨组织进行评价顺序优化算法的性能，将骨组织的初始位置设为(15，10，10，15，12，8)。进行50次实验。分别记录每一次的配准结果的值，使用如下公式进行评价：

其中xi为每次的配准误差，n为实验次数，为配准误差的均值，s标准差。

利用公式(15)(16)对50次配准实际值与理想值的误差进行计算，理想值为(15，10，10，15，12，8)，实际值为每次实验的得到的配准结果，将两个值相减即为误差值。

得到tx，ty，tz的误差均值(标准差)为-0.37，-0.45，0.16(±0.15，±0.16，±0.25)，rx，ry，rz的误差均值(标准差)为-0.15，-0.18，-0.05(±0.72，±0.87，±0.15)。配准时间均值(标准差)为28.32(±9.45)s。可以看出本优化配准方法具备更好的配准结果。

本发明的实施例公布的是其较好的实施方式，但并不限于此。本领域的技术人员能够容易地根据上述实施例，理解本发明的核心思想，只要不脱离本发明的技术方案的基础的变形或替换，都在本发明的保护范围之内。