CN104867104B - 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法 - Google Patents

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

Info

Publication number
CN104867104B
CN104867104B CN201510259870.5A CN201510259870A CN104867104B CN 104867104 B CN104867104 B CN 104867104B CN 201510259870 A CN201510259870 A CN 201510259870A CN 104867104 B CN104867104 B CN 104867104B
Authority
CN
China
Prior art keywords
mrow
msub
mouse
point set
target
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201510259870.5A
Other languages
English (en)
Other versions
CN104867104A (zh
Inventor
高峰
万文博
赵会娟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN201510259870.5A priority Critical patent/CN104867104B/zh
Publication of CN104867104A publication Critical patent/CN104867104A/zh
Application granted granted Critical
Publication of CN104867104B publication Critical patent/CN104867104B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G06T3/147—Transformations for image registration, e.g. adjusting or mapping for alignment of images using affine transformations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/12—Edge-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Image Processing (AREA)

Abstract

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

Description

基于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 diffuseoptical tomography holds promise for sensitivity enhancement of fluorescencemolecular tomography,"Biomedical Optics Express 10,3640-3659(2014)。
[2]L.-H.Wu,H.-J.Zhao and X.Wang et al,"Enhancement of fluorescencemolecular 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 DigitalMouse 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-DDigital Mouse Phantom for Molecular Imaging Research,"Molecular Imaging andBiology 6(3),149–159(2004)。
[5]B.Dogdas,D.Stout and A.F.Chatziioannou et al,"Digimouse:a 3D wholebody mouse atlas from CT and cryosection data,"PHYSICS IN MEDICINE ANDBIOLOGY 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 withMultilevel B-Splines,"IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTERGRAPHICS3(3),228-244(1997)。
发明内容
针对现有技术中存在的问题，本发明提出一种基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法，本发明方法将Digimouse模型[4]选定为小鼠的标准解剖结构，并通过非刚度图像配准算法，将Digimouse模型对应的XCT图像配准至目标鼠XCT图像，从而实现对目标小鼠各个组织器官的标定。
为了解决上述技术问题，本发明提出的一种基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法是：首先，将Digimouse模型设定为参考鼠解剖结构图谱，同时将Digimouse对应的XCT图像设定为参考图像；其次，对目标鼠进行XCT成像获取目标图像并进行预处理；然后，利用非刚度图像配准技术构建参考图像至目标图像的配准映射矩阵；最终，将配准映射矩阵作用于参考鼠解剖结构图谱，构建出目标鼠解剖结构图谱。
进一步讲，本发明基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法的具体步骤如下：
步骤一、设定参考鼠解剖结构图谱和参考图像：
将Digimouse结构设定为参考鼠解剖结构图谱Ar，将所述Digimouse的XCT图像设定为参考图像Ir；
步骤二、目标鼠XCT图像的获取与预处理：
利用XCT设备对目标鼠进行全身成像，获取目标鼠的XCT图像；对目标鼠XCT图像进行仿射变换，将目标鼠XCT图像中小鼠的头部与背部的朝向变换至与参考图像Ir相同；将变换后的XCT图像作为目标图像It；
步骤三、构建初步配准映射矩阵Mc和初步配准图像Ic：
利用图像分割算法分别分割出参考图像Ir中的小鼠骨骼区域与表皮区域，及目标图像It中的小鼠骨骼区域与表皮区域；
利用边缘检测算法分别提取上述参考图像Ir中小鼠骨骼区域和表皮区域及目标图像It中的小鼠骨骼区域和表皮区域共计四个区域的边界轮廓，对所述四个区域的边界轮廓进行等概率采样，从而计算出参考鼠骨骼特征点集Lrb、参考鼠表皮特征点集Lrs、目标鼠骨骼特征点集Ltb以及目标鼠表皮特征点集Lts；
利用TPS-RPM(Thin-plate Spline Robust Point Matching，薄板样条鲁棒点配准)算法[6]的开源程序，计算参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb；在计算时，将目标鼠骨骼特征点集Ltb设定为目标点集，将参考鼠骨骼特征点集Lrb设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代计算出骨骼特征点对应矩阵Cb；：
利用TPS-RPM算法的开源程序，计算参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs；在计算时，将目标鼠表皮特征点集Lts设定为目标点集，将参考鼠表皮特征点集Lrs设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代计算出表皮特征点对应矩阵Cs；
分别将求得的骨骼特征点对应矩阵Cb及表皮特征点对应矩阵Cs作用于参考鼠骨骼特征点集Lrb及参考鼠表皮特征点集Lrs，得到初步配准后的骨骼特征点集Lcb及初步配准后的表皮特征点集Lcs；
利用上述求得的初步特征点配准的结果，构建初步配准局部位移矩阵P：
将上述矩阵P转化为三组四维数据点集Px＝{(x,y,z,△x)}，Py＝{(x,y,z,△y)}，Pz＝{(x,y,z,△z)}，将Δx、Δy及Δz分别视为点(x,y,z)的函数值，即△x＝G1(x,y,z)，△y＝G2(x,y,z)，△z＝G3(x,y,z)；
采用多层次B样条拟合算法分别对所述三组四维数据点集Px＝{(x,y,z,△x)}，Py＝{(x,y,z,△y)}，Pz＝{(x,y,z,△z)}进行拟合，用于计算所述参考图像Ir中各个像素沿着x、y、z三个方向上的位移量；
利用所述多层次B样条拟合方法拟合数据点集Px＝{(x,y,z,△x)}的过程详细描述如下：
所述多层次是指，利用覆于参考图像Ir上的一组逐渐加密的立方控制网格Φ0,Φ1,...,Φk,...,Φh依次对迭代更新的四维数据点集进行B样条拟合，并将所求h层拟合函数之和作为最终多层次B样条拟合函数；其中，△0ξ＝△x，△k+1ξ＝△kξ-gk(x,y,z)，gk(x,y,z)为第k层B样条拟合结果；
所述第k层B样条拟合过程描述如下：
假设第k层控制网格Φk尺寸为Kx×Ky×Kz，则第k层B样条拟合函数如下式所示：
式(4)中， φk,(l+i,m+j,n+k)为位于控制网格Φk中坐标为(l+i,m+j,n+k)的控制节点值；l,m,n∈{0,1,2,3}；Bl、Bm及Bn分别为第l,m,n阶B样条基函数，其中0至3阶B样条函数的表达式描述如下:
所述式(4)中，控制网格Φk中各个控制节点的取值通过以下两步计算：
(a)计算中每个数据点对控制网格Φk中各个控制节点取值的影响量：
以中的一个数据点p＝(xp,yp,zp,△kξp)为例进行说明如下：
数据点p＝(xp,yp,zp,△kξp)对控制网格Φk中各个控制节点的影响矩阵表现为一个尺寸为Kx×Ky×Kz的矩阵Ψp；为计算简便，定义与矩阵Ψp尺寸相同的两个矩阵Γp与Ωp；所述矩阵Ψp、Γp与Ωp中坐标为(l+i,m+j,n+k)的元素分别通过式(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样条拟合函数gk(x,y,z)确立；综合各层次拟合函数，计算出多层次B样条拟合函数利用所述求取的多层次B样条拟合函数g(x,y,z)，计算出参考图像Ir中各个像素沿x轴的位移；
同理，利用多层次B样条拟合四维数据点集Py＝{(x,y,z,△y)}和Pz＝{(x,y,z,△z)}，从而计算出参考图像Ir中各个像素沿y、z轴向中的位移量，由此构建参考图像Ir的初步配准映射矩阵Mc；
利用所述初步配准映射矩阵Mc，反向求解出初步配准图像Ic；在构建初步配准图像时，灰度的赋值采用三线性插值方法；
步骤四、构建精细配准映射矩阵Mf：
利用图像分割算法，提取出初步配准图像Ic中的小鼠表皮区域及骨骼区域，并利用边缘检测算法分别提取出所述小鼠表皮区域及骨骼区域的轮廓；将小鼠表皮及骨骼的轮廓相互叠加，并使用矩形网格进行采样，由此得到一组初步配准图像特征点集L'c；
利用块匹配方法求取初步配准图像特征点集L'c中每个特征点在目标图像上的对应位置，以L'c中的任意一点pl＝(xp,yp,zp)为例，将所述过程描述如下：
在初步配准图像Ic中以坐标(xp,yp,zp)为中心选取尺寸为N1×N1×N1的立方邻域T，在目标图像It中以坐标(xp,yp,zp)为中心选取N2×N2×N2的立方邻域S，其中N2>N1；将T作为模板，S作为搜索域，在S区域中搜索与T具有最大相似度的子区域s1，并将s1的中心点pl′作为点pl的对应位置；
以此类推，依次寻找出初步配准图像特征点集L'c中各点在目标图像It上的对应位置，由此构建出精细配准特征点集Lf；
利用初步配准图像特征点集L'c及精细配准特征点集Lf构建精细配准局部位移矩阵Q：
Q＝[L'c,L'c-Lf]＝[x,y,z,Δx,Δy,Δz] (8)
将上述矩阵Q转化为三组四维数据点集Qx＝{(x,y,z,△x)}，Qy＝{(x,y,z,△y)}，Qz＝{(x,y,z,△z)}；与步骤三中多层次B样条拟合过程一致，分别对所述三组四维数据点集Qx＝{(x,y,z,△x)}，Qy＝{(x,y,z,△y)}进行多层次B样条拟合，从而分别计算出初步配准图像Ic中各个像素沿x、y、z三个轴向的位移量，由此构建初步配准图像Ic的精细配准映射矩阵Mf；
步骤五、构建目标鼠解剖结构图谱：
将所述参考鼠解剖结构图谱Ar投影至像素坐标系，得到像素坐标系下的参考鼠解剖结构图谱按顺序依次将初步配准映射矩阵Mc及精细配准映射矩阵Mf作用于所述像素坐标系下的参考鼠解剖结构图谱使所述像素坐标系下的参考鼠解剖结构图谱产生与步骤三初步配准及步骤四精细配准过程相同的变形，获得像素坐标系下的配准鼠解剖结构图谱将所述像素坐标系下的配准鼠解剖结构图谱投影至物理坐标系下，得到物理坐标系下的配准鼠解剖结构图谱Af，所述物理坐标系下的配准鼠解剖结构图谱Af即为目标鼠解剖结构图谱。
与现有技术相比，本发明的有益效果是：
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]设定为参考鼠解剖结构图谱Ar，将所述Digimouse的XCT图像设定为参考图像Ir；
至此，数字化的小鼠参考模型设定完毕；
步骤二、目标鼠XCT图像的获取与预处理：
利用XCT设备对目标鼠进行全身成像，获取目标鼠的XCT图像；对目标鼠XCT图像进行仿射变换，将目标鼠XCT图像中小鼠的头部与背部的朝向变换至与参考图像Ir相同；将变换后的XCT图像作为目标图像It；
至此，源于目标鼠的目标图像It获取完毕；
步骤三、构建初步配准映射矩阵Mc和初步配准图像Ic：
由于目标鼠与参考鼠的姿势、尺寸及体内空腔位置存在着较大的差别，因此本发明首先进行初步图像配准，以调整两幅图像中小鼠尺寸与体态的差异。所述步骤三的主要过程在于对参考图像Ir及目标图像It进行初步配准，构建出控制参考图像Ir产生非刚度形变的初步配准映射矩阵Mc，将Mc作用于参考图像Ir，构建出初步配准图像Ic。具体实施过程详细描述如下：
(3-1)分别提取参考图像Ir及目标图像It的特征点集
利用图像分割算法分别分割出参考图像Ir中的小鼠骨骼区域与表皮区域，及目标图像It中的小鼠骨骼区域与表皮区域；
利用边缘检测算法分别提取上述参考图像Ir中小鼠骨骼区域和表皮区域及目标图像It中的小鼠骨骼区域和表皮区域共计四个区域的边界轮廓，对所述四个区域的边界轮廓进行等概率采样，从而计算出参考鼠骨骼特征点集Lrb、参考鼠表皮特征点集Lrs、目标鼠骨骼特征点集Ltb以及目标鼠表皮特征点集Lts；
(3-2)基于TPS-RPM算法的特征点配准
利用Haili Chui提出的TPS-RPM算法[6]的开源程序，计算参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb、参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs；
所述参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb及参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs均为模糊对应矩阵，在模糊对应矩阵中各个元素均为[0,1]之间的浮点数，用以描述两点之间对应程度的强弱；
其中，在计算参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb时，将目标鼠骨骼特征点集Ltb设定为目标点集，将参考鼠骨骼特征点集Lrb设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代计算下列能量方程的极值点：
式(1)中，cb,ij为对应矩阵Cb中的元素；Nrb为参考鼠骨骼特征点集Lrb所包含特征点的个数，Ntb为目标鼠骨骼特征点集Ltb所包含点特征的个数；lrb,i为参考鼠骨骼特征点集Lrb中的第i个特征点的坐标，ltb,j为目标鼠骨骼特征点集Ltb中的第j个特征点的坐标；f为薄板样条函数；||Lf||2为对f的平滑性约束，其中L为微分算子；T为模拟退火过程中用于控制模糊对应程度的温度系数；λ为先验平滑度权值；γ为鲁棒性控制权值；随着温度系数T逐渐减小，对应矩阵Cb与薄板样条函数f交替迭代；对应矩阵Cb的最终迭代结果即为所求骨骼特征点对应矩阵Cb；
在计算参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs时，将目标鼠表皮特征点集Lts设定为目标点集，将参考鼠表皮特征点集Lrs设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代计算下列能量方程的极值点：
式(2)中，cs,ij为对应矩阵Cs中的元素；Nrs为参考鼠表皮特征点集Lrs所包含特征点的个数，Nts为目标鼠表皮特征点集Lts所包含点特征的个数；lrs,i为参考鼠表皮特征点集Lrs中的第i个点的坐标，lts,j为目标鼠表皮特征点集Lts中的第j个点的坐标；其余变量含义与式(1)相同；随着温度系数T逐渐减小，对应矩阵Cs与薄板样条函数f交替迭代；对应矩阵Cs的最终迭代结果即为所求表皮特征点对应矩阵Cs；
分别将求得的骨骼特征点对应矩阵Cb及表皮特征点对应矩阵Cs作用于参考鼠骨骼特征点集Lrb及参考鼠表皮特征点集Lrs，得到初步配准后的骨骼特征点集Lcb及初步配准后的表皮特征点集Lcs；
(3-3)初步配准局部位移矩阵P的构建
利用上述求得的初步特征点配准的结果，构建初步配准局部位移矩阵P：
(3-4)初步配准映射矩阵Mc的构建
将上述矩阵P转化为三组四维数据点集Px＝{(x,y,z,△x)}，Py＝{(x,y,z,△y)}，Pz＝{(x,y,z,△z)}，将Δx、Δy及Δz分别视为点(x,y,z)的函数值，即△x＝G1(x,y,z)，△y＝G2(x,y,z)，△z＝G3(x,y,z)；
采用多层次B样条拟合算法分别对所述三组四维数据点集Px＝{(x,y,z,△x)}，Py＝{(x,y,z,△y)}，Pz＝{(x,y,z,△z)}进行拟合，用于计算所述参考图像Ir中各个像素沿着x、y、z三个方向上的位移量；
本发明中多层次B样条拟合算法是Seungyong Lee提出的三维数据拟合算法[7]在四维空间上的扩展。利用所述多层次B样条拟合方法拟合数据点集Px＝{(x,y,z,△x)}的过程详细描述如下：
所述多层次是指，利用覆于参考图像Ir上的一组逐渐加密的立方控制网格Φ0,Φ1,...,Φk,...,Φh依次对迭代更新的四维数据点集进行B样条拟合，并将所求h层拟合函数之和作为最终多层次B样条拟合函数；其中，△0ξ＝△x，△k+1ξ＝△kξ-gk(x,y,z)，gk(x,y,z)为第k层B样条拟合结果；
所述第k层B样条拟合过程描述如下：
假设第k层控制网格Φk尺寸为Kx×Ky×Kz，则第k层B样条拟合函数如下式所示：
式(4)中， φk,(l+i,m+j,n+k)为位于控制网格Φk中坐标为(l+i,m+j,n+k)的控制节点值；l,m,n∈{0,1,2,3}；Bl、Bm及Bn分别为第l,m,n阶B样条基函数，其中0至3阶B样条函数的表达式描述如下:
所述式(4)中，控制网格Φk中各个控制节点的取值通过以下两步计算：
(a)计算中每个数据点对控制网格Φk中各个控制节点取值的影响量：
以中的一个数据点p＝(xp,yp,zp,△kξp)为例进行说明如下：
数据点p＝(xp,yp,zp,△kξp)对控制网格Φk中各个控制节点的影响矩阵表现为一个尺寸为Kx×Ky×Kz的矩阵Ψp；为计算简便，定义与矩阵Ψp尺寸相同的两个矩阵Γp与Ωp；所述矩阵Ψp、Γp与Ωp中坐标为(l+i,m+j,n+k)的元素分别通过式(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样条拟合函数gk(x,y,z)确立；综合各层次拟合函数，计算出多层次B样条拟合函数利用所述求取的多层次B样条拟合函数g(x,y,z)，计算出参考图像Ir中各个像素沿x轴的位移；
同理，利用多层次B样条拟合四维数据点集Py＝{(x,y,z,△y)}和Pz＝{(x,y,z,△z)}，从而计算出参考图像Ir中各个像素沿y、z轴向中的位移量，由此构建参考图像Ir的初步配准映射矩阵Mc；
(3-5)初步配准图像Ic的构建
利用所述初步配准映射矩阵Mc，反向求解出初步配准图像Ic；在构建初步配准图像时，灰度的赋值采用三线性插值方法；
至此，初步配准映射矩阵Mc，与初步配准图像Ic获取完毕。
步骤四、构建精细配准映射矩阵Mf及精细配准图像If：
经过初步配准之后，初步配准图像Ic与目标图像It仍具有较大差别，为了提高配准精度，需要进行更为精细的配准。所述步骤四的主要过程在于对初步配准图像Ic及目标图像It进行精细配准，构建出控制初步配准图像Ic产生非刚度形变的精细配准映射矩阵Mf，将Mf作用于初步配准图像Ic，构建出精细配准图像Icf。具体实施过程详细描述如下：
(4-1)重新提取初步配准图像特征点集
利用图像分割算法，提取出初步配准图像Ic中的小鼠表皮区域及骨骼区域，并利用边缘检测算法分别提取出所述小鼠表皮区域及骨骼区域的轮廓；将小鼠表皮及骨骼的轮廓相互叠加，并使用矩形网格进行采样，由此得到一组初步配准图像特征点集L'c；
(4-2)基于块匹配方法的特征点配准
利用块匹配方法求取初步配准图像特征点集L'c中每个特征点在目标图像上的对应位置，以L'c中的任意一点pl＝(xp,yp,zp)为例，将所述过程描述如下：
在初步配准图像Ic中以坐标(xp,yp,zp)为中心选取尺寸为N1×N1×N1的立方邻域T，在目标图像It中以坐标(xp,yp,zp)为中心选取N2×N2×N2的立方邻域S，其中N2>N1；将T作为模板，S作为搜索域，在S区域中搜索与T具有最大相似度的子区域s1，并将s1的中心点pl′作为点pl的对应位置；
以此类推，依次寻找出初步配准图像特征点集L'c中各点在目标图像It上的对应位置，由此构建出精细配准特征点集Lf；
(4-3)精细配准局部位移矩阵Q的构建
利用初步配准图像特征点集L'c及精细配准特征点集Lf构建精细配准局部位移矩阵Q：
Q＝[L'c,L'c-Lf]＝[x,y,z,Δx,Δy,Δz] (8)
(4-4)精细配准映射矩阵Mf的构建
将上述矩阵Q转化为三组四维数据点集Qx＝{(x,y,z,△x)}，Qy＝{(x,y,z,△y)}，Qz＝{(x,y,z,△z)}；与步骤三中多层次B样条拟合过程一致，分别对所述三组四维数据点集Qx＝{(x,y,z,△x)}，Qy＝{(x,y,z,△y)}进行多层次B样条拟合，从而分别计算出初步配准图像Ic中各个像素沿x、y、z三个轴向的位移量，由此构建初步配准图像Ic的精细配准映射矩阵Mf；
至此，精细配准映射矩阵Mf获取完毕。
(4-5)精细配准图像If的构建(可选)
此过程之后，可以选择将精细配准映射矩阵Mf作用于初步配准图像Ic，构建出精细配准图像If；在构建初步配准图像时，灰度的赋值同样采用三线性插值方法；通过判断最终配准图像(即精细配准图像If)与目标图像It的相似度以评定配准的有效程度；
步骤五、构建目标鼠解剖结构图谱：
在非刚度图像配准过程中，通过初步配准映射矩阵Mc及精细配准映射矩阵Mf的控制作用，参考图像Ir变形至精细配准图像If，精细配准图像If与目标图像It具有较高相似度。因此，将初步配准映射矩阵Mc及精细配准映射矩阵Mf依次作用于参考鼠的解剖结构图谱Ar，则可以近似构建目标鼠的解剖结构图谱。然而，参考解剖结构图谱建立于物理坐标系下(以mm为单位)，而上述两个配准映射矩阵建立于像素坐标系下(以像素为单位)，因此在此过程中需要对Ar进行坐标变换。具体实施过程详细描述如下：
(5-1)将所述参考鼠解剖结构图谱Ar投影至像素坐标系，得到像素坐标系下的参考鼠解剖结构图谱
(5-2)按顺序依次将初步配准映射矩阵Mc及精细配准映射矩阵Mf作用于所述像素坐标系下的参考鼠解剖结构图谱使所述像素坐标系下的参考鼠解剖结构图谱产生与步骤三初步配准及步骤四精细配准过程相同的变形，获得像素坐标系下的配准鼠解剖结构图谱
(5-3)将所述像素坐标系下的配准鼠解剖结构图谱投影至物理坐标系下，得到物理坐标系下的配准鼠解剖结构图谱Af，所述物理坐标系下的配准鼠解剖结构图谱Af即为目标鼠解剖结构图谱。
尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (1)

1.一种基于XCT图像非刚度配准的目标鼠解剖结构图谱获取方法，其特征在于，其基本步骤是：首先，将Digimouse模型设定为参考鼠解剖结构图谱，同时将Digimouse对应的XCT图像设定为参考图像；其次，对目标鼠进行XCT成像获取目标图像并进行预处理；然后，利用非刚度图像配准技术构建参考图像至目标图像的配准映射矩阵；最终，将配准映射矩阵作用于参考鼠解剖结构图谱，构建出目标鼠解剖结构图谱；具体步骤如下：
步骤一、设定参考鼠解剖结构图谱和参考图像：
将Digimouse结构设定为参考鼠解剖结构图谱Ar，将所述Digimouse的XCT图像设定为参考图像Ir；
步骤二、目标鼠XCT图像的获取与预处理：
利用XCT设备对目标鼠进行全身成像，获取目标鼠的XCT图像；对目标鼠XCT图像进行仿射变换，将目标鼠XCT图像中小鼠的头部与背部的朝向变换至与参考图像Ir相同；将变换后的XCT图像作为目标图像It；
步骤三、构建初步配准映射矩阵Mc和初步配准图像Ic：
利用图像分割算法分别分割出参考图像Ir中的小鼠骨骼区域与表皮区域，及目标图像It中的小鼠骨骼区域与表皮区域；
利用边缘检测算法分别提取上述参考图像Ir中小鼠骨骼区域和表皮区域及目标图像It中的小鼠骨骼区域和表皮区域共计四个区域的边界轮廓，对所述四个区域的边界轮廓进行等概率采样，从而计算出参考鼠骨骼特征点集Lrb、参考鼠表皮特征点集Lrs、目标鼠骨骼特征点集Ltb以及目标鼠表皮特征点集Lts；
利用TPS-RPM算法的开源程序，计算参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb；在计算时，将目标鼠骨骼特征点集Ltb设定为目标点集，将参考鼠骨骼特征点集Lrb设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代计算出骨骼特征点对应矩阵Cb；
利用TPS-RPM算法的开源程序，计算参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs；在计算时，将目标鼠表皮特征点集Lts设定为目标点集，将参考鼠表皮特征点集Lrs设定为待配准点集，并设定TPS-RPM所用模拟退火算法中的初始温度系数及迭代终止条件，通过迭代计算出表皮特征点对应矩阵Cs；
利用TPS-RPM算法的开源程序求取骨骼特征点对应矩阵Cb与表皮特征点对应矩阵Cs的具体过程如下：
所述参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb及参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs均为模糊对应矩阵，在模糊对应矩阵中各个元素均为[0,1]之间的浮点数，用以描述两点之间对应程度的强弱；
所述TPS-RPM算法在计算参考鼠骨骼特征点集Lrb与目标鼠骨骼特征点集Ltb之间的骨骼特征点对应矩阵Cb时，利用模拟退火算法计算如下能量方程的极小值点：
<mrow> <mi>E</mi> <mrow> <mo>(</mo> <msub> <mi>C</mi> <mi>b</mi> </msub> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>b</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>l</mi> <mrow> <mi>r</mi> <mi>b</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>-</mo> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>l</mi> <mrow> <mi>t</mi> <mi>b</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>L</mi> <mi>f</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>T</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>b</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mi>log</mi> <mi> </mi> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mi>&amp;gamma;</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>b</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
式(1)中，cb,ij为对应矩阵Cb中的元素；Nrb为参考鼠骨骼特征点集Lrb所包含特征点的个数，Ntb为目标鼠骨骼特征点集Ltb所包含点特征的个数；lrb,i为参考鼠骨骼特征点集Lrb中的第i个特征点的坐标，ltb,j为目标鼠骨骼特征点集Ltb中的第j个特征点的坐标；f为薄板样条函数；||Lf||2为对f的平滑性约束，其中L为微分算子；T为模拟退火过程中用于控制模糊对应程度的温度系数；λ为先验平滑度权值；γ为鲁棒性控制权值；随着温度系数T逐渐减小，对应矩阵Cb与薄板样条函数f交替迭代；对应矩阵Cb的最终迭代结果即为所求骨骼特征点对应矩阵Cb；
所述TPS-RPM算法在计算参考鼠表皮特征点集Lrs与目标鼠表皮特征点集Lts之间的表皮特征点对应矩阵Cs时，利用模拟退火算法计算如下能量方程的极小值点：
<mrow> <mi>E</mi> <mrow> <mo>(</mo> <msub> <mi>C</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>s</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>l</mi> <mrow> <mi>r</mi> <mi>s</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>-</mo> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>l</mi> <mrow> <mi>t</mi> <mi>s</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>L</mi> <mi>f</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>T</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>s</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mi>log</mi> <mi> </mi> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mi>&amp;gamma;</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>s</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
式(2)中，cs,ij为对应矩阵Cs中的元素；Nrs为参考鼠表皮特征点集Lrs所包含特征点的个数，Nts为目标鼠表皮特征点集Lts所包含点特征的个数；lrs,i为参考鼠表皮特征点集Lrs中的第i个点的坐标，lts,j为目标鼠表皮特征点集Lts中的第j个点的坐标；其余变量含义与式(1)中相同；随着温度系数T逐渐减小，对应矩阵Cs与薄板样条函数f交替迭代；对应矩阵Cs的最终迭代结果即为所求表皮特征点对应矩阵Cs；
分别将求得的骨骼特征点对应矩阵Cb及表皮特征点对应矩阵Cs作用于参考鼠骨骼特征点集Lrb及参考鼠表皮特征点集Lrs，得到初步配准后的骨骼特征点集Lcb及初步配准后的表皮特征点集Lcs；
利用上述求得的初步特征点配准的结果，构建初步配准局部位移矩阵P：
<mrow> <mi>P</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </mtd> <mtd> <mrow> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>L</mi> <mrow> <mi>c</mi> <mi>b</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </mtd> <mtd> <mrow> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>L</mi> <mrow> <mi>c</mi> <mi>s</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mo>&amp;lsqb;</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>&amp;Delta;</mi> <mi>x</mi> <mo>,</mo> <mi>&amp;Delta;</mi> <mi>y</mi> <mo>,</mo> <mi>&amp;Delta;</mi> <mi>z</mi> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
将上述矩阵P转化为三组四维数据点集Px＝{(x,y,z,Δx)}，Py＝{(x,y,z,Δy)}，Pz＝{(x,y,z,Δz)}，将Δx、Δy及Δz分别视为点(x,y,z)的函数值，即Δx＝G1(x,y,z)，Δy＝G2(x,y,z)，Δz＝G3(x,y,z)；
采用多层次B样条拟合算法分别对所述三组四维数据点集Px＝{(x,y,z,Δx)}，Py＝{(x,y,z,Δy)}，Pz＝{(x,y,z,Δz)}进行拟合，用于计算所述参考图像Ir中各个像素沿着x、y、z三个方向上的位移量；
利用所述多层次B样条拟合方法拟合数据点集Px＝{(x,y,z,Δx)}的过程详细描述如下：
所述多层次是指，利用覆于参考图像Ir上的一组逐渐加密的立方控制网格Φ0,Φ1,...,Φk,...,Φh依次对迭代更新的四维数据点集进行B样条拟合，并将所求h层拟合函数之和作为最终多层次B样条拟合函数；其中，Δ0ξ＝Δx，Δk+1ξ＝Δkξ-gk(x,y,z)，gk(x,y,z)为第k层B样条拟合结果；
所述第k层B样条拟合过程描述如下：
假设第k层控制网格Φk尺寸为Kx×Ky×Kz，则第k层B样条拟合函数如下式所示：
<mrow> <msub> <mi>g</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <msub> <mi>B</mi> <mi>l</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>z</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&amp;phi;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
式(4)中， φk,(l+i,m+j,n+k)为位于控制网格Φk中坐标为(l+i,m+j,n+k)的控制节点值；l,m,n∈{0,1,2,3}；Bl、Bm及Bn分别为第l,m,n阶B样条基函数，其中0至3阶B样条函数的表达式描述如下:
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>(</mo> <mn>3</mn> <msup> <mi>&amp;delta;</mi> <mn>3</mn> </msup> <mo>-</mo> <mn>6</mn> <msup> <mi>&amp;delta;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mo>)</mo> </mrow> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>(</mo> <mo>-</mo> <mn>3</mn> <msup> <mi>&amp;delta;</mi> <mn>3</mn> </msup> <mo>+</mo> <mn>3</mn> <msup> <mi>&amp;delta;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>3</mn> <mi>&amp;delta;</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&amp;delta;</mi> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
所述式(4)中，控制网格Φk中各个控制节点的取值通过以下两步计算：
(a)计算中每个数据点对控制网格Φk中各个控制节点取值的影响量：
以中的一个数据点p＝(xp,yp,zp,Δkξp)为例进行说明如下：
数据点p＝(xp,yp,zp,Δkξp)对控制网格Φk中各个控制节点的影响矩阵表现为一个尺寸为Kx×Ky×Kz的矩阵Ψp；为计算简便，定义与矩阵Ψp尺寸相同的两个矩阵Γp与Ωp；所述矩阵Ψp、Γp与Ωp中坐标为(l+i,m+j,n+k)的元素分别通过式(6)计算：
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;psi;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <msup> <mi>&amp;Delta;</mi> <mi>k</mi> </msup> <msub> <mi>&amp;xi;</mi> <mi>p</mi> </msub> </mrow> <msub> <mi>&amp;omega;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <msub> <mi>B</mi> <mi>l</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>x</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>y</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>z</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;omega;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </msubsup> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </msubsup> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>l</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>x</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>y</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>z</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
式(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样条拟合函数gk(x,y,z)确立；综合各层次拟合函数，计算出多层次B样条拟合函数利用所述求取的多层次B样条拟合函数g(x,y,z)，计算出参考图像Ir中各个像素沿x轴的位移；
同理，利用多层次B样条拟合四维数据点集Py＝{(x,y,z,Δy)}和Pz＝{(x,y,z,Δz)}，从而计算出参考图像Ir中各个像素沿y、z轴向中的位移量，由此构建参考图像Ir的初步配准映射矩阵Mc；
利用所述初步配准映射矩阵Mc，反向求解出初步配准图像Ic；在构建初步配准图像时，灰度的赋值采用三线性插值方法；
步骤四、构建精细配准映射矩阵Mf：
利用图像分割算法，提取出初步配准图像Ic中的小鼠表皮区域及骨骼区域，并利用边缘检测算法分别提取出所述小鼠表皮区域及骨骼区域的轮廓；将小鼠表皮及骨骼的轮廓相互叠加，并使用矩形网格进行采样，由此得到一组初步配准图像特征点集L'c；
利用块匹配方法求取初步配准图像特征点集L'c中每个特征点在目标图像上的对应位置，以L'c中的任意一点pl＝(xp,yp,zp)为例，将所述过程描述如下：
在初步配准图像Ic中以坐标(xp,yp,zp)为中心选取尺寸为N1×N1×N1的立方邻域T，在目标图像It中以坐标(xp,yp,zp)为中心选取N2×N2×N2的立方邻域S，其中N2>N1；将T作为模板，S作为搜索域，在S区域中搜索与T具有最大相似度的子区域s1，并将s1的中心点pl′作为点pl的对应位置；
以此类推，依次寻找出初步配准图像特征点集L'c中各点在目标图像It上的对应位置，由此构建出精细配准特征点集Lf；
利用初步配准图像特征点集L'c及精细配准特征点集Lf构建精细配准局部位移矩阵Q：
Q＝[L'c,L'c-Lf]＝[x,y,z,Δx,Δy,Δz] (8)
将上述矩阵Q转化为三组四维数据点集Qx＝{(x,y,z,Δx)}，Qy＝{(x,y,z,Δy)}，Qz＝{(x,y,z,Δz)}；与步骤三中多层次B样条拟合过程一致，分别对所述三组四维数据点集Qx＝{(x,y,z,Δx)}，Qy＝{(x,y,z,Δy)}进行多层次B样条拟合，从而分别计算出初步配准图像Ic中各个像素沿x、y、z三个轴向的位移量，由此构建初步配准图像Ic的精细配准映射矩阵Mf；
步骤五、构建目标鼠解剖结构图谱：
将所述参考鼠解剖结构图谱Ar投影至像素坐标系，得到像素坐标系下的参考鼠解剖结构图谱按顺序依次将初步配准映射矩阵Mc及精细配准映射矩阵Mf作用于所述像素坐标系下的参考鼠解剖结构图谱使所述像素坐标系下的参考鼠解剖结构图谱产生与步骤三初步配准及步骤四精细配准过程相同的变形，获得像素坐标系下的配准鼠解剖结构图谱将所述像素坐标系下的配准鼠解剖结构图谱投影至物理坐标系下，得到物理坐标系下的配准鼠解剖结构图谱Af，所述物理坐标系下的配准鼠解剖结构图谱Af即为目标鼠解剖结构图谱。
CN201510259870.5A 2015-05-20 2015-05-20 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法 Expired - Fee Related CN104867104B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510259870.5A CN104867104B (zh) 2015-05-20 2015-05-20 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510259870.5A CN104867104B (zh) 2015-05-20 2015-05-20 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法

Publications (2)

Publication Number Publication Date
CN104867104A CN104867104A (zh) 2015-08-26
CN104867104B true CN104867104B (zh) 2017-12-15

Family

ID=53912921

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510259870.5A Expired - Fee Related CN104867104B (zh) 2015-05-20 2015-05-20 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法

Country Status (1)

Country Link
CN (1) CN104867104B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105212936A (zh) * 2015-09-09 2016-01-06 首都医科大学宣武医院 一种脑模板的生成方法
CN105427237B (zh) * 2015-11-25 2018-04-24 长春乙天科技有限公司 一种大幅面光学测量***的钢网图像配准与检测方法
CN109754396B (zh) * 2018-12-29 2021-02-19 上海联影智能医疗科技有限公司 图像的配准方法、装置、计算机设备和存储介质
CN110322491B (zh) * 2019-06-11 2022-03-04 大连理工大学 一种可变形小鼠全身图谱与小鼠影像配准的算法
CN111091567B (zh) * 2020-03-23 2020-06-23 南京景三医疗科技有限公司 医学图像配准方法、医疗设备及存储介质
CN112116642A (zh) * 2020-09-27 2020-12-22 上海联影医疗科技股份有限公司 医学图像的配准方法、电子装置和存储介质
CN112995528B (zh) * 2021-05-06 2021-09-21 中国工程物理研究院流体物理研究所 一种光电分幅相机通道间图像配准方法
CN113298856B (zh) * 2021-05-28 2023-10-20 上海联影医疗科技股份有限公司 一种图像配准方法、装置、设备和介质
CN114742869B (zh) * 2022-06-15 2022-08-16 西安交通大学医学院第一附属医院 基于图形识别的脑部神经外科配准方法及电子设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103077522A (zh) * 2013-01-22 2013-05-01 南方医科大学 一种改进的局部形变明显图像的配准方法
CN104599268A (zh) * 2015-01-06 2015-05-06 广州医科大学附属肿瘤医院 一种联合点配准的局部区域精确变形配准算法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9524552B2 (en) * 2011-08-03 2016-12-20 The Regents Of The University Of California 2D/3D registration of a digital mouse atlas with X-ray projection images and optical camera photos

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103077522A (zh) * 2013-01-22 2013-05-01 南方医科大学 一种改进的局部形变明显图像的配准方法
CN104599268A (zh) * 2015-01-06 2015-05-06 广州医科大学附属肿瘤医院 一种联合点配准的局部区域精确变形配准算法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"一种多模态医学图像数据融合方法与应用";高峰等;《中国医疗设备》;20131231;第28卷(第5期);第164-167页 *
"基于区域标识的‘粗粒度’扩散光学层析成像方法稳态测量模式";赵会娟、高峰等;《光学学报》;20130228;第33卷(第2期);第1-8页 *
"小鼠荧光层析成像***及数据提取扩展方法";赵会娟等;《光子学报》;20141030;第43卷(第10期);第1-7页 *

Also Published As

Publication number Publication date
CN104867104A (zh) 2015-08-26

Similar Documents

Publication Publication Date Title
CN104867104B (zh) 基于xct图像非刚度配准的目标鼠解剖结构图谱获取方法
Cao et al. Deformable image registration using a cue-aware deep regression network
CN106920234B (zh) 一种复合式自动放疗计划的方法
Claes et al. Computerized craniofacial reconstruction: conceptual framework and review
CN101339670B (zh) 一种计算机辅助的三维颅面复原方法
CN110459301A (zh) 基于热力图和面部关键点的脑部神经外科导航配准方法
CN106204561A (zh) 基于混合模型的***多模态图像非刚性配准方法
CN102903103B (zh) 基于迁移活动轮廓模型的胃部ct序列图像分割方法
CN113570627B (zh) 深度学习分割网络的训练方法及医学图像分割方法
CN107507195A (zh) 基于超图模型的pet‑ct多模态鼻咽癌图像分割方法
CN109118455B (zh) 一种基于现代人软组织分布的古人类头骨颅面交互复原方法
Xie et al. Feature‐based rectal contour propagation from planning CT to cone beam CT
CN107220965A (zh) 一种图像分割方法及***
Desvignes et al. 3D semi-landmarks based statistical face reconstruction
CN114792326A (zh) 一种基于结构光的手术导航点云分割与配准方法
Wang et al. Deformable torso phantoms of Chinese adults for personalized anatomy modelling
CN106709867A (zh) 一种基于改进surf和改进互信息的医学图像配准方法
Markel et al. A 4D biomechanical lung phantom for joint segmentation/registration evaluation
Wang et al. Using optimal transport to improve spherical harmonic quantification of complex biological shapes
CN112598669B (zh) 一种基于数字人技术的肺叶分割方法
CN111798500B (zh) 一种基于层次邻域谱特征的微分同胚非刚性配准算法
CN110322491B (zh) 一种可变形小鼠全身图谱与小鼠影像配准的算法
CN111862320B (zh) 一种spect三维重建图像到标准视图的自动转向方法
Bond et al. Series of 4D adult XCAT phantoms for imaging research and dosimetry
CN114913149B (zh) 一种基于ct影像的头部可变形统计图谱构建方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
EXSB Decision made by sipo to initiate substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171215

Termination date: 20200520