CN103337065B - 小鼠三维ct图像的非刚性配准方法 - Google Patents

Abstract

小鼠三维CT图像的非刚性配准方法，所述方法包括：A对小鼠三维CT图像的源图像的方向位置进行校正；B用阈值分割算法分割小鼠图像，得到源图像和目标图像的小鼠骨骼，分别提取骨骼切片连通域的质心作为特征点，得到源图像特征点与目标图像特征点；C对步骤B中源图像特征点与目标图像特征点进行匹配，由匹配后的特征点计算薄板样条变换的矩阵，并采用分层变换的方式对小鼠源图像进行薄板样条变换；D对薄板样条变换后的小鼠源图像与目标图像进行精细配准。本发明提取图像特征点时减少了特征点的数量并且采用分层变换的方式，克服了配准时间过长的缺点，在配准图像时基本能够实现自动配准，配准三维小鼠图像更加方便。

Description

小鼠三维CT图像的非刚性配准方法

技术领域

本发明涉及医学图像处理领域，更进一步涉及医学图像的配准方法。本发明可以用于小鼠三维CT图像的配准。

背景技术

医学图像配准在医学图像分析和临床实验中都有着重要的作用，可以用于图像分割、图像融合、手术导航等技术。在医学实验中，小鼠常作为研究的对象，在分析和处理小鼠三维图像前，图像配准是很有必要的。但是小鼠的身体结构比较复杂，不同的图像之间经常有不规则的非刚性形变，而且软组织与骨骼的形变程度是不一样的，在配准的时候需要将软组织与骨骼区别对待，否则不能得到很好的配准结果。传统的刚性配准方法并不能很好的完成这种配准任务，必须采用非刚性配准方法，但是这样会造成配准时间太长，实用性不强等问题。

XiaLi等人在“AutomaticnonrigidregistrationofwholebodyCTmiceimages.MedicalPhysics35(4)20081507–1520.”中提出一种小鼠整体三维CT图像的配准方法。该方法首先利用基于仿射变换的配准方法对小鼠在图像中的相对位置进行校正，然后利用点配准算法对小鼠图像进行粗配准，最后利用一种非刚性配准算法对小鼠图像进行精细配准。该方法配准后的结果有非常高的配准精度，并且能够实现不同小鼠个体图像的配准。该方法的不足之处有二：其一是配准耗费的时间太长，没有太大的实用性；其二是不能用于小鼠身体形变太大的情况。

MartinBaiker等人在“Atlas-basedwhole-bodysegmentationofmicefromlow-contrastMicro-CTdata.”中提出用小鼠三维CT图像与小鼠图谱配准的方法。该方法首先对小鼠图谱中的骨骼按照小鼠的解剖学结构手工分割每个骨头关节，在配准是先将小鼠CT图像的骨骼与图谱的骨骼按每个关节进行配准，然后配准小鼠CT图像与图谱的皮肤表面，得到的变换应用于整个小鼠的软组织。该方法的配准结果有很高的配准精度，并且能够应用于小鼠器官的精确分割。该方法的不足之处：小鼠图谱骨骼的手工分割一般难以实现，配准的过程复杂。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种三维CT图像非刚性配准方法，用于实现小鼠三维CT图像的配准。

为实现上述目的，本发明1.小鼠三维CT图像的非刚性配准方法，所述方法包括：

A对小鼠三维CT图像的源图像的方向位置等进行校正;

B用阈值分割算法分割小鼠图像，得到源图像和目标图像的小鼠骨骼，分别提取骨骼切片连通域的质心作为特征点，得到源图像特征点与目标图像特征点；

C对步骤B中源图像特征点与目标图像特征点进行匹配，由匹配后的特征点计算薄板样条变换的矩阵，并采用分层变换的方式对小鼠源图像进行薄板样条变换。

D最后对薄板样条变换后的小鼠源图像与目标图像进行精细配准。

在上述技术方案的基础上，所述步骤C中匹配的算法采用TPS-RPM点匹配算法

在上述技术方案的基础上，所述步骤C中匹配的算法采用ICP点匹配算法。

本发明还提供小鼠三维CT图像的非刚性配准方法，包括如下步骤：

1）源图像刚性校准步骤，其包括：

对源图像进行刚性变换，使源图像中小鼠的头尾方向，小鼠仰卧姿势和小鼠的主轴线与目标图像一致；

（2）特征点提取步骤，其包括：

2a）采用阈值分割方法分割源图像和目标图像，得到源骨骼图像和目标骨骼图像；

2b）计算目标骨骼图像中每个与小鼠的主轴线垂直切面中连通域的面积，计算面积大于设定阈值的连通区域的质心，得到目标图像特征点；

2c）计算源骨骼图像中每个与小鼠主轴线垂直切面中连通域的面积，计算面积大于设定阈值的连通区域的质心，得到源图像特征点；

（3）特征点匹配步骤，其包括：

匹配源图像特征点与目标图像特征点，得到匹配特征点；

（4）薄板样条变换步骤，其包括：

4a）对源图像进行采样间隔为4*4*4的下采样，得到下采样图像；

4b）利用源图像特征点和匹配特征点，对下采样图像进行薄板样条变换，得到配准采样图像；

4c）计算配准采样图像在薄板样条变换时每个坐标的偏移量，得到薄板样条偏移场；

4d）采用三次线性插值方法对薄板样条偏移场进行采样间隔为4*4*4的上采样插值，得到配准偏移场；

4e）利用配准偏移场对源图像进行变换，得到点配准图像；

（5）精细配准步骤，其包括：

对点配准图像和目标图像进行精细配准，得到灰度配准图像，配准过程结束。

在上述技术方案的基础上，步骤（3）中匹配源图像特征点与目标图像特征点采用TPS-RPM点匹配算法，其按如下步骤进行：

第一步：选取模拟退火算法的退火系数，源图像特征点作为待匹配特征点；

第二步：利用下式计算待匹配特征点与目标图像特征点的对应矩阵：

$M=\frac{1}{t}\exp (-\frac{{(x-y)}^T(x-y)}{2t})$

其中，M表示待匹配特征点与目标图像特征点的匹配矩阵；t表示退火系数；exp表示指数运算；x表示目标图像特征点的坐标；y表示待匹配特征点的坐标；T表示矩阵的转置运算；

第三步：计算待匹配特征点与目标图像特征点之间薄板样条变换的核矩阵，非刚性变换矩阵和刚性变换矩阵；

第四步：按照下式重新计算待匹配特征点的坐标：

X=x*d+K*c

其中，X表示重新计算后的待匹配特征点的坐标；x表示重新计算前的待匹配特征点的坐标；*表示矩阵相乘运算；d表示薄板样条变换的刚性变换矩阵；+表示矩阵相加运算；K表示薄板样条变换的核矩阵；c表示薄板样条变换的非刚性变换矩阵；

第五步：判断退火系数是否低于设定的阈值；如果是，待匹配特征点即为匹配特征点，结束匹配过程，否则按下式计算退火系数：

t=0.93*t₀

其中，t表示退火系数；t₀表示前一次迭代的退火系数；并返回第二步。

在上述技术方案的基础上，所述步骤（5）利用B样条非刚性配准算法对点配准图像和目标图像进行精细配准，其具体步骤如下：

第一步：将点配准图像作为灰度预配准图像；

第二步：对点配准图像和目标图像进行随机采样，每幅图像得到2000-5000个采样点，计算两幅图像采样点之间的归一化互信息；

第三步：利用灰度预配准图像和目标图像采样点之间的归一化互信息，采用自适应梯度下降优化算法对灰度预配准图像与目标图像之间B样条变换参数进行梯度优化，得到最优B样条变换参数；

第四步：利用最优B样条变换参数对灰度预配准图像进行三次均匀B样条变换，采用一阶B样条差值方法对变换后的图像进行插值，将插值后的图像作为灰度预配准图像；

第五步：判断梯度优化时步长是否小于人工设定的最小步长，如果是，灰度预配准图像即是灰度配准图像，非刚性配准结束；如果否，继续返回第二步。

与现有技术相比较，本发明在提取图像特征点时通过采样减少了源图像的特征点的数量，在薄板样条变换时采用分层变换的方式，克服了现有技术中配准时间过长的缺点，使得本发明在保证配准精度的前提下，大大缩减了配准的时间，使用本发明配准三维图像更加快捷；

且在配准图像时基本能够实现自动配准，克服了现有技术中预处理复杂的缺点，使本发明易于实现，使用本发明配准三维图像更加方便。

附图说明

图1为本发明小鼠三维CT图像的非刚性配准方法流程图；

图2本发明Micro-CT扫描的目标图像；

图3为本发明Micro-CT扫描的需要配准的源图像；

图4为本发明目标图像阈值分割的小鼠骨骼图像；

图5为本发明源图像阈值分割的小鼠骨骼图像；

图6为本发明处理后的配准图像；

图7为本发明为配准结果图像阈值分割的小鼠骨骼图像。

具体实施方式

下面结合图1对本发明的具体实施步骤做进一步的详细描述。

步骤1.对源图像刚性校准

在刚性校准时，首先校正图像中小鼠的头部和尾部方向一致，然后校正小鼠仰卧姿势一致，最后小鼠的中轴线一致对齐。

步骤2.提取源图像和目标图像的特征点

（1）选取阈值，采用阈值分割方法提取源图像和目标图像中小鼠的骨骼，得到源骨骼图像和目标骨骼图像，在本实施例中，阈值的大小是105；

（2）计算目标骨骼图像中每个与小鼠的主轴线垂直切面连通域的面积，计算面积大于10的连通区域的质心，得到目标图像特征点；

（3）计算源骨骼图像中每个与小鼠的主轴线垂直切面连通域的面积，计算面积大于10的连通区域的质心，得到源图像特征点。

步骤3.匹配源图像和目标图像的特征点

将源图像特征点与目标图像特征点进行匹配，得到匹配特征点；其算法可采用TPS-RPM点匹配算法和ICP点匹配算法，其中，采用TPS-RPM点匹配算其具体步骤如下：

第一步：选取模拟退火算法的退火系数，源图像特征点作为待匹配特征点；

第二步：利用下式计算待匹配特征点与目标图像特征点的对应矩阵：

$M=\frac{1}{t}\exp (-\frac{{(x-y)}^T(x-y)}{2t})$

其中，M表示待匹配特征点与目标图像特征点的匹配矩阵；t表示退火系数；exp表示指数运算；x表示目标图像特征点的坐标；y表示待匹配特征点的坐标；T表示矩阵的转置运算；

第三步：计算待匹配特征点与目标图像特征点之间薄板样条变换的核矩阵，非刚性变换矩阵和刚性变换矩阵；

第四步：按照下式重新计算待匹配特征点的坐标：

X=x*d+K*c

其中，X表示重新计算后的待匹配特征点的坐标；x表示重新计算前的待匹配特征点的坐标；*表示矩阵相乘运算；d表示薄板样条变换的刚性变换矩阵；+表示矩阵相加运算；K表示薄板样条变换的核矩阵；c表示薄板样条变换的非刚性变换矩阵；

第五步：判断退火系数是否低于设定的阈值；如果是，待匹配特征点即为匹配特征点，结束匹配过程，否则按下式计算退火系数：

t=0.93*t₀

其中，t表示退火系数；t₀表示前一次迭代的退火系数；并返回第二步。

步骤4.对源图像进行薄板样条变换

第一步：对源图像进行采样间隔为4*4*4的下采样，得到下采样图像；

第二步：利用源图像特征点和匹配特征点，对下采样图像进行薄板样条变换，得到配准采样图像；

第三步：计算下采样图像在薄板样条变换时每个坐标的偏移量，得到薄板样条偏移场；

第四步：采用三线性插值对薄板样条偏移场进行采样间隔为4*4*4的上采样，得到配准偏移场；

第五步：利用配准偏移场对源图像进行变换，得到点配准图像。

步骤5.对薄板样条变换后的源图像进行精细配准

利用改进的B样条非刚性配准算法对点配准图像和目标图像进行精细配准。其具体步骤如下：

第一步：将点配准图像作为灰度预配准图像；

第二步：对点配准图像和目标图像进行随机采样，每幅图像得到2000-5000个采样点，计算两幅图像采样点之间的归一化互信息；

第三步：利用灰度预配准图像和目标图像采样点之间的归一化互信息，采用自适应梯度下降优化算法对灰度预配准图像与目标图像之间B样条变换参数进行梯度优化，得到最优B样条变换参数；

第四步：利用最优B样条变换参数对灰度预配准图像进行三次均匀B样条变换，采用一阶B样条差值方法对变换后的图像进行插值，将插值后的图像作为灰度预配准图像；

第五步：判断梯度优化时步长是否小于设定的最小步长，如果是，灰度预配准图像即是灰度配准图像，非刚性配准结束；如果否，继续返回第二步。

下面结合附图2至附图7对本发明的配准结果做进一步的描述。

1.仿真条件

在CPU为Intelcorei5、3.1GHZ、内存8G、Window7操作系统上使用MATLAB2009a，VS2010，3DMed4.1进行仿真

2.仿真内容

同一只小鼠由于在Micro-CT扫描时的时间不同，在扫描时的小鼠的姿势出现了变化，致使小鼠的骨骼和内脏器官出现了不可预知的变形。通过本发明对出现变形的图像与另一幅图像进行配准，使两幅图像能够共同分析和处理。

图2为Micro-CT扫描的目标图像；图3为Micro-CT扫描的需要配准的源图像。图4为目标图像阈值分割的小鼠骨骼；图5为源图像阈值分割的小鼠骨骼。图6为本发明处理后的配准图像；图7为配准图像阈值分割的小鼠骨骼。

3.仿真结果

根据图7可以看出本发明对小鼠发生姿势形变时的配准结果是很准确的，下面我们通过计算配准前后源图像和目标图像小鼠骨骼之间的相关系数和Dice系数来定量分析本发明的结果：

	配准前	配准后
			相关系数	0.179	0.958
Dice系数	0.172	0.893

本发明的另一个优点是配准三维图像的时间大大的缩减，在仿真时所用的时间是：

点匹配	薄板样条变换	B样条非刚性配准	总时间
				40.6s	206.33s	72s	318.93s

以上仿真结果可以看出，采用三维CT图像非刚性配准方法，能够很好的对三维CT小鼠图像非刚性变形进行配准，在保证一定配准精度的前提下，缩短了配准的时间。

Claims (3)

1.小鼠三维CT图像的非刚性配准方法，其特征在于，包括如下步骤：

(1)源图像刚性校准步骤，其包括：

对源图像进行刚性变换，使源图像中小鼠的头尾方向，小鼠仰卧姿势和小鼠的主轴线与目标图像一致；

(2)特征点提取步骤，其包括：

2a)采用阈值分割方法分割源图像和目标图像，得到源骨骼图像和目标骨骼图像；

2b)计算目标骨骼图像中每个与小鼠的主轴线垂直切面中连通域的面积，提取面积大于设定阈值的连通区域的质心，得到目标图像特征点；

2c)计算源骨骼图像中每个与小鼠主轴线垂直切面中连通域的面积，提取面积大于设定阈值的连通区域的质心，得到源图像特征点；

(3)特征点匹配步骤，其包括：

匹配源图像特征点与目标图像特征点，得到匹配特征点；

(4)薄板样条变换步骤，其包括：

4a)对源图像进行采样间隔为4*4*4的下采样，得到下采样图像；

4b)利用源图像特征点和匹配特征点，对下采样图像进行薄板样条变换，得到配准采样图像；

4c)计算配准采样图像在薄板样条变换时每个坐标的偏移量，得到薄板样条偏移场；

4d)采用三次线性插值方法对薄板样条偏移场进行采样间隔为4*4*4的上采样插值，得到配准偏移场；

4e)利用配准偏移场对源图像进行变换，得到点配准图像；

(5)精细配准步骤，其包括：

对点配准图像和目标图像进行精细配准，得到灰度配准图像，配准过程结束。

2.如权利要求1中所述的小鼠三维CT图像的非刚性配准方法，其特征在于：步骤(3)中匹配源图像特征点与目标图像特征点采用TPS-RPM点匹配算法，其按如下步骤进行：

第一步：选取模拟退火算法的退火系数，将源图像特征点作为待匹配特征点；

第二步：利用下式计算待匹配特征点与目标图像特征点的对应矩阵：

$M=\frac{1}{t}\exp (-\frac{{(x-y)}^T(x-y)}{2t})$

其中，M表示待匹配特征点与目标图像特征点的匹配矩阵；t表示退火系数；exp表示指数运算；x表示目标图像特征点的坐标；y表示待匹配特征点的坐标；T表示矩阵的转置运算；

第三步：计算待匹配特征点与目标图像特征点之间薄板样条变换的核矩阵，非刚性变换矩阵和刚性变换矩阵；

第四步：按照下式重新计算待匹配特征点的坐标：

X＝x*d+K*c

其中，X表示重新计算后的待匹配特征点的坐标；x表示重新计算前的待匹配特征点的坐标；*表示矩阵相乘运算；d表示薄板样条变换的刚性变换矩阵；+表示矩阵相加运算；K表示薄板样条变换的核矩阵；c表示薄板样条变换的非刚性变换矩阵；

第五步：判断退火系数是否低于设定的阈值；如果是，待匹配特征点即为匹配特征点，结束匹配过程，否则按下式计算退火系数：

t＝0.93*t₀

其中，t表示退火系数；t₀表示前一次迭代的退火系数；并返回第二步。

3.如权利要求1中所述的小鼠三维CT图像的非刚性配准方法，其特征在于：所述步骤(5)利用B样条非刚性配准算法对点配准图像和目标图像进行精细配准，其具体步骤如下：

第一步：将点配准图像作为灰度预配准图像；

第二步：对点配准图像和目标图像进行随机采样，每幅图像得到2000-5000个采样点，计算两幅图像采样点之间的归一化互信息；

第三步：利用灰度预配准图像和目标图像采样点之间的归一化互信息，采用自适应梯度下降优化算法对灰度预配准图像与目标图像之间B样条变换参数进行梯度优化，得到最优B样条变换参数；

第四步：利用最优B样条变换参数对灰度预配准图像进行三次均匀B样条变换，采用一阶B样条差值方法对变换后的图像进行插值，将插值后的图像作为灰度预配准图像；

第五步：判断梯度优化时步长是否小于设定的最小步长，如果是，灰度预配准图像即是灰度配准图像，非刚性配准结束；如果否，继续返回第二步。