基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法

摘要

本发明公开了一种基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法，其基本步骤是：首先，将Digimouse模型设定为参考鼠解剖结构图谱，同时将Digimouse对应的XCT图像设定为参考图像；其次，对目标鼠进行XCT成像获取目标图像并进行预处理；然后，利用非刚度图像配准技术构建参考图像至目标图像的配准映射矩阵；最终，将配准映射矩阵作用于参考鼠解剖结构图谱，构建出目标鼠解剖结构图谱，从而实现对目标鼠各个组织器官的标定。本发明所使用配准方法精度高，能够通过图像配准的方法简单易行的实现小鼠组织器官在小鼠XCT图像上的标识；本发明提出方法同样适用于其他医学应用领域如人脑结构学研究，即可以利用标准人脑解剖学图像实现目标人脑的解剖结构的获取。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及非刚度图像配准与目标鼠组织解剖结构的有 限元标定。

背景技术

目前荧光分子层析成像(Fluorescence Molecular Tomography,FMT)方法通常使用的均匀 光学结构背景将会引入光子输运建模上的显著误差，有效的光学结构先验信息对FMT重建 精度和灵敏度的提升具有重要意义。光学结构的建立与解剖结构信息的获取密切相关：一 方面它是赋予了解剖结构中各区域相关光学参数特性的物理几何信息，另一方面它是各区 域光学参数在体获取的先决条件^[1,2]。

常见的解剖成像模态用于光学结构获取均存在一定的局限性。高场小动物核磁共振成 像(Micro Magnetic Resonance Imaging,μMRI)对软组织具有高灰度分辨率，能够利用图像方 法获取各个软组织器官所在区域。Dhenain等人利用该技术对多个小鼠胚胎进行成像，获取 了在不同发育时期小鼠胚胎的解剖结构图谱(Atlas)^[3]。Segars等人利用多组成年小鼠核磁共 振数据建立四维MOBY模型，模拟小鼠在心跳呼吸等生理过程中的动态解剖结构图谱^[4]。然 而，μMRI成本极高，成像耗时较长，限制了其在小鼠FMT实验中的应用。X射线计算机断 层成像(X-ray Computed Tomography,XCT)作为一种常用的解剖学成像模式，其成像速度较 快且造价适中。然而X射线对软组织的分辨率较低，利用XCT图像分割软组织具有较大难度， 若辅以其它成像模态，则可更为精确的获取生物组织体解剖结构。Dogdas等人利用XCT、 正电子发射计算机断层成像及冷切片技术对小鼠进行多模态成像，并以此为依据建立 Digimouse模型，获取了小鼠的精确解剖结构图谱^[5]。此方法能够准确地获取小鼠的解剖结 构信息，为相关研究提供了重要的参考意义，但该方法所用成像系统复杂，实验过程繁琐， 成本较高，同样不利于在常规的小鼠FMT实验中采用。

[参考文献]

[1]L.-H.Wu,W.-B.Wan and X.Wang et al,"Shape-parameterized diffuse optical  tomography holds promise for sensitivity enhancement of fluorescence molecular tomography," Biomedical Optics Express 10,3640-3659(2014)。

[2]L.-H.Wu,H.-J.Zhao and X.Wang et al,"Enhancement of fluorescence molecular  tomography with structural-prior-based diffuse optical tomography:combating optical  background uncertainty,"APPLIED OPTICS 53(30),6970-6982(2014)。

[3]M.Dhenain,S.W.Ruffins,and R.E.Jacobs,"Three-Dimensional Digital Mouse Atlas  Using High-Resolution MRI,"Developmental Biology 232,458-470(2001)。

[4]W.P.Segars,B.M.W.Tsui,and E.C.Frey et al,"Development of a 4-D Digital Mouse  Phantom for Molecular Imaging Research,"Molecular Imaging and Biology 6(3), 149–159(2004)。

[5]B.Dogdas,D.Stout and A.F.Chatziioannou et al,"Digimouse:a 3D whole body mouse  atlas from CT and cryosection data,"PHYSICS IN MEDICINE AND BIOLOGY 52(3), 577-587(2007)。

[6]H.Chui and A.Rangarajan,"A new point matching algorithm for non-rigid registration," Computer Vision and Image Understanding 89,114-141(2003)。

[7]S.Lee,G.Wolberg,and S.Y.Shin,"Scattered Data Interpolation with Multilevel  B-Splines,"IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS 3(3),228-244(1997)。

发明内容

针对现有技术中存在的问题，本发明提出一种基于XCT图像非刚度配准的目标鼠解剖 结构图谱获取方法，本发明方法将Digimouse模型^[4]选定为小鼠的标准解剖结构，并通过非 刚度图像配准算法，将Digimouse模型对应的XCT图像配准至目标鼠XCT图像，从而实现对 目标小鼠各个组织器官的标定。

为了解决上述技术问题，本发明提出的一种基于XCT图像非刚度配准的目标鼠解剖结 构图谱获取方法是：首先，将Digimouse模型设定为参考鼠解剖结构图谱，同时将Digimouse 对应的XCT图像设定为参考图像；其次，对目标鼠进行XCT成像获取目标图像并进行预 处理；然后，利用非刚度图像配准技术构建参考图像至目标图像的配准映射矩阵；最终， 将配准映射矩阵作用于参考鼠解剖结构图谱，构建出目标鼠解剖结构图谱。

进一步讲，本发明基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法的具体步 骤如下：

步骤一、设定参考鼠解剖结构图谱和参考图像：

将Digimouse结构设定为参考鼠解剖结构图谱A_r，将所述Digimouse的XCT图像设定 为参考图像I_r；

步骤二、目标鼠XCT图像的获取与预处理：

利用XCT设备对目标鼠进行全身成像，获取目标鼠的XCT图像；对目标鼠XCT图像 进行仿射变换，将目标鼠XCT图像中小鼠的头部与背部的朝向变换至与参考图像I_r相同； 将变换后的XCT图像作为目标图像I_t；

步骤三、构建初步配准映射矩阵M_c和初步配准图像I_c：

利用图像分割算法分别分割出参考图像I_r中的小鼠骨骼区域与表皮区域，及目标图像I_t中的小鼠骨骼区域与表皮区域；

利用边缘检测算法分别提取上述参考图像I_r中小鼠骨骼区域和表皮区域及目标图像I_t中的小鼠骨骼区域和表皮区域共计四个区域的边界轮廓，对所述四个区域的边界轮廓进行 等概率采样，从而计算出参考鼠骨骼特征点集L_rb、参考鼠表皮特征点集L_rs、目标鼠骨骼 特征点集L_tb以及目标鼠表皮特征点集L_ts；

利用TPS-RPM(Thin-plate Spline Robust Point Matching，薄板样条鲁棒点配准)算法^[6]的开源程序，计算参考鼠骨骼特征点集L_rb与目标鼠骨骼特征点集L_tb之间的骨骼特征点对 应矩阵C_b；在计算时，将目标鼠骨骼特征点集L_tb设定为目标点集，将参考鼠骨骼特征点集 L_rb设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止 条件，通过迭代计算出骨骼特征点对应矩阵C_b；：

利用TPS-RPM算法的开源程序，计算参考鼠表皮特征点集L_rs与目标鼠表皮特征点集 L_ts之间的表皮特征点对应矩阵C_s；在计算时，将目标鼠表皮特征点集L_ts设定为目标点集， 将参考鼠表皮特征点集L_rs设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初 始温度系数及迭代终止条件，通过迭代计算出表皮特征点对应矩阵C_s；

分别将求得的骨骼特征点对应矩阵C_b及表皮特征点对应矩阵C_s作用于参考鼠骨骼特 征点集L_rb及参考鼠表皮特征点集L_rs，得到初步配准后的骨骼特征点集L_cb及初步配准后的 表皮特征点集L_cs；

利用上述求得的初步特征点配准的结果，构建初步配准局部位移矩阵P：

$P=\left(\begin{array}{cc}{L}_{\mathrm{rb}}& L_{\mathrm{rb}}-L_{\mathrm{cb}}\\ L_{\mathrm{rs}}& L_{\mathrm{rs}}-L_{\mathrm{cs}}\end{array}\right)=[x,y,z,\mathrm{\Delta x},\mathrm{\Delta y},\mathrm{\Delta z}]---\left(3\right)$

将上述矩阵P转化为三组四维数据点集P_x＝{(x,y,z,△x)}，P_y＝{(x,y,z,△y)}， P_z＝{(x,y,z,△z)}，将Δx、Δy及Δz分别视为点(x,y,z)的函数值，即△x＝G₁(x,y,z)， △y＝G₂(x,y,z)，△z＝G₃(x,y,z)；

采用多层次B样条拟合算法分别对所述三组四维数据点集P_x＝{(x,y,z,△x)}， P_y＝{(x,y,z,△y)}，P_z＝{(x,y,z,△z)}进行拟合，用于计算所述参考图像I_r中各个像素沿着x、 y、z三个方向上的位移量；

利用所述多层次B样条拟合方法拟合数据点集P_x＝{(x,y,z,△x)}的过程详细描述如下：

所述多层次是指，利用覆于参考图像I_r上的一组逐渐加密的立方控制网格 Φ₀,Φ₁,...,Φ_k,...,Φ_h依次对迭代更新的四维数据点集进行B样条拟合，并将 所求h层拟合函数之和作为最终多层次B样条拟合函数；其中， △⁰ξ＝△x，△^k+1ξ＝△^kξ-g_k(x,y,z)，g_k(x,y,z)为第k层B样条拟合结果；

所述第k层B样条拟合过程描述如下：

假设第k层控制网格Φ_k尺寸为K_x×K_y×K_z，则第k层B样条拟合函数如下式所示：

$g_k(x,y,z)=\underset{l=0}{\overset{3}{\Sigma }}\underset{m=0}{\overset{3}{\Sigma }}\underset{n=0}{\overset{3}{\Sigma }}{B}_l\left(d_x\right)B_m\left(d_y\right)B_n\left(d_z\right){\phi }_{k,(l+i,m+j,n+k)}---\left(4\right)$

式(4)中，φ_{k,(l+i,m+j,n+k)}为位于控制网格Φ_k中坐标为 (l+i,m+j,n+k)的控制节点值；l,m,n∈{0,1,2,3}；B_l、B_m及B_n分别为第l,m,n阶B样条基 函数，其中0至3阶B样条函数的表达式描述如下:

$\left(\begin{array}{c}{B}_0\left(\delta \right)={(1-\delta )}^3/6\\ B_1\left(\delta \right)=(3{\delta }^3-6{\delta }^2+4)/6\\ B_2\left(\delta \right)=(-3{\delta }^3+3{\delta }^2+3\delta +1)/6\\ B_3\left(\delta \right)={\delta }^3/6\end{array}\right)---\left(5\right)$

所述式(4)中，控制网格Φ_k中各个控制节点的取值通过以下两步计算：

(a)计算中每个数据点对控制网格Φ_k中各个控制节点取值的影响量：

以中的一个数据点p＝(x_p,y_p,z_p,△^kξ_p)为例进行说明如下：

数据点p＝(x_p,y_p,z_p,△^kξ_p)对控制网格Φ_k中各个控制节点的影响矩阵表现为一个尺寸 为K_x×K_y×K_z的矩阵Ψ_p；为计算简便，定义与矩阵Ψ_p尺寸相同的两个矩阵Γ_p与Ω_p；所 述矩阵Ψ_p、Γ_p与Ω_p中坐标为(l+i,m+j,n+k)的元素分别通过式(6)计算：

$\left(\begin{array}{c}{\psi }_{p,(l+i,m+j,n+k)}=\frac{{\gamma }_{p,(l+i,m+j,n+k)}{\Delta }^k{\xi }_p}{{\omega }_{p,(l+i,m+j,n+k)}}\\ {\gamma }_{p,(l+i,m+j,n+k)}=B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right)\\ {\omega }_{p,(l+i,m+j,n+k)}={\Sigma }_{l=0}^3{\Sigma }_{m=0}^3{\Sigma }_{n=0}^3[B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right){]}^2\end{array}\right)---\left(6\right)$

式(6)中，l,m,n∈{0,1,2,3}；

在矩阵Ψ_p、Γ_p与Ω_p中，除所述坐标为(l+i,m+j,n+k)共计64个元素以外其余位置， ψ_p、γ_p与ω_p均为0；

(b)求取控制网格Φ_k中每个控制节点的值

综合中每个数据点对控制网格Φ_k中各个控制节点取值的影响，求取 制网格Φ_k中每个控制节点的值；控制网格Φ_k中坐标为(a,b,c)的控制节点φ_k,(a,b,c)取值为：

式(7)中，γ_p,(a,b,c)、ω_p,(a,b,c)、ψ_p,(a,b,c)分别为Ψ_p、Γ_p与Ω_p中坐标为(a,b,c)的元素值；

至此，第k层B样条拟合函数g_k(x,y,z)确立；综合各层次拟合函数，计算出多层次B 样条拟合函数利用所述求取的多层次B样条拟合函数g(x,y,z)， 计算出参考图像I_r中各个像素沿x轴的位移；

同理，利用多层次B样条拟合四维数据点集P_y＝{(x,y,z,△y)}和P_z＝{(x,y,z,△z)}，从而计算 出参考图像I_r中各个像素沿y、z轴向中的位移量，由此构建参考图像I_r的初步配准映射矩 阵M_c；

利用所述初步配准映射矩阵M_c，反向求解出初步配准图像I_c；在构建初步配准图像时， 灰度的赋值采用三线性插值方法；

步骤四、构建精细配准映射矩阵M_f：

利用图像分割算法，提取出初步配准图像I_c中的小鼠表皮区域及骨骼区域，并利用边 缘检测算法分别提取出所述小鼠表皮区域及骨骼区域的轮廓；将小鼠表皮及骨骼的轮廓相 互叠加，并使用矩形网格进行采样，由此得到一组初步配准图像特征点集L'_c；

利用块匹配方法求取初步配准图像特征点集L'_c中每个特征点在目标图像上的对应位 置，以L'_c中的任意一点p_l＝(x_p,y_p,z_p)为例，将所述过程描述如下：

在初步配准图像I_c中以坐标(x_p,y_p,z_p)为中心选取尺寸为N₁×N₁×N₁的立方邻域T，在 目标图像I_t中以坐标(x_p,y_p,z_p)为中心选取N₂×N₂×N₂的立方邻域S，其中N₂>N₁；将T作 为模板，S作为搜索域，在S区域中搜索与T具有最大相似度的子区域s₁，并将s₁的中心点 p_l′作为点p_l的对应位置；

以此类推，依次寻找出初步配准图像特征点集L'_c中各点在目标图像I_t上的对应位置， 由此构建出精细配准特征点集L_f；

利用初步配准图像特征点集L'_c及精细配准特征点集L_f构建精细配准局部位移矩阵Q：

Q＝[L'_c,L'_c-L_f]＝[x,y,z,Δx,Δy,Δz]      (8)

将上述矩阵Q转化为三组四维数据点集Q_x＝{(x,y,z,△x)}，Q_y＝{(x,y,z,△y)}， Q_z＝{(x,y,z,△z)}；与步骤三中多层次B样条拟合过程一致，分别对所述三组四维数据点集 Q_x＝{(x,y,z,△x)}，Q_y＝{(x,y,z,△y)}进行多层次B样条拟合，从而分别计算出初步配准图像 I_c中各个像素沿x、y、z三个轴向的位移量，由此构建初步配准图像I_c的精细配准映射矩阵 M_f；

步骤五、构建目标鼠解剖结构图谱：

将所述参考鼠解剖结构图谱A_r投影至像素坐标系，得到像素坐标系下的参考鼠解剖结 构图谱按顺序依次将初步配准映射矩阵M_c及精细配准映射矩阵M_f作用于所述像素坐 标系下的参考鼠解剖结构图谱使所述像素坐标系下的参考鼠解剖结构图谱产生与步 骤三初步配准及步骤四精细配准过程相同的变形，获得像素坐标系下的配准鼠解剖结构图 谱将所述像素坐标系下的配准鼠解剖结构图谱投影至物理坐标系下，得到物理坐 标系下的配准鼠解剖结构图谱A_f，所述物理坐标系下的配准鼠解剖结构图谱A_f即为目标 鼠解剖结构图谱。

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

1.本发明仅使用XCT单模态成像方式对小鼠成像，实验简单，成本较低；

2.本发明所提出的方法能够有效避免XCT图像单模态成像方法导致的软组织辨识难度， 简单易行的实现小鼠组织器官在XCT图像上的标识问题；

3.本发明使用的两步配准能够实现较好的配准精度，从而为小鼠解剖结构的正确标识提 供条件；

4.本发明在具体实施过程中亦可使用其他小鼠解剖结构图谱作为参考，以期实现处于不 同体态或发育时期下目标鼠的解剖结构图谱获取；

5.本发明的主要思想同样适用于其他医学应用领域如人脑结构学研究，即利用标准人脑 解剖学图像及本发明所述方法，实现目标人脑的解剖结构的获取。

附图说明

图1是基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法框图；

图2是本发明提出的基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法的具体 实施流程图；

图3是本发明在具体实施过程中的数据流图。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施 例仅仅对本发明进行解释说明，并不用以限制本发明。

图1示出了本发明基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法的基本步 骤是：

首先，将Digimouse结构设定为参考鼠解剖结构图谱，同时将Digimouse对应的XCT 图像设定为参考图像；

其次，对目标鼠进行XCT成像获取目标图像并进行预处理；

然后，利用非刚度图像配准技术构建参考图像至目标图像的配准映射矩阵；

最终，将配准映射矩阵作用于参考鼠解剖结构图谱，构建出目标鼠解剖结构图谱；可 选择将配准映射矩阵作用于参考图像，构建出配准图像，用以评定图像配准精度。

由于不同小鼠在成像时存在较大差异，所述配准方法首先进行初步配准调整小鼠大致 结构，其次使用精细配准调整图像细节部分。因此，本发明基于XCT图像非刚度配准的目 标鼠解剖结构图谱获取方法在实施过程中可细化为5个步骤。所述实施过程的流程图如图2 所示，实施过程中的数据流如图3所示。对本发明基于XCT图像非刚度配准的目标鼠解剖 结构图谱获取方法具体实施过程中的5个步骤详细说明如下：

步骤一、设定参考鼠解剖结构图谱和参考图像：

将Dogdas等人文献中所提出的Digimouse模型^[5]设定为参考鼠解剖结构图谱A_r，将所 述Digimouse的XCT图像设定为参考图像I_r；

至此，数字化的小鼠参考模型设定完毕；

步骤二、目标鼠XCT图像的获取与预处理：

利用XCT设备对目标鼠进行全身成像，获取目标鼠的XCT图像；对目标鼠XCT图像 进行仿射变换，将目标鼠XCT图像中小鼠的头部与背部的朝向变换至与参考图像I_r相同； 将变换后的XCT图像作为目标图像I_t；

至此，源于目标鼠的目标图像I_t获取完毕；

步骤三、构建初步配准映射矩阵M_c和初步配准图像I_c：

由于目标鼠与参考鼠的姿势、尺寸及体内空腔位置存在着较大的差别，因此本发明首 先进行初步图像配准，以调整两幅图像中小鼠尺寸与体态的差异。所述步骤三的主要过程 在于对参考图像I_r及目标图像I_t进行初步配准，构建出控制参考图像I_r产生非刚度形变的初 步配准映射矩阵M_c，将M_c作用于参考图像I_r，构建出初步配准图像I_c。具体实施过程详 细描述如下：

(3-1)分别提取参考图像I_r及目标图像I_t的特征点集

利用图像分割算法分别分割出参考图像I_r中的小鼠骨骼区域与表皮区域，及目标图像I_t中的小鼠骨骼区域与表皮区域；

利用边缘检测算法分别提取上述参考图像I_r中小鼠骨骼区域和表皮区域及目标图像I_t中的小鼠骨骼区域和表皮区域共计四个区域的边界轮廓，对所述四个区域的边界轮廓进行 等概率采样，从而计算出参考鼠骨骼特征点集L_rb、参考鼠表皮特征点集L_rs、目标鼠骨骼 特征点集L_tb以及目标鼠表皮特征点集L_ts；

(3-2)基于TPS-RPM算法的特征点配准

利用Haili Chui提出的TPS-RPM算法^[6]的开源程序，计算参考鼠骨骼特征点集L_rb与目 标鼠骨骼特征点集L_tb之间的骨骼特征点对应矩阵C_b、参考鼠表皮特征点集L_rs与目标鼠表 皮特征点集L_ts之间的表皮特征点对应矩阵C_s；

所述参考鼠骨骼特征点集L_rb与目标鼠骨骼特征点集L_tb之间的骨骼特征点对应矩阵C_b及参考鼠表皮特征点集L_rs与目标鼠表皮特征点集L_ts之间的表皮特征点对应矩阵C_s均为模 糊对应矩阵，在模糊对应矩阵中各个元素均为[0,1]之间的浮点数，用以描述两点之间对应程 度的强弱；

其中，在计算参考鼠骨骼特征点集L_rb与目标鼠骨骼特征点集L_tb之间的骨骼特征点对 应矩阵C_b时，将目标鼠骨骼特征点集L_tb设定为目标点集，将参考鼠骨骼特征点集L_rb设定 为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通 过迭代计算下列能量方程的极值点：

$E(C_b,f)=\underset{i=1}{\overset{N_{\mathrm{rb}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\Sigma }}{c}_{b,\mathrm{ij}}{\left|\right|l_{\mathrm{rb},i}-f\left(l_{\mathrm{tb},j}\right)\left|\right|}^2+\lambda {\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N_{\mathrm{rb}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\Sigma }}{c}_{b,\mathrm{ij}}\log c_{b,\mathrm{ij}}-\gamma \underset{i=1}{\overset{N_{\mathrm{rb}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\Sigma }}{c}_{b,\mathrm{ij}}---\left(1\right)$

式(1)中，c_b,ij为对应矩阵C_b中的元素；N_rb为参考鼠骨骼特征点集L_rb所包含特征点 的个数，N_tb为目标鼠骨骼特征点集L_tb所包含点特征的个数；l_rb,i为参考鼠骨骼特征点集L_rb中的第i个特征点的坐标，l_tb,j为目标鼠骨骼特征点集L_tb中的第j个特征点的坐标；f为薄 板样条函数；||Lf||²为对f的平滑性约束，其中L为微分算子；T为模拟退火过程中用于控 制模糊对应程度的温度系数；λ为先验平滑度权值；γ为鲁棒性控制权值；随着温度系数T 逐渐减小，对应矩阵C_b与薄板样条函数f交替迭代；对应矩阵C_b的最终迭代结果即为所求 骨骼特征点对应矩阵C_b；

在计算参考鼠表皮特征点集L_rs与目标鼠表皮特征点集L_ts之间的表皮特征点对应矩阵 C_s时，将目标鼠表皮特征点集L_ts设定为目标点集，将参考鼠表皮特征点集L_rs设定为待配 准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代 计算下列能量方程的极值点：

$E(C_s,f)=\underset{i=1}{\overset{N_{\mathrm{rs}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\Sigma }}{c}_{s,\mathrm{ij}}{\left|\right|l_{\mathrm{rs},i}-f\left(l_{\mathrm{ts},j}\right)\left|\right|}^2+\lambda {\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N_{\mathrm{rs}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\Sigma }}{c}_{s,\mathrm{ij}}\log c_{s,\mathrm{ij}}-\gamma \underset{i=1}{\overset{N_{\mathrm{rs}}}{\Sigma }}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\Sigma }}{c}_{s,\mathrm{ij}}---\left(2\right)$

式(2)中，c_s,ij为对应矩阵C_s中的元素；N_rs为参考鼠表皮特征点集L_rs所包含特征点 的个数，N_ts为目标鼠表皮特征点集L_ts所包含点特征的个数；l_rs,i为参考鼠表皮特征点集L_rs中的第i个点的坐标，l_ts,j为目标鼠表皮特征点集L_ts中的第j个点的坐标；其余变量含义与 式(1)相同；随着温度系数T逐渐减小，对应矩阵C_s与薄板样条函数f交替迭代；对应矩 阵C_s的最终迭代结果即为所求表皮特征点对应矩阵C_s；

分别将求得的骨骼特征点对应矩阵C_b及表皮特征点对应矩阵C_s作用于参考鼠骨骼特 征点集L_rb及参考鼠表皮特征点集L_rs，得到初步配准后的骨骼特征点集L_cb及初步配准后的 表皮特征点集L_cs；

(3-3)初步配准局部位移矩阵P的构建

利用上述求得的初步特征点配准的结果，构建初步配准局部位移矩阵P：

$P=\left(\begin{array}{cc}{L}_{\mathrm{rb}}& L_{\mathrm{rb}}-L_{\mathrm{cb}}\\ L_{\mathrm{rs}}& L_{\mathrm{rs}}-L_{\mathrm{cs}}\end{array}\right)=[x,y,z,\mathrm{\Delta x},\mathrm{\Delta y},\mathrm{\Delta z}]---\left(3\right)$

(3-4)初步配准映射矩阵M_c的构建

将上述矩阵P转化为三组四维数据点集P_x＝{(x,y,z,△x)}，P_y＝{(x,y,z,△y)}， P_z＝{(x,y,z,△z)}，将Δx、Δy及Δz分别视为点(x,y,z)的函数值，即△x＝G₁(x,y,z)， △y＝G₂(x,y,z)，△z＝G₃(x,y,z)；

采用多层次B样条拟合算法分别对所述三组四维数据点集P_x＝{(x,y,z,△x)}， P_y＝{(x,y,z,△y)}，P_z＝{(x,y,z,△z)}进行拟合，用于计算所述参考图像I_r中各个像素沿着x、 y、z三个方向上的位移量；

本发明中多层次B样条拟合算法是Seungyong Lee提出的三维数据拟合算法^[7]在四维空 间上的扩展。利用所述多层次B样条拟合方法拟合数据点集P_x＝{(x,y,z,△x)}的过程详细描 述如下：

所述多层次是指，利用覆于参考图像I_r上的一组逐渐加密的立方控制网格 Φ₀,Φ₁,...,Φ_k,...,Φ_h依次对迭代更新的四维数据点集进行B样条拟合，并将 所求h层拟合函数之和作为最终多层次B样条拟合函数；其中， △⁰ξ＝△x，△^k+1ξ＝△^kξ-g_k(x,y,z)，g_k(x,y,z)为第k层B样条拟合结果；

所述第k层B样条拟合过程描述如下：

假设第k层控制网格Φ_k尺寸为K_x×K_y×K_z，则第k层B样条拟合函数如下式所示：

$g_k(x,y,z)=\underset{l=0}{\overset{3}{\Sigma }}\underset{m=0}{\overset{3}{\Sigma }}\underset{n=0}{\overset{3}{\Sigma }}{B}_l\left(d_x\right)B_m\left(d_y\right)B_n\left(d_z\right){\phi }_{k,(l+i,m+j,n+k)}---\left(4\right)$

式(4)中，φ_{k,(l+i,m+j,n+k)}为位于控制网格Φ_k中坐标为 (l+i,m+j,n+k)的控制节点值；l,m,n∈{0,1,2,3}；B_l、B_m及B_n分别为第l,m,n阶B样条基 函数，其中0至3阶B样条函数的表达式描述如下:

$\left(\begin{array}{c}{B}_0\left(\delta \right)={(1-\delta )}^3/6\\ B_1\left(\delta \right)=(3{\delta }^3-6{\delta }^2+4)/6\\ B_2\left(\delta \right)=(-3{\delta }^3+3{\delta }^2+3\delta +1)/6\\ B_3\left(\delta \right)={\delta }^3/6\end{array}\right)---\left(5\right)$

所述式(4)中，控制网格Φ_k中各个控制节点的取值通过以下两步计算：

(a)计算中每个数据点对控制网格Φ_k中各个控制节点取值的影响量：

以中的一个数据点p＝(x_p,y_p,z_p,△^kξ_p)为例进行说明如下：

数据点p＝(x_p,y_p,z_p,△^kξ_p)对控制网格Φ_k中各个控制节点的影响矩阵表现为一个尺寸 为K_x×K_y×K_z的矩阵Ψ_p；为计算简便，定义与矩阵Ψ_p尺寸相同的两个矩阵Γ_p与Ω_p；所 述矩阵Ψ_p、Γ_p与Ω_p中坐标为(l+i,m+j,n+k)的元素分别通过式(6)计算：

$\left(\begin{array}{c}{\psi }_{p,(l+i,m+j,n+k)}=\frac{{\gamma }_{p,(l+i,m+j,n+k)}{\Delta }^k{\xi }_p}{{\omega }_{p,(l+i,m+j,n+k)}}\\ {\gamma }_{p,(l+i,m+j,n+k)}=B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right)\\ {\omega }_{p,(l+i,m+j,n+k)}={\Sigma }_{l=0}^3{\Sigma }_{m=0}^3{\Sigma }_{n=0}^3[B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right){]}^2\end{array}\right)---\left(6\right)$

式(6)中，l,m,n∈{0,1,2,3}；

在矩阵Ψ_p、Γ_p与Ω_p中，除所述坐标为(l+i,m+j,n+k)共计64个元素以外其余位置， ψ_p、γ_p与ω_p均为0；

(b)求取控制网格Φ_k中每个控制节点的值

综合中每个数据点对控制网格Φ_k中各个控制节点取值的影响，求取 制网格Φ_k中每个控制节点的值；控制网格Φ_k中坐标为(a,b,c)的控制节点φ_k,(a,b,c)取值为：

式(7)中，γ_p,(a,b,c)、ω_p,(a,b,c)、ψ_p,(a,b,c)分别为Ψ_p、Γ_p与Ω_p中坐标为(a,b,c)的元素值；

至此，第k层B样条拟合函数g_k(x,y,z)确立；综合各层次拟合函数，计算出多层次B 样条拟合函数利用所述求取的多层次B样条拟合函数g(x,y,z)， 计算出参考图像I_r中各个像素沿x轴的位移；

同理，利用多层次B样条拟合四维数据点集P_y＝{(x,y,z,△y)}和P_z＝{(x,y,z,△z)}，从而计算 出参考图像I_r中各个像素沿y、z轴向中的位移量，由此构建参考图像I_r的初步配准映射矩 阵M_c；

(3-5)初步配准图像I_c的构建

利用所述初步配准映射矩阵M_c，反向求解出初步配准图像I_c；在构建初步配准图像时， 灰度的赋值采用三线性插值方法；

至此，初步配准映射矩阵M_c，与初步配准图像I_c获取完毕。

步骤四、构建精细配准映射矩阵M_f及精细配准图像I_f：

经过初步配准之后，初步配准图像I_c与目标图像I_t仍具有较大差别，为了提高配准精度， 需要进行更为精细的配准。所述步骤四的主要过程在于对初步配准图像I_c及目标图像I_t进行 精细配准，构建出控制初步配准图像I_c产生非刚度形变的精细配准映射矩阵M_f，将M_f作 用于初步配准图像I_c，构建出精细配准图像I_cf。具体实施过程详细描述如下：

(4-1)重新提取初步配准图像特征点集

利用图像分割算法，提取出初步配准图像I_c中的小鼠表皮区域及骨骼区域，并利用边 缘检测算法分别提取出所述小鼠表皮区域及骨骼区域的轮廓；将小鼠表皮及骨骼的轮廓相 互叠加，并使用矩形网格进行采样，由此得到一组初步配准图像特征点集L'_c；

(4-2)基于块匹配方法的特征点配准

利用块匹配方法求取初步配准图像特征点集L'_c中每个特征点在目标图像上的对应位 置，以L'_c中的任意一点p_l＝(x_p,y_p,z_p)为例，将所述过程描述如下：

在初步配准图像I_c中以坐标(x_p,y_p,z_p)为中心选取尺寸为N₁×N₁×N₁的立方邻域T，在 目标图像I_t中以坐标(x_p,y_p,z_p)为中心选取N₂×N₂×N₂的立方邻域S，其中N₂>N₁；将T作 为模板，S作为搜索域，在S区域中搜索与T具有最大相似度的子区域s₁，并将s₁的中心点 p_l′作为点p_l的对应位置；

以此类推，依次寻找出初步配准图像特征点集L'_c中各点在目标图像I_t上的对应位置， 由此构建出精细配准特征点集L_f；

(4-3)精细配准局部位移矩阵Q的构建

利用初步配准图像特征点集L'_c及精细配准特征点集L_f构建精细配准局部位移矩阵Q：

Q＝[L'_c,L'_c-L_f]＝[x,y,z,Δx,Δy,Δz]      (8)

(4-4)精细配准映射矩阵M_f的构建

将上述矩阵Q转化为三组四维数据点集Q_x＝{(x,y,z,△x)}，Q_y＝{(x,y,z,△y)}， Q_z＝{(x,y,z,△z)}；与步骤三中多层次B样条拟合过程一致，分别对所述三组四维数据点集 Q_x＝{(x,y,z,△x)}，Q_y＝{(x,y,z,△y)}进行多层次B样条拟合，从而分别计算出初步配准图像 I_c中各个像素沿x、y、z三个轴向的位移量，由此构建初步配准图像I_c的精细配准映射矩阵 M_f；

至此，精细配准映射矩阵M_f获取完毕。

(4-5)精细配准图像I_f的构建(可选)

此过程之后，可以选择将精细配准映射矩阵M_f作用于初步配准图像I_c，构建出精细配 准图像I_f；在构建初步配准图像时，灰度的赋值同样采用三线性插值方法；通过判断最终 配准图像(即精细配准图像I_f)与目标图像I_t的相似度以评定配准的有效程度；

步骤五、构建目标鼠解剖结构图谱：

在非刚度图像配准过程中，通过初步配准映射矩阵M_c及精细配准映射矩阵M_f的控制 作用，参考图像I_r变形至精细配准图像I_f，精细配准图像I_f与目标图像I_t具有较高相似度。 因此，将初步配准映射矩阵M_c及精细配准映射矩阵M_f依次作用于参考鼠的解剖结构图谱 A_r，则可以近似构建目标鼠的解剖结构图谱。然而，参考解剖结构图谱建立于物理坐标系 下(以mm为单位)，而上述两个配准映射矩阵建立于像素坐标系下(以像素为单位)，因 此在此过程中需要对A_r进行坐标变换。具体实施过程详细描述如下：

(5-1)将所述参考鼠解剖结构图谱A_r投影至像素坐标系，得到像素坐标系下的参考鼠 解剖结构图谱

(5-2)按顺序依次将初步配准映射矩阵M_c及精细配准映射矩阵M_f作用于所述像素坐 标系下的参考鼠解剖结构图谱使所述像素坐标系下的参考鼠解剖结构图谱产生与步 骤三初步配准及步骤四精细配准过程相同的变形，获得像素坐标系下的配准鼠解剖结构图 谱

(5-3)将所述像素坐标系下的配准鼠解剖结构图谱投影至物理坐标系下，得到物理 坐标系下的配准鼠解剖结构图谱A_f，所述物理坐标系下的配准鼠解剖结构图谱A_f即为目 标鼠解剖结构图谱。

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式， 上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明 的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保 护之内。