CN103006263A - 一种基于线性扫描的医学超声三维成像的位置标定方法 - Google Patents

Abstract

本发明公开了一种基于线性扫描的医学超声三维成像的位置标定方法。使用已知尺寸的长方体作为扫描对象，在假定原始超声图像平面与线性滑轨之间不存在角度变换关系的情况下对二维图像进行三维重建，则可计算出原始超声图像平面与线性滑轨之间实际存在的角度与重建长方体形变量之间的数学关系；实际标定过程中，以线性扫描方式扫描长方体，对超声图像序列进行三维重建，将重建长方体的形变量代入数学关系式计算出原始超声图像平面（即超声探头的位置）与线性滑轨之间的旋转角度，即为线性扫描方式下探头位置标定结果。本发明简单易懂、可操作性强，试验简单、标定方便准确，可用于基于自由臂线性扫描或者机械驱动线性扫描方式的三维超声成像系统。

Description

一种基于线性扫描的医学超声三维成像的位置标定方法

技术领域

本发明属于医学超声波的图像处理和图像显示技术领域，特别涉及一种基于线性扫描的医学超声三维成像的位置标定方法。

背景技术

医学超声诊断是现代影像医学的重要组成部分，在人身体很多组织（比如心脏、肝脏、神经血管等）的病理诊断中有很大的价值。其中三维超声成像技术因其图像显示直观、诊断精确性高等特点在超声诊断中得到越来越多的应用。目前，医学三维超声成像中常用的扫描方法有：机械定位平行扫描、扇形扫描、旋转扫描和自由臂(Free-Hand)扫描四种方式。对原始超声图像序列进行三维重建之前，需要确定超声探头与原始二维超声图像之间的位置关系，才能准确的重建出三维对象，因此必须对超声探头做一个位置标定。在自由臂扫描方法中已经出现比较多的位置标定方法，比如交叉线扫描方法、单个反射面扫描方法等。而对于线性扫描方式的位置标定，目前大多是通过扫描一排均匀排列小球来做位置标定。制作这样的标定仿体较为复杂，因此，设计一种既能够减小医学超声三维成像中线性扫描位置标定的复杂度，同时还能保证三维成像准确性的位置标定方法极具实际意义。

发明内容

本发明的主要目的在于克服现有技术的缺点与不足，提供一种基于线性扫描的医学超声三维成像的位置标定方法，该方法既能够减小医学超声三维成像中线性扫描位置标定的复杂度，同时还能够提高三维成像的准确性。

本发明的目的通过以下的技术方案实现：一种基于线性扫描的医学超声三维成像的位置标定方法，首先为线性扫描用的线性滑轨建立一个线性滑轨坐标系，然后将一个测试用的、已知尺寸大小的长方体置于该线性滑轨坐标系中，并保证线性滑轨和长方体的一条边相互平行放置，以长方体一个顶点为原点建立世界坐标系，确定该世界坐标系下每个顶点的坐标；假定超声探头与线性滑轨之间存在一个旋转角度，将世界坐标系下的每个顶点坐标转换到原始超声图像平面坐标系中；再假定超声探头与线性滑轨之间不存在坐标转换，再将原始超声图像平面坐标系下的顶点坐标转换到线性滑轨坐标系中得到这些顶点新的坐标；根据各个顶点的转换坐标计算出每两条边的夹角与超声探头和线性滑轨之间的旋转角度的关系，得出线性扫描方式下探头位置标定公式；然后线性扫描并测量重建长方体形变量，根据标定公式计算出超声探头与线性滑轨之间的旋转角度，即为位置标定结果。

具体包括以下步骤：

（1）初始化：为线性扫描用的线性滑轨建立一个线性滑轨坐标系，一个位置传感器与线性滑轨固定相连，并可以通过机械驱动或手动产生移动；将测试用的、已知尺寸的长方体置于该线性滑轨坐标系中，并保证线性滑轨和长方体的一条边相互平行放置；根据长方体的位置，建立世界坐标系，根据其尺寸标注每一个顶点在世界坐标系中的坐标；线性滑轨坐标系与世界坐标系之间仅有平移关系，而没有旋转关系。

（2）推导长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式：将超声探头与位置传感器绑定，一起连接在线性滑轨之上，根据超声探头的位置，建立原始超声图像平面坐标系；假定超声探头产生的原始超声图像平面与线性滑轨的坐标系之间存在旋转角度，使用坐标转换公式将长方体各个顶点坐标从世界坐标系中反推到原始超声图像平面坐标系中；然后在假定原始超声图像平面与线性滑轨之间不存在旋转角度的情况下，对原始超声图像进行三维重建，即将原始超声图像平面坐标系中长方体各顶点坐标转换到线性滑轨坐标系中，因为假定原始超声图像平面与线性滑轨坐标系之间不存在旋转角度的变换通常是不成立的，因此此时计算得到各个顶点在线性滑轨坐标系中的坐标与其在世界坐标系中的坐标点之间产生扭曲，这种扭曲使重建后的长方体的相邻边之间不再保持90度直角关系，这包含了超声探头与线性滑轨之间的旋转角度信息；根据这些转换得到的线性滑轨坐标系下各个顶点坐标重新计算长方体相邻两条边的角度，则可以在理论上推导得到长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式；

（3）线性扫描并测量重建长方体形变量：采用线性扫描方式对标准长方体进行扫描，得到原始二维超声图像序列，假定超声探头与线性滑轨之间不存在坐标转换关系，对原始超声图像进行三维重建，在重建后的三维图像中测量长方体各个顶点在线性滑轨坐标系下的重建坐标，进而计算长方体相邻边之间夹角变化，即长方体形变量。超声探头与线性滑轨之间实际上是存在旋转角度的，实际上线性扫描得到的原始二维超声图像序列也是在探头与线性滑轨之间存在旋转角度的情况下得到的，而在步骤（2）中，推导得到的数学关系式，是在假定线性滑轨与超声探头之间不存在坐标转换的情况下得到的，因此这一步中在对二维超声图像序列进行重建时假定线性滑轨与探头不存在坐标转换关系，计算出重建长方体的形变量，代入步骤（2）推导出的公式，可以最终确定超声探头与线性滑轨的旋转角度。

（4）标定结果计算：将步骤（3）中测量得到的长方体形变量代入到步骤（2）中推导出的长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式，计算出超声探头与线性滑轨之间的旋转角度，即为系统的位置标定结果。

优选的，所述步骤（1）中，采用平行度测试仪来保证线性滑轨扫描方向与长方体上表面的一条边相互平行。

作为另一种优选方案，所述步骤（1）中，为保证线性滑轨扫描方向与长方体上表面的一条边相互平行，采用如下方法：将长方体和线性滑轨置于水平桌面上，然后调节线性滑轨的高度使其与长方体等高，接着调整线性滑轨的方向使其平行于长方体上表面的一条边，最后再调节线性滑轨的高度到合适高度。在实际应用中，任何其他可行的方法都可以用来测量这两者的平行度。

优选的，所述步骤（1）中，长方体的某个顶点位于建立的世界坐标系的原点位置，从而根据长方体已知的尺寸，确定长方体其他顶点在世界坐标系中的坐标，设长方体ABCDHEFG的顶点A位于世界坐标系的原点位置，长方体的三条边AE、AD和AB分别位于坐标系的三条坐标轴上，三条边的边长分别为：Xe、Yd、Zb，则长方体中A、B、D、E四个顶点的坐标分别为：

^cXa=(0,0,0)，^cXb=(0,0,Zb)，^cXd=(0,Yd,0)，^cXe=(Xe,0,0)。     （1）

优选的，所述步骤（2）中，首先将长方体各个顶点坐标从世界坐标系中转换到原始超声图像平面坐标系中，具体步骤如下：

将超声探头连接在线性滑轨之上，假定超声探头与线性滑轨之间存在一个旋转角度，根据超声探头的位置，建立原始超声图像平面坐标系，然后将世界坐标系下长方体各个顶点坐标转换到原始超声图像平面坐标系中，坐标转换公式如下：

^cX=^cT_l ^lT_p ^pX       （2）

将公式（2）做变换得：

^pX=(^lT_p)^-1(^cT_l)^-1(^cX)；       （3）

其中：^pX、^cX分别表示原始超声图像平面坐标系和世界坐标系中长方体各个顶点的顶点位置；^lT_p表示从原始超声图像平面坐标系到线性滑轨坐标系的转换表达式；^cT_l表示从线性滑轨坐标系到世界坐标系的转换表达式。坐标转换矩阵表达式如下：

$T_I^J(x,y,z,\mathrm{\α},\mathrm{\β},\mathrm{\γ})=\left[\begin{array}{cccc}\cos \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& x\\ \sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ}& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& y\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\β}\cos \mathrm{\γ}& z\\ 0& 0& 0& 1\end{array}\right];---\left(4\right)$

公式中，x、y、z表示待转换的两个坐标系的相对位移，α、β、γ表示待转换的两个坐标系在三个坐标轴X、Y、Z上的旋转角度。

对于线性滑轨坐标系与原始超声图像平面坐标系，在公式（4）中，因为超声探头设置在线性滑轨上，所以线性滑轨坐标系与原始超声图像平面坐标系之间的相对位移为：x=0，y＝0，z=Pz，将这些值代入公式（4）得到：

$T_p^l=\left[\begin{array}{cccc}\cos \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& 0\\ \sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ}& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& 0\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\β}\cos \mathrm{\γ}& \mathrm{Pz}\\ 0& 0& 0& 1\end{array}\right];---\left(5\right)$

对于线性滑轨坐标系与世界坐标系，根据步骤（1）中长方体和线性滑轨的位置摆放关系，即线性滑轨扫描方向与长方体上表面的一条边相互平行，则可知世界坐标系与线性滑轨坐标系的旋转角度为零，即α=0、β=0、γ=0。两个坐标系之间存在位移关系，设相对位移为：x=H、y＝W、z=0，代入坐标转换矩阵（4），得：

$T_l^c=\left[\begin{array}{cccc}1& 0& 0& H\\ 0& 1& 0& W\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right];---\left(6\right)$

将坐标转换矩阵（5）和（6）代入坐标转换公式（3）中，得：

$\left[\begin{array}{c}{P}_x\\ P_y\\ 0\\ 1\end{array}\right]={\left[\begin{array}{cccc}\cos \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& 0\\ \sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ}& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& 0\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\β}\cos \mathrm{\γ}& \mathrm{Pz}\\ 0& 0& 0& 1\end{array}\right]}^{-1}*{\left[\begin{array}{cccc}1& 0& 0& H\\ 0& 1& 0& W\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]}^{-1}*\left[\begin{array}{c}{V}_x\\ V_y\\ V_z\\ 1\end{array}\right]---\left(7\right)$

将长方体四个顶点在世界坐标系中坐标^cXa=(0,0,0),^cXb=(0,0,Zb),^cXd＝(0,Yd,0)和^cXe=(Xe,0,0)分别代入公式（7）中的(Vx,Vy,Vz)得出原始超声图像平面坐标系下的长方体各个顶点坐标：

^pXa=(P_zasinβ-Wsinαcosβ-Hcosαcosβ，                    （8）

-P_zacosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

^pXb=(-Zbsinβ+P_zbsinβ-Wsinαcosβ-Hcosαcosβ，Zbcosβsinγ

-P_zbcosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

^pXd=(Ydsinαcosβ+P_zcsinβ-Wsinαcosβ-Hcosαcosβ，Ydsinαsinβsinγ+Ydcosαsinγ

-P_zccosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

^pXe=(Xecosαcosβ+P_zesinβ-Wsinαcosβ-Hcosαcosβ,Xecosαsinβsinγ-Xesinαcosγ

-P_zecosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

其中Pza、Pzb、Pzd和Pze分别表示长方体四个顶点A、B、D和E经过坐标转换，在原始超声图像平面坐标系中的对应点相对于线性滑轨坐标系Z轴方向的位置，这些未知量在后面的计算中将被约掉。

优选的，所述步骤（2）中，将原始超声图像平面坐标系中各顶点坐标转换到线性滑轨坐标系中，具体包括以下步骤：

设线性滑轨坐标系与原始超声图像平面坐标系之间在三个坐标轴X、Y、Z上的旋转角度为α、β、γ，且均为0，则原始超声图像与线性滑轨之间的坐标转换关系式为：

^cX′=^lT_p′^pX；     （9）

其中^cX'表示线性滑轨坐标系下长方体各顶点的坐标，^pX表示原始超声图像平面坐标系中长方体各个顶点的顶点位置；^lT_p′表示从原始超声图像平面坐标系到线性滑轨坐标系的转换表达式，其具体表达式为：

$X_{\′}^c=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& \mathrm{Pz}\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{Px}\\ \mathrm{Py}\\ 0\\ 1\end{array}\right];---\left(10\right)$

将公式（8）中的原始超声图像平面坐标系下长方体的顶点坐标代入坐标转换公式（10）中的(Px,Py)，而四个顶点A、B、D、E在线性滑轨坐标系Z轴方向的位置分别为Pza、Pzb、Pzd和Pze，最后得到四个顶点在线性滑轨坐标系中的对应点a'、b'、d'、e'，以及它们的坐标：

^cXa′=((Wcosα-Hsinα)tanβtanγ-(Hcosα+Wsinα)/cosβ，

-(Wcosα-Hsinα)/cosγ,                           （11）

(Wcosαtanγ-Hsinαtanγ)/cosβ-(Hcosαtanβ+Wsinαtanβ))

^cXb'=((Wcosα-Hsinα)tanβtanγ-(Hcosα+Wsinα)/cosβ，

-(Wcosα-Hsinα)/cosγ,

Zb+(Wcosαtanγ-Hsinαtanγ)/cosβ-(Hcosαtanβ+Wsinαtanβ))

^cXd'=((Wcosα-Hsinα-Ydcosα)tanβtanγ-(Hcosα+Wsinα-Ydsinα)/cosβ，

-(Wcosα-Hsinα-Ydcosα)/cosγ,

(Wcosα-Hsinα-Ydcosα)*tanγ/cosβ-(Hcosα+Wsinα-Ydsinα)*tanβ)

^cXe'=((Wcosα-Hsinα+Xesinα)tanβtanγ-(Hcosα+Wsinα-Xecosα)/cosβ，

-(Wcosα-Hsinα+Xesinα)/cosγ,。

(Wcosα-Hsinα+Xesinα)*tanγ/cosβ-(Hcosα+Wsinα-Xecosα)*tanβ)

优选的，所述步骤（2）中，根据公式（11）中各个顶点在线性滑轨坐标系中的坐标，得出长方体的相邻两条边的边长为：

a'b'＝Zb；              （12）

$a^{\′}{d}^{\′}=Y_d*\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2{\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}/\left(\cos \mathrm{\β}\cos \mathrm{\γ}\right)$

$a^{\′}{e}^{\′}=X_e*\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}/\left(\cos \mathrm{\β}\cos \mathrm{\γ}\right);$

将边长代入到余弦定理公式：

$\∠A=\arccos \frac{{\mathrm{AB}}^2+{\mathrm{AC}}^2-{\mathrm{BC}}^2}{2*\mathrm{AB}*\mathrm{AC}};---\left(13\right)$

计算出长方体顶点A经过两次坐标转换之后，在线性滑轨坐标系下的对应点a'处每两条边的夹角：

$\∠b^{\′}{a}^{\′}{d}^{\′}=\arccos \left(\frac{\sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)---\left(14\right)$

$\∠b^{\′}{a}^{\′}{e}^{\′}=\arccos \left(\frac{\cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\sin \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{\left({2-\cos }^2\mathrm{\β}\right)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-{2\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)$

$\∠d^{\′}{a}^{\′}{e}^{\′}=\arccos (\frac{-\sin \left(2\mathrm{\α}\right)({\cos }^2\mathrm{\β}{\cos }^2\mathrm{\γ}-\cos \left(2\mathrm{\γ}\right)-\cos \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}};$

$*\frac{1}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}})$

以上三个公式即为线性滑轨坐标系中长方体形变量（用重建长方体中顶点a'处每两条边的夹角表示）与超声探头和线性滑轨之间的旋转角度的数学关系式。

优选的，所述步骤（3）中，采用线性扫描方式扫描已知长方体，并对采集得到的原始超声图像序列进行三维重建，测量重建长方体中四个顶点A'、B'、D'、E'的坐标，并使用余弦定理计算顶点A'处每两条边的夹角∠B'A'D′、∠B'A'E'和∠D'A'E'。

优选的，所述步骤（4）中，计算标定结果具体步骤是：

将步骤（3）中测量得到的长方体形变量，即相邻边之间夹角∠B'A'D'、∠B'A'E'和∠D'A'E'代入到步骤（2）中推导出的长方体形变量与超声探头和线性滑轨之间旋转角度的数学关系式（14）中，得到如下表达式：

$\∠B^{\′}{A}^{\′}{D}^{\′}=\arccos \left(\frac{\sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)$

$\∠B^{\′}{A}^{\′}{E}^{\′}=\arccos \left(\frac{\cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\sin \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{\left({2-\cos }^2\mathrm{\β}\right)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2-2\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)---\left(15\right)$

$\∠D^{\′}{A}^{\′}{E}^{\′}=\arccos (\frac{-\sin \left(2\mathrm{\α}\right)({\cos }^2\mathrm{\β}{\cos }^2\mathrm{\γ}-\cos \left(2\mathrm{\γ}\right)-\cos \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}},$

$*\frac{1}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}})$

进而求解得到的解向量即超声探头和线性滑轨之间的旋转角度的数学关系。

为了求解出α、β、γ的值，将等号左边式子移到等号右边得到三个表达式，分别用

和

表示，其中

为待求解向量，如下所示：

$f_1\stackrel{\→}{\left(x\right)}=\∠B^{\′}{A}^{\′}{D}^{\′}-\arccos (\frac{\sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2{\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

$f_2\stackrel{\→}{\left(x\right)}=\∠B^{\′}{A}^{\′}{D}^{\′}-\arccos (\frac{\cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\sin \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

$f_3\stackrel{\→}{\left(x\right)}=\∠D^{\′}{A}^{\′}{E}^{\′}=\arccos (\frac{-\sin \left(2\mathrm{\α}\right)({\cos }^2\mathrm{\β}{\cos }^2\mathrm{\γ}-\cos \left(2\mathrm{\γ}\right)-\cos \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

$*\frac{1}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

最后使用非线性方程求解方法求解以上非线性方程，本发明采用最小二乘法计算非线性方程组的根。matlab中的fsolve函数即使用了最小二乘法，该函数的具体表达式如下：

$\stackrel{\→}{x}=\mathrm{fsolve}\left(\right[f_1\stackrel{\→}{\left(x\right)},f_2\stackrel{\→}{\left(x\right)},f_3\stackrel{\→}{\left(x\right)}],\stackrel{\→}{x_0})$

其中

为初始解向量，将初始解向量设为(0,0,0)。最小二乘法通过计算 $F\stackrel{\→}{\left(x\right)}=\sqrt{{f_1}^2\stackrel{\→}{\left(x\right)}+{f_2}^2\stackrel{\→}{\left(x\right)}+{f_3}^2\stackrel{\→}{\left(x\right)}}$ 的值，通过不断调整

向量使得该目标函数值最小，从而得到的解向量

（即为最优解），即超声探头和线性滑轨之间的旋转角度。

本发明的工作过程：

本发明提出在线性扫描方式下沿着一个方向扫描一个标准长方体对象。通过几何关系变换，本发明推导并建立了重建长方体形变量与超声探头和线性滑轨之间的旋转角度的数学关系式，然后通过线性扫描方式下的三维重建实验测量重建长方体的形变量，代入数学关系式得到超声探头与线性滑轨之间的旋转角度，以此可以实现对一基于线性扫描方式的医学超声三维成像的位置标定。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明创造性地提出一种简单易懂、可操作性强的医学超声三维成像线性扫描方式下位置标定方法，该标定方法是针对一个已知尺寸的长方体进行试验。首先在滑轨坐标系与原始超声图像平面坐标系存在旋转角度的情况下，将世界坐标系中长方体各个顶点坐标转换到原始超声图像平面坐标系中，接着在假定超声探头和线性滑轨之间不存在旋转角度转换的情况下，对原始超声图像进行三维重建，这一过程将长方体各个顶点坐标从原始超声图像平面坐标系转换到线性滑轨坐标系中。因为假定超声探头与线性滑轨之间不存在旋转变换通常不是严格成立的，因此重建后的长方体会产生扭曲（形变）。通过推导长方体形变量与超声探头和线性滑轨旋转角度之间的数学关系式可以确定原始超声图像平面与线性滑轨之间的旋转角度关系；通过实际线性扫描方式下的三维重建实验测量重建长方体形变量，代入数学关系式，即可对超声探头相对线性滑轨的坐标变换关系进行标定。

2、本发明方法对是否采集到包含关键信息的二维图像帧依赖性不高，只需一次扫描即可，因此试验简单、标定方便准确，可用于基于自由臂线性扫描或者机械驱动线性扫描方式的三维超声成像系统。

3、本发明提出的已知尺寸的标准长方体进行位置标定的方法，相比以往的扫描关键标记点方法（比如线性排列小球），制作方法更加简单，可以通过直接购置各种标准长方体获得，无需专门定做，廉价快捷，易于实现。

附图说明

图1是从世界坐标系到原始超声图像平面直角坐标系的坐标转换示意图；

图2是重建长方体在世界坐标系中的示意图；

图3是本发明中位置标定方法的流程图；

图4是本实施例1中扫描空间中交于一点的三条轴重建所得效果图；

图5是本实施例1中实验测量重建物体坐标示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1

实施本发明首先需要从理论上推导出长方体形变量与超声探头和线性滑轨之间旋转角度的数学关系式，然后使用超声设备仪，将采集到的二维超声图像序列保存到计算机，使用三维成像软件进行三维重建并显示，然后测量重建物体位置，通过坐标转换计算出超声探头相对线性滑轨的位置，便能较好地实施本发明。

如图3所示，一种基于线性扫描的医学超声三维成像的位置标定方法，首先为线性滑轨建立线性滑轨坐标系，然后将一个测试用的、已知尺寸大小的长方体置于该线性滑轨坐标系中，并保证线性滑轨和长方体的一条边相互平行放置，以长方体一个顶点为原点建立世界坐标系，确定该世界坐标系下每个顶点的坐标；假定超声探头与线性滑轨之间存在一个旋转角度，将世界坐标系下的每个顶点坐标转换到原始超声图像平面坐标系中；再假定超声探头与线性滑轨之间不存在坐标转换，再将原始超声图像平面坐标系下的顶点坐标转换到线性滑轨坐标系中得到这些顶点新的坐标；根据各个顶点的转换坐标计算出每两条边的夹角与超声探头和线性滑轨之间的旋转角度的关系，得出线性扫描方式下探头位置标定公式；然后线性扫描并测量重建长方体形变量，根据标定公式计算出超声探头与线性滑轨之间的旋转角度，即为位置标定结果。

具体包括以下步骤：

（1）初始化：为线性扫描用的线性滑轨建立一个线性滑轨坐标系，一个位置传感器与线性滑轨固定相连，并可以通过机械驱动或手动产生移动；将测试用的、已知尺寸的长方体置于该线性滑轨坐标系中，并保证线性滑轨和长方体的一条边相互平行放置；根据长方体的位置，建立世界坐标系，根据其尺寸标注每一个顶点在世界坐标系中的坐标；

（2）推导长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式：将超声探头与位置传感器绑定，一起连接在线性滑轨之上，根据超声探头的位置，建立原始超声图像平面坐标系；假定超声探头产生的原始超声图像平面与线性滑轨的坐标系之间存在旋转角度，使用坐标转换公式将长方体各个顶点坐标从世界坐标系中反推到原始超声图像平面坐标系中；然后在假定原始超声图像平面与线性滑轨之间不存在旋转角度的情况下，对原始超声图像进行三维重建，即将原始超声图像平面坐标系中长方体各顶点坐标转换到线性滑轨坐标系中，因为假定原始超声图像平面与线性滑轨坐标系之间不存在旋转角度的变换通常是不成立的，因此此时计算得到各个顶点在线性滑轨坐标系中的坐标与其在世界坐标系中的坐标点之间产生扭曲，这种扭曲使重建后的长方体的相邻边之间不再保持90度直角关系，这包含了超声探头与线性滑轨之间的旋转角度信息；根据这些转换得到的线性滑轨坐标系下各个顶点坐标重新计算长方体相邻两条边的角度，则可以在理论上推导得到长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式；

（3）线性扫描并测量重建长方体形变量：采用线性扫描方式对标准长方体进行扫描，得到原始二维超声图像序列，假定超声探头与线性滑轨之间不存在坐标转换关系，对原始超声图像进行三维重建，在重建后的三维图像中测量长方体各个顶点在线性滑轨坐标系下的重建坐标，进而计算长方体相邻边之间夹角变化，即长方体形变量；

（4）标定结果计算：将步骤（3）中测量得到的长方体形变量代入到步骤（2）中推导出的长方体形变量与超声探头和线性滑轨之间旋转角度的数学关系式，计算出超声探头与线性滑轨之间的旋转角度，即为系统的位置标定结果。

本实施例中，采用平行度测试仪来保证线性滑轨扫描方向与长方体上表面的一条边相互平行。

为保证线性滑轨扫描方向与长方体上表面的一条边相互平行，采用如下方法：将长方体和线性滑轨置于水平桌面上，然后调节线性滑轨的高度使其与长方体等高，接着调整线性滑轨的方向使其平行于长方体上表面的一条边，最后再调节线性滑轨的高度到合适高度。在实际应用中，任何其他可行的方法都可以用来测量这两者的平行度。

本实施例中，长方体的某个顶点位于建立的世界坐标系的原点位置，从而根据长方体已知的尺寸，确定长方体其他顶点在世界坐标系中的坐标，设长方体ABCDHEFG的顶点A位于世界坐标系的原点位置，长方体的三条边AE、AD和AB分别位于坐标系的三条坐标轴上，三条边的边长分别为：Xe、Yd、Zb，则长方体中A、B、D、E四个顶点的坐标分别为：

^cXa=(0,0,0)，^cXb=(0,0,Zb)，^cXd=(0,Yd,0)，^cXe=(Xe,0,0)。        （1）

本实施例中，首先将长方体各个顶点坐标从世界坐标系中转换到原始超声图像平面坐标系中，具体步骤如下：

将超声探头连接在线性滑轨之上，假定超声探头与线性滑轨之间存在一个旋转角度，根据超声探头的位置，建立原始超声图像平面坐标系，然后将世界坐标系下长方体各个顶点坐标转换到原始超声图像平面坐标系中，坐标转换公式如下：

^cX=^cT_l ^lT_p ^pX          （2）

将公式（2）做变换得：

^pX=(^lT_p)^-1(^cT_l)^-1(^cX)；          （3）

其中：^pX、^cX分别表示原始超声图像平面坐标系和世界坐标系中长方体各个顶点的顶点位置；^lT_p表示从原始超声图像平面坐标系到线性滑轨坐标系的转换表达式；^cT_l表示从线性滑轨坐标系到世界坐标系的转换表达式。坐标转换矩阵表达式如下：

$T_I^J(x,y,z,\mathrm{\α},\mathrm{\β},\mathrm{\γ})=\left[\begin{array}{cccc}\cos \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& x\\ \sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ}& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& y\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\β}\cos \mathrm{\γ}& z\\ 0& 0& 0& 1\end{array}\right];---\left(4\right)$

公式中，x、y、z表示待转换的两个坐标系的相对位移，α、β、γ表示待转换的两个坐标系在三个坐标轴X、Y、Z上的旋转角度。

对于线性滑轨坐标系与原始超声图像平面坐标系，在公式（4）中，因为超声探头设置在线性滑轨上，所以线性滑轨坐标系与原始超声图像平面坐标系之间的相对位移为：x=0，y＝0，z=Pz，将这些值代入公式（4）得到：

$T_p^l=\left[\begin{array}{cccc}\cos \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& 0\\ \sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ}& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& 0\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\β}\cos \mathrm{\γ}& \mathrm{Pz}\\ 0& 0& 0& 1\end{array}\right];---\left(5\right)$

对于线性滑轨坐标系与世界坐标系，根据步骤（1）中长方体和线性滑轨的位置摆放关系，即线性滑轨扫描方向与长方体上表面的一条边相互平行，则可知世界坐标系与线性滑轨坐标系的旋转角度为零，即α=0、β=0、γ=0。两个坐标系之间存在位移关系，设相对位移为：x=H、y＝W、z=0，代入坐标转换矩阵（4），得：

$T_l^c=\left[\begin{array}{cccc}1& 0& 0& H\\ 0& 1& 0& W\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right];---\left(6\right)$

将坐标转换矩阵（5）和（6）代入坐标转换公式（3）中，得：

$\left[\begin{array}{c}{P}_x\\ P_y\\ 0\\ 1\end{array}\right]={\left[\begin{array}{cccc}\cos \mathrm{\α}\cos \mathrm{\β}& \cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}-\sin \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\γ}& 0\\ \sin \mathrm{\α}\cos \mathrm{\β}& \sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}+\cos \mathrm{\α}\cos \mathrm{\γ}& \sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}& 0\\ -\sin \mathrm{\β}& \cos \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\β}\cos \mathrm{\γ}& \mathrm{Pz}\\ 0& 0& 0& 1\end{array}\right]}^{-1}*{\left[\begin{array}{cccc}1& 0& 0& H\\ 0& 1& 0& W\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]}^{-1}*\left[\begin{array}{c}{V}_x\\ V_y\\ V_z\\ 1\end{array}\right]---\left(7\right)$

将长方体四个顶点在世界坐标系中坐标^cXa=(0,0,0),^cXb=(0,0,Zb),^cXd＝(0,Yd,0)和^cXe=(Xe,0,0)分别代入公式（7）中的(Vx,Vy,Vz)得出原始超声图像平面坐标系下的长方体各个顶点坐标：

^pXa=(P_zasinβ-Wsinαcosβ-Hcosαcosβ，                    （8）

-P_zacosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

^pXb=(-Zbsinβ+P_zbsinβ-Wsinαcosβ-Hcosαcosβ，Zbcosβsinγ

-P_zbcosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

^pXd=(Ydsinαcosβ+P_zcsinβ-Wsinαcosβ-Hcosαcosβ，Ydsinαsinβsinγ+Ydcosαsinγ

-P_zccosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

^pXe=(Xecosαcosβ+P_zesinβ-Wsinαcosβ-Hcosαcosβ,Xecosαsinβsinγ-Xesinαcosγ

-P_zecosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ)

其中Pza、Pzb、Pzd和Pze分别表示长方体四个顶点A、B、D和E经过坐标转换，在原始超声图像平面坐标系中的对应点相对于线性滑轨坐标系Z轴方向的位置，这些未知量在后面的计算中将被约掉。

本实施例中，将原始超声图像平面坐标系中各顶点坐标转换到线性滑轨坐标系中，具体包括以下步骤：

设线性滑轨坐标系与原始超声图像平面坐标系之间在三个坐标轴X、Y、Z上的旋转角度为α、β、γ，且均为0，则原始超声图像与线性滑轨之间的坐标转换关系式为：

^cX′=^lT_p′^pX；          （9）

其中^cX'表示线性滑轨坐标系下长方体各顶点的坐标，^pX表示原始超声图像平面坐标系中长方体各个顶点的顶点位置；^lT_p′表示从原始超声图像平面坐标系到线性滑轨坐标系的转换表达式，其具体表达式为：

$X_{\′}^c=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& \mathrm{Pz}\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{Px}\\ \mathrm{Py}\\ 0\\ 1\end{array}\right];---\left(10\right)$

将公式（8）中的原始超声图像平面坐标系下长方体的顶点坐标代入坐标转换公式（10）中的(Px,Py)，而四个顶点A、B、D、E在线性滑轨坐标系Z轴方向的位置分别为Pza、Pzb、Pzd和Pze，最后得到四个顶点在线性滑轨坐标系中的对应点a'、b'、d'、e′，以及它们的坐标：

^cXa'=((Wcosα-Hsinα)tanβtanγ-(Hcosα+Wsinα)/cosβ，

-(Wcosα-Hsinα)/cosγ,                            （11）

(Wcosαtanγ-Hsinαtanγ)/cosβ-(Hcosαtanβ+Wsinαtanβ))

^cXb′=((Wcosα-Hsinα)tanβtanγ-(Hcosα+Wsinα)/cosβ，

-(Wcosα-Hsinα)/cosγ,

Zb+(Wcosαtanγ-Hsinαtanγ)/cosβ-(Hcosαtanβ+Wsinαtanβ))

^cXd'=((Wcosα-Hsinα-Ydcosα)tanβtanγ-(Hcosα+Wsinα-Ydsinα)/cosβ，

-(Wcosα-Hsinα-Ydcosα)/cosγ,

(Wcosα-Hsinα-Ydcosα)*tanγ/cosβ-(Hcosα+Wsinα-Ydsinα)*tanβ)

^cXe'=((Wcosα-Hsinα+Xesinα)tanβtanγ-(Hcosα+Wsinα-Xecosα)/cosβ，

-(Wcosα-Hsinα+Xesinα)/cosγ,。

(Wcosα-Hsinα+Xesinα)*tanγ/cosβ-(Hcosα+Wsinα-Xecosα)*tanβ)

本实施例中，根据公式（11）中各个顶点在线性滑轨坐标系中的坐标，得出长方体的相邻两条边的边长为：

a'b'=Zb；               （12）

$a^{\′}{d}^{\′}=Y_d*\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2{\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}/\left(\cos \mathrm{\β}\cos \mathrm{\γ}\right)$

$a^{\′}{e}^{\′}=X_e*\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}/\left(\cos \mathrm{\β}\cos \mathrm{\γ}\right);$

将边长代入到余弦定理公式：

$\∠A=\arccos \frac{{\mathrm{AB}}^2+{\mathrm{AC}}^2-{\mathrm{BC}}^2}{2*\mathrm{AB}*\mathrm{AC}};---\left(13\right)$

计算出长方体顶点A经过两次坐标转换之后，在线性滑轨坐标系下的对应点a'处每两条边的夹角：

$\∠b^{\′}{a}^{\′}{d}^{\′}=\arccos \left(\frac{\sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)---\left(14\right)$

$\∠b^{\′}{a}^{\′}{e}^{\′}=\arccos \left(\frac{\cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\sin \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{\left({2-\cos }^2\mathrm{\β}\right)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-{2\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)$

$\∠d^{\′}{a}^{\′}{e}^{\′}=\arccos (\frac{-\sin \left(2\mathrm{\α}\right)({\cos }^2\mathrm{\β}{\cos }^2\mathrm{\γ}-\cos \left(2\mathrm{\γ}\right)-\cos \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}};$

$*\frac{1}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}})$

以上三个公式即为线性滑轨坐标系中长方体形变量（用重建长方体中顶点a'处每两条边的夹角表示）与超声探头和线性滑轨之间的旋转角度的数学关系式。

本实施例中，采用线性扫描方式扫描已知长方体，并对采集得到的原始超声图像序列进行三维重建，测量重建长方体中四个顶点A'、B'、D'、E'的坐标，并使用余弦定理计算顶点A'处每两条边的夹角∠B'A'D′、∠B'A'E'和∠D'A'E′。

本实施例中，将步骤（3）中测量得到的长方体形变量，即相邻边之间夹角∠B'A'D'、∠B'A'E'和∠D'A'E'代入到步骤（2）中推导出的长方体形变量与超声探头和线性滑轨之间旋转角度的数学关系式（14）中，得到如下表达式：

$\∠B^{\′}{A}^{\′}{D}^{\′}=\arccos \left(\frac{\sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)$

$\∠B^{\′}{A}^{\′}{E}^{\′}=\arccos \left(\frac{\cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\sin \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{\left({2-\cos }^2\mathrm{\β}\right)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2-2\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}\right)---\left(15\right)$

$\∠D^{\′}{A}^{\′}{E}^{\′}=\arccos (\frac{-\sin \left(2\mathrm{\α}\right)({\cos }^2\mathrm{\β}{\cos }^2\mathrm{\γ}-\cos \left(2\mathrm{\γ}\right)-\cos \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}},$

$*\frac{1}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}})$

进而求解得到的解向量

即超声探头和线性滑轨之间的旋转角度的数学关系。

为了求解出α、β、γ的值，将等号左边式子移到等号右边得到三个表达式，分别用

和

表示，其中

为待求解向量，如下所示：

$f_1\stackrel{\→}{\left(x\right)}=\∠B^{\′}{A}^{\′}{D}^{\′}-\arccos (\frac{\sin \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\cos \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2{\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

$f_2\stackrel{\→}{\left(x\right)}=\∠B^{\′}{A}^{\′}{D}^{\′}-\arccos (\frac{\cos \mathrm{\α}\sin \mathrm{\β}\cos \mathrm{\γ}-\sin \mathrm{\α}\sin \mathrm{\γ}}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

$f_3\stackrel{\→}{\left(x\right)}=\∠D^{\′}{A}^{\′}{E}^{\′}=\arccos (\frac{-\sin \left(2\mathrm{\α}\right)({\cos }^2\mathrm{\β}{\cos }^2\mathrm{\γ}-\cos \left(2\mathrm{\γ}\right)-\cos \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}{\sqrt{({\cos }^2\mathrm{\β}-2)*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+{2\cos }^2\mathrm{\α}-\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

$*\frac{1}{\sqrt{(2-{\cos }^2\mathrm{\β})*{\cos }^2\mathrm{\γ}*\cos \left(2\mathrm{\α}\right)+2-2{\cos }^2\mathrm{\α}+\sin \left(2\mathrm{\α}\right)\sin \mathrm{\β}\sin \left(2\mathrm{\γ}\right)}}$

最后使用非线性方程求解方法求解以上非线性方程，本发明采用最小二乘法计算非线性方程组的根。matlab中的fsolve函数即使用了最小二乘法，该函数的具体表达式如下：

$\stackrel{\→}{x}=\mathrm{fsolve}\left(\right[f_1\stackrel{\→}{\left(x\right)},f_2\stackrel{\→}{\left(x\right)},f_3\stackrel{\→}{\left(x\right)}],\stackrel{\→}{x_0})$

其中

为初始解向量，将初始解向量设为(0,0,0)。最小二乘法通过计算 $F\stackrel{\→}{\left(x\right)}=\sqrt{{f_1}^2\stackrel{\→}{\left(x\right)}+{f_2}^2\stackrel{\→}{\left(x\right)}+{f_3}^2\stackrel{\→}{\left(x\right)}}$ 的值，通过不断调整

向量使得该目标函数值最小，从而得到的解向量（即为最优解），即超声探头和线性滑轨之间的旋转角度。

下面的实验是采用本实施例所述的方法标定实验室所用的线性扫描滑轨上探头位置。

本实验采用线性扫描方式扫描三维空间中由垂直相交于一点的三条轴构成的直角支架，已知每条轴长度分别为21.0mm、24.0mm、58.0mm，用三维重建软件对采集到的二维超声图像进行重建，三维重建结果如图4所示，测量重建三维图形中每条轴上端点坐标，如图5所示，然后计算每两条轴的角度，最后根据每两条轴的角度计算超声探头与线性滑轨之间的旋转角度，三次实验结果如下：

	∠B'A'D'	∠B'A'E'	∠D'A'E'	α	β	γ
							实验1	91.65	87.86	90.12	1.31	2.08	1.83
实验2	91.77	88.62	90.15	1.44	1.95	1.84
							实验3	91.53	88.40	90.14	1.27	2.02	1.84

从表中可以看出，测量重建三维图像中每两条轴的夹角，可以确定超声探头与线性滑轨之间的旋转角度，即线性位置传感器标定结果，验证了本发明的有效性。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，首先为线性扫描用的线性滑轨建立线性滑轨坐标系，然后将一个测试用的、已知尺寸大小的长方体置于该线性滑轨坐标系中，并保证线性滑轨和长方体的一条边相互平行放置，以长方体一个顶点为原点建立世界坐标系，确定该世界坐标系下每个顶点的坐标；假定超声探头与线性滑轨之间存在一个旋转角度，将世界坐标系下的每个顶点坐标转换到原始超声图像平面坐标系中；再假定超声探头与线性滑轨之间不存在坐标转换，再将原始超声图像平面坐标系下的顶点坐标转换到线性滑轨坐标系中得到这些顶点新的坐标；根据各个顶点的转换坐标计算出每两条边的夹角与超声探头和线性滑轨之间的旋转角度的关系，得出线性扫描方式下探头位置标定公式；然后线性扫描并测量重建长方体形变量，根据标定公式计算出超声探头与线性滑轨之间的旋转角度，即为位置标定结果。 

2.根据权利要求1所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，包括以下步骤： 

（1）初始化：为线性扫描用的线性滑轨建立一个线性滑轨坐标系，一个位置传感器与线性滑轨固定相连，并可以通过机械驱动或手动产生移动；将测试用的、已知尺寸的长方体置于该线性滑轨坐标系中，并保证线性滑轨和长方体的一条边相互平行放置；根据长方体的位置，建立世界坐标系，根据其尺寸标注每一个顶点在世界坐标系中的坐标； 

（2）推导长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式：将超声探头与位置传感器绑定，一起连接在线性滑轨之上，根据超声探头的位置，建立原始超声图像平面坐标系；假定超声探头产生的原始超声图像平面与线性滑轨的坐标系之间存在旋转角度，使用坐标转换公式将长方体各个顶点坐标从世界坐标系中反推到原始超声图像平面坐标系中；然后在假定原始超声图像平面与线性滑轨之间不存在旋转角度的情况下，对原始超声图像进行三维重建，即将原始超声图像平面坐标系中长方体各顶点坐标转换到线性滑轨坐标系中；根据转换得到的线性滑轨坐标系下各个顶点坐标重新计算长方体相邻两条边的角度，则在理论上推导得到长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式； 

（3）线性扫描并测量重建长方体形变量：采用线性扫描方式对标准长方体 进行扫描，得到原始二维超声图像序列，假定超声探头与线性滑轨之间不存在坐标转换关系，对原始超声图像进行三维重建，在重建后的三维图像中测量长方体各个顶点在线性滑轨坐标系下的重建坐标，进而计算长方体相邻边之间夹角变化，即长方体形变量； 

（4）标定结果计算：将步骤（3）中测量得到的长方体形变量代入到步骤（2）中推导出的长方体形变量与超声探头与线性滑轨之间旋转角度的数学关系式，计算出超声探头与线性滑轨之间的旋转角度，即为系统的位置标定结果。 

3.根据权利要求2所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（1）中，采用平行度测试仪来保证线性滑轨扫描方向与长方体上表面的一条边相互平行。 

4.根据权利要求2所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（1）中，为保证线性滑轨扫描方向与长方体上表面的一条边相互平行，采用如下方法：将长方体和线性滑轨置于水平桌面上，然后调节线性滑轨的高度使其与长方体等高，接着调整线性滑轨的方向使其平行于长方体上表面的一条边，最后再调节线性滑轨的高度到合适高度。 

5.根据权利要求3或4所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（1）中，长方体的某个顶点位于建立的世界坐标系的原点位置，从而根据长方体已知的尺寸，确定长方体其他顶点在世界坐标系中的坐标，设长方体ABCDHEFG的顶点A位于世界坐标系的原点位置，长方体的三条边AE、AD和AB分别位于坐标系的三条坐标轴上，三条边的边长分别为：Xe、Yd、Zb，则长方体中A、B、D、E四个顶点的坐标分别为： 

^cXa=(0,0,0)，^cXb=(0,0,Zb)，^cXd=(0,Yd,0)，^cXe=(Xe,0,0)。      （1） 

6.根据权利要求5所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（2）中，首先将长方体各个顶点坐标从世界坐标系中转换到原始超声图像平面坐标系中，具体步骤如下：将超声探头连接在线性滑轨之上，假定超声探头与线性滑轨之间存在一个旋转角度，根据超声探头的位置，建立原始超声图像平面坐标系，然后将世界坐标系下长方体各个顶点坐标转换到原始超声图像平面坐标系中，坐标转换公式如下： 

^cX=^cT_l ^lT_p ^pX         （2） 

将公式（2）做变换得： 

^pX=(^lT_p)^-1(^cT_l)^-1(^cX)；（3） 

其中：^pX、^cX分别表示原始超声图像平面坐标系和世界坐标系中长方体各个顶点的顶点位置；^lT_p表示从原始超声图像平面坐标系到线性滑轨坐标系的转换表达式；^cT_l表示从线性滑轨坐标系到世界坐标系的转换表达式，其中坐标转换矩阵表达式如下： 

公式（4）中，x、y、z表示待转换的两个坐标系的相对位移，α、β、γ表示待转换的两个坐标系在三个坐标轴X、Y、Z上的旋转角度； 

对于线性滑轨坐标系与原始超声图像平面坐标系，在公式（4）中，因为超声探头设置在线性滑轨上，所以线性滑轨坐标系与原始超声图像平面坐标系之间的相对位移为：x＝0，y=0，z＝Pz，将这些值代入公式（4）得到： 

对于线性滑轨坐标系与世界坐标系来说，根据步骤（1）中假定的长方体和线性滑轨的位置摆放关系，即线性滑轨扫描方向与长方体上表面的一条边相互平行，则可知世界坐标系与线性滑轨坐标系的旋转角度为零，即α＝0，β=0，γ＝0，而且两个坐标系之间存在位移关系，假设相对位移为：x＝H、y＝W、z＝0，代入坐标转换矩阵（4），得： 

将坐标转换矩阵（5）和（6）代入坐标转换公式（3）中，得： 

将长方体四个顶点在世界坐标系中坐标 ^cXa=(0,0,0),^cXb=(0,0,Zb),^cXd＝(0,Yd,0)和^cXe=(Xe,0,0)分别代入公式（7）中的(Vx,Vy,Vz)得出原始超声图像平面坐标系下的长方体各个顶点坐标： 

^pXa=(P_zasinβ-Wsinαcosβ-Hcosαcosβ，                      （8） 

-P_zacosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ) 

^pXb=(-Zbsinβ+P_zbsinβ-Wsinαcosβ-Hcosαcosβ，Zbcosβsinγ 

-P_zbcosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ) 

^pXd=(Ydsinαcosβ+P_zcsinβ-Wsinαcosβ-Hcosαcosβ，Ydsinαsinβsinγ+Ydcosαsinγ 

-P_zccosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ) 

^pXe=(Xecosαcosβ+P_zesinβ-Wsinαcosβ-Hcosαcosβ,Xecosαsinβsinγ-Xesinαcosγ 

-P_zecosβsinγ-Wsinαsinβsinγ-Hcosαsinβsinγ+Hsinαcosγ-Wcosαcosγ) 

其中Pza、Pzb、Pzd和Pze分别表示长方体四个顶点A、B、D和E经过坐标转换，在原始超声图像平面坐标系中的对应点相对于线性滑轨坐标系Z轴方向的位置。 

7.根据权利要求6所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（2）中，设线性滑轨坐标系与原始超声图像平面坐标系之间在三个坐标轴X、Y、Z上的旋转角度为α、β、γ，且均为0，使用坐标转换公式将原始超声图像平面坐标系下长方体各个顶点坐标转换到线性滑轨坐标系中，原始超声图像与线性滑轨之间的坐标转换关系式为： 

^cX′=^lT_p′^pX；      （9） 

其中^cX'表示线性滑轨坐标系下长方体各顶点的坐标，^pX表示原始超声图像平面坐标系中长方体各个顶点的顶点位置；^lT_p′表示从原始超声图像平面坐标系到线性滑轨坐标系的转换表达式，其具体表达式为： 

将公式（8）中的原始超声图像平面坐标系下长方体的顶点坐标代入坐标转换公式（10）中的(Px,Py)，而四个顶点A、B、D、E的Pz分别为Pza、Pzb、Pzd和Pze，最后得到四个顶点在线性滑轨坐标系中的对应点a'、b'、d'、e'，以及它 们的坐标： 

^cXa'=((Wcosα-Hsinα)tanβtanγ-(Hcosα+Wsinα)/cosβ， 

-(Wcosα-Hsinα)/cosγ,                           （11） 

(Wcosαtanγ-Hsinαtanγ)/cosβ-(Hcosαtanβ+Wsinαtanβ)) 

^cXb′=((Wcosα-Hsinα)tanβtanγ-(Hcosα+Wsinα)/cosβ， 

-(Wcosα-Hsinα)/cosγ, 

Zb+(Wcosαtanγ-Hsinαtanγ)/cosβ-(Hcosαtanβ+Wsinαtanβ)) 

^cXd'=((Wcosα-Hsinα-Ydcosα)tanβtanγ-(Hcosα+Wsinα-Ydsinα)/cosβ， 

-(Wcosα-Hsinα-Ydcosα)/cosγ, 

(Wcosα-Hsinα-Ydcosα)*tanγ/cosβ-(Hcosα+Wsinα-Ydsinα)*tanβ) 

^cXe'=((Wcosα-Hsinα+Xesinα)tanβtanγ-(Hcosα+Wsinα-Xecosα)/cosβ， 

-(Wcosα-Hsinα+Xesinα)/cosγ,。 

(Wcosα-Hsinα+Xesinα)*tanγ/cosβ-(Hcosα+Wsinα-Xecosα)*tanβ)。 

8.根据权利要求7所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（2）中，根据公式（11）中各个顶点在线性滑轨坐标系中的坐标得出长方体的相邻两条边的边长为： 

a'b'=Zb；                     （12） 

将边长代入到余弦定理公式： 

计算出长方体顶点A经过两次坐标转换之后，在线性滑轨坐标系下的对应点a'处每两条边的夹角： 

以上三个公式即为线性滑轨坐标系中长方体形变量与超声探头与线性滑轨之间的旋转角度α、β、γ的数学关系式。 

9.根据权利要求1所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（3）中，采用线性扫描方式扫描已知长方体，并对采集得到的原始超声图像序列进行三维重建，测量重建长方体中四个顶点A'、B'、D'、E'坐标，并使用余弦定理计算顶点A'处每两条边的夹角∠B'A'D'、∠B'A'E'和∠D'A'E'。 

10.根据权利要求8所述的基于线性扫描的医学超声三维成像的位置标定方法，其特征在于，所述步骤（4）中，将步骤（3）中测量得到的长方体形变量，即相邻边之间夹角∠B'A'D'、∠B'A'E'和∠D'A'E'代入到步骤（2）中推导出的长方体形变量与线性滑轨和超声探头之间旋转角度的数学关系式（14）中，得到如下表达式： 

对上述方程组求解得到最优解向量

即超声探头和线性滑轨之间的旋转角度。 

Images