患者的骨幾何模型是個體化生物力學建模與分析、手術規劃方案制定、手術導板設計、人工關節定制的重要基礎。本研究針對股骨提出一種基于統計形狀模型的三維快速重建方法，結合患者X線平片數據，實現了個體化患者股骨幾何的快速三維建模。模型三維重建的平均誤差在1.597~1.842 mm之間，均方根誤差在1.453~2.341 mm之間。重建模型的股骨頭直徑、頸干角、偏置距離和前傾角的平均誤差分別為0.597 mm、1.163°、1.389 mm和1.354°。相比傳統建模方法，新方法能夠在更短時間內較為準確地實現股骨三維快速重建。本文為骨幾何三維快速建模提供了新技術，有助于推動針對患者的快速生物力學分析，為骨科植入物的選擇與定制植入物假體的快速研發提供新的思路。
引用本文： 張志偉, 陳瑱賢, 張志峰, 王彩梅, 靳忠民. 基于統計形狀模型的股骨快速建模方法. 生物醫學工程學雜志, 2022, 39(5): 862-869. doi: 10.7507/1001-5515.202111005 復制
引言
在骨科臨床，骨性關節炎、骨折、骨缺損和骨腫瘤等骨、關節疾病一直困擾著患者，為了達到精準治療的目的，建立骨骼或關節病區的三維幾何模型，有利于幫助醫生觀察診斷、手術規劃、器械定制以及開展生物力學機制研究。骨幾何三維模型已在假體設計[1-2]、生物力學分析[3-4]、術前規劃[5]和臨床診斷[6]中得到了廣泛的應用。傳統獲得骨幾何三維模型主要是通過計算機斷層掃描（computed tomography，CT），然后基于Mimics等建模軟件進行提取分割建模。然而，CT掃描費用較高、輻射劑量大，傳統Mimics提取分割建模方法時間長，不利于骨三維幾何模型的快速獲取。因此，降低骨幾何影像學數據獲取成本以及縮短骨幾何三維重建時間是骨幾何三維建模的發展趨勢。
針對現有獲得骨骼三維模型方法存在的缺陷，研究人員提出利用二維醫學影像數據（如X線平片和超聲圖像）獲得患者特定區域解剖結構三維模型的方法[7]。在該方法中，統計形狀模型（statistical shape model，SSM）[8-9]對引導重建過程至關重要。前人根據患者的CT或核磁共振數據，分別針對顱骨[10]、肱骨[11]、盆骨[12]、股骨[13]、下肢骨骼[14-15]、膝關節[16-17]、髖關節[18]、踝關節[19]建立了統計形狀模型，并將其應用于骨修復、骨幾何重建和假體設計等領域。基于相應骨骼部位的統計形狀模型，患者不再需要采集CT或磁共振影像學數據，通過對患者X線平片中的解剖特征信息的識別，就可以快速重建出與之相匹配的骨三維幾何模型。Schumann等[20]通過建立股骨近端統計形狀模型，實現了利用X線圖像對股骨近端的三維重建。相比CT等建模方法，基于統計形狀模型的二維到三維快速建模的方法其成像成本與輻射劑量顯著減小[21]。近年來定制化髖關節假體的興起，人們已不滿足于股骨柄的髓內定制。多項研究發現，股骨柄的髓外設計對于髖關節力也具有明顯的影響[22-24]，因此股骨柄髓外幾何參數的定制和選擇也受到了人們的重視。所以，如果將基于統計形狀模型的二維到三維快速建模的方法應用于生物力學分析、髖關節假體選擇與參數化定制，將大大縮短時間成本。
數字重建射線影像（digitally reconstructed radiograph，DRR）[25-26]和特征先驗信息重建[27-28]是基于統計形狀模型實現二維到三維重建的兩種主要方案。DRR是利用三維模型與二維影像之間的映射關系，通過不斷迭代確定統計形狀模型形變域，使得統計形狀模型的二維投影與患者X線平片差別最小，進而以較小的誤差實現三維重建[29-31]。但該方案的迭代性質有可能使得結果陷入局部最優解，同時為了保持較高精度的三維重建，需要生成成百上千個DRR，導致耗時較長[7,32]。特征先驗信息重建是將統計形狀模型與患者二維醫學影像中的特征信息之間建立聯系，以如輪廓線、曲面片以及外部觀測信息等[33]特征信息作為先驗信息，使用多元回歸方法實現模型的三維重建。該方案的優勢在于可以利用少量訓練數據和稀疏的先驗信息實現較高精度的三維重建[34]。在滿足建模需求的情況下，這種方法將大大減少時間成本。然而，目前國內還沒有利用特征先驗信息重建開展基于統計形狀模型二維到三維快速建模的研究報道。
本文利用迭代最近點（iterative closest point，ICP）算法建立統計形狀模型，將統計形狀模型與患者的二維醫學影像特征信息之間建立聯系，研究個體化患者X線平片到三維骨幾何快速建模的方法，為個體化生物力學建模與分析、手術方案規劃和醫療器械的快速設計提供技術基礎。
1 骨幾何數據
本文從醫學圖像SMIR數據庫（
表1
訓練集分組情況
Table1.
Grouping of training sets
2 基于ICP的點云配準
點分布模型（point distribution model，PDM）[35]是將三維模型離散化，用點坐標的改變表示模型形狀的變化。基于股骨的點分布模型建立股骨統計形狀模型，其中最關鍵的環節是建立不同骨模型特征點之間的對應關系。如圖1a所示，由于股骨三維幾何模型是來自不同患者、不同設備、不同位姿情況下采集的數據，其三維模型并不處于同一坐標系下。未經過處理的股骨模型，模型之間相互交錯，這種情況無法建立模型點之間的對應關系。所以，在建立統計形狀模型之前，需要將所有樣本統一到同一坐標系下。
圖1
數據預處理
a. 未處理數據集；b. 三維模型離散化
Figure1. Data preprocessinga. unprocessed data set; b. three-dimensional model discretization
ICP算法是Besl等[36]提出的用于點分布模型配準的方法，其核心是尋找空間變換矩陣，使得兩個模型之間差異最小化，實現對不同坐標系下的模型進行配準。由于ICP算法配準精度高且操作簡單，成為了點分布模型配準的主流算法[37]，在治療骨缺損和骨幾何重建領域得到了廣泛的應用[38-39]。由于ICP算法是針對點分布模型，所以在配準之前需要將三維模型離散化，如圖1b所示。
從數據集中選取一個模型作為參考模型，利用ICP算法，將數據集中其他股骨模型轉換到參考模型坐標系下，如圖2所示為點云模型的配準過程及模型配準效果。
圖2
點云配準過程
Figure2.
Point cloud registration process
3 股骨統計形狀建模
從上述離散的股骨幾何數據集中選取一個點云作為參考點云，將其他點云統一到該點云坐標系下。使用K-D樹算法[39]建立所有點云點之間的對應關系，構建統計形狀模型的訓練數據點集，記為R
，n為數據集中包含的點云個數，每個點云中包含m個三維空間點
，則每個股骨形狀點云可由
維向量表示為：
 |
本文采用高斯過程形變模型（Gaussian process morphable models，GPMM）[40]來建立股骨統計形狀模型。首先確定參考形狀，設 
為參考形狀點云中所有點構成的集合，然后將其他形狀視為參考形狀的形變。設 
為參考形狀點云 
中一個點所在的空間坐標，
為點 
在X、Y、Z三個方向形變量的大小，則任意點云中與 
相對應點的位置為 
，如圖3所示。本文假設點云中點的形變量遵循多元正態分布，于是形變量 
表達式如下：
圖3
單點形變示意圖
Figure3.
Schematic diagram of single point deformation
 |
式中 
、
、
分別代表 
在X、Y、Z三個方向的形變量，
、
、
分別代表 
在X、Y、Z三個方向的形變量的均值函數，
、
、
分別代表 
在X、Y、Z三個方向的形變量的協方差函數。
根據上述內容，可以將該形變模型推廣到包含任意數量點的點云模型中。得到參考形狀點云 
到任意形狀點云的高斯形變過程表達式為：
 |
其中，
為數據集中點云形變域的均值函數，
為數據集中點云形變域的協方差函數。
高斯過程形變模型就是通過給定數據集確定形變的均值函數和協方差函數，建立多元正態分布模型確定形變域，進而根據形變域建立統計形狀模型。高斯過程形變模型認為對于任意給定形狀均可以通過參考形狀點云 
形變得到，形變域為 
，新得到模型中的點與參考形狀點云中的點必須一一對應，所以任意形狀點云q可以表示為：
 '/> |
其中，V為參考形狀點云中所有點形變量構成的形變域，
為參考形狀點云。
在訓練集中，任意形狀點云 
到參考形狀點云 
中每個點形變量所構成的形變域 
，可以通過式(5)得到：
 '/> |
其中
代表數據集中點云數量。
根據上式計算得到數據集中每個點云相對于參考形狀點云的形變域，進而可以得到數據集中點云相對于參考形狀點云形變域的均值函數和協方差函數，其計算過程如式(6)～(7)所示：
 |
 |
其中，
為數據集中任意點云形變域。
當點云中點個數過多時會導致計算時間過長，所以本文采用主成分分析法（principal component analysis，PCA）[40]對數據進行降維，以減少運算量。經過降維處理后的高斯形變過程模型簡化為：
 |
其中k表示主要模式個數，
和 
分別表示協方差矩陣由大到小排列的第i個特征值和特征向量。
服從標準正態分布，決定形變域的變化范圍。
4 二維到三維的骨幾何快速重建
基于所建立的股骨統計形狀模型，本文利用高斯過程回歸（Gaussian process regression，GPR）[35]統計決策模型建立了二維到三維的骨幾何快速重建方法，可以利用1張正位X線平片快速重建完整股骨個體化三維模型。高斯過程回歸的核心是建立數據集中多個變量之間的聯合概率分布。
針對股骨三維重建，設 
為數據集中所有點云形變域構成的集合，其中n為數據集中點云的數量，
為與數據集中點云對應X線平片中所能夠觀察到信息的集合。其中 
表示點云中每一點沿X、Y、Z三個方向形變量所構成的形變域，m表示點云中所包含點的數目。
為每張二維醫學影像中所包含的信息，
為二維醫學影像中所包含信息的個數。假設 
和 
服從正態分布，則 
與 
的聯合概率分布亦為正態分布，其表達式如(9)所示[41]：
 |
其中 
為點云形變域的均值函數；
為觀察值的均值函數。
在得到多元聯合正態分布之后，取
的邊緣分布，根據貝葉斯推斷，建立后驗概率模型如式(10)所示：
 |
其中
為均值函數，
為協方差函數，其計算過過程如式(11)～(12)所示：
 |
 |
其中 
、
、
和 
通過式(13)～(16)計算得到：
 |
 |
 |
 |
綜上所述，本文定義了高斯過程回歸模型：
。這是一種通過控制參數，利用統計推斷確定形變域的模型，但是由于在實際觀測過程中可能會遇到噪音干擾，根據高斯過程回歸定義，通常要在計算過程中添加一個噪音項，因此高斯過程回歸均值函數 
和協方差函數 
轉變為以下形式：
 |
 |
其中 
服從標準正態分布，I為單位矩陣。
通過上述方法得到高斯過程回歸模型，該模型可以利用稀疏信息重建復雜三維模型，因此我們可以利用從二維醫學影像中獲得的數據，確定點云形變域，進而完成股骨三維模型的重建。
基于上述方法，本文將患者正位X線片中頸干角、內外側髁之間寬度（w）、股骨頭直徑和股骨頭沿股骨解剖軸與膝關節內側髁之間的最大距離（L）作為控制參數（如圖4所示），判斷L屬于哪個長度區間，將數據代入到相應長度區間所建立的統計形狀模型中進行計算，確定股骨點云形變域，進而重建股骨三維模型。
圖4
X線片觀察信息
通過正位X線片所觀察到5例測試集樣本的L、w、頸干角和股骨頭直徑等信息如表2所示。
表2
患者X線片中所包含的信息
Table2.
Information contained in patients’ x-rays
5 模型驗證及誤差分析
將上述測試集骨骼數據代入到建好的股骨統計形狀模型中，得到三維重建股骨模型。之后將重建股骨模型與測試集股骨進行對比，結果如圖5所示，從左到右依次為編號1到5測試集模型的誤差云圖。通過觀察誤差云圖發現，重建模型誤差較大的區域主要集中在股骨近端和股骨遠端，可能是由于這兩個區域形狀復雜。同時由正態分布規律可知，樣本分布較為集中且預測值在均值附近時，統計模型預測精度較高，反之預測精度較低，所以統計形狀模型在形狀差異較小的區域有較好的預測精度，在差異較大區域預測精度則較低，而不同個體在股骨近端和股骨遠端會存在較大差異，所以這也是導致統計形狀模型對股骨近端和股骨遠端預測誤差較大的原因。
圖5
重建模型誤差云圖（mm）
Figure5.
Reconstruction model error (mm)
在本研究中，將45個股骨模型離散化后作為訓練集，并根據股骨長度進行分組，建立股骨的統計形狀模型。根據統計學原理，當患者股骨的特征參數與訓練集中模型的特征參數差別過大時就會導致重建模型在該區域有較大誤差，并且當組內樣本較少時，重建精度較差，所以要使用盡可能多的樣本作為訓練集以提高重建精度。
本文使用平均誤差、均方根誤差和最大誤差作為量化指標，對比分析了重建模型與樣本之間的誤差，如表3所示。重建模型的平均誤差范圍在1.597～1.842 mm之間，均方根誤差范圍在1.453～2.341 mm之間，最大誤差范圍在5.830～10.761 mm之間。其中5號模型的平均誤差、均方根誤差和最大誤差均是五個預測模型中最大的，分析其原因可能是股骨頭沿股骨解剖軸與膝關節內側髁之間的最大距離在373～413 mm內的訓練樣本過少所導致的。總體上本文所建立的股骨統計形狀模型的重建誤差均小于2 mm。而且，在輸入患者特征數據之后，模型重建時間可控制在5 s之內，極大地縮短了模型重建時長。
表3
統計形狀模型誤差
Table3.
Statistical shape model error
然而由圖5可知，雖然重建模型平均誤差在2 mm以內，但是在股骨近端和遠端會有較大的誤差，沒有辦法直接應用到臨床當中。所以在得到重建模型之后，我們測量了股骨近端的幾何參數，包括股骨頭直徑、頸干角、偏置距離和前傾角。這些數據是選擇股骨柄和確定股骨柄近端形態的重要指標。我們將重建模型的幾何數據與原模型的數據進行對比，計算得到其偏差的絕對值，如表4所示。股骨近端幾何參數對于髖關節假體選擇有著重要的指導意義。再加上近年來定制化髖關節假體的興起，股骨近端幾何參數為髖關節股骨柄髓外部分定制提供了設計依據。Schumann等[20]的重建模型主要對比了股骨近端幾何參數的重建精度，其重建模型的偏置距離和頸干角的平均誤差分別為1.56 mm和2.21°。而根據本文方法重建得到的三維模型，偏置距離和頸干角的平均誤差分別為1.389 mm和1.163°。除股骨頭直徑外本文重建模型的近端幾何參數精度都有所提高，并且股骨頭直徑誤差也保持在較低水平。根據表4的計算結果，股骨頭直徑、頸干角、偏置距離和前傾角的重建精度都保持在較高水平，說明重建模型可以滿足髖關節假體的定制化使用需求。
表4
重建模型部分幾何參數誤差值
Table4.
The error values of some geometric parameters of the reconstructed model
之前有關股骨統計形狀模型的研究主要集中在重建股骨的局部區域，并且沒有考慮患者之間的差異性，無法實現針對患者的個性化重建，顯然是不太合理的。本文基于高斯過程形變以及高斯過程回歸方法，將患者股骨頭直徑、頸干角以及股骨長度等信息作為三維重建的參考依據，實現了在較短時間內針對患者股骨的個性化重建。根據建模結果誤差分析可知，較少的訓練樣本數量會造成較大的建模誤差。所以后續研究中可以通過增加訓練集樣本數量，提高模型的重建精度。
本文還存在一定的局限性。首先由于本文采用的數據集來源于公開數據庫，其中并沒有每位志愿者的個人信息，所以沒有考慮年齡和性別對于模型重建的影響。其次本文的重建模型精度還有待提高，尤其在股骨近端和遠端，在后續研究中可以收集更多的數據來提高模型重建精度。同時，本文也具有一定的優勢。首先本文利用較少的外部數據實現了股骨三維模型的重構，為骨骼模型的三維重建提供了一個新的思路。其次本文考慮到患者之間的差異性，實現了利用不同患者數據，重構股骨三維模型，為髖關節假體的選擇和股骨柄髓外幾何參數的定制提供了新的方法。
6 結論
本文基于統計形狀模型，建立了利用患者二維醫學影像數據實現股骨快速三維重建的方法。該方法可以根據患者X線平片中包含的參數信息，5 s之內實現對患者股骨較高精度的三維重建，模型表面平均誤差在2 mm以內，重建模型的股骨頭直徑、頸干角、偏置距離和前傾角的平均誤差分別為0.597 mm、1.163 °、1.389 mm和1.354 °，極大地節省了獲得三維骨幾何模型所需要的時間。本文建立的骨幾何快速三維重建方法可以應用到骨科手術規劃、髖關節假體選擇及生物力學研究當中，以縮短建模時間，降低建模成本，快速獲得所需要的骨幾何三維模型。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：張志偉和陳瑱賢構思并設計了研究方案，張志偉實施方案并處理數據，張志偉和張志峰撰寫論文，王彩梅、陳瑱賢和靳忠民討論并修改了論文。
引言
在骨科臨床，骨性關節炎、骨折、骨缺損和骨腫瘤等骨、關節疾病一直困擾著患者，為了達到精準治療的目的，建立骨骼或關節病區的三維幾何模型，有利于幫助醫生觀察診斷、手術規劃、器械定制以及開展生物力學機制研究。骨幾何三維模型已在假體設計[1-2]、生物力學分析[3-4]、術前規劃[5]和臨床診斷[6]中得到了廣泛的應用。傳統獲得骨幾何三維模型主要是通過計算機斷層掃描（computed tomography，CT），然后基于Mimics等建模軟件進行提取分割建模。然而，CT掃描費用較高、輻射劑量大，傳統Mimics提取分割建模方法時間長，不利于骨三維幾何模型的快速獲取。因此，降低骨幾何影像學數據獲取成本以及縮短骨幾何三維重建時間是骨幾何三維建模的發展趨勢。
針對現有獲得骨骼三維模型方法存在的缺陷，研究人員提出利用二維醫學影像數據（如X線平片和超聲圖像）獲得患者特定區域解剖結構三維模型的方法[7]。在該方法中，統計形狀模型（statistical shape model，SSM）[8-9]對引導重建過程至關重要。前人根據患者的CT或核磁共振數據，分別針對顱骨[10]、肱骨[11]、盆骨[12]、股骨[13]、下肢骨骼[14-15]、膝關節[16-17]、髖關節[18]、踝關節[19]建立了統計形狀模型，并將其應用于骨修復、骨幾何重建和假體設計等領域。基于相應骨骼部位的統計形狀模型，患者不再需要采集CT或磁共振影像學數據，通過對患者X線平片中的解剖特征信息的識別，就可以快速重建出與之相匹配的骨三維幾何模型。Schumann等[20]通過建立股骨近端統計形狀模型，實現了利用X線圖像對股骨近端的三維重建。相比CT等建模方法，基于統計形狀模型的二維到三維快速建模的方法其成像成本與輻射劑量顯著減小[21]。近年來定制化髖關節假體的興起，人們已不滿足于股骨柄的髓內定制。多項研究發現，股骨柄的髓外設計對于髖關節力也具有明顯的影響[22-24]，因此股骨柄髓外幾何參數的定制和選擇也受到了人們的重視。所以，如果將基于統計形狀模型的二維到三維快速建模的方法應用于生物力學分析、髖關節假體選擇與參數化定制，將大大縮短時間成本。
數字重建射線影像（digitally reconstructed radiograph，DRR）[25-26]和特征先驗信息重建[27-28]是基于統計形狀模型實現二維到三維重建的兩種主要方案。DRR是利用三維模型與二維影像之間的映射關系，通過不斷迭代確定統計形狀模型形變域，使得統計形狀模型的二維投影與患者X線平片差別最小，進而以較小的誤差實現三維重建[29-31]。但該方案的迭代性質有可能使得結果陷入局部最優解，同時為了保持較高精度的三維重建，需要生成成百上千個DRR，導致耗時較長[7,32]。特征先驗信息重建是將統計形狀模型與患者二維醫學影像中的特征信息之間建立聯系，以如輪廓線、曲面片以及外部觀測信息等[33]特征信息作為先驗信息，使用多元回歸方法實現模型的三維重建。該方案的優勢在于可以利用少量訓練數據和稀疏的先驗信息實現較高精度的三維重建[34]。在滿足建模需求的情況下，這種方法將大大減少時間成本。然而，目前國內還沒有利用特征先驗信息重建開展基于統計形狀模型二維到三維快速建模的研究報道。
本文利用迭代最近點（iterative closest point，ICP）算法建立統計形狀模型，將統計形狀模型與患者的二維醫學影像特征信息之間建立聯系，研究個體化患者X線平片到三維骨幾何快速建模的方法，為個體化生物力學建模與分析、手術方案規劃和醫療器械的快速設計提供技術基礎。
1 骨幾何數據
本文從醫學圖像SMIR數據庫（
表1
訓練集分組情況
Table1.
Grouping of training sets
2 基于ICP的點云配準
點分布模型（point distribution model，PDM）[35]是將三維模型離散化，用點坐標的改變表示模型形狀的變化。基于股骨的點分布模型建立股骨統計形狀模型，其中最關鍵的環節是建立不同骨模型特征點之間的對應關系。如圖1a所示，由于股骨三維幾何模型是來自不同患者、不同設備、不同位姿情況下采集的數據，其三維模型并不處于同一坐標系下。未經過處理的股骨模型，模型之間相互交錯，這種情況無法建立模型點之間的對應關系。所以，在建立統計形狀模型之前，需要將所有樣本統一到同一坐標系下。
圖1
數據預處理
a. 未處理數據集；b. 三維模型離散化
Figure1. Data preprocessinga. unprocessed data set; b. three-dimensional model discretization
ICP算法是Besl等[36]提出的用于點分布模型配準的方法，其核心是尋找空間變換矩陣，使得兩個模型之間差異最小化，實現對不同坐標系下的模型進行配準。由于ICP算法配準精度高且操作簡單，成為了點分布模型配準的主流算法[37]，在治療骨缺損和骨幾何重建領域得到了廣泛的應用[38-39]。由于ICP算法是針對點分布模型，所以在配準之前需要將三維模型離散化，如圖1b所示。
從數據集中選取一個模型作為參考模型，利用ICP算法，將數據集中其他股骨模型轉換到參考模型坐標系下，如圖2所示為點云模型的配準過程及模型配準效果。
圖2
點云配準過程
Figure2.
Point cloud registration process
3 股骨統計形狀建模
從上述離散的股骨幾何數據集中選取一個點云作為參考點云，將其他點云統一到該點云坐標系下。使用K-D樹算法[39]建立所有點云點之間的對應關系，構建統計形狀模型的訓練數據點集，記為R，n為數據集中包含的點云個數，每個點云中包含m個三維空間點，則每個股骨形狀點云可由維向量表示為：
|
本文采用高斯過程形變模型（Gaussian process morphable models，GPMM）[40]來建立股骨統計形狀模型。首先確定參考形狀，設 為參考形狀點云中所有點構成的集合，然后將其他形狀視為參考形狀的形變。設 為參考形狀點云 中一個點所在的空間坐標， 為點 在X、Y、Z三個方向形變量的大小，則任意點云中與 相對應點的位置為 ，如圖3所示。本文假設點云中點的形變量遵循多元正態分布，于是形變量 表達式如下：
圖3
單點形變示意圖
Figure3.
Schematic diagram of single point deformation
|
式中 、、 分別代表 在X、Y、Z三個方向的形變量，、、 分別代表 在X、Y、Z三個方向的形變量的均值函數，、、 分別代表 在X、Y、Z三個方向的形變量的協方差函數。
根據上述內容，可以將該形變模型推廣到包含任意數量點的點云模型中。得到參考形狀點云 到任意形狀點云的高斯形變過程表達式為：
|
其中，為數據集中點云形變域的均值函數，為數據集中點云形變域的協方差函數。
高斯過程形變模型就是通過給定數據集確定形變的均值函數和協方差函數，建立多元正態分布模型確定形變域，進而根據形變域建立統計形狀模型。高斯過程形變模型認為對于任意給定形狀均可以通過參考形狀點云 形變得到，形變域為 ，新得到模型中的點與參考形狀點云中的點必須一一對應，所以任意形狀點云q可以表示為：
| '/> |
其中，V為參考形狀點云中所有點形變量構成的形變域， 為參考形狀點云。
在訓練集中，任意形狀點云 到參考形狀點云 中每個點形變量所構成的形變域 ，可以通過式(5)得到：
| '/> |
其中代表數據集中點云數量。
根據上式計算得到數據集中每個點云相對于參考形狀點云的形變域，進而可以得到數據集中點云相對于參考形狀點云形變域的均值函數和協方差函數，其計算過程如式(6)～(7)所示：
|
|
其中，為數據集中任意點云形變域。
當點云中點個數過多時會導致計算時間過長，所以本文采用主成分分析法（principal component analysis，PCA）[40]對數據進行降維，以減少運算量。經過降維處理后的高斯形變過程模型簡化為：
|
其中k表示主要模式個數， 和 分別表示協方差矩陣由大到小排列的第i個特征值和特征向量。 服從標準正態分布，決定形變域的變化范圍。
4 二維到三維的骨幾何快速重建
基于所建立的股骨統計形狀模型，本文利用高斯過程回歸（Gaussian process regression，GPR）[35]統計決策模型建立了二維到三維的骨幾何快速重建方法，可以利用1張正位X線平片快速重建完整股骨個體化三維模型。高斯過程回歸的核心是建立數據集中多個變量之間的聯合概率分布。
針對股骨三維重建，設 為數據集中所有點云形變域構成的集合，其中n為數據集中點云的數量， 為與數據集中點云對應X線平片中所能夠觀察到信息的集合。其中 表示點云中每一點沿X、Y、Z三個方向形變量所構成的形變域，m表示點云中所包含點的數目。 為每張二維醫學影像中所包含的信息， 為二維醫學影像中所包含信息的個數。假設 和 服從正態分布，則 與 的聯合概率分布亦為正態分布，其表達式如(9)所示[41]：
|
其中 為點云形變域的均值函數； 為觀察值的均值函數。
在得到多元聯合正態分布之后，取的邊緣分布，根據貝葉斯推斷，建立后驗概率模型如式(10)所示：
|
其中為均值函數，為協方差函數，其計算過過程如式(11)～(12)所示：
|
|
其中 、、 和 通過式(13)～(16)計算得到：
|
|
|
|
綜上所述，本文定義了高斯過程回歸模型：。這是一種通過控制參數，利用統計推斷確定形變域的模型，但是由于在實際觀測過程中可能會遇到噪音干擾，根據高斯過程回歸定義，通常要在計算過程中添加一個噪音項，因此高斯過程回歸均值函數 和協方差函數 轉變為以下形式：
|
|
其中 服從標準正態分布，I為單位矩陣。
通過上述方法得到高斯過程回歸模型，該模型可以利用稀疏信息重建復雜三維模型，因此我們可以利用從二維醫學影像中獲得的數據，確定點云形變域，進而完成股骨三維模型的重建。
基于上述方法，本文將患者正位X線片中頸干角、內外側髁之間寬度（w）、股骨頭直徑和股骨頭沿股骨解剖軸與膝關節內側髁之間的最大距離（L）作為控制參數（如圖4所示），判斷L屬于哪個長度區間，將數據代入到相應長度區間所建立的統計形狀模型中進行計算，確定股骨點云形變域，進而重建股骨三維模型。
圖4
X線片觀察信息
通過正位X線片所觀察到5例測試集樣本的L、w、頸干角和股骨頭直徑等信息如表2所示。
表2
患者X線片中所包含的信息
Table2.
Information contained in patients’ x-rays
5 模型驗證及誤差分析
將上述測試集骨骼數據代入到建好的股骨統計形狀模型中，得到三維重建股骨模型。之后將重建股骨模型與測試集股骨進行對比，結果如圖5所示，從左到右依次為編號1到5測試集模型的誤差云圖。通過觀察誤差云圖發現，重建模型誤差較大的區域主要集中在股骨近端和股骨遠端，可能是由于這兩個區域形狀復雜。同時由正態分布規律可知，樣本分布較為集中且預測值在均值附近時，統計模型預測精度較高，反之預測精度較低，所以統計形狀模型在形狀差異較小的區域有較好的預測精度，在差異較大區域預測精度則較低，而不同個體在股骨近端和股骨遠端會存在較大差異，所以這也是導致統計形狀模型對股骨近端和股骨遠端預測誤差較大的原因。
圖5
重建模型誤差云圖（mm）
Figure5.
Reconstruction model error (mm)
在本研究中，將45個股骨模型離散化后作為訓練集，并根據股骨長度進行分組，建立股骨的統計形狀模型。根據統計學原理，當患者股骨的特征參數與訓練集中模型的特征參數差別過大時就會導致重建模型在該區域有較大誤差，并且當組內樣本較少時，重建精度較差，所以要使用盡可能多的樣本作為訓練集以提高重建精度。
本文使用平均誤差、均方根誤差和最大誤差作為量化指標，對比分析了重建模型與樣本之間的誤差，如表3所示。重建模型的平均誤差范圍在1.597～1.842 mm之間，均方根誤差范圍在1.453～2.341 mm之間，最大誤差范圍在5.830～10.761 mm之間。其中5號模型的平均誤差、均方根誤差和最大誤差均是五個預測模型中最大的，分析其原因可能是股骨頭沿股骨解剖軸與膝關節內側髁之間的最大距離在373～413 mm內的訓練樣本過少所導致的。總體上本文所建立的股骨統計形狀模型的重建誤差均小于2 mm。而且，在輸入患者特征數據之后，模型重建時間可控制在5 s之內，極大地縮短了模型重建時長。
表3
統計形狀模型誤差
Table3.
Statistical shape model error
然而由圖5可知，雖然重建模型平均誤差在2 mm以內，但是在股骨近端和遠端會有較大的誤差，沒有辦法直接應用到臨床當中。所以在得到重建模型之后，我們測量了股骨近端的幾何參數，包括股骨頭直徑、頸干角、偏置距離和前傾角。這些數據是選擇股骨柄和確定股骨柄近端形態的重要指標。我們將重建模型的幾何數據與原模型的數據進行對比，計算得到其偏差的絕對值，如表4所示。股骨近端幾何參數對于髖關節假體選擇有著重要的指導意義。再加上近年來定制化髖關節假體的興起，股骨近端幾何參數為髖關節股骨柄髓外部分定制提供了設計依據。Schumann等[20]的重建模型主要對比了股骨近端幾何參數的重建精度，其重建模型的偏置距離和頸干角的平均誤差分別為1.56 mm和2.21°。而根據本文方法重建得到的三維模型，偏置距離和頸干角的平均誤差分別為1.389 mm和1.163°。除股骨頭直徑外本文重建模型的近端幾何參數精度都有所提高，并且股骨頭直徑誤差也保持在較低水平。根據表4的計算結果，股骨頭直徑、頸干角、偏置距離和前傾角的重建精度都保持在較高水平，說明重建模型可以滿足髖關節假體的定制化使用需求。
表4
重建模型部分幾何參數誤差值
Table4.
The error values of some geometric parameters of the reconstructed model
之前有關股骨統計形狀模型的研究主要集中在重建股骨的局部區域，并且沒有考慮患者之間的差異性，無法實現針對患者的個性化重建，顯然是不太合理的。本文基于高斯過程形變以及高斯過程回歸方法，將患者股骨頭直徑、頸干角以及股骨長度等信息作為三維重建的參考依據，實現了在較短時間內針對患者股骨的個性化重建。根據建模結果誤差分析可知，較少的訓練樣本數量會造成較大的建模誤差。所以后續研究中可以通過增加訓練集樣本數量，提高模型的重建精度。
本文還存在一定的局限性。首先由于本文采用的數據集來源于公開數據庫，其中并沒有每位志愿者的個人信息，所以沒有考慮年齡和性別對于模型重建的影響。其次本文的重建模型精度還有待提高，尤其在股骨近端和遠端，在后續研究中可以收集更多的數據來提高模型重建精度。同時，本文也具有一定的優勢。首先本文利用較少的外部數據實現了股骨三維模型的重構，為骨骼模型的三維重建提供了一個新的思路。其次本文考慮到患者之間的差異性，實現了利用不同患者數據，重構股骨三維模型，為髖關節假體的選擇和股骨柄髓外幾何參數的定制提供了新的方法。
6 結論
本文基于統計形狀模型，建立了利用患者二維醫學影像數據實現股骨快速三維重建的方法。該方法可以根據患者X線平片中包含的參數信息，5 s之內實現對患者股骨較高精度的三維重建，模型表面平均誤差在2 mm以內，重建模型的股骨頭直徑、頸干角、偏置距離和前傾角的平均誤差分別為0.597 mm、1.163 °、1.389 mm和1.354 °，極大地節省了獲得三維骨幾何模型所需要的時間。本文建立的骨幾何快速三維重建方法可以應用到骨科手術規劃、髖關節假體選擇及生物力學研究當中，以縮短建模時間，降低建模成本，快速獲得所需要的骨幾何三維模型。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：張志偉和陳瑱賢構思并設計了研究方案，張志偉實施方案并處理數據，張志偉和張志峰撰寫論文，王彩梅、陳瑱賢和靳忠民討論并修改了論文。
