CN102609979A - 一种基于傅里叶梅林域的二维三维图像配准方法 - Google Patents

Abstract

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

Description

一种基于傅里叶梅林域的二维三维图像配准方法

技术领域

本发明属于医学图像配准领域，特别是涉及二维三维图像配准。 

背景技术

计算机辅助的导航手术是一个集医学和计算机技术等诸多学科为一体的新型交叉研究领域。其目的是借助计算机以及跟踪设备模拟手术中涉及的各个过程包括手术规划，手术导航等来指导医生实现高精度的外科手术。 

近年来，计算机辅助手术规划逐渐替代了传统的基于X光片的手术规划。目前市场上流行的手术规划软件有很多如西门子公司的PreOPlan，Stryker公司的CasePlan，iCRco公司的Orthocrat等。这些软件借助自带的标定功能可以识别二维投影图像中像素的大小。利用人机交互可以测量骨骼长度，选择合适的植入器械，定义骨骼修正后的位置以及植入器械的固定方法。在创伤手术的流程中，术前计划一般都是基于二维的投影图像。因为三维扫描过程比较长，需要病人的移动较多，常常会给创伤部位带来不必要的附加伤害。然而，术中采集的数据因为要作为手术导航中的基础数据通常都是三维的。这就决定了要用二维三维配准的方法来匹配术前投影数据和术中采集的三维数据。目前在手术导航系统中，术前规划还只是以示意图的形式作为一种指导手术的途径。 

随着手术导航系统在临床中的广泛应用，二维三维图像配准已经成为系统开发中的一个瓶颈。目前主流的算法是通过对仿真投影DRR和投影图像相似度的比较，优化三维图像和二维投影之间的空间位置关系。由于该算法是一个在高维空间上的优化过程。因此它的复杂度高，并且鲁棒性低。 

发明内容

本发明基于傅里叶梅林变换提出了二维三维配准方法。该方法包含如下步骤： 

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 Heightx 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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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优化结果 

E02令绕B₂轴旋转角度R2的值为 

在绕B₁轴旋转角度R1的第Count次优化空间 

内进行间隔为 

等间隔采样，其中 

为 $({R1}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{R1}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ N是采样数量其取值范围是15≤N≤25，N∈整数； 

E03在绕B₁轴旋转角度R1的每个采样值上，对绕B₂轴旋转角度R2的第Count次优化空间 为  $({R2}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{R2}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ 进行间隔为 

的等间隔采样，其中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)，其中  $F_1(u,v)=\frac{1}{P^2}\underset{x=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{P-1}{\mathrm{\Σ}}}{J}_1{(x,y)e}^{-j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ $F_2(u,v)=\frac{1}{P^2}\underset{x=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{P-1}{\mathrm{\Σ}}}{J}_2{(x,y)e}^{-j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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₂(ρ，θ)，其中 $\mathrm{\ρ}=0.5\log (u^2+v^2)\·\frac{P}{\log (0.5\·P)},\mathrm{\θ}=\frac{(\mathrm{\π}+\arctan \left(\frac{v}{u}\right))\·P}{2\mathrm{\π}};$

E0505对二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)再次进行DFT变换得到二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)；其中 $G_1(u,v)=\frac{1}{P^2}\underset{\mathrm{\ρ}=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{\mathrm{\θ}=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{FM}}_1{(\mathrm{\ρ},\mathrm{\θ})e}^{-j2\mathrm{\π}\frac{\mathrm{u\ρ}+\mathrm{v\θ}}{P}},$ $G_2(u,v)=\frac{1}{P^2}\underset{\mathrm{\ρ}=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{\mathrm{\θ}=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{FM}}_2{(\mathrm{\ρ},\mathrm{\θ})e}^{-j2\mathrm{\π}\frac{\mathrm{u\ρ}+\mathrm{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，其中  $\mathrm{Sim}=\frac{\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}\left|G_1(u,v)\right|\·\left|G_2(u,v)\right|}{\sqrt{\left(\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{G}_1{(u,v)}^2\right)\·\left(\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{G}_2{(u,v)}^2\right)}}$ 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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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)=\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{R(u,v)e}^{j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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的优化空间 

为 $({D1}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{D1}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ 进行间隔为 

的 采样，其中D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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)=\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{R(u,v)e}^{j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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的最优值，D2C_urr是投影模型坐标系坐标原点到虚拟投影平面的初始距离，其值为0.5*D，D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N是采样数量其取值范围是15≤N≤25，N∈整数，D2_Est是投影模型坐标系坐标原点到虚拟投影平面距离D2的修正值； 

H06令优化次数Count＝0，第Count-1次优化投影模型坐标系坐标原点到虚拟投影平面距离D2的初始值 

H07在参数优化范围 $({D2}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{0.2*D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{D2}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{0.2*D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}})$ 内对 投影模型坐标系坐标原点到虚拟投影平面距离D2进行间隔为 进行等间隔采样，其中D是三维图像对角线长度，  $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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)=\frac{F_1(u,v)F_2^{*}(u,v)}{\left|F_1(u,v)F_2^{*}(u,v)\right|},$

是二维仿真投影图像的频谱F₂(u，v)的共轭； 

H0903对二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)进行IDFT变换得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)，其中 $q(x,y)=\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{R(u,v)e}^{j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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=\underset{r_k=0}{\overset{\mathrm{MaxVal}}{\mathrm{\Σ}}}-\mathrm{histoN}\left(r_k\right)*\log \left(\mathrm{histoN}\left(r_k\right)\right),$ 其中MaxVal是二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)中最大灰度值，E代表着真实投影图像和仿真投影图像的缩放比例一致性； 

H10求优化空间中最小值，其对应的投影坐标系坐标原点到虚拟投影平面距离D2的值记为第Count次优化的投影坐标系坐标原点到虚拟投影平面距离D2的最优值 

H11令优化次数Count＝Count+1，重复步骤H07-H10，直到Count≥2记优化后参数值为D2_opt； 

I：计算平面内平移，步骤有： 

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_opt

J：更新初始化参数，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝R3_opt X射线源到投影模型坐标系坐标原点距离D1＝D1_opt，投影模型坐标系坐标原点到虚拟投影平面距离D2＝D2_opt，沿B1轴平移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的变化量。 

本发明利用傅里叶梅林变换的尺度，旋转，平移不变性将高维空间中的优化过程转化为一系列的低维空间中的优化过程。增强了配准的鲁棒性。 

附图说明

图1二维三维图像配准系统框图； 

图2二维三维图像配准中应用的投影模型； 

图3二维三维图像配准的流程图； 

图4投影图像与仿真投影图像相似度的比较流程； 

图5Sobel滤波器模板，A模板坐标定义；B为横向模板；C为纵向模板； 

图6平面外旋转参数Ru，Rv的优化流程； 

图7X射线源-三维图像距离D1的优化流程； 

图8三维图像-接收器距离D2的优化流程； 

图9根据当前接收器位置估算实际接收器位置示意图； 

图10配准精确度，A配准后的平移误差；B配准后的旋转误差。 

具体实施方式

本发明涉及一种基于傅里叶梅林域的二维三维图像配准方法。用于术前规划的二维投影图像可以由X光设备来采集。这类设备的一个实例是Siemens Mobilet。用于手术导航的三维数据可以由平板C形臂采集。这类设备的一个实例是Siemens Artis Zeego。 

如图1所示，采集的三维图像和二维投影图像首先被传送到DICOM服务器上，在配准算法运行之前再将数据由DICOM服务器传输到配准PC上。配准PC得到的三维图像记为V，其大小为Width x Height x Depth，其中Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，配准PC得到的二维真实投影图像记为I₁(x，y)，其大小为 PxP。 

如图2所示，用于配准的投影模型由三维图像和投影结构两部分组成，其中直角坐标系A₁A₂A₃的原点位于三维图像中心且坐标轴A₁，A₂，A₃分别正交于三维图像相对应的外平面；直角坐标系B₁B₂B₃和直角坐标系A₁A₂A₃重合且B₁，B₂，B₃坐标轴分别与A₁，A₂，A₃坐标轴同向，记B₁B₂B₃坐标原点为ISO，X光源位于B₃轴正半轴上。X光源到ISO中心的距离是D1，投影平面垂直于B₃轴，投影平面和B₃轴的交点位于B₃轴的负半轴，交点位于投影平面的中心，交点和ISO中心的距离是D2。 

B₁B₂B₃在坐标系A₁A₂A₃下的空间变化可以由以下变化组成： 

围绕B₂轴的旋转R2→围绕B₁轴的旋转R1→围绕B₃轴的旋转R3→沿着B₁轴平移T1→沿着B₂轴平移T2。D1的变化影响了仿真投影DRR中的投影失真。D2的变化引起DRR中物体的缩放比例。对于该投影模型，待优化的参数有：绕B₁轴旋转参数R1，绕B₂轴的旋转参数R2，绕B₃轴的旋转参数R3，沿B₁轴的平移T1，沿B₂轴的平移T2，虚拟X射线源到投影模型坐标系坐标原点距离D1，投影模型坐标系坐标原点ISO到虚拟投影平面距离D2。虚拟X射线源到投影模型坐标系坐标原点距离D1的改变只会影响到投影图像中由于投影失真产生的形变。投影模型坐标系坐标原点ISO到虚拟投影平面距离D2的改变只会影响到投影图像中物体的大小。因此这种新的投影模型可以有效地将投影失真和物体大小的估计分离，增加了配准的鲁棒性。 

如图3所示，二维三维图像配准方法包括如下主要步骤： 

1设定参数初始值。令绕B₁轴的旋转角度R1＝0，绕B₂轴的旋转角度R2＝0，绕B₃轴的旋转角度R3＝0；虚拟X射线源到投影模型坐标系坐标原点距离D1＝2*D其中D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽；投影模型坐标系坐标原点到虚拟投影平面距离D2＝0.5*D；沿B₁轴平移T1＝0，沿B₂轴平移T2＝0。 

2对绕B₁轴的旋转角度R1以及绕B₂轴的旋转角度R2进行优化。这两个参数优化的流程由图6示出。对这两个参数的优化过程实际上是一个在二维参数优化空间中找最优参数R1_opt，R2_opt使其对应的仿真投影DRR和真实投影得相似度最大的过程。为了加速优化过程，本发明使用了由粗略到精细的优化结构。首先对优化空间进行粗略采样，在得到粗略 最优值之后在其邻域进行精细采样和优化，直到两次迭代步骤间优化参数改变量小于精度要求。在优化之前，首先对参数进行初始化，包括初始化优化次数变量Count＝0，第Count-1次对绕B₁轴旋转角度R1优化结果 

第Count-1次对绕B₂轴旋转角度R2优化结果 

优化过程包括如下步骤：第一步，在绕B₁轴旋转角度R1的第Count次优化空间 

内进行间隔为 

等间隔采样，其中 

为  $({R1}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{R1}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ N＝20是采样数量。第二步，在绕B₁轴旋转角度R1的每个采样值上，对绕B₂轴旋转角度R2的第Count次优化空间 

为  $({R2}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{R2}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ 进行间隔为 

的等间隔采样，其中N＝20是采样数量。在完成对R2的优化范围采样后，即完成了对R1，R2优化空间的采样。第三步，在每一个采样点上生成仿真投影DRR，记为I₂(x，y)，其中x是图像的横坐标y是图像的纵坐标且满足0≤x≤P-1；0≤y≤P-1；x，y∈整数，P是二维图像的单边长度。如图6所示，为了对采样空间中的采样点进行遍历移入了临时变量i和j。第四步，计算每一个DRR和真实投影间的相似度。如图4所示，相似度的计算可以由如下步骤完成： 

201应用Sobel滤波器对二维真实投影图像I₁(x，y)和二维仿真投影图像I₂(x，y)分别进行滤波。首先用Sobel算子的横向模板B和纵向模板C(参见图5)分别对二维真实投影图像I₁(x，y)和二维仿真投影图像I₂(x，y)进行空间滤波得到二维真实投影图像的水平增强图像I_1B(x，y)，二维真实投影图像的竖直增强图像I_1C(x，y)，二维仿真投影图像的水平增强图像I_2B(x，y)和二维仿真投影图像的竖直增强图像I_2C(x，y)，参照图5A中模板脚标的定义，滤波的过程可以表示为： 

$I_{1B}(x,y)=\underset{m=-1}{\overset{1}{\mathrm{\Σ}}}\underset{n=-1}{\overset{1}{\mathrm{\Σ}}}B(m,n)I_1(x+m,y+n),$ $I_{1C}(x,y)=\underset{m=-1}{\overset{1}{\mathrm{\Σ}}}\underset{n=-1}{\overset{1}{\mathrm{\Σ}}}C(m,n)I_1(x+m,y+n),$

$I_{2B}(x,y)=\underset{m=-1}{\overset{1}{\mathrm{\Σ}}}\underset{n=-1}{\overset{1}{\mathrm{\Σ}}}B(m,n)I_2(x+m,y+n),$ $I_{2C}(x,y)=\underset{m=-1}{\overset{1}{\mathrm{\Σ}}}\underset{n=-1}{\overset{1}{\mathrm{\Σ}}}C(m,n)I_2(x+m,y+n)$

计算二维真实投影图像的增强图像J₁(x，y)和二维仿真投影图像的增强图像J₂(x，y)，其中 $J_1(x,y)=\sqrt{I_{1B}{(x,y)}^2+I_{1C}{(x,y)}^2},$ $J_2(x,y)=\sqrt{I_{2B}{(x,y)}^2+I_{2C}{(x,y)}^2};$

202对二维真实投影图像的增强图像J₁(x，y)和二维仿真投影图像的增强图像 J₂(x，y)进行DFT变换得到F₁(u，v)和F₂(u，v)，其中 $F_1(u,v)=\frac{1}{P^2}\underset{x=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{P-1}{\mathrm{\Σ}}}{J}_1{(x,y)e}^{-j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ $F_2(u,v)=\frac{1}{P^2}\underset{x=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{P-1}{\mathrm{\Σ}}}{J}_2{(x,y)e}^{-j2\mathrm{\π}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ P为二维图像的单边长度，x，y分别为二维真实投影图像的增强图像J₁(x，y)和二维仿真投影图像的增强图像J₂(x，y)的横、纵坐标，u，v分别为二维真实投影图像频谱F₁(u，v)和二维仿真投影图像频谱F₂(u，v)的横、纵坐标； 

203计算二维真实投影图像频谱F₁(u，v)的幅度谱|F₁(u，v)|和二维仿真投影图像频谱F₂(u，v)的幅度谱|F₂(u，v)|，其中 

分别代表二维真实投影图像的频谱F₁(u，v)的实部和虚部， 

分别代表二维仿真投影图像的频谱F₂(u，v)的实部和虚部； 

204对二维真实投影图像频谱的幅度谱|F₁(u，v)|和二维仿真投影图像频谱的幅度谱|F₂(u，v)|进行梅林变换得到二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)，其中  $\mathrm{\ρ}=0.5\log (u^2+v^2)\·\frac{P}{\log (0.5\·P)},$ $\mathrm{\θ}=\frac{(\mathrm{\π}+\arctan \left(\frac{v}{u}\right))\·P}{2\mathrm{\π}},$ P为二维图像的单边长度。 

205对二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)再次进行DFT变换得到二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)；其中  $G_1(u,v)=\frac{1}{P^2}\underset{\mathrm{\ρ}=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{\mathrm{\θ}=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{FM}}_1{(\mathrm{\ρ},\mathrm{\θ})e}^{-j2\mathrm{\π}\frac{\mathrm{u\ρ}+\mathrm{v\θ}}{P}},$ $G_2(u,v)=\frac{1}{P^2}\underset{\mathrm{\ρ}=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{\mathrm{\θ}=0}{\overset{P-1}{\mathrm{\Σ}}}{\mathrm{FM}}_2{(\mathrm{\ρ},\mathrm{\θ})e}^{-j2\mathrm{\π}\frac{\mathrm{u\ρ}+\mathrm{v\θ}}{P}},$ P为二维图像的单边长度，ρ，θ分别为二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)的横、纵座标，u，v分别为二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的横、纵座标。 

206计算二维真实投影图像的傅立叶梅林频谱G₁(u，v)的幅度谱|G₁(u，v)|和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的幅度谱|G₂(u，v)|，其中 

分别代表二维真实投影图像的傅立叶梅林频谱G₁(u，v)的实部和虚部， 

分别代表二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的实部和虚部。 

207计算二维真实投影图像的傅立叶梅林频谱G₁(u，v)的幅度谱|G₁(u，v)|和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的幅度谱|G₂(u，v)|的相似度Sim，其中  $\mathrm{Sim}=\frac{\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}\left|G_1(u,v)\right|\·\left|G_2(u,v)\right|}{\sqrt{\left(\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{G}_1{(u,v)}^2\right)\·\left(\underset{u=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{P-1}{\mathrm{\Σ}}}{G}_2{(u,v)}^2\right)}}$ P为二维图像的单边长度。 

 第五步，在对所有的采样点进行相似度计算后得到优化空间，在优化空间中寻找全局最大值，其对应的坐标分别是第Count次优化后绕B₁轴旋转角度最优值 

和第Count次优化后绕B₂轴旋转角度最优值 

其中参数Count是优化的次数。令Count＝Count+1重复第一步到第五步，直到 $|{R1}_{\mathrm{opt}}^{\mathrm{Count}}-{R1}_{\mathrm{opt}}^{\mathrm{Count}-1}|<{\mathrm{\&epsiv;}}_{R1}$ 并且 $|{R2}_{\mathrm{opt}}^{\mathrm{Count}}-{R2}_{\mathrm{opt}}^{\mathrm{Count}-1}|<{\mathrm{\&epsiv;}}_{R2},$ 其中ε_R1是绕B₁轴旋转角度R1要求达到的精确度，ε_R2是绕B₂轴旋转角度R2要求达到的精确度，记绕B₁轴旋转角度最终优化值为R1_opt，绕B₂轴旋转角度最终优化值为R2_opt。 

3以最新估计的参数为基准，包括绕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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，R1_opt是绕B₁轴旋转角度R1的最优值，R2_opt是绕B₂轴旋转角度R2的最优值。重复步骤201-205，得到二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)。对二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)进行相位相关运算得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)，其中  $R(u,v)=\frac{G_1(u,v)G_2^{*}(u,v)}{\left|G_1(u,v)G_2^{*}(u,v)\right|},$ 其中 

是二维仿真投影图像的傅立叶梅林频谱G₂(u，v)的共轭。对二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相 位相关谱R(u，v)进行IDFT变换得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)，其中 $r(x,y)=\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{R(u,v)e}^{j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ P是二维图像的单边长度，u，v是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)的横、纵坐标值，x，y是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)的横、纵坐标值。在二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)上查找最大值点，其对应的纵轴上的坐标值乘以系数 

得到绕B₃旋转角度R3的最优值R3_opt。 

4对虚拟X射线源到投影模型坐标系坐标原点距离D1进行优化。按照由粗略到精细的优化方法对该参数进行优化。如图7所示，首先初始化参数，令优化次数CounT＝0，第Count-1次优化的虚拟X射线源到投影模型坐标系坐标原点距离D1的最优值 优化过程包括如下步骤：第一步，对虚拟X射线源到投影模型坐标系坐标原点D1的优化空间 

为 $({D1}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{D1}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ 进行间隔为 

的采样，其中D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N＝20是采样数量。第二步，以最新优化的参数为基准，包括绕B₁轴旋转角度R1＝R1_opt，绕B₂轴旋转角度R2＝R2_opt，绕B₃轴旋转角度R3＝R3_opt投影模型坐标系坐标原点到虚拟投影平面距离D2＝0.5*D，沿B₁轴平移T1＝0，沿B₂轴平移T2＝0，对所有采样点生成仿真投影DRR。如图7所示，为了对采样空间中每一个采样点进行遍历引入了临时变量i。第三步，对每一幅仿真投影DRR重复步骤201-207计算DRR与真实投影图像的相似度，得到相似度分布。第四步，在相似度分布中找到最大相似度并记相应的虚拟X射线源到投影坐标系坐标原点距离D1值为第Count次优化的虚拟X射线源到投影模型坐标系坐标原点距离D1的最优值 

第五步，令Count＝Count+1重复步骤一至步骤五直到Count≥2，记优化后的参数值为D1_opt。 

5对三维图像-接收器距离D2进行优化。如图8所示，首先根据仿真投影图像和真实 投影图像之间的尺度差粗略估计三维图像-接收器距离D2，其步骤包括：第一步，首先以最新估计的参数为基准，包括绕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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽。第二步，重复步骤201-205，得到二维真实投影图像傅立叶梅林频谱G₁(u，v)和二维仿真投影图像傅立叶梅林频谱G₂(u，v)。第三步，将二维真实投影图像傅立叶梅林频谱G₁(u，v)和二维仿真投影图像傅立叶梅林频谱G₂(u，v)进行相位相关运算得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)，其中 $R(u,v)=\frac{G_1(u,v)G_2^{*}(u,v)}{\left|G_1(u,v)G_2^{*}(u,v)\right|},$

是二维仿真投影图像傅立叶梅林频谱G₂(u，v)的共轭。第四步，对二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)进行IDFT变换得到二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)，其中  $r(x,y)=\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{R(u,v)e}^{j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ P是二维图像的单边长度，u，v是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关谱R(u，v)的横、纵坐标值，x，y是二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)的横、纵坐标值。第五步，在二维真实投影图像的傅立叶梅林频谱和二维仿真投影图像的傅立叶梅林频谱的相位相关图像r(x，y)上查找最大值点，找到全局最大值，其对应的横轴上的值乘以系数 

得到尺度参数，通过a可以计算出三维图像-接收器之间的近似距离D2_Est，求解步骤如下：如图9所示，接收器真实位置上的投影和接收器当前位置上的投影存在缩放比例 

其中l是二维真实投影图像中物体的大小，l′是二维仿真投影图像中物体的大小。可以推导出 

整理后得到D2_Est＝a*(D1_opt+D2_Curr)-D1_opt，式中D1_opt是虚拟X射线源到投影模型坐标系坐标原点 距离D1的最优值，D2_Curr是投影模型坐标系坐标原点到虚拟投影平面的初始距离，其值为0.5*D，D是三维图像对角线长度，D2_Est是投影模型坐标系坐标原点到虚拟投影平面距离D2的修正值。 

在得到三维图像-接收器距离D2的修正值D2_Est后，按照由粗略到精细的优化方法对该参数进行精细优化。首先初始化优化参数，令优化次数Count＝0，第Count-1次优化投影模型坐标系坐标原点到虚拟投影平面距离D2的初始值 

对该参数的优化过程包含如下步骤：第一步，在参数优化范围  $({D2}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{0.2*D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{D2}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{0.2*D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}})$ 内对投影模型坐标系坐标原点到虚拟投影平面距离D2进行间隔为 

进行等间隔采样，其中D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N是采样数量其取值范围是15≤N≤25，N∈整数。第二步，对每一个采样点生成仿真投影DRR。如图8所示，为了遍历采样空间中的采样点引入临时变量i。第三步，对每一个仿真投影DRR求其与真实投影缩放比例的一致性，步骤包括： 

501重复步骤201-202得到二维真实投影图像的频谱F₁(u，v)和二维仿真投影图像的频谱F₂(u，v)。 

502对二维真实投影图像的频谱F₁(u，v)和二维仿真投影图像的频谱F₂(u，v)进行相位相关运算得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)，其中 $Q(u,v)=\frac{F_1(u,v)F_2^{*}(u,v)}{\left|F_1(u,v)F_2^{*}(u,v)\right|},$

是二维仿真投影图像的频谱F₂(u，v)的共轭。 

503对二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)进行IDFT变换得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)，其中 $q(x,y)=\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{R(u,v)e}^{j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ P是二维图像的单边长度，x，y分别为二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)的横、纵座标，u，v分别为二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)的横、纵座标。 

504求二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像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=\underset{r_k=0}{\overset{\mathrm{MaxVal}}{\mathrm{\&Sigma;}}}-\mathrm{histoN}\left(r_k\right)*\log \left(\mathrm{histoN}\left(r_k\right)\right),$ 其中MaxVal是二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)中最大灰度值，E代表真实投影图像和仿真投影图像的缩放比例一致性。 

第四步，求优化空间中最小值，其对应的投影坐标系坐标原点到虚拟投影平面距离D2的值记为第Count次优化的投影坐标系坐标原点到虚拟投影平面距离D2的最优值 

第五步，令优化次数Count＝Count+1，重复第一步-第四步，直到Count≥2记优化后参数值为D2_opt； 

6计算沿B₁轴平移T1和沿B₂轴平移T2。首先根据当前的优化参数值，包括绕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。重复步骤501-503得到得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)。查找图像q(x，y)中的最大值，其对应的横坐标记为沿B₁轴平移T1的最优值T1_opt，纵坐标记为沿B₂轴平移T2的最优值T2_opt

7更新初始化参数，包括绕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的变化量。 

实验结果 

在对20幅二维投影和一幅三维图像的配准结果显示，本方法具有较高的鲁棒性和较高的配准精度。如图10-A所示，本方法配准后的平均平移误差随着优化迭代次数的增加而减少，在三次优化之后平均平移误差为1.42毫米。如图10-B所示，本方法配准后的平均旋转误差随着优化迭代次数的增加而减少，在三次优化之后平均旋转误差为0.63度。 

Claims (1)

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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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优化结果 ${R2}_{\mathrm{opt}}^{\mathrm{Count}-1}=0$

E02令绕B₂轴旋转角度R2的值为

在绕B₁轴旋转角度R1的第Count次优化空间

内进行间隔为等间隔采样，其中

为 $({R1}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{R1}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ N是采样数量其取值范围是15≤N≤25，N∈整数；

E03在绕B₁轴旋转角度R1的每个采样值上，对绕B₂轴旋转角度R2的第Count次优化空间为 $({R2}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{R2}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{180*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ 进行间隔为的等间隔采样，其中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)，其中

$F_1(u,v)=\frac{1}{P^2}\underset{x=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{J}_1{(x,y)e}^{-j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$

$F_2(u,v)=\frac{1}{P^2}\underset{x=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{J}_2{(x,y)e}^{-j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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₂(ρ，θ)，其中 $\mathrm{\&rho;}=0.5\log (u^2+v^2)\&CenterDot;\frac{P}{\log (0.5\&CenterDot;P)},$ $\mathrm{\&theta;}=\frac{(\mathrm{\&pi;}+\arctan \left(\frac{v}{u}\right))\&CenterDot;P}{2\mathrm{\&pi;}};$

E0505对二维真实投影图像的傅立叶梅林图像FM₁(ρ，θ)和二维仿真投影图像的傅立叶梅林图像FM₂(ρ，θ)再次进行DFT变换得到二维真实投影图像的傅立叶梅林频谱G₁(u，v)和二维仿真投影图像的傅立叶梅林频谱G₂(u，v)；其中 $G_1(u,v)=\frac{1}{P^2}\underset{\mathrm{\&rho;}=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&theta;}=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{\mathrm{FM}}_1{(\mathrm{\&rho;},\mathrm{\&theta;})e}^{-j2\mathrm{\&pi;}\frac{\mathrm{u\&rho;}+\mathrm{v\&theta;}}{P}},$ $G_2(u,v)=\frac{1}{P^2}\underset{\mathrm{\&rho;}=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{\mathrm{\&theta;}=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{\mathrm{FM}}_2{(\mathrm{\&rho;},\mathrm{\&theta;})e}^{-j2\mathrm{\&pi;}\frac{\mathrm{u\&rho;}+\mathrm{v\&theta;}}{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，其中 $\mathrm{Sim}=\frac{\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\left|G_1(u,v)\right|\&CenterDot;\left|G_2(u,v)\right|}{\sqrt{\left(\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{G}_1{(u,v)}^2\right)\&CenterDot;\left(\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{G}_2{(u,v)}^2\right)}}$ 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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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)=\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{R(u,v)e}^{j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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的优化空间为 $({D1}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{D1}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}}),$ 进行间隔为

的采样，其中D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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)=\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{R(u,v)e}^{j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ Width是三维图像的长，Height是三维图像的高，Depth是三维图像的宽，其中N是采样数量其取值范围是15≤N≤25，N∈整数，D2_Est是投影模型坐标系坐标原点到虚拟投影平面距离D2的修正值；

H06令优化次数Count＝0，第Count-1次优化投影模型坐标系坐标原点到虚拟投影平面距离D2的初始值

H07在参数优化范围 $({D2}_{\mathrm{opt}}^{\mathrm{Count}-1}-\frac{{0.2*D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}},{D2}_{\mathrm{opt}}^{\mathrm{Count}-1}+\frac{{0.2*D*2}^{\mathrm{Count}}}{N^{\mathrm{Count}}})$ 内对投影模型坐标系坐标原点到虚拟投影平面距离D2进行间隔为进行等间隔采样，其中D是三维图像对角线长度， $D=\sqrt{{\mathrm{Width}}^2+{\mathrm{Height}}^2+{\mathrm{Depth}}^2},$ 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)=\frac{F_1(u,v)F_2^{*}(u,v)}{\left|F_1(u,v)F_2^{*}(u,v)\right|},$

是二维仿真投影图像的频谱F₂(u，v)的共轭；

H0903对二维真实投影图像的频谱和二维仿真投影图像的频谱的相关谱Q(u，v)进行IDFT变换得到二维真实投影图像的频谱和二维仿真投影图像的频谱的相关图像q(x，y)，其中 $q(x,y)=\underset{u=0}{\overset{P-1}{\mathrm{\&Sigma;}}}\underset{v=0}{\overset{P-1}{\mathrm{\&Sigma;}}}{R(u,v)e}^{j2\mathrm{\&pi;}\frac{\mathrm{ux}+\mathrm{vy}}{P}},$ 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=\underset{r_k=0}{\overset{\mathrm{MaxVal}}{\mathrm{\&Sigma;}}}-\mathrm{histoN}\left(r_k\right)*\log \left(\mathrm{histoN}\left(r_k\right)\right),$ 其中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_opt

J：更新初始化参数，包括绕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的变化量。

Images