本發(fā)明涉及一種基于顱面稠密對應點云的顱面形態(tài)分析及面貌復原方法，屬于計算機圖形學、數字幾何處理、人體測量學技術領域。
背景技術：
人體頭面部幾何形態(tài)復雜，顱骨形態(tài)決定了面貌的基本形態(tài)，顱骨和面貌形態(tài)間存在相關性，對其關系進行研究已經成為法醫(yī)人類學領域研究的熱點。傳統(tǒng)研究方法大多針對尸源進行軟組織厚度測量、頭面部測量指標統(tǒng)計分析，然而該方法存在研究對象樣本有限，測量指標數量少的局限。隨著醫(yī)學影像設備的發(fā)展，CT、MRI等影像設備已經被廣泛應用于活體樣本的數據采集和三維建模，為新技術和新方法的開展提供了數據基礎。
目前，針對顱面形態(tài)關系的研究方法主要包括三大類：1，軟組織厚度測量和統(tǒng)計。法醫(yī)人類學家利用計算機軟件在影像數據或者三維模型上標定特征點，然后統(tǒng)計不同年齡、性別、種族人群特征點處的軟組織厚度。目前，各國學者已經測量了澳大利亞、葡萄牙、埃及、高加索、捷克、土耳其、匈牙利、斯洛伐克、印度、美國、法國、日本、巴西、韓國、巴基斯坦、中國多個國家或民族的軟組織厚度，并通過t-test、多元回歸分析方法研究面部軟組織的分布，借此發(fā)現顱骨和面貌之間的相關性。為了保證軟組織測量的一致性，Stephan基于已經發(fā)布的軟組織厚度數據集，統(tǒng)計分析了6700余個樣本的軟組織厚度。2，頭面部幾何測量指標的統(tǒng)計。法醫(yī)人類學家利用計算機軟件在影像數據或者三維模型上標定顱骨和對應的面貌特征點，然后分別測量顱骨幾何測量項和對應的面貌幾何測量項，包括距離、角度、面積的指標，進而分析顱骨和面貌形態(tài)之間的關系，如顱骨形狀和臉型的關系、顱骨五官和面貌五官之間的關系。3，基于稠密點云的顱面形態(tài)表示。Berar等在2006年針對顱面稠密對應點云首次提出基于線性計算模型的顱面形態(tài)關系方法，該方法將顱骨點云和對應的面貌點云看作一個向量，通過最小二乘擬合方法實現顱面形態(tài)關系的表示?？紤]到顱骨五官與面貌五官間更為復雜的關系，Zhang在Berar方法的基礎上將顱骨和面貌模型按照五官分區(qū)，將每部分的顱骨點云和對應的面貌點云看作一個向量。為了顯示地表示顱面形態(tài)關系，針對顱骨和面貌模型，利用偏最小二乘回歸方法建立顱骨點云坐標和面貌點云坐標間的相互關系，然而由于顱骨和面貌點云數量大，給表示模型計算造成了極大困難。因此，依據統(tǒng)計形變模型(statistic morphable model)的思想，現有方法大多將顱骨和面貌模型分別表示為稠密對應點云，進而利用主成分分析方法實現顱骨和面貌模型的降維表示，最后通過機器學習方法學習顱骨樣本主成分系數和面貌樣本主成分系數間的關系，實現顱面形態(tài)關系的定量表示。該類方法的核心是利用主成分系數表示顱骨和面貌的幾何形態(tài)。目前，支持向量機、特征根回歸、偏最小二乘回歸等機器學習方法已經被用于顱面形態(tài)關系的定量表示。段福慶將顱骨和面貌模型按照五官分區(qū)，然后利用偏最小二乘回歸方法實現顱面形態(tài)關系的學習。鄧擎瓊建立了考慮年齡、性別、身體質量指數的回歸模型，提高了表示的準確性。
但上述方法仍不能有效解釋顱骨主成分和面貌主成分的幾何意義，不能有效解釋基于主成分和線性回歸的顱面形態(tài)分析和表示方法的可行性。另一方面，由于顱面形態(tài)關系復雜且面部軟組織在各區(qū)域的分布并不一致，在頭部整體數據集上采用主成分分析技術進行表示將降低顱面形態(tài)關系表示的準確性。因此需要依靠面部軟組織分布，建立一種更有效的顱面分區(qū)方法，使得各分區(qū)內形態(tài)關系盡可能一致，從而提高顱面形態(tài)關系表示的準確性。此外，顱骨稠密對應點云和面貌稠密對應點云的建立方法還亟待加強研究。
技術實現要素：
本發(fā)明的目的是針對現有技術存在的問題，提供一種基于顱面稠密對應點云的顱面形態(tài)分析及面貌復原方法，本發(fā)明通過特征自動提取和基于隱式函數和頂點能量約束最優(yōu)結合的非剛性配準算法，實現顱骨點云間和面貌點云間的稠密對應。通過可視和定量分析結合的方法，解決基于主成分分析的顱面形態(tài)關系幾何意義難以解釋的問題。通過基于軟組織分區(qū)的顱面形態(tài)關系線性學習方法，解決復雜顱面形態(tài)關系難以準確表示的問題，提高無身源顱骨面貌復原的準確性。本發(fā)明應用于刑事案件、考古學中未知身源顱骨的面貌復原，并能推廣到應用于醫(yī)學領域中人體器官、腫瘤等形狀的幾何分析，本發(fā)明包括以下步驟：
步驟1：顱骨稠密對應點云；
步驟2：面貌稠密對應點云；
步驟3：顱面形態(tài)關系可視分析；
步驟4：基于軟組織分區(qū)的顱面形態(tài)關系表示；
步驟5：未知身源顱骨的面貌復原。
進一步地，所述步驟1包括以下具體步驟：
步驟1.1：針對顱骨三維網格模型，定義內點和邊界點，實現顱骨孔洞邊緣的自動識別；并通過孔洞邊界長度、孔洞形狀和位置，分析確定孔洞之間的對應關系；
步驟1.2：基于高斯映射和動態(tài)區(qū)域增長算法，實現顱骨幾何形狀變化較大區(qū)域包括上頜、下頜、顴骨等的自動分割；
步驟1.3：從顱骨模型數據集中選擇選擇兩個模型，一個作為參考顱骨，另一個作為目標顱骨，依據步驟1.1和步驟1.2產生的兩個顱骨模型特征，利用最近點迭代算法實現參考顱骨向目標顱骨的剛性配準，進一步提出通過基于隱式函數和頂點能量約束最優(yōu)結合的方法，實現兩個顱骨模型的非剛性配準；
步驟1.4：計算目標模型的每個頂點與變形后的參考模型的最近點，記錄其頂點序號，建立目標模型與參考模型頂點的對應關系；
步驟1.5：依據步驟1.4計算獲得的頂點序號，針對原始參考顱骨模型，生成與目標顱骨點云對應的參考顱骨頂點坐標；
步驟1.6：從顱骨模型數據集中選擇其他模型作為參考顱骨，重復步驟1.3-步驟1.5，直到遍歷所有模型為止，從而建立顱骨模型間的頂點對應關系，即所有顱骨模型具有相同的頂點個數且對應頂點具有近似的解剖學位置。
所述步驟1.1進一步包括：
步驟1.1.1：定義顱骨三維網格模型為skull＝{P,E}，其中P＝{p1,p2,...,pn},pi＝(xi,yi,zi)∈R3表示n個顱骨頂點，E＝{ek＝(pi,pj),k＝1,2,3,...,m}表示m條邊。孔洞識別前首選對模型進行檢驗，判斷其是否滿足流行，并記錄不滿足條件的邊序號?？锥醋R別過程中，首先遍歷模型中各個頂點pi，確定其對應的鄰接頂點集合Adjpi，如果Adjpi中的點通過模型的邊ek直接連成封閉的多邊形，則pi為內點，否則為邊界點。然后，以任意邊界點pj為起點，根據ek尋找下一個邊界點，直到遍歷完所有邊界點。最終則會獲得c條封閉邊界，記為boundaryi＝{pj},i＝1,2,3,...,c。
步驟1.1.2：計算每條邊界的長度、中心，通過分析邊界長度和中心位置即能識別左眼眶、右眼眶和鼻骨三個輪廓。通過統(tǒng)計訓練集中樣本的左眼眶、右眼眶、鼻骨邊界輪廓的長度和中心坐標確定閾值。
所述步驟1.2進一步包括：
步驟1.2.1：計算顱骨模型各頂點pi的法矢ni＝(xi,yi,zi)和高斯曲率gaussi，通過高斯映射將顱骨頂點映射到單位球上，映射后的頂點坐標為
步驟1.2.2：顱骨模型的上頜、下頜、顴骨區(qū)域的幾何形狀復雜，曲面幾何形狀變化大，通過分析頂點的曲率和法矢法向完成上述區(qū)域的分割。選擇高斯曲率gaussi最大的頂點作為種子點，采用動態(tài)區(qū)域增長算法合并鄰接頂點，直到遍歷所有頂點后停止，區(qū)域增長的條件決定了特征提取的準確性，其區(qū)域合并條件定義如下：①||p'i-p'j||＜δ，高斯球上兩個相鄰頂點p'i和p'j間的距離描述了兩個法矢方向的差異，其閾值為δ；②||gaussi-gaussj||＜ε，兩個相鄰頂點的高斯曲率gaussi和gaussj的差異，其閾值為ε；閾值均由統(tǒng)計反饋動態(tài)生成。
所述步驟1.3進一步包括：
步驟1.3.1：針對由步驟1.1和步驟1.2獲得的參考顱骨模型的特征點集，其包括顱骨上頜、下頜和顴骨，以及眼眶和鼻骨的孔洞邊緣點集兩部分，和步驟1.1和步驟1.2獲得的目標顱骨模型特征點集，采用最近點迭代算法實現兩個顱骨模型的剛性配準，配準過程中采用隨機抽樣一致算法去除特征點集的誤對應，建立頂點間的準確對應關系，提高配準結果的準確性，參考顱骨模型的特征點集記為S＝{si},si＝(xi,yi,zi)∈R3，目標顱骨對應的特征點集記為Q＝{qi},qi＝(xi,yi,zi)∈R3；
步驟1.3.2：為了快速建立參考顱骨模型和目標顱骨模型頂點間的對應關系，采用徑向基函數和帶有緊支撐的徑向基函數，實現參考顱骨模型與目標顱骨模型的非剛性配準，進而將最近點作為對應點。參考顱骨模型點集記為V＝{vi},vi＝(xi,yi,zi)∈R3，目標顱骨對應點記為U＝{ui},ui＝(xi,yi,zi)∈R3，該對應將作為步驟1.3.3中對應點誤差能量項中的初始對應關系；
步驟1.3.3：定義具有保剛性的能量函數，計算參考顱骨模型每個頂點pi對應的仿射變換進而根據仿射變換矩陣X＝[x1 x2 … xn]T實現參考顱骨和目標顱骨的非剛性配準。該能量函數由對應點誤差能量項Ed(X)、特征誤差項El(X)和局部剛性能量項Es(X)三部分構成：E(X)＝Ed(X)+αEs(X)+βEl(X)，其中α和β為權值。具體而言，對應點誤差能量項為wi為權值，參考模型與目標模型的初始對應關系由步驟1.3.2獲得，此后每次代過程中通過搜尋最近點確定對應關系(vi,ui)；局部剛性能量項為其中F表示由目標模型的頂點-邊組成的鄰接矩陣,為克羅內克乘積算子，G＝diag(1，1，1，1)為對角矩陣；特征誤差能量項為即由步驟1.3.1產生的對應點集。
進一步地，所述步驟2包括以下具體步驟：
步驟2.1：面貌模型表面特征線提取，計算面貌三維模型每個頂點的法線和高斯曲率，統(tǒng)計高斯曲率值局部最大且鄰接頂點法矢夾角差異較大的頂點，該頂點作為特征點；面貌五官模型耳朵、鼻子、嘴、眼睛自動分割，以頂點高斯曲率局部最大的頂點作為種子頂點，基于高斯映射和動態(tài)區(qū)域增長算法實現分割；
步驟2.2：從面貌模型數據集中選擇兩個模型，一個作為參考面貌，另一個作為目標面貌，針對步驟2.1生成的兩個面貌模型特征，利用最近點迭代算法實現參考面貌向目標面貌的剛性配準，進一步提出通過基于隱式函數和頂點能量約束最優(yōu)結合的方法，實現兩個面貌模型的非剛性配準；
步驟2.3：變形后的參考模型的每個頂點與目標模型的最近點即為當前頂點的對應點，記錄其頂點序號。針對原始參考面貌模型，生成與目標面貌點云對應的參考面貌頂點坐標；
步驟2.4：從面貌數據集中選擇其他模型作為參考面貌，重復步驟2.1-步驟2.3，直到所有面貌均已建立對應關系時停止。
所述步驟2.2進一步包括：
步驟2.2.1：針對由步驟2.1獲得的參考面貌模型特征點集和目標面貌模型特征點集，采用最近點迭代算法實現兩個面貌模型的剛性配準，配準過程中采用隨機抽樣一致算法去除特征點集的誤對應，建立頂點間的準確對應關系，提高配準結果的準確性，參考面貌模型的特征點集記為S＝{si},si＝(xi,yi,zi)∈R3，目標面貌對應的特征點集記為Q＝{qi},qi＝(xi,yi,zi)∈R3；
步驟2.2.2：為了快速建立參考面貌模型和目標面貌模型頂點間的對應關系，采用徑向基函數和帶有緊支撐的徑向基函數，實現參考面貌模型與目標面貌模型的非剛性配準，進而將最近點作為對應點，參考面貌模型的點集記為V＝{vi},vi＝(xi,yi,zi)∈R3，目標面貌的對應點記為U＝{ui},ui＝(xi,yi,zi)∈R3，該對應作為步驟2.2.3中對應點誤差能量項的初始對應關系；
步驟2.2.3：定義具有保剛性的能量函數，計算參考面貌模型每個頂點pi對應的仿射變換進而根據仿射變換矩陣X＝[x1 x2 … xn]T實現參考面貌和目標面貌的非剛性配準。該能量函數由對應點誤差能量項Ed(X)、特征誤差項El(X)和局部剛性能量項Es(X)三部分構成：E(X)＝Ed(X)+αEs(X)+βEl(X)，其中α和β為權值。具體而言，對應點誤差能量項為wi為權值，參考模型與目前模型的初始對應關系由步驟2.2.2獲得，此后每次代過程中通過搜尋最近點確定對應關系(vi,ui)；局部剛性能量項為其中F表示由目標模型的頂點-邊組成的鄰接矩陣,為克羅內克乘積算子，G＝diag(1,1,1,1)為對角矩陣；特征誤差能量項為即由步驟2.2.1產生的對應點集。
進一步地，所述步驟3包括以下具體步驟：
步驟3.1：利用主成分分析方法，降維表示顱骨稠密點云，計算特征值和特征向量。利用主成分分析方法，降維表示面貌稠密點云，計算特征值和特征向量。
步驟3.2：為了觀察各主成分對模型幾何形狀的影響，針對顱骨模型，從第一個主成分開始，將其對應的主成分系數設置為給定值value＝3·λ1δ1，其中λ1＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ1為該主成分系數的方差，同理，針對面貌模型，從第一個主成分開始，將其對應的主成分系數設置為給定值value＝3·λ2δ2，其中λ1＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ2為該主成分系數的方差；
步驟3.3：為了觀察基于主成分表示的顱骨和面貌模型間的相互關系。每次分別取相同貢獻率的顱骨主成分和面貌主成分，將當前選擇的顱骨主成分系數和對應的面貌主成分系數分別設置為給定值如-1.0、-0.8、-0.6、-0.4、-0.2、0.2、0.4、0.6、0.8、1.0，其他主成分系數值為0，顯示顱骨模型和對應的面貌模型；
步驟3.4：計算顱骨主成分系數和對應的面貌主成分系數的相關性，判斷其是否近似滿足線性關系；如果滿足線性相關，則采用步驟4.4中的最小二乘法進行形態(tài)關系學習。
進一步地，所述步驟4包括以下具體步驟：
步驟4.1：針對數據集中的每個樣本，計算每個頂點的軟組織厚度值；針對所有樣本，計算每個頂點的軟組織厚度均值和方差，利用改進的K均值聚類算法，按軟組織厚度進行聚類如分為四類；聚類過程中首先均勻采樣設定聚類中心和較小粒度的聚類條件，然后完成初始分類并建立各分類間的鄰接關系無向圖，最后以含有頂點數較多的分類為中心通過合并鄰接分類完成指定數量的聚類；
步驟4.2：針對目標顱骨和目標面貌模型，依據每個頂點對應的軟組織厚度分類，將顱骨頂點和面貌頂點進行分區(qū)，實現基于軟組織厚度的顱面分區(qū)；
步驟4.3：針對顱骨各個分區(qū)點云，計算主成分系數和特征向量；針對面貌各個分區(qū)點云，計算主成分系數和特征向量；
步驟4.4：針對每個分區(qū)數據集，設Skulll×p＝[α1,p,α2,p,...,αl,p]和Facel×q＝[b1,q,b2,q,...,bl,q]分別為該分區(qū)中每個樣本顱骨的主成分和對應面貌的主成分，則顱骨和面貌之間的形態(tài)關系M＝{Mi,i＝1,2,L,k}能表示為Mi＝argmin||Skull×Mi-Face||2+λ2||Mi||2。利用最小二乘法求解M得Mi＝(SkullT·Skull+λI)-1·SkullT·Face，其中λ為權值，I為單位矩陣，則顱骨和面貌間的形態(tài)關系表示為M＝{(M1,M2,L,Mk)}。
進一步地，所述步驟5包括以下具體步驟：
步驟5.1：針對待復原的未知身源顱骨，利用步驟1.1-步驟1.5實現未知身源顱骨模型與目標顱骨模型的非剛性配準，建立頂點間的對應，實現未知身源顱骨的分區(qū)和頂點對應；
步驟5.2：針對未知身源顱骨的每個分區(qū)，計算每個分區(qū)對應的主成分(skull1,skull2,L,skullk)。依據顱面形態(tài)關系M＝{(M1,M2,L,Mk)}，計算每個顱骨分區(qū)主成分對應的面貌分區(qū)主成分(Face1,Face2,L,Facek)。根據計算獲得每個面貌分區(qū)主成分及其對應的特征向量，計算面貌分區(qū)復原結果。
步驟5.3：建立邊緣約束的能量方程，求解各分區(qū)中每個頂點對應的仿射變換，實現面貌分區(qū)復原結果的光滑融合；
所述步驟5.3進一步包括：
步驟5.3.1：對每個分區(qū)點云進行三角剖分，建立頂點間的鄰接關系。統(tǒng)計每個分區(qū)的頂點數量，并將頂點數量最多的分區(qū)作為目標模型，其他分區(qū)作為參考模型并向該分區(qū)變形；
步驟5.3.2：提取各分區(qū)復原結果的邊緣輪廓點集。計算目標模型與每個參考模型邊緣輪廓間的對應關系，對應的參考模型的邊緣點集記為V＝{vi},vi＝(xi,yi,zi)∈R3，目標模型對應的點集記為U＝{ui},ui＝(xi,yi,zi)∈R3；
步驟5.3.3：定義具有保剛性的能量函數，計算參考模型每個頂點pi對應的仿射變換進而根據仿射變換矩陣X＝[x1 x2 … xn]T實現各分區(qū)融合。該能量函數由邊緣誤差項Eedge(X)和局部剛性能量項Es(X)兩部分構成：E(X)＝Eedge(X)+αEs(X)，其中α為權值。具體而言，邊緣誤差能量項為wi為權值，其對應關系在步驟5.3.2中生成；局部剛性能量項為其中F表示由目標模型的頂點-邊組成的鄰接矩陣,為克羅內克乘積算子，G＝diag(1,1,1,1)為對角矩陣。
本發(fā)明相對于現有技術具有如下的優(yōu)點和積極效果：
1本發(fā)明通過基于高斯映射和區(qū)域增長結合的方法，自動提取顱骨形狀變化大的區(qū)域作為特征，將自動提取的對應輪廓邊緣和特征作為特征點集，建立了隱式函數和頂點能量約束最優(yōu)結合的非剛性配準方法，提高了顱骨模型配準的準確性。本發(fā)明由于能夠自動提取特征點集，能有效克服配準過程中需要人工參與標定特征點的不足。通過建立非剛性能量優(yōu)化和頂點對應關系自動尋找的迭代機制，在非剛性配準過程中實現頂點對應關系的建立，提高了顱骨模型頂點間對應關系尋找的準確性。
2本發(fā)明能夠自動提取面貌特征點集，通過建立非剛性能量優(yōu)化和頂點對應關系自動確定的迭代機制，提高了面貌頂點間對應關系尋找的準確性，有效克服配準過程中需要人工參與標定特征點的不足。
3本發(fā)明通過合理選擇和設置顱骨主成分和面貌主成分，能夠可視地顯示顱骨幾何形態(tài)、面貌幾何形態(tài)與主成分之間的關系，發(fā)現了顱骨主成分和面貌主成分的幾何意義。通過主成分之間的相關性檢測，進一步解釋了基于主成分和線性回歸的顱面形態(tài)關系分析方法的幾何意義。
4本發(fā)明通過統(tǒng)計面部稠密頂點的軟組織厚度，實現了基于軟組織分布的顱面分區(qū)。相比較于基于位置的五官分區(qū)方法，使得各分區(qū)內的顱面形態(tài)關系更為一致，提高了利用最小二乘回歸方法進行顱面形態(tài)關系定量表示的準確性，克服了過擬合和擬合不足的問題。
5本發(fā)明通過建立邊緣約束的能量方程，提高了分區(qū)復原結果融合的準確性。
本發(fā)明提供一種基于顱面稠密對應點云的顱面形態(tài)分析及面貌復原方法，是對現有顱面形態(tài)分析方法進行了創(chuàng)新和擴展。通過自動確定對應特征和隱式函數和頂點能量約束最優(yōu)結合的非剛性配準方法，提高了顱骨點云和面貌點云稠密對應的準確性；通過分析面部稠密軟組織分布，建立了基于軟組織分布的顱面分區(qū)模型，提高了顱面形態(tài)關系表示的準確性。通過建立邊緣約束的能量方程，提高了分區(qū)復原結果融合的準確性。最終應用于無身源顱骨的面貌復原。
附圖說明
圖1本發(fā)明所述一種基于顱面稠密對應點云的顱面形態(tài)分析及面貌復原方法的流程圖。
具體實施方式
以下結合附圖對本發(fā)明的實施例進行詳細描述，但本發(fā)明并不僅僅限于這些實施例。本發(fā)明涵蓋任何在本發(fā)明的精髓和范圍上做的替代、修改、等效方法以及方案。
如圖1所示，本發(fā)明一種基于顱面稠密對應點云的顱面形態(tài)分析及面貌復原方法包括如下步驟：
1、顱骨稠密對應點云；
2、面貌稠密對應點云；
3、顱面形態(tài)關系可視分析；
4、基于軟組織分區(qū)的顱面形態(tài)關系表示；
5、未知身源顱骨的面貌復原。
如圖1所示，在所述步驟1顱骨稠密對應點云和步驟2面貌稠密點云完成后，根據顱面形態(tài)分析的方式進行步驟3顱面形態(tài)關系可視分析和步驟4基于軟組織分區(qū)的顱面形態(tài)關系表示。最后依據顱面形態(tài)關系完成步驟5未知身源顱骨的面貌復原。
所述步驟1包括以下具體步驟：
步驟1.1：針對顱骨三維網格模型，定義內點和邊界點，實現顱骨孔洞邊緣的自動識別；并通過孔洞邊界長度、孔洞形狀和位置，分析確定孔洞之間的對應關系；
步驟1.2：基于高斯映射和動態(tài)區(qū)域增長算法，實現顱骨幾何形狀變化較大區(qū)域包括上頜、下頜、顴骨等的自動分割；
步驟1.3：從顱骨模型數據集中選擇選擇兩個模型，一個作為參考顱骨，另一個作為目標顱骨，依據步驟1.1和步驟1.2產生的兩個顱骨模型特征，利用最近點迭代算法實現參考顱骨向目標顱骨的剛性配準，進一步提出通過基于隱式函數和頂點能量約束最優(yōu)結合的方法，實現兩個顱骨模型的非剛性配準；
步驟1.4：計算目標模型的每個頂點與變形后的參考模型的最近點，記錄其頂點序號，建立目標模型與參考模型頂點的對應關系；
步驟1.5：依據步驟1.4計算獲得的頂點序號，針對原始參考顱骨模型，生成與目標顱骨點云對應的參考顱骨頂點坐標；
步驟1.6：從顱骨模型數據集中選擇其他模型作為參考顱骨，重復步驟1.3-步驟1.5，直到遍歷所有模型為止，從而建立顱骨模型間的頂點對應關系，即所有顱骨模型具有相同的頂點個數且對應頂點具有近似的解剖學位置。
所述步驟1.1進一步包括：
步驟1.1.1：定義顱骨三維網格模型為skull＝{P,E}，其中P＝{p1,p2,...,pn},pi＝(xi,yi,zi)∈R3表示n個顱骨頂點，E＝{ek＝(pi,pj),k＝1,2,3,...,m}表示m條邊?？锥醋R別前首選對模型進行檢驗，判斷其是否滿足流行，并記錄不滿足條件的邊序號?？锥醋R別過程中，首先遍歷模型中各個頂點pi，確定其對應的鄰接頂點集合Adjpi，如果Adjpi中的點通過模型的邊ek直接連成封閉的多邊形，則pi為內點，否則為邊界點。然后，以任意邊界點pj為起點，根據ek尋找下一個邊界點，直到遍歷完所有邊界點。最終則會獲得c條封閉邊界，記為boundaryi＝{pj},i＝1,2,3,...,c。
步驟1.1.2：計算每條邊界的長度、中心，通過分析邊界長度和中心位置即能識別左眼眶、右眼眶和鼻骨三個輪廓。能通過統(tǒng)計訓練集中樣本的左眼眶、右眼眶、鼻骨邊界輪廓的長度和中心坐標確定閾值。
所述步驟1.2進一步包括：
步驟1.2.1：計算顱骨模型各頂點pi的法矢ni＝(xi,yi,zi)和高斯曲率gaussi，通過高斯映射將顱骨頂點映射到單位球上，映射后的頂點坐標為
步驟1.2.2：顱骨模型的上頜、下頜、顴骨區(qū)域的幾何形狀復雜，曲面幾何形狀變化大，通過分析頂點的曲率和法矢法向完成上述區(qū)域的分割。選擇高斯曲率gaussi最大的頂點作為種子點，采用動態(tài)區(qū)域增長算法合并鄰接頂點，直到遍歷所有頂點后停止，區(qū)域增長的條件決定了特征提取的準確性，其區(qū)域合并條件定義如下：①||p'i-p'j||＜δ，高斯球上兩個相鄰頂點p'i和p'j間的距離描述了兩個法矢方向的差異，其閾值為δ；②||gaussi-gaussj||＜ε，兩個相鄰頂點的高斯曲率gaussi和gaussj的差異，其閾值為ε；閾值均由統(tǒng)計反饋動態(tài)生成。
所述步驟1.3進一步包括：
步驟1.3.1：針對由步驟1.1和步驟1.2獲得的參考顱骨模型的特征點集，其包括顱骨上頜、下頜和顴骨，以及眼眶和鼻骨的孔洞邊緣點集兩部分，和步驟1.1和步驟1.2獲得的目標顱骨模型特征點集，采用最近點迭代算法實現兩個顱骨模型的剛性配準，配準過程中采用隨機抽樣一致算法去除特征點集的誤對應，建立頂點間的準確對應關系，提高配準結果的準確性，參考顱骨模型的特征點集記為S＝{si},si＝(xi,yi,zi)∈R3，目標顱骨對應的特征點集記為Q＝{qi},qi＝(xi,yi,zi)∈R3；
步驟1.3.2：為了快速建立參考顱骨模型和目標顱骨模型頂點間的對應關系，采用徑向基函數和帶有緊支撐的徑向基函數，實現參考顱骨模型與目標顱骨模型的非剛性配準，進而將最近點作為對應點，參考顱骨模型點集記為V＝{vi},vi＝(xi,yi,zi)∈R3，目標顱骨對應點記為U＝{ui},ui＝(xi,yi,zi)∈R3，該對應將作為步驟1.3.3中對應點誤差能量項中的初始對應關系；
步驟1.3.3：定義具有保剛性的能量函數，計算參考顱骨模型每個頂點pi對應的仿射變換進而根據仿射變換矩陣X＝[x1 x2 … xn]T實現參考顱骨和目標顱骨的非剛性配準。該能量函數由對應點誤差能量項Ed(X)、特征誤差項El(X)和局部剛性能量項Es(X)三部分構成：E(X)＝Ed(X)+αEs(X)+βEl(X)，其中α和β為權值。具體而言，對應點誤差能量項為wi為權值，參考模型與目標模型的初始對應關系由步驟1.3.2獲得，此后每次代過程中通過搜尋最近點確定對應關系(vi,ui)；局部剛性能量項為其中F表示由目標模型的頂點-邊組成的鄰接矩陣,為克羅內克乘積算子，G＝diag(1,1,1,1)為對角矩陣；特征誤差能量項為即由步驟1.3.1產生的對應點集。
進一步地，所述步驟2包括以下具體步驟：
步驟2.1：面貌模型表面特征線提取，計算面貌三維模型每個頂點的法線和高斯曲率，統(tǒng)計高斯曲率值局部最大且鄰接頂點法矢夾角差異較大的頂點，該頂點作為特征點；面貌五官模型耳朵、鼻子、嘴、眼睛自動分割，以頂點高斯曲率局部最大的頂點作為種子頂點，基于高斯映射和動態(tài)區(qū)域增長算法實現分割；
步驟2.2：從面貌模型數據集中選擇兩個模型，一個作為參考面貌，另一個作為目標面貌，針對步驟2.1生成的兩個面貌模型特征，利用最近點迭代算法實現參考面貌向目標面貌的剛性配準，進一步提出通過基于隱式函數和頂點能量約束最優(yōu)結合的方法，實現兩個面貌模型的非剛性配準；
步驟2.3：變形后的參考模型的每個頂點與目標模型的最近點即為當前頂點的對應點，記錄其頂點序號。針對原始參考面貌模型，生成與目標面貌點云對應的參考面貌頂點坐標；
步驟2.4：從面貌數據集中選擇其他模型作為參考面貌，重復步驟2.1-步驟2.3，直到所有面貌均已建立對應關系時停止。
所述步驟2.2進一步包括：
步驟2.2.1：針對由步驟2.1獲得的參考面貌模型特征點集和目標面貌模型特征點集，采用最近點迭代算法實現兩個面貌模型的剛性配準。配準過程中采用隨機抽樣一致算法去除特征點集的誤對應，建立頂點間的準確對應關系，提高配準結果的準確性。參考面貌模型的特征點集記為S＝{si},si＝(xi,yi,zi)∈R3，目標面貌對應的特征點集記為Q＝{qi},qi＝(xi,yi,zi)∈R3。
步驟2.2.2：為了快速建立參考面貌模型和目標面貌模型頂點間的對應關系，采用徑向基函數和帶有緊支撐的徑向基函數，實現參考面貌模型與目標面貌模型的非剛性配準，進而將最近點作為對應點。參考面貌模型的點集記為V＝{vi},vi＝(xi,yi,zi)∈R3，目標面貌的對應點記為U＝{ui},ui＝(xi,yi,zi)∈R3，該對應作為步驟2.2.3中對應點誤差能量項的初始對應關系。
步驟2.2.3：定義具有保剛性的能量函數，計算參考面貌模型每個頂點pi對應的仿射變換進而根據仿射變換矩陣X＝[x1 x2 … xn]T實現參考面貌和目標面貌的非剛性配準。該能量函數由對應點誤差能量項Ed(X)、特征誤差項El(X)和局部剛性能量項Es(X)三部分構成：E(X)＝Ed(X)+αEs(X)+βEl(X)，其中α和β為權值。具體而言，對應點誤差能量項為wi為權值，參考模型與目前模型的初始對應關系由步驟2.2.2獲得，此后每次代過程中通過搜尋最近點確定對應關系(vi,ui)；局部剛性能量項為其中F表示由目標模型的頂點-邊組成的鄰接矩陣,為克羅內克乘積算子，G＝diag(1,1,1,1)為對角矩陣；特征誤差能量項為即由步驟2.2.1產生的對應點集。
進一步地，所述步驟3包括以下具體步驟：
步驟3.1：利用主成分分析方法，降維表示顱骨稠密點云，計算特征值和特征向量。利用主成分分析方法，降維表示面貌稠密點云，計算特征值和特征向量。
步驟3.2：為了觀察各主成分對模型幾何形狀的影響，針對顱骨模型，從第一個主成分開始，將其對應的主成分系數設置為給定值value＝3·λ1δ1，其中λ1＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ1為該主成分系數的方差，同理，針對面貌模型，從第一個主成分開始，將其對應的主成分系數設置為給定值value＝3·λ2δ2，其中λ1＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ2為該主成分系數的方差；
步驟3.3：為了觀察基于主成分表示的顱骨和面貌模型間的相互關系。每次分別取相同貢獻率的顱骨主成分和面貌主成分，將當前選擇的顱骨主成分系數和對應的面貌主成分系數分別設置為給定值如-1.0、-0.8、-0.6、-0.4、-0.2、0.2、0.4、0.6、0.8、1.0，其他主成分系數值為0，顯示顱骨模型和對應的面貌模型；
步驟3.4：計算顱骨主成分系數和對應的面貌主成分系數的相關性，判斷其是否近似滿足線性關系。如果滿足線性相關，則采用步驟4.4中的最小二乘法進行形態(tài)關系學習。
進一步地，所述步驟4包括以下具體步驟：
步驟4.1：針對數據集中的每個樣本，計算每個頂點的軟組織厚度值；針對所有樣本，計算每個頂點的軟組織厚度均值和方差，利用改進的K均值聚類算法，按軟組織厚度進行聚類如分為四類；聚類過程中首先均勻采樣設定聚類中心和較小粒度的聚類條件，然后完成初始分類并建立各分類間的鄰接關系無向圖，最后以含有頂點數較多的分類為中心通過合并鄰接分類完成指定數量的聚類；
步驟4.2：針對目標顱骨和目標面貌模型，依據每個頂點對應的軟組織厚度分類，將顱骨頂點和面貌頂點進行分區(qū)，實現基于軟組織厚度的顱面分區(qū)；
步驟4.3：針對顱骨各個分區(qū)點云，計算主成分系數和特征向量；針對面貌各個分區(qū)點云，計算主成分系數和特征向量；
步驟4.4：針對每個分區(qū)數據集，設Skulll×p＝[α1,p,α2,p,...,αl,p]和Facel×q＝[b1,q,b2,q,...,bl,q]分別為該分區(qū)中每個樣本顱骨的主成分和對應面貌的主成分，則顱骨和面貌之間的形態(tài)關系M＝{Mi,i＝1,2,L,k}能表示為Mi＝argmin||Skull×Mi-Face||2+λ2||Mi||2。利用最小二乘法求解M得Mi＝(SkullT·Skull+λI)-1·SkullT·Face，其中λ為權值，I為單位矩陣，則顱骨和面貌間的形態(tài)關系表示為M＝{(M1,M2,L,Mk)}。
進一步地，所述步驟5包括以下具體步驟：
步驟5.1：針對待復原的未知身源顱骨，利用步驟1.1-步驟1.5實現未知身源顱骨模型與目標顱骨模型的非剛性配準，建立頂點間的對應，實現未知身源顱骨的分區(qū)和頂點對應；
步驟5.2：針對未知身源顱骨的每個分區(qū)，計算每個分區(qū)對應的主成分(skull1,skull2,L,skullk)。依據顱面形態(tài)關系M＝{(M1,M2,L,Mk)}，計算每個顱骨分區(qū)主成分對應的面貌分區(qū)主成分(Face1,Face2,L,Facek)。根據計算獲得每個面貌分區(qū)主成分及其對應的特征向量，計算面貌分區(qū)復原結果。
步驟5.3：建立邊緣約束的能量方程，求解各分區(qū)中每個頂點對應的仿射變換，實現面貌分區(qū)復原結果的光滑融合；
所述步驟5.3進一步包括：
步驟5.3.1：對每個分區(qū)點云進行三角剖分，建立頂點間的鄰接關系。統(tǒng)計每個分區(qū)的頂點數量，并將頂點數量最多的分區(qū)作為目標模型，其他分區(qū)作為參考模型并向該分區(qū)變形；
步驟5.3.2：提取各分區(qū)復原結果的邊緣輪廓點集。計算目標模型與每個參考模型邊緣輪廓間的對應關系，對應的參考模型的邊緣點集記為V＝{vi},vi＝(xi,yi,zi)∈R3，目標模型對應的點集記為U＝{ui},ui＝(xi,yi,zi)∈R3；
步驟5.3.3：定義具有保剛性的能量函數，計算參考模型每個頂點pi對應的仿射變換進而根據仿射變換矩陣X＝[x1 x2 … xn]T實現各分區(qū)融合。該能量函數由邊緣誤差項Eedge(X)和局部剛性能量項Es(X)兩部分構成：E(X)＝Eedge(X)+αEs(X)，其中α為權值。具體而言，邊緣誤差能量項為wi為權值，其對應關系在步驟5.3.2中生成；局部剛性能量項為其中F表示由目標模型的頂點-邊組成的鄰接矩陣,為克羅內克乘積算子，G＝diag(1,1,1,1)為對角矩陣。
在本發(fā)明的實施例中：首先，針對人體頭部CT影像數據，通過圖像分割和邊緣輪廓跟蹤技術，提取顱骨和面貌的外表面，進而將獲得的顱骨點云和面貌點云模型變換到法蘭克福坐標系下。然后，通過提取顱骨點云模型的特征、孔洞邊緣，基于能量優(yōu)化的非剛性配準方法，建立顱骨稠密對應點云。特征自動提取過程中，為了消除特征誤對應的問題，需要考慮曲率一致性、法線一致性、距離的約束條件并結合RANSAC方法，提高特征對應的準確性。面貌稠密點云對應過程中也需要確保特征對應的一致性，消除誤匹配。最后，在顱面形態(tài)關系可視分析過程中，為了保證分析的有效性，需要選擇具有相同貢獻率的顱骨主成分和面貌主成分?？梢暦治鲞^程中，用戶應該同時設置顱骨主成分系數和面貌主成分系數值，比較模型的長度、寬度和角度的變化。
基于軟組織分區(qū)的顱面形態(tài)關系定量分析過程中必須針對已經建立好稠密對應關系的顱骨點云和面貌點云。首先計算每個頂點的軟組織厚度和整個樣本集的軟組織厚度均值，每個分區(qū)的軟組織均值可以通過統(tǒng)計臉頰處、額頭、下頜、顴骨的軟組織厚度均值獲得，從而提高分區(qū)結果的有效性。如果聚類后的分區(qū)過小，需要與相鄰分區(qū)進行合并。針對每個分區(qū)利用最小二乘擬合顱面形態(tài)關系時，需要引入擾動項，克服過擬合和擬合不足的問題。
以上所述，僅為本發(fā)明的具體實施方式，但本發(fā)明的保護范圍并不局限于此，任何熟悉本技術領域的技術人員在本發(fā)明公開的范圍內，能夠輕易想到的變化或替換，都應涵蓋在本發(fā)明權利要求的保護范圍內。