一种基于局部和全局联合约束的三维生物图像形变方法

本发明涉及图像形变和图像配准，尤其是一种基于局部和全局联合约束的三维生物图像形变方法。

背景技术：

1、在图像形变与图像配准领域，典型的任务是将一个图像的形态位置变形为另一个图像的形态位置，或者将一个图像与一个参考图像进行形变配准。薄板样条插值形变(tps)是一种常用的方法。tps形变的基本原理是通过在图像上定义一组控制点，并通过对这些控制点进行变形来实现整体的图像形变。tps形变时，需要计算每个像素点相对于控制点的变形位移，然后，根据这些变形位移，可以对图像进行插值，以生成最终的形变图像。

2、在tps形变的基础上可以通过增加约束方案使tps形变效果更为平滑(smooth)，即stps。stps形变可以通过添加更多的控制点来增加形变的精度，但其使用整体的控制点来为stps形变增加平滑约束的时候，会使得局部区域的控制点受到其他区域控制点的平滑约束的影响，导致计算出的位移场产生误差。而且较大尺寸的三维图像，需要更多的控制点集来进行位移场计算，会导致形变时长较长，下采样形变位移场是一个典型的方案，便又产生了下采样位移场的插值问题。为了提高3d图像变形性能与效率，屈等人设计了一种超快图像翘曲工具(littlequickwarp)方案，其采用的方案为应用待配准生物图像与目标生物图像生成一组匹配点集，利用tps算法计算下采样位移场，之后为下采样位移场块进行三线性或b样条插值，插值过程中为了达到整体效果平滑，采用局部缺失位移场块插入大量其他区块位移场值并做平滑操作的处理，该方案在提升形变速度的同时达到整体形变平滑的效果。但是其在计算控制点集之间的平滑约束时，没有有效利用生物医学图像的局部解剖信息，使用整体控制点集计算的约束会导致生物医学图像局部生物组织的位移场出现误差，且该方法在下采样位移场插值过程中采用的方案在区块位移场中引入了大量其他区块位移场值，此操作也会引入形变误差。

技术实现思路

1、为解决三维生物图像形变时局部形变效果精度低的问题，本发明的目的在于提供一种能够有效提高生物图像形变及配准过程中局部形变的精度，同时做到整体形变效果平滑的基于局部和全局联合约束的三维生物图像形变方法。

2、为实现上述目的，本发明采用了以下技术方案：一种基于局部和全局联合约束的三维生物图像形变方法，该方法包括下列顺序的步骤：

3、(1)获取整体匹配点集：基于目标生物图像解剖图谱获取控制点集tar_markers，复制控制点集tar_markers并将其映射到待形变生物图像上，获取待形变图像点集sub_markers，控制点集tar_markers与待形变图像点集sub_markers形成整体匹配点集；

4、(2)根据标签索引信息ind计算局部约束：提取出标签索引信息ind对应的局部的匹配点集，计算标签索引信息ind对应区域的stps形变的仿射项参数矩阵dind以及非仿射项参数矩阵cind；

5、(3)计算全局约束：使用整体匹配点集来计算stps形变的仿射项参数矩阵dall以及非仿射项参数矩阵call；

6、(4)计算基于局部约束与全局约束的下采样位移场df：根据目标生物图像解剖图谱，使用对应的局部约束来计算局部体素位置处的位移场；对于背景区域，使用全局约束计算背景体素位置处的位移场；

7、(5)计算三次b样条基函数矩阵；

8、(6)遍历下采样位移场df，在进行位移场区块初步插值后，通过移动三次b样条基函数矩阵影响区域与初步插值之后的位移场块进行平滑计算，得到最终插值之后的目标图像尺寸的位移场即最终生成的位移场df_whole；

9、(7)根据最终生成的位移场df_whole进行三维生物图像形变，得到形变后的三维生物图像。

10、所述步骤(1)具体包括以下步骤：

11、(1a)在目标生物图像解剖图谱的每个局部区域边界处，利用pcl点云库生成一组点集，并添加对应局部区域的标签索引信息ind，将这组点集整体记为控制点集tar_markers；

12、(1b)复制控制点集tar_markers并将其映射到待形变生物图像上，通过对映射后的点集进行调整，使待形变生物图像中映射后的点集与目标生物图像中的控制点集在对应位置处实现精准匹配，精准匹配后的点集记为待形变图像点集sub_markers；

13、(1c)控制点集tar_markers与待形变图像点集sub_markers形成整体匹配点集；控制点集tar_markers与待形变图像点集sub_markers的尺寸均为ncpt×3，其中，行数ncpt为匹配点对的数量，3列分别存储每个点的x、y、z坐标。

14、所述步骤(2)具体包括以下步骤：

15、(2a)按照标签索引信息ind进行索引，从控制点集tar_markers与待形变图像点集sub_markers中提取出对应标签索引信息ind区域的匹配点集分别存储进局部控制点集cpt_targetind与局部待形变图像点集cpt_subjectind，局部控制点集cpt_targetind与局部待形变图像点集cpt_subjectind的尺寸均为ncptind×3，其中，行数ncptind为标签索引信息ind对应区域匹配点对的数量，3列分别存储每个点的x、y、z坐标；

16、(2b)定义用于存放局部控制点集cpt_targetind每两点之间欧氏距离的矩阵为xnxn_kind，xnxn_kind的尺寸为ncptind×ncptind；计算局部控制点集cpt_targetind每两点之间欧氏距离dlocal_ij，并将欧氏距离结果取反后存储进矩阵xnxn_kind，计算方式如下：

17、dlocal_ij_x＝cpt_targetind[i][0]-cpt_targetind[j][0]；

18、dlocal_ij_y＝cpt_targetind[i][1]-cpt_targetind[j][1]；

19、dlocal_ij_z＝cpt_targetind[i][2]-cpt_targetind[j][2]；

20、

21、其中，cpt_targetind[i][0]为局部控制点集cpt_targetind中第i个点的x坐标，cpt_targetind[i][1]为局部控制点集cpt_targetind中第i个点的y坐标,cpt_targetind[i][2]为局部控制点集cpt_targetind中第i个点的z坐标，cpt_targetind[j][0]为局部控制点集cpt_targetind中第j个点的x坐标，cpt_targetind[j][1]为局部控制点集cpt_targetind中第j个点的y坐标,cpt_targetind[j][2]为局部控制点集cpt_targetind中第j个点的z坐标，dlocal_ij_x为cpt_targetind[i][0]与cpt_targetind[j][0]之间的距离，dlocal_ij_y为cpt_targetind[i][1]与cpt_targetind[j][1]之间的距离，dlocal_ij_z为cpt_targetind[i][2]与cpt_targetind[j][2]之间的距离；

22、(2c)分配两个尺寸为ncptind×4的矩阵xind与yind，分别用于存放局部控制点集cpt_targetind与局部待形变图像点集cpt_subjectind的位置信息，其中，行数ncptind为标签索引信息ind对应区域匹配点对的数量，xind中的4列第一列设值为1，其余三列分别存储局部控制点集cpt_targetind每个点的x、y、z坐标；yind中的4列第一列设值为1，其余三列分别存储局部待形变图像点集cpt_subjectind每个点的x、y、z坐标；

23、(2d)分配一个尺寸为ncptind×ncptind的矩阵qind，矩阵qind前4列与矩阵xind赋值相同；使用qr分解将矩阵qind进一步分解为一个正交矩阵和一个上三角矩阵rind，并将得到的正交矩阵重新赋值给矩阵qind，使用gram-schmidt正交化过程将求得的矩阵qind扩展为具有正交列的矩阵qind；

24、(2e)分配矩阵q1ind、q2ind、rind，将所求得的qind第1到第4列存入矩阵q1ind,将qind第5到第ncptind列存入矩阵q2ind，将所求得的矩阵rind的第1到第4行和第1到第4列存入rind；

25、(2f)定义标签索引信息ind对应区域的stps形变非仿射项参数矩阵为cind,定义标签索引信息ind对应区域的stps形变仿射项参数矩阵为dind,定义cind的尺寸为ncptind×4，dind的尺寸为4×4，cind与dind计算过程如下：

26、alocal＝q2indt·xnxn_kind·q2ind+i·0.2

27、cind＝q2ind·(alocal-1·q2indt·yind)

28、dind＝rind-1·q1indt·(yind-xnxn_kind·cind)

29、其中，i是(ncptind-4)阶的单位矩阵；alocal为局部tps形变的核矩阵；

30、(2g)使用标签索引信息ind重复执行步骤(2a)至步骤(2f)，直到对所有分区的局部约束求解完成。

31、所述步骤(3)具体包括以下步骤：

32、(3a)定义用于存放控制点集tar_markers每两点之间欧氏距离的矩阵为xnxn_k，xnxn_k的尺寸为ncpt×ncpt；计算控制点集tar_markers每两点之间欧氏距离dglobal_ij，并将欧氏距离结果取反后存储进xnxn_k，计算方式如下：

33、dglobal_ij_x＝tar_markers[i][0]-tar_markers[j][0]；

34、dglobal_ij_y＝tar_markers[i][1]-tar_markers[j][1]；

35、dglobal_ij_z＝tar_markers[i][2]-tar_markers[j][2]；

36、

37、其中，tar_markers[i][0]为控制点集tar_markers中第i个点的x坐标，tar_markers[i][1]为控制点集tar_markers中第i个点的y坐标,tar_markers[i][2]为控制点集tar_markers中第i个点的z坐标，tar_markers[j][0]为控制点集tar_markers中第j个点的x坐标，tar_markers[j][1]为控制点集tar_markers中第j个点的y坐标,tar_markers[j][2]为控制点集tar_markers中第j个点的z坐标，dglobal_ij_x为tar_markers0]与tar_markers[j][0]之间的距离，dglobal_ij_y为tar_markers[i][1]与tar_markers[j][1]之间的距离，dglobal_ij_z为tar_markers[i][2]与tar_markers[j][2]之间的距离；

38、(3b)定义两个尺寸为ncpt×4的矩阵x与y，分别用于存放控制点集tar_markers与待形变图像点集sub_markers的位置信息，其中，行数ncpt为匹配点对的数量，x中的4列第一列设值为1，其余三列分别存储tar_markers每个点的x、y、z坐标；y中的4列第一列设值为1，其余三列分别存储sub_markers每个点的x、y、z坐标；

39、(3c)定义一个尺寸为ncpt×ncpt的矩阵q，q前4列与x赋值相同，使用qr分解将矩阵q进一步分解为一个正交矩阵和一个上三角矩阵r，并将得到的正交矩阵重新赋值给q，使用gram-schmidt正交化过程将求得的矩阵q扩展为具有正交列的矩阵q；

40、(3d)定义矩阵q1、q2、r，将所求得的q第1到第4列存入矩阵q1,将q第5到第ncpt列存入矩阵q2，将所求得的矩阵r的第1到第4行和第1到第4列存入r；

41、(3e)定义整体stps形变非仿射项参数矩阵为call,仿射项参数矩阵为dall,定义call的尺寸为ncpt×4，dall的尺寸为4×4；call与dall计算过程如下：

42、aglobal＝q2t·xnxn_k·q2+ia·0.2

43、call＝q2·(aglobal-1·q2 t·y)

44、dall＝r-1·q1 t·(y-xnxn_k·call)

45、其中，ia是(ncpt-4)阶的单位矩阵，aglobal为全局tps形变的核矩阵。

46、所述步骤(4)具体包括以下步骤：

47、(4a)分配一个下采样位移场df，其尺寸按照目标图像尺寸的4倍下采样尺寸确定，并且每个位移场元素包含x、y、z三个方向上的位移场信息；

48、(4b)遍历下采样位移场df的空间位置，记遍历索引为(xdfi,ydfi,zdfi)，用下采样位移场的遍历索引乘以下采样倍数得到真实的当前体素位置(xori,yori,zori)：

49、xori＝xdfi×4

50、yori＝ydfi×4

51、zori＝zdfi×4

52、分配一个用于存储当前体素位置信息的矩阵xori，其尺寸为1×4，其中第一列设值为1，第二至第四列分别存储xori、yori、zori；

53、(4c)选取当前体素位置处用于计算当前控制点集cpt_target_ref以及对应的约束，所述约束包括仿射项参数d_ref以及非仿射项参数c_ref：判断当前体素位置在目标图像解剖图谱上所属的区域，如果当前体素位置属于标签索引信息ind所对应的区域，则选取局部控制点集cpt_targetind作为cpt_target_ref，选取dind作为d_ref，选取cind作为c_ref，以计算该区域的对应位置位移场；如果当前体素位置属于背景区域，则选取控制点集tar_markers作为cpt_target_ref，选取仿射项参数矩阵dall作为d_ref，选取非仿射项参数矩阵call作为c_ref，以计算背景区域的对应位置位移场；

54、(4d)分配用于计算当前体素位置与cpt_target_ref点集欧式距离的矩阵xmxn_k，尺寸为1×ncpt_ref，其中，ncpt_ref为当前控制点集cpt_target_ref中的点的个数，计算当前体素位置与第i个点的欧式距离计算过程如下：

55、dppx_i＝xori-cpt_target_ref[i][0]；

56、dppy_i＝yori-cpt_target_ref[i][1]；

57、dppz_i＝zori-cpt_target_ref[i][2]；

58、

59、其中，cpt_target_ref[i][0]为当前控制点集cpt_target_ref中第i个点的x坐标，cpt_target_ref[i][1]为当前控制点集cpt_target_ref中第i个点的y坐标,cpt_target_ref[i][2]为当前控制点集cpt_target_ref中第i个点的z坐标，dppx_i为当前体素位置的x坐标与cpt_target_ref[i][0]的距离，dppy_i为当前体素位置的y坐标与cpt_target_ref[i][1]的距离,dppz_i为当前体素位置的z坐标与cpt_target_ref[i][2]的距离；

60、(4e)分配用于计算下采样位移场df的stps位移场的矩阵xstps,其尺寸为1×4，并利用当前体素位置处对应约束d_ref、c_ref计算当前体素位置的stps位移场矩阵：

61、xstps＝xori·d_ref+xmxn_k·c_ref

62、(4f)将当前体素位置处计算的stps位移场的矩阵xstps整理到下采样位移场df当前遍历索引(xdfi,ydfi,zdfi)处：

63、df[xdfi][ydfi][zdfi].sx＝xstps(1,2)-xori

64、df[xdfi][ydfi][zdfi].sy＝xstps(1,3)-yori

65、df[xdfi][ydfi][zdfi].sz＝xstps(1,4)-zori

66、其中df[xdif][ydfi][zdfi].sx为下采样位移场df在(xdfi,ydfi,zdfi)索引处的x方向的位移场值，y方向、z方向同理；

67、(4g)通过遍历下采样位移场df的空间位置，重复上述(4a)至(4f)步骤，得到最终基于局部约束与全局约束的下采样位移场df。

68、所述步骤(5)具体包括以下步骤：

69、(5a)分配尺寸为4×4的矩阵b，作为三次b样条基函数的系数矩阵，其值如下：

70、

71、(5b)分配尺寸为4×4的矩阵t，以存放三次b样条基函数的数据：

72、

73、其中t1＝0、

74、(5c)计算尺寸为4×4的初步三次b样条基函数矩阵tb：

75、tb＝t·b

76、(5d)计算tb和tb的kronecker积，记结果矩阵为bxb，其尺寸为42×42：

77、

78、(5e)按照步骤(5d)过程中kronecker积的计算原理，将矩阵bxb和矩阵tb做kronecker积运算，得到尺寸为43×43的用于三维插值的三次b样条基函数矩阵xbspline。

79、所述步骤(6)具体包括以下步骤：

80、(6a)遍历下采样位移场df，遍历索引位置处(xdfi,ydfi,zdfi)对应的真实体素位置(xdfi*4,ydfi*4,zdfi*4)定为startp；

81、(6b)以startp为起点向后找8个存在下采样位移场df的体素位置，这8个体素位置处的位移场在下采样位移场df中索引值分别为：(xdfi,ydfi,zdfi)、(xdfi+1,ydfi,zdfi)、(xdfi,ydfi+1,zdfi)、(xdfi,ydfi,zdfi+1)、(xdfi+1,ydfi+1,zdfi)、(xdfi+1,ydfi,zdfi+1)、(xdfi,ydfi+1,zdfi+1)、(xdfi+1,ydfi+1,zdfi+1)；由于下采样位移场df尺寸为目标图像尺寸四倍下采样的尺寸，上面找到的8个体素位置，每一个后面均存在一个缺乏位移场的4×4×4的体素空间，合起来总共是一个8×8×8的缺乏位移场的体素空间；

82、(6c)对步骤(6b)中找到的8个缺乏位移场的4×4×4的体素空间，每一个缺乏位移场的4×4×4的体素空间均按照体素在目标解剖图谱中的位置进行初步位移场插值：每个缺乏位移场的4×4×4的体素空间只有空间8个顶点处存在下采样位移场df值，先根据目标解剖图谱判断8个顶点所属的区域，再按照目标解剖图谱判断4×4×4的体素空间中的体素与哪个顶点处于同一区域，将对应体素位置处缺乏的位移场值插值为与此体素处于同一区域的顶点处的位移场值；如果体素位置处与8个顶点所处区域均不相同，则将此体素位置处的位移场值插值为与此体素位置最近的顶点处的位移场值；

83、(6d)遍历以startp为起点的4×4×4的体素空间，通过移动三次b样条基函数矩阵xbspline的作用区域对此空间内初步插值之后的位移场块进行平滑计算，具体过程为：遍历的当前体素位置定为currentp，以currentp为起点向后找寻一个经过步骤(6c)初步插值的尺寸为4×4×4的位移场块,将此位移场块命名为df_smallblock，将df_smallblock按位移场值整理尺寸为64×3的矩阵df_smallblock_m，其中64行表示4×4×4位移场块中的64个位移场，3列分别表示每个位置位移场在x、y、z三个方向的位移场分量；再使用三次b样条基函数矩阵xbspline对df_smallblock_m进行平滑计算，将结果存储进矩阵df_b_block：

84、df_b_block＝xbspline·df_smallblock_m

85、将df_b_block结果中的第一行作为currentp体素位置处的最终插值的位移场，其中第一行的三列的值分别作为currentp体素位置处位移场x、y、z方向的分量；

86、(6e)重复步骤(6a)至(6d)直至遍历完下采样位移场df，同时控制取值不超出边界，将得到完整的经过三次b样条基函数插值的且局部与全局联合约束的最终生成的位移场df_whole。

87、所述步骤(7)具体包括以下步骤：

88、(7a)分配一个目标图像尺寸的空白空间img_warp，用以存放图像形变的结果；

89、(7b)遍历空白空间img_warp中的体素位置，使用该体素位置pos_tar的x、y、z的坐标加上最终生成的位移场df_whole对应位置处位移场的x、y、z分量，得到位置pos_sub，将待形变图像中的位置pos_sub处的体素经过三线性插值，赋值到img_warp中pos_tar的位置；

90、(7c)经过步骤(7b)的操作，最终得到基于局部约束与全局约束形变后的三维生物图像img_warp。

91、由上述技术方案可知，本发明的有益效果为：第一，本发明利用生物图像解剖图谱进行局部与全局联合约束的下采样位移场计算，以及通过滑动三次b样条基函数矩阵的影响区域对进行下采样位移场块进行插值及平滑操作，得到最终位移场并进行三维生物图像形变；第二，本发明生成局部与全局相联系的平滑位移场，能够有效提高生物图像形变及配准过程中局部形变的精度，同时做到整体形变效果平滑。