[发明专利]基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法

权利要求书

1.一种基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法，其特征在于，其基本步骤是：首先，将Digimouse模型设定为参考鼠解剖结构图谱，同时将Digimouse对应的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算法的开源程序，计算参考鼠骨骼特征点集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；

利用TPS-RPM算法的开源程序求取骨骼特征点对应矩阵C_b与表皮特征点对应矩阵C_s的具体过程如下：

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

所述TPS-RPM算法在计算参考鼠骨骼特征点集L_rb与目标鼠骨骼特征点集L_tb之间的骨骼特征点对应矩阵C_b时，利用模拟退火算法计算如下能量方程的极小值点：

E(Cb,f)=Σi=1NrbΣj=1Ntbcb,ij||lrb,i-f(ltb,j)||2+λ||Lf||2+TΣi=1NrbΣj=1Ntbcb,ijlog cb,ij-γΣi=1NrbΣj=1Ntbcb,ij---(1)]]>

式(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；

所述TPS-RPM算法在计算参考鼠表皮特征点集L_rs与目标鼠表皮特征点集L_ts之间的表皮特征点对应矩阵C_s时，利用模拟退火算法计算如下能量方程的极小值点：

E(Cs,f)=Σi=1NrsΣj=1Ntscs,ij||lrs,i-f(lts,j)||2+λ||Lf||2+TΣi=1NrsΣj=1Ntscs,ijlog cs,ij-γΣi=1NrsΣj=1Ntscs,ij---(2)]]>

式(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；

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

P=LrbLrb-LcbLrsLrs-Lcs=[x,y,z,Δx,Δy,Δz]---(3)]]>

将上述矩阵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样条拟合函数如下式所示：

gk(x,y,z)=Σl=03Σm=03Σn=03Bl(dx)Bm(dy)Bn(dz)φk,(l+i,m+j,n+k)---(4)]]>

式(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样条函数的表达式描述如下:

B0(δ)=(1-δ)3/6B1(δ)=(3δ3-6δ2+4)/6B2(δ)=(-3δ3+3δ2+3δ+1)/6B3(δ)=δ3/6---(5)]]>

所述式(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)计算：

ψp,(l+i,m+j,n+k)=γp,(1+i,m+j,n+k)Δkξpωp,(l+i,m+j,n+k)γp,(l+i,m+j,n+k)=Bl(dxp)Bm(dyp)Bn(dzp)ωp,(l+i,m+j,n+k)=Σl=03Σm=03Σn=03[Bl(dxp)Bm(dyp)Bn(dzp)]2---(6)]]>

式(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即为目标鼠解剖结构图谱。