一种基于图像特征增强数字体图像相关法的变形测量方法

1.本发明涉及变形测量的技术领域，尤其是指一种基于图像特征增强数字体图像相关法的变形测量方法。


背景技术：

2.准确的位移测量一直都是实验力学中最基本的任务，对通过对不同阶段的试样进行测量，得出两个阶段之间的位移数值。目前工程中广泛应用的位移测量方法主要是以应变电测为代表的接触式测量方法。在使用应变电测时，通过将应变计附着到待测试样表面上，当试样受到外界影响发生变形时，应变计一同产生变形，此时应形计中的电阻值发生相应的变化，通过电阻应变仪测量并转换应变计中的电阻值，即可得到最终的变形结果。然而这种方法仅能用于单点测量，因此在测量前需要进行精确的定位放置。同时该方式仅能测量试样表面位移，而对试样内部位移束手无策。
3.随着诸如磁共振成像、光学相干断层扫描和共聚焦激光扫描显微镜等三维成像技术的发展，变形测量从原来的试样表面变形测量扩展到了试样内部变形测量，然而，现有的成像技术也并非十分完美，因为其成像过程十分耗时，这就导致待测量的变形是大而复杂的内部变形。在这种情况下容易出现脱相关现象，即配准算法可能无法在变形前后获得的体图像中定位到相应的体块。
4.数字体图像相关法是一种非接触式光测力学方法，其基本原理是通过采集和分析变形前后的试样体图像得到试样内部三维全场位移情况，从而避免了传统接触式测量的繁琐操作，同时在不破坏待测试样的情况下对试样内部变形进行测量。brisk算法是一种特征提取描述算法，其通过构建尺度空间搜索关键点，并利用周围像素点进行二进制描述，这些特征点往往具有尺度不变和旋转不变性，常用以计算机视觉的匹配任务中。


技术实现要素：

5.本发明的目的在于克服现有技术的缺点与不足，提出了一种基于图像特征增强数字体图像相关法的变形测量方法，在图像特征提取描述算法的基础之上进行扩展，并与数字体图像相关法进行结合，扩展的图像特征提取描述算法的复杂度将远小于数字体图像相关法的复杂度。在保证计算效率的同时，进一步处理大而复杂的变形情况。
6.为实现上述目的，本发明所提供的技术方案为：一种基于图像特征增强数字体图像相关法的变形测量方法，包括以下步骤：
7.s1：取变形前三维图像为参考图像，变形后三维图像为目标图像，分别构造参考图像的尺度空间和目标图像的尺度空间，所述参考图像的尺度空间和目标图像的尺度空间均包含2n层，其中n为正整数；
8.s2：取步骤s1中生成的参考图像的尺度空间和目标图像的尺度空间，对两种图像的尺度空间中的每个像素点进行检测，将通过检测的像素点作为待定关键点，计算待定关键点得分，利用非极大值抑制方法对待定关键点进行筛选，得到两种图像的最终可用关键
点；
9.s3：取步骤s2中提取到的两种图像的可用关键点，对于每一个关键点，利用其周围像素点信息计算出旋转矩阵，用于表示其主方向；
10.s4：取步骤s2中提取到的两种图像的可用关键点及步骤s3中计算的旋转矩阵，对于每一个关键点，使用其对应的旋转矩阵将其周围的像素点进行旋转，利用旋转后的像素点计算出二进制描述符，用于描述该关键点，形成三维图像的图像特征；
11.s5：取步骤s2中提取到的两种图像的可用关键点及步骤s4中计算的二进制描述符，枚举计算参考图像中提取到的关键点以及目标图像中提取到的关键点对应描述符之间的汉明距离，以汉明距离最小为最佳匹配结果，形成关键点对；
12.s6：取步骤s5中计算的关键点对，选用三维一阶形函数对兴趣点进行描述，使用兴趣点周围的关键点对对三维一阶形函数进行拟合，形成变形初值，再选用零均值归一化平方误差和作为ic-gn算法的目标函数，将变形初值作为ic-gn算法的初值进一步迭代得到最终的变形结果。
13.进一步，在步骤s1中，三维图像是指由磁共振成像、光学相干断层扫描和共聚焦激光扫描显微镜三维成像技术描扫形成的图像。
14.进一步，在步骤s1中，尺度空间是指由多个不同尺度的相同图像所组成的图像层scalespace，表示为：
15.scalespace＝{c0,d0,c1,d1,...,c
n-1
,d
n-1
}
16.式中，c0,c1,...,c
n-1
是指n层octave层，d0,d1,...,d
n-1
是指n层intra-octave层，共有2n层，其中octave层表示为：
[0017][0018]
式中，ci表示为第i层octave层，image表示构造尺度空间的三维图像，dosampling(c
i-1
)表示将第i-1层octave层进行2倍下采样操作；
[0019]
而intra-octave层则表示为：
[0020][0021]
式中，di表示为第i层intra-octave层，dowsamping
1.5
(image)表示将构造尺度空间的三维图像进行1.5倍下采样操作，dosampling(d
i-1
)表示将第i-1层intra-octave层进行2倍下采样操作；
[0022]
尺度空间中的各层有其相应的尺度，表示为：
[0023]
t(ci)＝2i,i＝0,1,...,n-1
[0024]
t(di)＝1.5
×2i
,i＝0,1,...,n-1
[0025]
式中，t(ci)表示第i层octave层的尺度，t(di)表示为第i层intra-octave层的尺度。
[0026]
进一步，在步骤s2中，像素点是三维图像中的坐标点，表示为：
[0027]
pi＝(xi,yi,zi,si),i＝0,1,...,m-1
[0028]
式中，pi是指第i个像素点，xi,yi,zi分别是pi在像素坐标系下x轴、y轴及z轴分量，
si为pi所在层的尺度，共有m个点，m为正整数；而检测则分为粗检测以及细检测两种检测，依次进行；粗检测首先取待检测点周围6个点，比较6个点与待检测点之间的灰度值差的绝对值是否超过阈值，统计满足的个数，若大于或等于5个点满足，则通过检测，否则不通过检测，表示为：
[0029][0030]
式中，test
coarse
(pi)表示对pi进行粗检测，其中1表示为通过测试，0表示为未通过测试，q(pi)表示x轴、y轴及z轴三个方向与pi距离为3的点集，qj表示该点集中的第j个点，valid(pi,qj)表示pi处的灰度值与qj处的灰度值之差的绝对值是否大于自设定的阈值t；其中，valid(pi,qj)表示为：
[0031][0032]
式中，i(pi)、i(qj)分别表示尺度空间中pi、qj处的灰度值；而q(pi)则表示为：
[0033]
q(pi)＝{(x
i-3,yi,zi),(xi 3,yi,zi),(xi,y
i-3,zi),(xi,yi 3,zi),(xi,yi,z
i-3),(xi,yi,zi 3)}
[0034]
细检测首先取与待检测点距离在2.5～3.5之间的整数点，共计98个点，比较该98个点与待检测点之间的灰度值差的绝对值是否超过阈值，统计满足的个数，若大于或等于50个点满足，且满足的点连通则通过检测，否则不通过检测，表示为：
[0035][0036]
式中，test
fine
(pi)表示对pi进行细检测，其中1表示为通过测试，0表示为未通过测试，s(pi)表示与pi距离在2.5～3.5之间的整点集中满足与pi之间的灰度值差的绝对值是否超过阈值的点集，connect(s(pi))表示点集s(pi)是否连通；其中，s(pi)表示为：
[0037][0038]
式中，表示全体整数组成的集合，δp表示x轴、y轴及z轴分量为整数的坐标点；
[0039]
而connect(s(pi))则表示为：
[0040]
使得||q
j-qk||＝1
[0041]
式中，表示对于s(pi)中任意一点，表示s(pi)中存在一点。
[0042]
进一步，所述步骤s2包括以下步骤：
[0043]
s21：取步骤s1中的参考图像的尺度空间和目标图像的尺度空间，对于尺度空间中的所有像素点pi，使用test
coarse
粗检测方法对pi进行检测，得到一个粗关键点集；
[0044]
s22：取步骤s21中的粗关键点集，对于其中的所有关键点，使用test
fine
细检测方法对关键点进行检测，得到待定关键点；
[0045]
s23：取步骤s22中的待定关键点集，对待定关键点中的所有关键点进行非极大值抑制，非极大值抑制中采用的得分表示为：
[0046][0047]
式中，socre(pi)表示pi与s(pi)之间的灰度值之差的绝对值之和；qj表示s(pi)中第j个点，i(pi)、i(qj)分别表示尺度空间中pi、qj处的灰度值；取关键点上、下、左、右、前及后相邻的6个像素点，分别计算它们的得分，判断该关键点的得分是否比这6个像素点的得分高，若全大于则为可用关键点，否则不是可用关键点。
[0048]
进一步，所述步骤s3包括以下步骤：
[0049]
s31：利用采样方式对关键点周围像素点进行选取，其中采样方式为以关键点为中心，设定4个半径r1,r2,r3,r4构成四个同心球面，利用等边三角形将球面进行三角化，以三角化后的顶点作为采样点，分别在四个球面上采样u1,u2,u3,u4个点，共计采样u＝u1 u2 u3 u4 1个点，即四个球面上的点及关键点所在的一个点，则对第i个可用关键点进行采样后的采样点集表示为sampling(kpi)；
[0050]
s32：取步骤s31中对第i个可用关键点进行采样后的采样点集sampling(kpi)，将其两两配对得到u
×
(u-1)
÷
2对采样点对，表示为：
[0051]
pair(kpi)＝{《rj,rk》|rj∈samping(kpi),rk∈samping(kpi)}
[0052]
式中，pair(kpi)表示第i个可用关键点所生成的采样点对，rj,rk为从samping(kpi)中拿到的采样点；计算采样点对之间的距离，大于阈值a的采样点对被称为长距离点对，形成长距离集，表示为：
[0053]
long(kpi)＝{《rj,rk》||||r
j-rk||＞a,《rj,rk》∈pair(kpi)}
[0054]
式中，long(kpi)表示第i个可用关键点所生成的长距离集，《rj,rk》为从pair(kpi)中拿到的采样点对；
[0055]
s33：取步骤s32中的长距离集，求长距离集中点对梯度的平均值，表示为：
[0056][0057]
式中，g(kpi)表示长距离集kpi的梯度平均值，g
x
,gy,gz表示梯度在x轴、y轴、z轴的分量，gradient(rj,rk)表示rj,rk之间的梯度值：
[0058][0059]
式中，scale为缩放系数，人为设定，sj,sk分别表示rj,rk中的尺度分量，g(rj,scale
×
sj)表示在rj处采用标准差为scale
×
sj的三维高斯模糊后的灰度值，g(rk,scale
×
sk)表示在rk处采用标准差为scale
×
sk的三维高斯模糊后的灰度值；
[0060]
利用梯度平均值构造结构张量，利用奇异值分解方法求取该结构张量的特征向量，表示为：
[0061][0062]
式中，为构造的结构张量，v1,v2,v3为该结构张量的特征向量；选取特征向量中的v1,v2及v1×
v2构成旋转矩阵rotate(kpi)作为关键点的主方向：
[0063]
rotate(kpi)＝[v
1 v
2 v1×
v2]。
[0064]
进一步，所述步骤s4包括以下步骤：
[0065]
s41：利用采样方式对关键点周围像素点进行选取，得到第i个可用关键点进行采样后的采样点集sampling(kpi)；
[0066]
步骤s42：取步骤s41中对第i个可用关键点进行采样后的采样点集sampling(kpi)，将其两两配对得到第i个可用关键点的采样点对pair(kpi)，计算采样点对之间的距离，小于阈值b的采样点对被称为短距离点对，形成短距离集，表示为：
[0067]
short(kpi)＝{《rj,rk》||||r
j-rk||＜b,《rj,rk》∈pair(kpi)}
[0068]
式中，short(kpi)表示第i个可用关键点所生成的短距离集，《rj,rk》为从pair(kpi)中拿到的采样点对；
[0069]
s43：取步骤s42中的短距离集short(kpi)及步骤s33中的旋转矩阵rotate(kpi)，将采样点进行旋转后，比较短距离点对之间的大小关系作为kpi的二进制描述符，表示为：
[0070][0071]
式中，descriptor(kpi)表示kpi的二进制描述符，l表示《rj,rk》是short(kpi)第l对点对，compare(rj,rk)表示比较rj,rk旋转后的灰度值大小；其中，compare(rj,rk)表示为：
[0072][0073]
式中，scale为缩放系数，人为设定；sj,sk分别表示rj,rk中的尺度分量；g(rotate(kpi)rj,scale＝sj)表示在rj处使用rotate(kpi)旋转后采用标准差为scale
×
sj的三维高斯模糊后的灰度值；g(rotate(kpi)rk,scale
×
sk)表示在rk处使用rotate(kpi)旋转后采用标准差为scale
×
sk的三维高斯模糊后的灰度值。
[0074]
进一步，所述步骤s6包括以下步骤：
[0075]
步骤s61：取步骤s5中关键点对，对兴趣点poi进行初值估计，其中兴趣点poi为参考图像中的像素点，表示为：
[0076]
poii＝(xi,yi,zi,1)
[0077]
式中，poii表示参考图像中的第i个兴趣点；xi,yi,zi分别表示poii在x轴、y轴及z轴上的分量，i取正整数；选用三维一阶形函数对变形描述：
[0078][0079]
式中，poii'表示为poii对应目标图像中的像素点；u,v,w分别为沿x轴、y轴及z轴方向的位移；u
x
,uy,uz分别表示x轴方向上的位移沿x轴、y轴及z轴方向的梯度分量；v
x
,vy,vz分别表示y轴方向上的位移沿x轴、y轴及z轴方向的梯度分量；w
x
,wy,wz分别表示z轴方向上的位移沿x轴、y轴及z轴方向的梯度分量；u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz为待定参数，由后续步骤拟合得到；以兴趣点poii为中心，划定subu×
subv×
subw大小的子体块，其中subu,subv,subw分别表示子体块在x轴、y轴及z轴方向上的大小，人为指定，获取子体块内的关键点对作为邻关键点对，代入到三维一阶形函数中进行拟合，得到该兴趣点poii的变形初值u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz；
[0080]
s62：取步骤s61中的兴趣点poii的变形初值u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz，通过三维一阶形函数计算得到poi
′i；以poii为中心，在参考图像上划定大小为subu×
subv×
subw的子体块以poii'为中心，在目标图像上划定大小为subu×
subv×
subw的子体块选用零均值归一化平方误差和计算子体块与子体块之间像素点的灰度值的相似度c
znssd
：
[0081][0082]
式中，为变形前子体块中第i个像素点的零均值化灰度值，为变形后子体块中第i个像素点的零均值化灰度值；零均值化灰度值表示为：
[0083][0084]
式中，fi,hi分别表示变形前后第i个像素点灰度值，n为子体块中像素点个数；
[0085]
将c
znssd
作为ic-gn迭代算法的目标函数，迭代更新计算u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz，直至目标函数小于迭代阈值λ或迭代次数大于迭代次数阈值γ，λ、γ人为设定；将迭代结果作为最终变形结果。
[0086]
本发明与现有技术相比，具有如下优点与有益效果：
[0087]
1、简单方便。利用试样变形前后的三维图像，可以在不破坏试样的情况下进行非接触式测量，更加方便在测量试样内部的三维位移情况。
[0088]
2、外界影响小，稳定性高。，该图像特征从尺度空间出发，构造主方向，针对图像的旋转、尺度以及亮度变化可以保持不变，因此其计算较为稳定。
[0089]
3、适用范围大。通过引入图像特征，从全局进行匹配，有效地进行解决了大而复杂的变形，同时也能应对小变形情况。
附图说明
[0090]
图1是本发明方法的流程图。
[0091]
图2是参考图像和目标图像的示意图。
具体实施方式
[0092]
下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
[0093]
如图1和图2所示，本实施例所提供的基于图像特征增强数字体图像相关法的变形测量方法，使肜了图像特征及数字体图像相关法，其包括以下步骤：
[0094]
s1：取变形前三维图像为参考图像，变形后三维图像为目标图像，分别构造参考图像的尺度空间和目标图像的尺度空间，所述参考图像的尺度空间和目标图像的尺度空间均包含2n层，其中n为正整数，实施中取4。
[0095]
三维图像是指由磁共振成像、光学相干断层扫描和共聚焦激光扫描显微镜三维成像技术描扫形成的图像。
[0096]
尺度空间是指由多个不同尺度的相同图像所组成的图像层scalespace，表示为：
[0097]
scalespace＝{c0,d0,c1,d1,...,c
n-1
,d
n-1
}
[0098]
式中，c0,c1,...,c
n-1
是指n层octave层，d0,d1,...,d
n-1
是指n层intra-octave层，共有2n层，其中octave层表示为：
[0099][0100]
式中，ci表示为第i层octave层，image表示构造尺度空间的三维图像，dosampling(c
i-1
)表示将第i-1层octave层进行2倍下采样操作；
[0101]
而intra-octave层则表示为：
[0102][0103]
式中，di表示为第i层intra-octave层，dowsamping
1.5
(image)表示将构造尺度空间的三维图像进行1.5倍下采样操作，dosampling(d
i-1
)表示将第i-1层intra-octave层进行2倍下采样操作；
[0104]
尺度空间中的各层有其相应的尺度，表示为：
[0105]
t(ci)＝2i,i＝0,1,...,n-1
[0106]
t(di)＝1.5
×2i
,i＝0,1,...,n-1
[0107]
式中，t(ci)表示第i层octave层的尺度，t(di)表示为第i层intra-octave层的尺度。
[0108]
s2：取步骤s1中生成的参考图像的尺度空间和目标图像的尺度空间，对两种图像的尺度空间中的每个像素点进行检测，将通过检测的像素点作为待定关键点，计算待定关键点得分，利用非极大值抑制方法对待定关键点进行筛选，得到两种图像的最终可用关键点；包括以下步骤：
[0109]
s21：取步骤s1中的参考图像的尺度空间和目标图像的尺度空间，对于尺度空间中的所有像素点pi，使用test
coarse
粗检测方法对pi进行检测，得到一个粗关键点集；
[0110]
s22：取步骤s21中的粗关键点集，对于其中的所有关键点，使用test
fine
细检测方
法对关键点进行检测，得到待定关键点；
[0111]
s23：取步骤s22中的待定关键点集，对待定关键点中的所有关键点进行非极大值抑制，非极大值抑制中采用的得分表示为：
[0112][0113]
socre(pi)表示pi与s(pi)之间的灰度值之差的绝对值之和；qj表示s(pi)中第j个点，i(pi)、i(qj)分别表示尺度空间中pi、qj处的灰度值；取关键点上、下、左、右、前及后相邻的6个像素点，分别计算它们的得分，判断该关键点的得分是否比这6个像素点的得分高，若全大于则为可用关键点，否则不是可用关键点。
[0114]
像素点是三维图像中的坐标点，表示为：
[0115]
pi＝(xi,yi,zi,si),i＝0,1,...,m-1
[0116]
式中，pi是指第i个像素点，xi,yi,zi分别是pi在像素坐标系下x轴、y轴及z轴分量，si为pi所在层的尺度，共有m个点，m为正整数；而检测则分为粗检测以及细检测两种检测，依次进行；粗检测首先取待检测点周围6个点，比较6个点与待检测点之间的灰度值差的绝对值是否超过阈值，统计满足的个数，若大于或等于5个点满足，则通过检测，否则不通过检测，表示为：
[0117][0118]
式中，test
coarse
(pi)表示对pi进行粗检测，其中1表示为通过测试，0表示为未通过测试，q(pi)表示x轴、y轴及z轴三个方向与pi距离为3的点集，qj表示该点集中的第j个点，valid(pi,qj)表示pi处的灰度值与qj处的灰度值之差的绝对值是否大于自设定的阈值t；其中，valid(pi,qj)表示为：
[0119][0120]
式中，i(pi)、i(qj)分别表示尺度空间中pi、qj处的灰度值；而q(pi)则表示为：
[0121]
q(pi)＝{(x
i-3,yi,zi),(xi 3,yi,zi),(xi,y
i-3,zi),(xi,yi 3,zi),(xi,yi,z
i-3),(xi,yi,zi 3)}
[0122]
细检测首先取与待检测点距离在2.5～3.5之间的整数点，共计98个点，比较该98个点与待检测点之间的灰度值差的绝对值是否超过阈值，统计满足的个数，若大于或等于50个点满足，且满足的点连通则通过检测，否则不通过检测，表示为：
[0123][0124]
式中，test
fine
(pi)表示对pi进行细检测，其中1表示为通过测试，0表示为未通过测试，s(pi)表示与pi距离在2.5～3.5之间的整点集中满足与pi之间的灰度值差的绝对值是否
超过阈值的点集，connect(s(pi))表示点集s(pi)是否连通；其中，s(pi)表示为：
[0125][0126]
式中，表示全体整数组成的集合，δp表示x轴、y轴及z轴分量为整数的坐标点；
[0127]
而connect(s(pi))则表示为
[0128]
使得||p
j-pk||＝1∧j≠k
[0129]
式中，表示对于s(pi)中任意一点，表示s(pi)中存在一点。
[0130]
s3：取步骤s2中提取到的两种图像的可用关键点，对于每一个关键点，利用其周围像素点信息计算出旋转矩阵，用于表示其主方向，包括以下步骤：
[0131]
s31：利用采样方式对关键点周围像素点进行选取，其中采样方式为以关键点为中心，设定4个半径r1,r2,r3,r4构成四个同心球面，利用等边三角形将球面进行三角化，以三角化后的顶点作为采样点，分别在四个球面上采样u1,u2,u3,u4个点，共计采样u＝u1 u2 u3 u4 1个点，即四个球面上的点及关键点所在的一个点，则对第i个可用关键点进行采样后的采样点集表示为sampling(kpi)；
[0132]
s32：取步骤s31中对第i个可用关键点进行采样后的采样点集sampling(kpi)，将其两两配对得到u
×
(u-1)
÷
2对采样点对，表示为：
[0133]
pair(kpi)＝{《rj,rk》|rj∈samping(kpi),rk∈samping(kpi)}
[0134]
式中，pair(kpi)表示第i个可用关键点所生成的采样点对，rj,rk为从samping(kpi)中拿到的采样点；计算采样点对之间的距离，大于阈值a的采样点对被称为长距离点对，形成长距离集，表示为：
[0135]
long(kpi)＝{《rj,rk》||||r
j-rk||＞a,《rj,rk》∈pair(kpi)}
[0136]
式中，long(kpi)表示第i个可用关键点所生成的长距离集，《rj,rk》为从pair(kpi)中拿到的采样点对；
[0137]
s33：取步骤s32中的长距离集，求长距离集中点对梯度的平均值，表示为：
[0138][0139]
式中，g(kpi)表示长距离集kpi的梯度平均值，g
x
,gy,gz表示梯度在x轴、y轴、z轴的分量，gradient(rj,rk)表示rj,rk之间的梯度值：
[0140][0141]
式中，scale为缩放系数，人为设定，sj,sk分别表示rj,rk中的尺度分量，g(rj,scale
×
sj)表示在rj处采用标准差为scale
×
sj的三维高斯模糊的灰度值，g(rk,scale
×
sk)表示在rk处采用标准差为scale
×
sk的三维高斯模糊后的灰度值；
[0142]
利用梯度平均值构造结构张量，利用奇异值分解方法求取该结构张量的特征向量，表示为：
[0143][0144]
式中，为构造的结构张量，v1,v2,v3为该结构张量的特征向量；选取特征向量中的v1,v2及v1×
v2构成旋转矩阵rotate(kpi)作为关键点的主方向：
[0145]
rotate(kpi)＝[v
1 v
2 v1×
v2]。
[0146]
s4：取步骤s2中提取到的两种图像的可用关键点及步骤s3中计算的旋转矩阵，对于每一个关键点，使用其对应的旋转矩阵将其周围的像素点进行旋转，利用旋转后的像素点计算出二进制描述符，用于描述该关键点，形成三维图像的图像特征，包括以下步骤：
[0147]
s41：利用采样方式对关键点周围像素点进行选取，得到第i个可用关键点进行采样后的采样点集sampling(kpi)；
[0148]
步骤s42：取步骤s41中对第i个可用关键点进行采样后的采样点集sampling(kpi)，将其两两配对得到第i个可用关键点的采样点对pair(kpi)，计算采样点对之间的距离，小于阈值b的采样点对被称为短距离点对，形成短距离集，表示为：
[0149]
short(kpi)＝{《rj,rk》||||r
j-rk||＜b,《rj,rk》∈pair(kpi)}
[0150]
式中，short(kpi)表示第i个可用关键点所生成的短距离集，《rj,rk》为从pair(kpi)中拿到的采样点对；
[0151]
s43：取步骤s42中的短距离集short(kpi)及步骤s33中的旋转矩阵rotate(kpi)，将采样点进行旋转后，比较短距离点对之间的大小关系作为kpi的二进制描述符，表示为：
[0152][0153]
式中，descriptor(kpi)表示kpi的二进制描述符，l表示《rj,rk》是short(kpi)第l对点对，compare(rj,rk)表示比较rj,rk旋转后的灰度值大小；其中，compare(rj,rk)表示为：
[0154][0155]
式中，scale为缩放系数，人为设定；sj,sk分别表示rj,rk中的尺度分量；g(rotate(kpi)rj,scale
×
sj)表示在rj处使用rotate(kpi)旋转后采用标准差为scale
×
sj的三维高斯模糊后的灰度值；g(rotate(kpi)rk,scale
×
sk)表示在rk处使用rotate(kpi)旋转后采用标准差为scale
×
sk的三维高斯模糊后的灰度值。
[0156]
s5：取步骤s2中提取到的两种图像的可用关键点及步骤s4中计算的二进制描述符，枚举计算参考图像中提取到的关键点以及目标图像中提取到的关键点对应描述符之间的汉明距离，以汉明距离最小为最佳匹配结果，形成关键点对。
[0157]
s6：取步骤s5中计算的关键点对，选用三维一阶形函数对兴趣点进行描述，使用兴趣点周围的关键点对对三维一阶形函数进行拟合，形成变形初值，再选用零均值归一化平
方误差和作为ic-gn算法的目标函数，将变形初值作为ic-gn算法的初值进一步迭代得到最终的变形结果；包括以下步骤：
[0158]
步骤s61：取步骤s5中关键点对，对兴趣点poi进行初值估计，其中兴趣点poi为参考图像中的像素点，表示为：
[0159]
poii＝(xi,yi,zi,1)
[0160]
式中，poii表示参考图像中的第i个兴趣点；xi,yi,zi分别表示poii在x轴、y轴及z轴上的分量，i取正整数；选用三维一阶形函数对变形描述：
[0161][0162]
式中，poii'表示为poii对应目标图像中的像素点；u,v,w分别为沿x轴、y轴及z轴方向的位移；u
x
,uy,uz分别表示x轴方向上的位移沿x轴、y轴及z轴方向的梯度分量；v
x
,vy,vz分别表示y轴方向上的位移沿x轴、y轴及z轴方向的梯度分量；w
x
,wy,wz分别表示z轴方向上的位移沿x轴、y轴及z轴方向的梯度分量；u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz为待定参数，由后续步骤拟合得到；以兴趣点poii为中心，划定subu×
subv×
subw大小的子体块，其中subu,subv,subw分别表示子体块在x轴、y轴及z轴方向上的大小，实施中取值分别为33，33，33获取子体块内的关键点对作为邻关键点对，代入到三维一阶形函数中形成拟合，得到该兴趣点poii的变形初值u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz；
[0163]
s62：取步骤s61中的兴趣点poii的变形初值u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz，通过三维一阶形函数计算得到poi
′i；以poii为中心，在参考图像上划定大小为subu×
subv×
subw的子体块以poi
′i为中心，在目标图像上划定大小为subu×
subv×
subw的子体块选用零均值归一化平方误差和计算子体块与子体块之间像素点的灰度值的相似度c
znssd
：
[0164][0165]
式中，为变形前子体块中第i个像素点的零均值化灰度值，为变形后子体块中第i个像素点的零均值化灰度值；零均值化灰度值表示为：
[0166][0167]
式中，fi,hi分别表示变形前后第i个像素点灰度值，n为子体块中像素点个数；
[0168]
将c
znssd
作为ic-gn迭代算法的目标函数，迭代更新计算u,v,w,u
x
,uy,uz,v
x
,vy,vz,w
x
,wy,wz，直至目标函数小于迭代阈值λ或迭代次数大于迭代次数阈值γ，λ、γ人为设定；将迭代结果作为最终变形结果。
[0169]
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。