患者病患部位的3D解剖结构模型^[1]在计算机辅助诊断的定量和定性分析^[2]、识别膝关节生物力学特征^[3]、虚拟重建面部缺损^[4]、手术规划和对齐术中导航^[5-7]等诸多医学图像处理领域都有着广泛的应用. 在医学图像的处理过程中，患者病患部位的3D解剖结构模型可以帮助医生更好地将相关信息转变成可视化数据^[8]，以便对患者的病情做出快速、准确的判断，降低误诊的风险，有效制订后续的治疗计划. 获取患者病患部位的3D解剖结构模型的常见方法是使用成像技术，例如计算机断层扫描技术^[9]（computed tomography，CT）和核磁共振图像技术^[10]（magnetic resonance imaging，MRI）. 然而，在一些偏远地区和资源有限的地区，医院通常无法使用这些技术获得患者病患部位的3D解剖结构模型^[11]. 因此，用非成像技术构建患者病患部位的3D解剖结构模型成为研究热点.

在建立可变模型的过程中，统计形状模型^[12]能够对解剖结构的几何变形进行参数化表示，有效地构建出可变模型. 越来越多的学者通过统计形状模型构建患者病患部位的3D解剖结构模型. Davies^[13]提出在统计形状模型的构建阶段手动标记特征点，在模板样本与目标样本间粗略地建立对应关系，这样不仅耗时长还容易造成操作错误，导致构建的3D模型性能较差. Kelemen等^[14]提出球谱函数描述法（spherical harmonics descriptors method，SPHARM）并用该方法构建出患者病患部位的3D解剖结构模型，但是构建模型的训练样本只适用于球体. Eck等^[15]采用基于球谐函数的3D强度模型分析异染色质焦点的形状和强度，将统计模型与图像强度模型拟合确定最终的模型参数. 与此同时，Davies等^[16]提出最小描述长度法（minimum description length，MDL），通过Euclidean距离计算3D模型与测试样本间的距离完成二者的拟合. 但是在构建模型的过程中需要对每个训练样本手动标记特征点，操作过程复杂且拟合效果较差，不能满足临床应用需求. Perez等^[17]提出标记点滑动法（landmark sliding method，SLIDE），在MDL手动标记特征点的基础上，将初始特征点沿切线方向移动. 但SLIDE需要不断计算每个特征点的切线方向，计算过程复杂并且仅适用在特征点数量较少的图像的情况下.

以上方法忽略了如何对解剖结构模型进行有效的评估. Plessers等^[18]提出通过肩胛骨的统计形状模型实现肩胛骨的重构. 该方法能够有效构建出肩胛骨的统计形状模型，但是肩胛骨的统计形状模型与测试肩胛骨的拟合误差相对较大，并且计算关节盂的倾斜度和倾斜角的过程比较复杂. Nolte等^[19]采用40组女性股骨的CT扫描图像构建统计形状模型. 该方法在构建模型的过程中需要根据三角形表面创建四面体网格，采用1.5的粗化因子将每个股骨用105 064个四面体重新表示. 该过程烦琐且耗时，不能满足临床应用需求.

针对传统建模方法不能有效构建拟合误差较小的3D模型，本研究提出基于点分布模型（point distribution model，PDM）的3D模型拟合方法（3D model fitting method based on PDM, 3DMFB-PDM）. 通过引入附加项和转换非线性方程组，使点分布模型与测试样本具有最小的拟合误差.

训练样本集是指感兴趣区域的3D模型集合. 本研究将感兴趣区域从原始的2D CT扫描图像中分割出来，再将每个2D的感兴趣区域图像用3D模型表示，所有感兴趣区域的3D模型集合构成训练样本集. 模板样本是指从训练样本集中任意选择的单个3D模型. 测试样本是指根据患者病患部位的几何数据形成的3D点云模型. 这些几何数据可以有效表征患者病患部位的3D形态，通过3D位置传感器将患者病患部位的几何数据形成2D图像，再对2D图像进行数据处理，将2D图像以3D点云形式表示. 获取测试样本实际上就是获取患者病患部位的几何数据，使用跟踪探针来数字化在手术过程中暴露的病患部位和解剖标志区域. 在这个过程中，通常将单个3D位置传感装置连接到跟踪探头上，并通过跟踪探头捕捉的几何数据形成患者病患部位的图像，这些几何数据定义了每个图像相对于患者所在坐标系的3D位置和方向. 点分布模型是指训练样本集中的3D模型经过对齐、建立对应关系、主成分分析等过程后最终生成的3D模型，即3D点云模型.

作为统计形状模型，点分布模型是基于对应点位置的变形模型. 模型的基本思想是从训练样本中提取关键点，用其坐标构成的特征向量描述目标样本，通过训练样本获得目标样本的统计信息.

假设 ${{{\boldsymbol{\varOmega }}}} = \left[ {{{{{\boldsymbol{S}}}}_{{1}}},{{{{\boldsymbol{S}}}}_{{2}}}, \cdots ,{{{{\boldsymbol{S}}}}_{{{N}}}}} \right]$， $N$为训练样本的个数， ${{{{\boldsymbol{P}}}}_i} = \left[ {{x_i},{y_i},{z_i}} \right]$， $i = 1,2, \cdots ,n$.这 $N$个训练样本是拓扑同胚的， ${{{{\boldsymbol{P}}}}_{{{i}}}}$为训练样本的三角形网格顶点， $n$为三角形网格顶点的个数，训练样本的点分布模型 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$可以由 $3n$个元素表示为

(1) $ {{{{\boldsymbol{S}}}}_{{{\rm{P}}}}} = {\left[ {{x_1},{y_1},{z_1}, \cdots ,{x_n},{y_n},{z_n}} \right]^{\text{T}}}. $

用正态分布 $\theta $表示训练样本的形状变化：

(2) $ {{{{\boldsymbol{S}}}}_{{{\rm{P}}}}} \sim \theta \left( {{{\bar {\boldsymbol{S}}}},{\boldsymbol{C}}} \right). $

式中： ${{\bar {\boldsymbol{S}}}}$为训练样本的平均形状， ${\boldsymbol{C}}$为协方差矩阵.

(3) $ {{\bar {\boldsymbol{S}}}} = {N^{ - 1}}\sum\nolimits_{i = 1}^N{{\boldsymbol{S}}_i}, $

(4) $ {\boldsymbol{C}} = \left( {{{N- 1}}} \right)^{-1}\sum\nolimits_{i = 1}^N {\left( {{{{{{\boldsymbol{S}}}}}_{{{\boldsymbol{i}}}}} - {{\bar {\boldsymbol{S}}}}} \right)} {\left( {{{{{\boldsymbol{S}}}}_{{{i}}}} - {{\bar {\boldsymbol{S}}}}} \right)^{\rm{T}}}. $

通过统计的方法生成 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$的一般表达式，将其表示为平均形状与主要变形模型之和：

(5) $ {{{{\boldsymbol{S}}}}_{{{\rm{P}}}}} = {{\bar {\boldsymbol{S}}}} + \sum\nolimits_{i = 1}^M {{\lambda _i}{{{{\boldsymbol{p}}}}_i}} . $

式中： $M$为主要变形模型的个数， ${{{{\boldsymbol{p}}}}_{{{i}}}}$为第 $i$个特征向量， ${\lambda _i}$为 ${{{{\boldsymbol{p}}}}_{{{i}}}}$所对应的权重. ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$的形状受 ${\lambda _i}$的影响，因此将 ${\lambda _i}$作为点分布模型的形状参数. 特征向量按降序排列，特征值满足：

(6) $ {\lambda _i} > {\lambda _{i + 1}}. $

当 ${\lambda _i}$快速下降时，点分布模型可以由前 $G$个变形模型精确近似为

(7) $ {{{{\boldsymbol{S}}}}_{{{\rm{P}}}}} = {{\bar {\boldsymbol{S}}}} + \sum\nolimits_{i = 1}^G {{\lambda _i}{{{{\boldsymbol{p}}}}_{{{i}}}}}, $

(8) $ 100\% {\sum\nolimits _{i = 1}^G} {\lambda _i}\Biggr / {\sum\nolimits _{i = 1}^M} {\lambda _i} \geqslant P . $

式中： $P$为用前 $G$个变形模型表示整个变形模型的占比.

在完成训练样本的点分布模型构建之后，须拟合点分布模型与测试样本，即在点分布模型中寻找与测试样本中特征点对应的点，并优化两者间的距离，使点分布模型逐渐适应患者的解剖结构. 在拟合过程中，先将点分布模型进行旋转和平移，使点分布模型与测试样本对齐；再用Euclidean距离计算测试样本特征点与点分布模型对应点间的距离，使间距最小.

假设点分布模型通过旋转变换 ${\boldsymbol{R}}$和平移变换 ${\boldsymbol{T}}$后与测试样本对齐，将对齐后的点分布模型记为 ${{{\boldsymbol{S}}}}_{{{\rm{p}}}}^{\mathbf{'}}$，计算公式为

(9) $ {{{\boldsymbol{S}}}}_{{{\rm{p}}}}^{\mathbf{'}} = {\boldsymbol{R}}\left( {{{\bar {\boldsymbol{S}}}} + \sum\nolimits_{i = 1}^M {{\lambda _i}{{{{\boldsymbol{p}}}}_{{{i}}}}} } \right) + {\boldsymbol{T}} . $

在点分布模型和测试样本拟合时，须对齐点分布模型和测试样本，以减小二者在形状和位置上的差异. 采用普氏分析（procrustes analysis，PA）^[20]进行对齐操作，其中 $R$、 $T$分别为PA中的旋转和尺度参数. 对齐过程实际上是2个模型的归一化过程，通过归一化处理使点分布模型尽可能地接近测试样本.

用目标函数 ${{{{f}}}_{1}}$表示测试样本与点分布模型之间的Euclidean距离：

(10) $\begin{split} {{{{f}}}_{1}} = {\sum\nolimits_{i = 1}^D {\sum\nolimits_{j = 1}^K {\left\| {{{{{\boldsymbol{d}}}}_{{{i}}}} - \left. {{{{\boldsymbol{m}}}}\left( {{{j}}} \right)} \right\|^2} \right.} } } ;\\ i = 1,2, \cdots ,D;j = 1,2, \cdots ,K. \end{split}$

式中： ${{{{\boldsymbol{d}}}}_{{{i}}}}$为测试样本上的第 $i$个点； $D$为测试样本中的特征点的个数； ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$为点分布模型上与 ${{{{\boldsymbol{d}}}}_{{{i}}}}$一一对应的点； $K$为点分布模型上与 ${{{{\boldsymbol{d}}}}_{{{i}}}}$对应的点的个数，满足 $D = K$. ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$是 ${{{\boldsymbol{S}}}}_{{{\boldsymbol{p}}}}^{\mathbf{'}}$上的点，满足：

(11) $ {{{\boldsymbol{m}}}}\left( {{{j}}} \right) = {\boldsymbol{R}}\left( {{{\bar {\boldsymbol{S}}}}\left( {{{j}}} \right) + \sum\nolimits_{i = 1}^M {{\lambda _i}{{{{\boldsymbol{p}}}}_{{{i}}}}\left( {{{j}}} \right)} } \right) + {\boldsymbol{T}} . $

式中： ${{\bar {\boldsymbol{S}}}}\left( {{{j}}} \right)$为平均形状 ${{\bar {\boldsymbol{S}}}}$上的第 $ j $个点， ${{{{\boldsymbol{p}}}}_{{{i}}}}\left( {{{j}}} \right)$为第 $i$个特征向量所对应的点， $ j $为点分布模型上的特征点. 将式（11）带入式（10），得到目标函数 ${{{{f}}}_{1}}$的一般表示式为

(12) $ \begin{split} {{{{f}}}_{1}} = \sum\nolimits_{i = 1}^D {\sum\nolimits_{j = 1}^K\left\| {{\boldsymbol{d}_i} - \left( {{\boldsymbol{R}}\left({ \bar {\boldsymbol{S}} \left( j \right) + \sum\nolimits_{i = 1}^M {{\lambda _i}{\boldsymbol{p}_i}\left( j \right)} } \right) + {\boldsymbol{T}}} \right)} \right\|^2}. \hfill \\ \end{split} $

最小化 ${{{{f}}}_{1}}$能够直接计算出 ${{{{\boldsymbol{d}}}}_{{{i}}}}$和 ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$的距离，但是在点分布模型与测试样本拟合的过程中， ${{{{f}}}_{1}}$计算出的距离并不是最短，原因如下。 1） ${{{{f}}}_{1}}$认为 ${{{{\boldsymbol{d}}}}_{{{i}}}}$和 ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$是2个不相关的点集，二者不服从任何分布. 但是在拟合的过程中， ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$服从正态分布^[21]，并且协方差 ${{{\boldsymbol{C}}}}$分解后产生的模型 ${{{{\boldsymbol{S}}}}_{{{{{\rm{P}}}}}}}$与患者的身高、体重、性别等信息有关. 2） ${{{{f}}}_{1}}$认为 ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$的形状固定，是没有变形的模型. 实际上， ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$是由平均形状和变形模型构成的，而且这些变形模型是相关的，因此变形模型通常可以由前几个主要模型表示.

为了解决 ${{{{\boldsymbol{S}}}}_{{{{\rm{P}}}}}}$与测试样本拟合的过程中 ${{{{f}}}_{1}}$忽略的问题，实现点分布模型与测试样本的最佳拟合，提出在计算 ${{{{\boldsymbol{d}}}}_{{{i}}}}$与 ${{{\boldsymbol{m}}}}\left( {j} \right)$间的距离时，引入Mahalanobis距离 ${{{{\boldsymbol{D}}}}_{{{\rm{m}}}}}$作为附加项，并将非线性的目标函数线性化，然后对线性目标函数求解，实现 ${{{{\boldsymbol{d}}}}_{{{i}}}}$与 ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$之间的距离最短. Mahalanobis距离表示点与单个分布间的距离，它能够有效地衡量二者的相似程度，并且考虑变形模型间的相关性. 令 ${{{{\boldsymbol{P}}}}_{{{{x}}_{{j}}}}}、{{{{\boldsymbol{P}}}}_{{{{y}}_{{j}}}}}、{{{{\boldsymbol{P}}}}_{{{{z}}_{{j}}}}}$分别为点分布模型上的点 $j$在 $x$、 $y$、 $z$方向上的基向量，这3个基向量的Mahalanobis距离可以表示为

(13) $ {{{{D}}}_{{{\rm{m}}}}}\left( {{{{{\boldsymbol{P}}}}_{{{k}}}}} \right) = {{{{\boldsymbol{P}}}}_{{{k}}}}^{{{\rm{T}}}}{{{{\boldsymbol{\varLambda}} }}^{{\mathbf{ - 1}}}}{{{{\boldsymbol{P}}}}_{{{k}}}} ,\;k \in \left\{ {{x_j},{y_j},{z_j}} \right\}. $

式中， ${{{\boldsymbol{\varLambda}} }}$为协方差矩阵 ${{{\boldsymbol{C}}}}$的特征值所构成的对角矩阵。

令目标函数 ${{{{f}}}_{2}} = {{{{D}}}_{{{\boldsymbol{m}}}}}\left( {{{{{\boldsymbol{P}}}}_{{{\boldsymbol{k}}}}}} \right)$，则有：

(14) $ {{{{f}}}_{2}} = \sum\nolimits_{i = 1}^M {\frac{{\lambda _i^2}}{{{{{{u}}}_{{{i}}}}}}} . $

式中： ${\lambda _i}$为式（5）中的形状参数， ${{{{u}}}_{{{i}}}}$为 ${{{\boldsymbol{\varLambda}} }}$中对角线上的元素，设目标函数 ${{{f}}}$为

(15) $ {{{f}}} = \mu {{{{f}}}_{1}} + \left( {1 - \mu } \right){{{{f}}}_{2}} . $

式中： $\;\mu $为放宽目标函数对点分布模型的约束参数. 将式（12）、（14）代入式（15），得到目标函数 ${{{f}}}$的最终表达式为

(16) $ \begin{split} {{{f}}} =& \mu \begin{gathered} \left\{ \sum\nolimits_{i = 1}^D {\sum\nolimits_{j = 1}^K {\left\| {{{{{\boldsymbol{d}}}}_{{{i}}}} - } \right.} } {\left. {\left( {{\boldsymbol{R}}\left( {{{\bar {\boldsymbol{S}}}}\left( {{{j}}} \right) + \sum\nolimits_{i = 1}^M {{\lambda _i}{{{{\boldsymbol{p}}}}_{{{i}}}}\left( {{{j}}} \right)} } \right) + {\boldsymbol{T}}} \right)} \right\|^2} \right\} \hfill \\ \end{gathered} + \hfill \\ & \left( {1 - \mu } \right)\sum\nolimits_{i = 1}^M {\frac{{\lambda _i^2}}{{{{{{u}}}_{{{i}}}}}}} {\kern 1pt} \;. \\[-10pt] \end{split} $

目标函数 ${{{{f}}}}$是非线性的，求解比较困难，因此将其转化为线性方程组来进行求解. 在求解模型的未知形状参数 ${\lambda _i}$时，将目标函数 ${{{{f}}}}$相对于 ${\lambda _i}$进行微分，并使等式为 $0$，则有：

(17) $ \begin{split} \frac{{\partial {{{f}}}}}{{\partial {\lambda _i}}} =&\mu \frac{\partial }{{\partial {\lambda _i}}} \begin{gathered} \left\{ \sum\nolimits_{i = 1}^D {\sum\nolimits_{j = 1}^K {\left\| {{{{{\boldsymbol{d}}}}_{{{\boldsymbol{i}}}}} - } \right.} } {\left. {\left( {{\boldsymbol{R}}\left( {{{\bar {\boldsymbol{S}}}}\left( {{{\boldsymbol{j}}}} \right) + \sum\nolimits_{i = 1}^M {{\lambda _i}{{{{\boldsymbol{p}}}}_{{{\boldsymbol{i}}}}}\left( {{{\boldsymbol{j}}}} \right)} } \right) + } \right.} \right.}\right. \end{gathered} \hfill \\ & {\left.{\boldsymbol{T}} \right)} {\left. \right\|^2} \left. \right\} +\left( {1 - \mu } \right)\sum\nolimits_{i = 1}^M {\frac{{2{\lambda _i}}}{{{{{{u}}}_{{{i}}}}}}} = 0\;. \hfill \\[-10pt] \end{split} $

结果可以由线性方程组的形式表示为

(18) $ {{{\boldsymbol{AX}}}} = {{{\boldsymbol{B}}}} , $

(19) $\begin{split}{ {{{\boldsymbol{A}}}} = \left( {\begin{array}{*{20}{c}} {{{{{\boldsymbol{p}}}}_{1}}\left( {{{{j}}}} \right){{{{\boldsymbol{p}}}}_{1}}\left( {{{{j}}}} \right) + \dfrac{{1 - \mu }}{\mu }\dfrac{1}{{{{{{u}}}_{1}}}}}&\cdots&{ {{{{\boldsymbol{p}}}}_{{{{M}}}{\text{ - }}{1}}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{1}}\left( {{{j}}} \right)}&{{{{{\boldsymbol{p}}}}_{{{M}}}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{1}}\left( {{{j}}} \right)} \\ {{{{\boldsymbol{p}}}}_{1}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{{{{M} - 1}}}}\left( {{{j}}} \right) &\cdots&{{{{\boldsymbol{p}}}}_{{{{M} - 1}}}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{{{{M} - 1}}}}\left( {{{j}}} \right) + \dfrac{{1 - \mu }}{\mu }\frac{1}{{{{{{u}}}_{{{{M} - 1}}}}}} &{{{{\boldsymbol{p}}}}_{{{M}}}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{{{{M} - 1}}}}\left( {{{j}}} \right) \hfill \\ {{{{\boldsymbol{p}}}}_{1}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{f{{M}}}}\left( {{{j}}} \right) \hfill &\cdots& {{{{\boldsymbol{p}}}}_{{{{M} - 1}}}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{{{M}}}}\left( {{{j}}} \right) &{{{{\boldsymbol{p}}}}_{{{M}}}}\left( {{{j}}} \right){{{{\boldsymbol{p}}}}_{{{M}}}}\left( {{{j}}} \right) + \dfrac{{1 - \mu }}{\mu }\dfrac{1}{{{{{{u}}}_{{{M}}}}}} \hfill \end{array}} \right) }_{3\times M}, \end{split}$

(20) $ \begin{split} \\ {{{\boldsymbol{B}}}} = \sum\nolimits_{i = 1}^D {\left[ \begin{gathered} \left( {{{{{\boldsymbol{d}}}}_{{{i}}}} - {{\bar {\boldsymbol{S}}}}\left( {{{j}}} \right)} \right){{{{\boldsymbol{p}}}}_{{1}}}\left( {{{j}}} \right) \hfill \\ \left( {{{{{\boldsymbol{d}}}}_{{{i}}}} - {{\bar {\boldsymbol{S}}}}\left( {{{j}}} \right)} \right){{{{\boldsymbol{p}}}}_{{{{M} - 1}}}}\left( {{{j}}} \right) \hfill \\ \left( {{{{{\boldsymbol{d}}}}_{{{i}}}} - {{\bar {\boldsymbol{S}}}}\left( {{{j}}} \right)} \right){{{{\boldsymbol{p}}}}_{{{M}}}}\left( {{{j}}} \right) \hfill \\ \end{gathered} \right]} , \end{split} $

(21) $ {{{\boldsymbol{X}}}} = {\left[ {{\lambda _1}, \cdots, {\lambda _{M - 1}},{\lambda _M}} \right]^{\text{T}}} . $

通过正交分解即可求得未知参数 $\left\{ {{\lambda _1}, \cdots, {\lambda _{M - 1}},{\lambda _M}} \right\}$. 由此可见，引入Mahalanobis距离作为附加项，并将非线性方程组转化为线性方程组求解的方法，能够估计形状参数 ${\lambda _i}$，同时计算出 ${{{{\boldsymbol{d}}}}_{{{i}}}}$和 ${{{\boldsymbol{m}}}}\left( {{{j}}} \right)$的最短距离，实现测试样本与点分布模型之间的拟合.

3DMFB-PDM的建模与拟合过程如图1所示， ${{N}_{{\rm{S}},1}}$表示第一个测试样本， $h = 1, \cdots ,{H}$为测试样本的个数， ${{\boldsymbol{S}}_{{\rm{P}},1}}\left( {{\lambda _i}} \right)$为与 $ {{N}_{{\rm{S}},1}}$进行拟合的点分布模型， ${\rm{Sim}}\left\langle {{{\boldsymbol{S}}_{{\rm{P}},1}},{{\textit{N}}_{{\rm{S}},1}}} \right\rangle$为点分布模型与 $ {{N}_{{\rm{S}},1}} $进行相似性度量. 由图1可知，在构建点分布模型时，首先通过普氏分析将训练样本集中的模板样本和目标样本对齐，然后在样本间建立对应关系，并根据已经建立的对应关系生成点分布模型. 在点分布模型与测试样本进行相似性度量时，应计算点分布模型与测试样本的距离，如果二者之间的距离最小，则进行拟合；反之，则引入附加项，调整点模型的参数 ${\lambda _i}$，使其与测试样本的距离最小，以达到最佳拟合效果.

用30组骨盆的CT扫描图像作为原始输入图像，实验验证3DMFB-PDM的通用性，CT分辨率在平面上为0.9 mm，在切片之间为1.5 mm.将骨盆的左侧股骨区域作为感兴趣区域，通过对骨盆的CT扫描图像进行数据处理，得到30组左侧股骨的3D模型并将其作为训练样本集，将30组左侧股骨的3D模型用CloudCompare^[22]在一个窗口中显示，结果如图2所示.

给定含有30组左侧股骨的训练样本集 ${{{\boldsymbol{\varOmega}} }} = \left[ {{{{{\boldsymbol{S}}}}_{1}},{{{{\boldsymbol{S}}}}_{2}}, \cdots ,{{{{\boldsymbol{S}}}}_{{{N}}}}} \right]$， $N = 30$. 训练样本集的人体测量学信息如表1所示，其中W、H、F、V分别为训练样本的体重、身高、三角形表面个数、三角形顶点个数. 任意选择训练样本集 ${{{\boldsymbol{\varOmega}} }}$中的3D模型（如011），将其作为模板样本 ${{{{\boldsymbol{S}}}}_{{{\rm{t}}}}}$并可视化，结果如图3所示. 在得到 ${{{{\boldsymbol{S}}}}_{{{\rm{t}}}}}$的可视化结果后用点云将其表示出来，并记录顶点的个数，如图4所示. 对图4分析可知，其顶点个数为57 224，将 ${{{\boldsymbol{\varOmega}} }}$中的所有股骨图像都用点云表示后的股骨图像的顶点个数为4 432~78 176. 将 ${{{{\boldsymbol{S}}}}_{{{{\rm{t}}}}}}$、 ${\mathbf{\Omega }}$中的剩余样本对齐，并通过自由变形（free-form deformation，FFD）^[23]的方式在训练样本间建立对应关系并生成训练样本的点分布模型，如图5所示. 当得到训练样本的点分布模型后，计算点分布模型与测试样本的距离，并根据结果调整点分布模型的参数，使点分布模型与测试样本有效地拟合. 当30组左侧股骨的点分布模型与测试样本进行拟合时，首先将30组左侧股骨 $\left\{ {{{{{\boldsymbol{S}}}}_{1}}, \cdots {{{{\boldsymbol{S}}}}_{{11}}}, \cdots ,{{{{\boldsymbol{S}}}}_{{30}}}} \right\}$用点云的形式表示出来，得到 $\left\{ {{{{\boldsymbol{S}}}}_{1}^{\mathbf{'}}, \cdots {{{\boldsymbol{S}}}}_{{11}}^{\mathbf{'}}, \cdots ,{{{\boldsymbol{S}}}}_{{30}}^{\mathbf{'}}} \right\}$，然后由 $\left\{ {{{{\boldsymbol{S}}}}_{1}^{\mathbf{'}}, \cdots {{{\boldsymbol{S}}}}_{{11}}^{\mathbf{'}}, \cdots ,{{{\boldsymbol{S}}}}_{{30}}^{\mathbf{'}}} \right\}$生成点分布模型 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$，并将 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$与测试样本拟合，通过调整 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$的形状参数使 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$与测试样本的距离最小，完成拟合.

由式（7）、（8）可知，特征向量的取值对应变形模型的取值，并且点分布模型可以由前 $G$个变形模型精确近似. 对于30组左侧股骨，计算前20个特征向量. 发现由前7个特征向量可以表征98.5%的变形模型，随着特征向量的继续增加，特征向量表征变形模型的百分比增加幅度很小，因此可以用前7个变形模型来近似变形模型. 图6为用部分变形模型表征的变形模型过程. 图中，EV为特征向量的数量，SM为累加的变形模型占整个变形模型的百分比. 由图可以看出，前7个变形模型可以近似表征整个变形模型.

在判断距离是否达到最小时，根据图1中的相似度量函数判断

(22) $ {\rm{Si}}{{\rm{m}}_{\rm{S}}}\left( {{\lambda _i}} \right) = \sum\limits_{h = 1}^H {{\rm{Sim}}\left\langle {{{\boldsymbol{S}}_{{\rm{P}},h}}\left( {{\lambda _i}} \right),{N_{{\rm{S}},h}}} \right\rangle } . $

给 ${\rm{Si}}{{\rm{m}}_{\rm{S}}}$设定阈值. 阈值的设定方式主要有2种：1）非经验性的阈值设定；2）根据不同的实验数据，反复尝试并观察实验结果，根据经验设定阈值. 虽然根据经验确定阈值存在误差，但相较方式1），此方式的计算量小且计算过程也相对简单，因此应用广泛，邹涛等^[24]通过经验确定文本分类模型的阈值参数，实现中文文档自动分类系统；Guha等^[25]通过经验设定相似性阈值，将输入数据点作为图的连通分量，极大程度地减少了算法运行的时间. 本研究阈值的设定根据经验确定，通过多次实验并观察实验结果是否最优，设定合理的阈值. 阈值大小的设定也会对实验结果产生影响. 通常来说，阈值大小的设定会影响点分布模型与测试样本的拟合误差以及形状参数 ${\lambda _i}$的迭代次数. 阈值设定过大，点分布模型与测试样本之间的距离越大，使二者间的拟合误差越大；阈值设定过小， ${\lambda _i}$的迭代次数越多，点分布模型与测试样本拟合所需的时间越长.

设定 ${\rm{Si}}{{\rm{m}}_{\rm{S}}}$的阈值后，计算点分布模型与测试样本之间的距离，当二者间的距离大于预先定义的阈值，说明间距还没有达到最小，引入附加项，对形状参数 ${\lambda _i}$进行迭代，这时形状参数由 ${\lambda _i}$变为 ${\lambda _{i + 1}}$. 在引入附加项后，点分布模型与测试样本间的拟合误差减小，但是不一定达到阈值，因此继续迭代点分布模型的形状参数，直到二者间的拟合误差小于阈值的时候，迭代停止，将点分布模型与测试样本进行拟合，输出拟合结果.

为了评估 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$与测试样本的拟合误差，比较3DMFB-PDM在未引入附加项时的拟合误差ED、附加项后的拟合误差（ED+MD）. 用最大误差Max、平均误差Mean和均方根误差RMS作为评价指标，结果如表2所示. 由表可知，引入附加项后的拟合误差小于未引入附加项的拟合误差. 表2为 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$与单个测试样本拟合的结果，为了进一步验证3DMFB-PDM的通用性，选取多个测试样本与 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$拟合，结果如表3所示. 由表可知，3DMFB-PDM在多个测试样本进行拟合时仍然能够满足要求.

为了进一步分析3DMFB-PDM的效率，将30组左侧股骨作为训练样本集，并对比分析3DMFB-PDM与MDL、SLIDE、文献[15]、文献[18]、文献[19]在构建左侧股骨的点分布模型所需的时间，结果如表4所示. 表中，t为每种方法运行的时间. 由表可知，3DMFB-PDM构建左侧股骨的点分布模型所需的时间最短. 文献[18]的方法略长于3DMFB-PDM；SLIDE通过滑动特征点，不断计算每个特征点的切线方向，因此耗时最长. 由此可见，相比于其他方法，3DMFB-PDM在构建点分布模型所需的时间上具有一定的优势.

为了比较 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$与测试样本在不同方法下的拟合误差，将文献[15]、18]、[19]中的方法与3DMFB-PDM进行对比，结果如表5所示. 文献[15]、[19]的最大拟合误差分别为15.06 、13.68 mm，大于3DMFB-PDM的12.27 mm，表明3DMFB-PDM具有较小的拟合误差. 3DMFB-PDM的平均误差和均方根误差分别为2.42 、3.11 mm，与文献[15]、[19]的相比也明显较小，表明3DMFB-PDM的拟合误差明显小于文献[15]、[19]的. 此外，文献[18]的最大拟合误差、平均误差和均方根误差分别为12.37、2.48、3.23 mm，小于文献[15]、[19]的，但略大于3DMFB-PDM的，这是由于文献[18]方法在计算拟合误差时采用欧氏距离计算点分布模型与测试样本间的最小距离. 由此可见，3DMFB-PDM的拟合误差在4种方法中最小，表明3DMFB-PDM能够有效地将 ${{{{\boldsymbol{S}}}}_{{{\rm{P}}}}}$与测试样本拟合. 将这4种方法在3个测试样本下的拟合误差进行对比，结果如表6所示. 当文献[15]、[19]方法与多个测试样本进行拟合时，均具有较大的最大误差、平均误差和均方根误差，因此它们的拟合误差较大，拟合效果较差；文献[18]方法与多个测试样本拟合时的误差较小，当3DMFB-PDM与多个测试样本拟合时的最大误差、平均误差和均方根误差均为4种方法中最小的，表明3DMFB-PDM具有最小的拟合误差. 因此，与文献[15]、[18]、[19]方法相比，3DMFB-PDM构建的点分布模型能够有效拟合测试样本.

为了进一步评估3DMFB-PDM构建的点分布模型的性能，将文献[15]、[18]、[19]方法与3DMFB-PDM进行对比验证，并将模型的紧凑性、特异性和通用性作为评价指标. 紧凑性^[26]（compactness）表示用尽可能少的主要变形模型表征整个点分布模型的变形能力，它被定义为对应于变形模型的特征值的累积. 特异性^[27]（specificity）表示模型生成的样本与训练样本集中的样本间的相似程度，它被定义为创建的形状样本与训练样本集中最相似的形状样本间的平均拟合误差. 通用性^[28]（generality）是模型生成未由训练样本集中明确提供的样本的能力，它被定义为遗漏训练样本集中的样本，用训练样本集中的剩余样本构建模型，计算构建的模型和遗漏的样本间的拟合误差. 通常来说，紧凑性、特异性和通用性的值越小，模型的性能越好. 4种方法构建点分布模型的性能对比结果如图7所示. 图中，DM为变形模型的个数，C、S、G分别为紧凑性、特异性和通用性表征的指标. 由图可知，文献[15]、[19]方法紧凑性明显差于文献[18]方法和3DMFB-PDM，文献[18]方法紧凑性则略差于3DMFB-PDM；3DMFB-PDM的特异性明显优于其他3种方法；文献[15]方法的通用性最差，文献[18]、[19]方法的通用性接近，3DMFB-PDM的通用性明显优于其他3种方法. 上述结果表明，3DMFB-PDM构建的点分布模型具有较好的性能.

将Mahalanobis距离作为附加项引入计算患者病患部位的3D解剖结构模型与测试样本的最小距离过程中，提出基于点分布模型的3D模型拟合方法. 该方法通过引入Mahalanobis距离，将非线性方程组转换为线性方程组求解，同时判断点分布模型与测试样本距离是否最小，调整点分布模型的形状参数，使患者病患部位的3D解剖结构模型与测试样本间的拟合误差最小，解决传统方法构建的患者病患部位3D解剖结构模型与测试样本拟合效果较差的问题. 实验结果表明，3DMFB-PDM能够有效地将3D模型与测试样本拟合. 在计算3D模型与测试样本的拟合误差时，本研究采用最大误差、平均误差、均方根误差和相似性系数衡量拟合程度，这在大程度上依赖于3D模型与测试样本的最小距离计算. 计算二者的最小距离的过程较为烦琐且耗时. 未来将会研究如何尽可能的简化计算3D模型与测试样本的最小距离的过程，在保证二者之间的距离最小的同时简化计算的复杂度.