[发明专利]一种基于傅里叶梅林域的二维三维图像配准方法

摘要

基于傅里叶梅林域的术前二维投影图像和术中三维图像配准方法，属医学图像配准领域。该方法可用于骨创伤手术导航系统中，将术前规划与术中三维图像关联达到指导手术的目的。方法利用傅里叶梅林域中图像的尺度，旋转和平移不变性分离待优化参数间的耦合度，按照估计平面外旋转参数，计算平面内旋转参数，估计X射线-三维图像距离，估计三维图像-接收器距离，计算平面内平移的顺序进行配准。该二维三维图像方法为术前数据在术中三维图像中的重建奠定了基础。

主权项

1.一种基于傅里叶梅林域的二维三维配准方法，其特征在于，该方法包含如下步骤：A：将二维真实投影图像I₁(x，y)由DICOM服务器调入配准PC，图像的大小是PxP，其中x是图像的横坐标y是图像的纵坐标且满足0≤x≤P-1；0≤y≤P-1；x，y∈整数，P是二维图像的单边长度。B：将三维图像V由DICOM服务器调入配准PC，三维图像的大小是Width x Height x Depth，其中Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽。C：建立投影模型，包括步骤：C01在三维图像中心建立三维图像直角坐标系A₁A₂A₃，且坐标系原点O位于三维图像中心，坐标轴A₁，A₂，A₃分别正交于三维图像相对应的外平面；C02建立投影模型直角坐标系B₁B₂B₃，投影模型直角坐标系B₁B₂B₃和三维图像直角坐标系A₁A₂A₃重合且B₁，B₂，B₃坐标轴分别与A₁，A₂，A₃坐标轴同向，记投影模型坐标系B₁B₂B₃的坐标原点为ISO；C03虚拟X射线源设定在B₃轴正半轴上，由虚拟X射线源到投影模型坐标系坐标原点ISO的距离记为D1；C04虚拟投影平面垂直于B₃轴，虚拟投影平面和B₃轴的交点位于B₃轴的负半轴，交点位于虚拟投影平面的中心，交点和投影模型坐标系坐标原点ISO的距离记为D2；C05将虚拟X射线源和虚拟投影平面相对三维图像的运动描述为投影模型坐标系B₁B₂B₃在三维图像坐标系A₁A₂A₃内的旋转和平移，包括：绕B₁轴旋转参数R1，绕B₂轴的旋转参数R2，绕B₃轴的旋转参数R3，沿B₁轴的平移T1，沿B₂轴的平移T2；C06投影模型自身的变化由虚拟X射线源到投影模型坐标系坐标原点ISO的距离D1和投影模型坐标系坐标原点ISO到虚拟投影平面距离D2描述；D：初始化待确定的参数，令绕B₁轴的旋转角度R1＝0，绕B₂轴的旋转角度R2＝0，绕B₃轴的旋转角度R3＝0；虚拟X射线源到投影模型坐标系坐标原点距离D1＝2*D其中D是三维图像对角线长度，D=Width2+Height2+Depth2,]]>Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽；投影模型坐标系坐标原点到虚拟投影平面距离D2＝0.5*D；沿B₁轴平移T1＝0，沿B₂轴平移T2＝0；E：将二维真实投影图像作为配准的对象，对绕B₁轴旋转角度R1和绕B₂轴旋转角度R2进行优化，设优化次数为Count，第Count次对绕B₁轴旋转角度R1优化结果为第Count-1次对绕B₁轴旋转角度R1优化结果为第Count次对绕B₁轴旋转角度R1优化空间为第Count次对绕B₂轴旋转角度R2优化结果为第Count-1次对绕B₂轴旋转角度R2优化结果为第Count次对绕B₂轴旋转角度R2优化空间为优化步骤包括：E01令优化次数Count＝0，第Count-1次对绕B₁轴旋转角度R1优化结果第Count-1次对绕B₂轴旋转角度R2优化结果R2optCount-1=0]]>E02令绕B₂轴旋转角度R2的值为在绕B₁轴旋转角度R1的第Count次优化空间内进行间隔为等间隔采样，其中为(R1optCount-1-180*2CountNCount,R1optCount-1+180*2CountNCount),]]>N是采样数量其取值范围是15≤N≤25，N∈整数；E03在绕B₁轴旋转角度R1的每个采样值上，对绕B₂轴旋转角度R2的第Count次优化空间为(R2optCount-1-180*2CountNCount,R2optCount-1+180*2CountNCount),]]>进行间隔为的等间隔采样，其中N是采样数量其取值范围是15≤N≤25，N∈整数；E04在每一个采样点上生成仿真投影DRR，记为I₂(x，y)，其中x是图像的横坐标y是图像的纵坐标且满足0≤x≤P-1；0≤y≤P-1；x，y∈整数，P是二维图像的单边长度；E05计算每一个DRR和真实投影间的相似度；相似度计算包含以下步骤：E0501应用Sobel滤波器对二维真实投影图像I₁(x，y)和二维仿真投影图像I₂(x，y)分别进行滤波得到二维真实投影图像的增强图像J₁(x，y)和二维仿真投影图像的增强图像J₂(x，y)。E0502对二维真实投影图像的增强图像J₁(x，y)和二维仿真投影图像的增强图像J₂(x，y)进行DFT变换得到二维真实投影图像的频谱F₁(u，v)和二维仿真投影图像的频谱F₂(u，v)，其中F1(u,v)=1P2Σx=0P-1Σy=0P-1J1(x,y)e-j2πux+vyP,]]>F2(u,v)=1P2Σx=0P-1Σy=0P-1J2(x,y)e-j2πux+vyP,]]>P为二维图像的单边长度，x，y分别为二维真实投影图像的增强图像J₁(x，y)和二维仿真投影图像的增强图像J₂(x，y)的横、纵座标，u，v分别为二维真实投影图像的频谱F₁(u，v)和二维仿真投影图像的频谱F₂(u，v)的横、纵坐标；E0503计算二维真实投影图像的频谱F₁(u，v)的幅度谱|F₁(u，v)|和二维仿真投影图像的频谱F₂(u，v)的幅度谱|F₂(u，v)|，其中分别代表二维真实投影图像的频谱F₁(u，v)的实部和虚部，分别代表二维仿真投影图像的频谱F₂(u，v)的实部和虚部；E0504对二维真实投影图像频谱的幅度谱|F₁(u，v)|和二维仿真投影图像频谱的幅度谱|F₂(u，v)|进行梅林变换得到二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)，其中ρ=0.5log(u2+v2)·Plog(0.5·P),]]>θ=(π+arctan(vu))·P2π;]]>E0505对二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)再次进行DFT变换得到二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)；其中G1(u,v)=1P2Σρ=0P-1Σθ=0P-1FM1(ρ,θ)e-j2πuρ+vθP,]]>G2(u,v)=1P2Σρ=0P-1Σθ=0P-1FM2(ρ,θ)e-j2πuρ+vθP,]]>P为二维图像的单边长度，ρ，θ分别为二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)的横、纵座标，u，v分别为二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的横、纵座标；E0506计算二维真实投影图像的傅立叶梅林频谱G₁(u，v)的幅度谱|G₁(u，v)|和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的幅度谱|G₂(u，v)|，其中分别代表二维真实投影图像的傅立叶梅林频谱G₁(u，v)的实部和虚部，分别代表二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的实部和虚部；E0507计算二维真实投影图像的傅立叶梅林频谱G₁(u，v)的幅度谱|G₁(u，v)|和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的幅度谱|G₂(u，v)|的相似度Sim，其中Sim=Σu=0P-1Σv=0P-1|G1(u,v)|·|G2(u,v)|(Σu=0P-1Σv=0P-1G1(u,v)2)·(Σu=0P-1Σv=0P-1G2(u,v)2)]]>P为二维图像的单边长度；E06在对所有的采样点进行相似度计算后得到优化空间，在优化空间中寻找全局最大值，其对应的坐标分别是第Count次优化后绕B₁轴旋转角度最优值和第Count次优化后绕B₂轴旋转角度最优值其中参数Count是优化的次数；E07令Count＝Count+1重复步骤E03-E05直到并且其中ε_R1是绕B₁轴旋转角度R1要求达到的精确度，ε_R2是绕B₂轴旋转角度R2要求达到的精确度，记绕B₁轴旋转角度最终优化值为R1_opt，绕B₂轴旋转角度最终优化值为R2_opt；F：计算绕B₃轴旋转角度R3，包括以下步骤：F01以最新估计的参数为基准，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝0，虚拟X射线源到投影模型坐标系坐标原点距离D1＝2*D，投影模型坐标系坐标原点到虚拟投影平面距离D2＝0.5*D，沿B₁轴平移T1＝0，沿B₂轴平移T2＝0生成仿真投影DRR，其中D是三维图像对角线长度，D=Width2+Height2+Depth2,]]>Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，R1_opt是绕B₁轴旋转角度R1的最优值，R2_opt是绕B₂轴旋转角度R2的最优值；F02重复步骤E0501-E0505，得到二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)；F03对二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)进行相位相关运算得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)，其中其中是二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的共轭；F04对二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)进行IDFT变换得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)，其中r(x,y)=Σu=0P-1Σv=0P-1R(u,v)ej2πux+vyP,]]>P是二维图像的单边长度，u，v是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)的横、纵坐标值，x，y是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)的横、纵坐标值；F05在二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)上查找最大值点，其对应的纵轴上的坐标值乘以系数得到绕B₃旋转角度R3的最优值R3_opt；G：对虚拟X射线源到投影模型坐标系坐标原点D1进行优化，包括以下步骤：G01令优化次数Count＝0，第Count-1次优化的虚拟X射线源到投影模型坐标系坐标原点距离D1的最优值G02对虚拟X射线源到投影模型坐标系坐标原点D1的优化空间为(D1optCount-1-D*2CountNCount,D1optCount-1+D*2CountNCount),]]>进行间隔为的采样，其中D是三维图像对角线长度，D=Width2+Height2+Depth2,]]>Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N是采样数量其取值范围是15≤N≤25，N∈整数；G03以最新优化的参数为基准，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝R3_opt投影模型坐标系坐标原点到虚拟投影平面距离D2＝0.5*D，沿B₁轴平移T1＝0，沿B₂轴平移T2＝0，对所有采样点生成仿真投影DRR；G04对每一幅仿真投影DRR重复步骤E0501-E0507计算DRR与真实投影图像的相似度，得到相似度分布；G05在相似度分布中找到最大相似度并记相应的虚拟X射线源到投影坐标系坐标原点距离D1值为第Count次优化的虚拟X射线源到投影模型坐标系坐标原点距离D1的最优值G06令Count＝Count+1重复G02-G05直到Count≥2，记优化后的虚拟X射线源到投影模型坐标系坐标原点距离D1的最优值为D1_opt；H：对投影模型坐标系坐标原点到虚拟投影平面距离D2进行优化，其步骤包括：H01首先以最新估计的参数为基准，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝R3_opt，虚拟X射线源到投影模型坐标系坐标原点距离D1＝D1_opt，投影模型坐标系坐标原点到虚拟投影平面距离D2＝0.5*D，沿B₁轴平移T1＝0，沿B₂轴平移T2＝0生成仿真投影DRR，其中D是三维图像对角线长度，D=Width2+Height2+Depth2,]]>Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽；H02重复步骤E0501-E0505，得到二维真实投影图像傅立叶梅林频谱G₁(u，v)和二维仿真投影图像傅立叶梅林频谱G₂(u，v)；H03将二维真实投影图像傅立叶梅林频谱G₁(u，v)和二维仿真投影图像傅立叶梅林频谱G₂(u，v)进行相位相关运算得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)，其中是二维仿真投影图像傅立叶梅林频谱G₂(u，v)的共轭；H04对二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)进行IDFT变换得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)，其中r(x,y)=Σu=0P-1Σv=0P-1R(u,v)ej2πux+vyP,]]>P是二维图像的单边长度，u，v是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)的横、纵坐标值，x，y是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)的横、纵坐标值；H05在二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)上查找最大值点，找到全局最大值，其对应的横轴上的值乘以系数得到尺度参数a，a代表着图像的缩放比例，通过尺度参数a可以计算出投影模型坐标系坐标原点到虚拟投影平面之间的近似距离D2_Est，D2_Est＝a*(D1_opt+D2_Curr)-D1_opt，式中D1_opt是虚拟X射线源到投影模型坐标系坐标原点距离D1的最优值，D2_Curr是投影模型坐标系坐标原点到虚拟投影平面的初始距离，其值为0.5*D，D是三维图像对角线长度，D=Width2+Height2+Depth2,]]>Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N是采样数量其取值范围是15≤N≤25，N∈整数，D2_Est是投影模型坐标系坐标原点到虚拟投影平面距离D2的修正值；H06令优化次数Count＝0，第Count-1次优化投影模型坐标系坐标原点到虚拟投影平面距离D2的初始值H07在参数优化范围(D2optCount-1-0.2*D*2CountNCount,D2optCount-1+0.2*D*2CountNCount)]]>内对投影模型坐标系坐标原点到虚拟投影平面距离D2进行间隔为进行等间隔采样，其中D是三维图像对角线长度，D=Width2+Height2+Depth2,]]>Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N是采样数量其取值范围是15≤N≤25，N∈整数；H08对每一个采样点生成仿真投影DRR；H09对每一个仿真投影DRR求其与真实投影图像缩放比例的一致性，步骤包括：H0901重复步骤E0501-E0502得到二维真实投影图像的频谱F₁(u，v)和二维仿真投影图像的频谱F₂(u，v)；H0902对二维真实投影图像的频谱F₁(u，v)和二维仿真投影图像的频谱F₂(u，v)进行相位相关运算得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)，其中Q(u,v)=F1(u,v)F2*(u,v)|F1(u,v)F2*(u,v)|,]]>是二维仿真投影图像的频谱F₂(u，v)的共轭；H0903对二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)进行IDFT变换得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)，其中q(x,y)=Σu=0P-1Σv=0P-1R(u,v)ej2πux+vyP,]]>P是二维图像的单边长度，x，y分别为二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)的横、纵座标，u，v分别为二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)的横、纵座标；H0904求二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)的直方图histo(r_k)＝n_k，其中r_k是第k级灰度，n_k是图像中灰度级为r_k的像素个数，对直方图进行能量归一化得到能量归一化直方图r_K＝0，1，...MaxVal，其中MaxVal表示二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)中最大灰度值，n表示二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)中像素总数，根据能量归一化直方图求二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)的信息熵E，其中E=Σrk=0MaxVal-histoN(rk)*log(histoN(rk)),]]>其中MaxVal是二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)中最大灰度值，E代表着真实投影图像和仿真投影图像的缩放比例一致性；H10求优化空间中最小值，其对应的投影坐标系坐标原点到虚拟投影平面距离D2的值记为第Count次优化的投影坐标系坐标原点到虚拟投影平面距离D2的最优值H11令优化次数Count＝Count+1，重复步骤H07-H10，直到Count≥2记优化后参数值为D2_opt；I：计算沿B₁轴平移T1和沿B₂轴平移T2，步骤有：I01根据当前的优化参数值，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝R3_opt X射线源到投影模型坐标系坐标原点距离D1＝D1_opt，投影模型坐标系坐标原点到虚拟投影平面距离D2＝D2_opt，沿B₁轴平移T1＝0，沿B₂轴平移T2＝0生成仿真投影DRR；I02重复步骤H0901-H0903得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)I05查找图像q(x，y)中的最大值，其对应的横坐标记为沿B₁轴平移T1的最优值T1_opt，纵坐标记为沿B₂轴平移T2的最优值T2_optJ：更新初始化参数，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝R3_opt X射线源到投影模型坐标系坐标原点距离D1＝D1_opt，投影模型坐标系坐标原点到虚拟投影平面距离D2＝D2_opt，沿B₁轴平移T1＝T1_opt，沿B₂轴平移T2＝T2_opt，重复步骤D-H，使得参数变化小于阈值，包括|ΔR1|＜1°，|ΔR2|＜1°，|ΔR3|＜1°，|ΔT1|＜2，|ΔT2|＜2，其中ΔR1是两次优化间绕B₁轴旋转角度R1的变化量，ΔR2是两次优化间绕B₂轴旋转角度R2的变化量，ΔR3是两次优化间绕B₃轴旋转角度R3的变化量，ΔT1是两次优化间沿B₁平移T1的变化量，ΔT2是两次优化间沿B₂平移T2的变化量。