完整牙齒三維模型為正畸醫師進行錯頜畸形診斷、治療方案規劃提供了必要的輔助信息。目前，牙齒三維模型主要通過口腔計算機斷層掃描（CT）圖像分割與重構得到。然而，利用 CT 圖像重構的模型精度較低，不適用于無托槽隱形矯治器的輔助設計；另一方面，在正畸治療的不同周期對患者進行重復口腔 CT 掃描會給患者帶來輻射傷害。為此，本文提出基于口腔 CT 圖像與激光掃描圖像融合重構牙齒三維模型的方法。該方法利用口腔 CT 圖像重構的牙根與激光掃描圖像重構的牙冠經配準融合后建立完整的牙齒三維模型，所建立的模型牙冠部位精度更高。同時，該方法只需在正畸治療前對患者進行一次口腔 CT 掃描，在治療過程中進行激光掃描即可得到患者各治療周期完整的牙列模型，從而通過減少 CT 掃描次數降低給患者帶來的輻射傷害。
引用本文： 張東霞, 甘陽洲, 熊璟, 夏澤洋. 基于口腔計算機斷層掃描圖像與激光掃描圖像融合的牙齒三維模型重構. 生物醫學工程學雜志, 2017, 34(1): 7-14. doi: 10.7507/1001-5515.201604014 復制
引言
錯頜畸形不僅影響患者頜面美觀，還影響患者頜面發育及口腔功能，給患者的身體、社交及心理行為帶來不可忽視的傷害。錯頜畸形主要通過臨床口腔正畸治療進行矯治[1-2]。臨床上，正畸醫師需要利用牙列的三維信息進行錯頜畸形診斷、治療方案規劃以及矯治器設計。傳統的物理石膏模型制作過程復雜、制作周期長而且不含牙根。同時，物理石膏模型存在不易保存、攜帶不便及反復使用易損壞等缺點。隨著醫學成像及計算機技術的發展，利用口腔醫學圖像重建患者牙列的數字化三維模型已逐步應用于臨床口腔正畸輔助治療。
口腔領域常用的三維成像技術包括激光掃描（物理牙模激光掃描或口內激光掃描）和口腔計算機斷層掃描（computed tomography，CT）。激光掃描圖像主要包含患者牙冠的三維信息，其具有十微米級的圖像分辨率，所重構的模型精度較高。利用激光掃描圖像重建的患者牙冠模型被廣泛應用于無托槽隱形矯治器輔助設計中[3-7]。正畸治療的目標是要矯治患者牙齒使牙冠與牙根均排列整齊。在正畸治療前及治療過程中，醫師需要患者完整牙齒的三維信息配合診斷以及治療方案規劃。利用激光掃描圖像建立的牙冠三維模型在臨床應用中具有一定的局限性。
目前，完整牙齒的三維模型主要通過口腔 CT 圖像分割重構得到[8-9]。利用口腔 CT 圖像重建的模型雖然包含了完整的牙齒三維信息，但其精度相對較低，不適用于無托槽隱形矯正器的輔助設計[10]。同時，為了建立各周期不同牙頜形態下的牙齒三維模型，在治療過程中需要多次反復地對患者進行口腔 CT 掃描，給患者帶來輻射傷害。
本文提出了一種基于口腔 CT 圖像與激光掃描圖像融合的牙齒三維模型重構方法。該方法利用口腔 CT 圖像重構的牙根與激光掃描圖像重構的牙冠經配準融合后，建立完整的牙齒三維模型。所建立的模型牙冠部位精度更高，可適用于無托槽隱形矯治器的輔助設計。同時，利用該方法，患者只需在正畸治療前進行一次口腔 CT 掃描，在各治療周期中進行激光掃描即可建立各周期不同牙頜形態下完整牙齒的三維模型，從而通過減少 CT 掃描次數來降低給患者帶來的輻射傷害。
1 材料和方法
1.1 測試數據
本研究已通過中國科學院深圳先進技術研究院倫理審查委員會的審查和批準，并獲得受試者簽署的知情同意書。本文所采用的口腔 CT 圖像和激光掃描圖像均在中山大學光華口腔醫院獲取。本文共使用 5 位錯頜畸形患者的圖像進行實驗測試。所有患者的牙冠無缺損且牙根完整；部分患者存在拔牙情況。采集圖像受試者包括 3 名男性成年患者和 2 名女性成年患者，平均年齡為（23±1.41）歲。所采用的激光掃描圖像通過口內掃描儀（Konica Minolta，Aichi，Japan，50 μm）對患者進行口內掃描得到。所采用的口腔 CT 圖像由 CT 掃描機（NewTom VG，Italy）在 120 kV 的掃描電壓，5 mA 的掃描電流和曝光時間 6 s 情況下，對患者進行口腔 CT 掃描得到。CT 圖像的層厚為 1 mm，像素分辨率為 0.30 mm。
1.2 測試平臺
本研究在 Windows 8.1 系統下 MATLAB R2014a 軟件中通過自編代碼完成。
1.3 方法
基于口腔 CT 圖像與激光掃描圖像融合重建三維牙齒模型的流程如圖 1 所示。首先從激光掃描圖像中分割重構出牙冠模型；從口腔 CT 圖像中分割重構出完整的牙齒模型。然后，將由激光掃描圖像重構的牙冠模型配準到由 CT 圖像重構的牙齒模型的牙冠區域。最后，在配準模型上將激光掃描圖像重構出的牙冠與 CT 圖像重構的牙根融合重建，得到最終的三維牙齒模型。
圖1
基于口腔 CT 圖像與激光掃描圖像融合的牙齒三維模型重構流程圖
Figure1.
Diagram of 3D tooth model reconstruction based on fusion of dental CT images and laser-scanned images
2 模型分割與重構
2.1 激光掃描圖像的分割重構
激光掃描圖像獲取的是點云數據，本文先采用移動立方體（marching cube，MC）算法將激光掃描圖像點云數據重構為三角網格模型，然后采用一種改進的基于交互標記的快速分水嶺網格分割方法進行分割，得到獨立的激光牙冠模型[11]。
傳統基于交互標記的快速分水嶺網格分割方法以網格模型中三角面片的彎曲程度來設計高度函數，并以高度函數確定分割邊界[12]。該方法分割效率較高但往往存在過分割現象。同時，由于相鄰牙齒邊界處三角面片彎曲程度變化緩慢，利用該方法也難以實現相鄰牙齒的分割。本文提出一種改進的網格分割方法，其分割流程如圖 2 所示。對傳統方法的改進主要包括兩方面：①在使用快速分水嶺算法之前，在每顆牙齒的咬合區域使用區域生長算法進行初始分割，避免傳統方法在咬合區域出現過分割現象[13]。②同時利用三角網格模型中三角網格面片的彎曲程度和三角面片的面積設計高度函數以實現相鄰牙齒的分割。
圖2
由激光掃描圖像分割重建獨立的牙冠模型流程圖
Figure2.
Diagram of separated crown segmentation and reconstruction from laser-scanned images
對于給定的三角面片f1，f1 與其一階鄰域面片f2間的高度函數H（f1，f2）定義如式（1）所示：
| $ H({f_1},{f_2}) = {w_1}\frac{{area({f_1}) \! + \! area({f_2})}}{{\max \{ area({f_i})\} }} \! + \! {w_2}C({f_1},{f_2}) $ |
其中：fi 表示f1 的所有一階鄰域面片，area（·）表示三角面片的面積，max{·}表示取最大值算子，C（f1，f2）表示三角面片的彎曲程度[12]。w1 和w2 分別表示三角面片面積和三角面片彎曲程度的權重系數，根據模型差異需要調整，本文選擇w1=0.5，w2=60。
與傳統方法相比，本文提出的改進方法對牙冠網格模型的分割結果更準確，部分對比結果如圖 3 所示。
圖3
牙冠網格模型分割結果對比
Figure3.
Result comparison using different methods for crown mesh model segmentation
2.2 口腔 CT 圖像的分割重構
口腔 CT 圖像的分割重構流程圖如圖 4 所示。用戶首先需要手動選擇一張牙頸處的切片作為初始切片，并在這張初始切片上為每顆牙齒標記初始種子點。接著，從這張切片開始分別沿著牙冠方向和牙根方向自動實現各切片中獨立牙齒輪廓的分割。輪廓自動分割過程中，前一張切片的分割結果將作為當前待分割切片中牙齒的初始輪廓，利用混合水平集方法完成牙齒輪廓分割[8]。分割完成后，對分割出的每顆牙齒區域設置不同的灰度值，然后利用 MC 算法設置相應的閾值實現各獨立牙齒三維模型重構。
圖4
口腔 CT 圖像分割重構流程圖
Figure4.
Diagram of segmentation and reconstruction of dental CT images
3 模型配準
激光掃描圖像重構的牙冠模型與口腔 CT 圖像重構的牙齒模型之間的配準分為 2 個階段：粗配準和精配準。粗配準采用主成分分析（principal component analysis，PCA）算法，以實現兩種模型的大致對齊；精配準采用迭代最近點（iterative closest point，ICP）算法，對大致對齊后的兩模型的位置進行精細調整[14-15]。在進行配準操作前，本文用一平面從口腔 CT 圖像重構的牙齒模型中截取牙冠部分，所截取的牙冠將與激光掃描圖像重構的牙冠模型進行配準。
利用 PCA 算法進行粗配準的流程為：①分別計算兩種模型所有頂點的協方差矩陣Covlaser 和CovCT ；②分別計算Covlaser 和CovCT 的三個相互正交的單位特征向量，并以此建立兩種模型的貼體坐標系，坐標原點在各模型所有頂點的質心位置。③通過坐標系變換使兩種模型的貼體坐標系大致對齊，即求取粗配準的旋轉矩陣R0 和平移矩陣T0。設Covlaser 和CovCT 的單位特征向量矩陣分別為Rlaser 和RCT ，旋轉矩陣R0 和平移矩陣T0 的求解公式為：
| $ {R}_0 = {R}_{CT}{\rm {inv}} ({R}_{laser}) $ |
| $ {T}_0 = center\left( p \right) - {{R}_0 }center\left( q \right) $ |
其中，inv（·）表示矩陣的逆操作算子，p 和q 分別表示兩種模型所有頂點坐標構成的矩陣（以各頂點坐標為列向量），center（·）表示求質心。
利用 ICP 算法進行精配準的過程，即搜索使兩種模型中對應點距離的平方和均值誤差最小時的旋轉矩陣R1 和平移矩陣T1。
設粗配準后激光掃描圖像重構的牙冠模型頂點集合為p’，ICP 算法實現精配準的流程如下：
①利用k 維樹結構化粗配準后兩種模型的頂點集合[16]；
②對于p’ 中的任一頂點p’i ，在 CT 圖像重構得到的牙冠模型頂點集合q 中搜索與p’i 距離最近的頂點qj ，并求所有的點對（p’i ，qj ）之間距離平方和均方根誤差e：
| $ e =\sqrt { \frac{{\sum\limits_{i = 1}^N {||{{p'}_i} - {q_j}|{|^2}} }}{{ N}}} $ |
其中，N 為p’ 中頂點個數，||·||表示取模算子。
③利用四元數法求兩模型配準矩陣R1k 和T1k[15]；
中間矩陣Q 的求解如下：
| $ {\rm{Q}} = \left[ {\begin{array}{*{20}{c}} {tr\left( {\sum {({{p'}_k},q)} } \right)} & {{\Delta ^{\rm {\rm T}}}}\\ \Delta & {\sum {({{p'}_k},q)} + {{\left( {\sum {({{p'}_k},q)} } \right)}^{\rm T}} - {\rm tr}\left( {\sum {({{p'}_k},q)} } \right){{\rm{I}}_3}} \end{array}} \right] $ |
其中，tr（·）表示矩陣的跡，（·）T表示求矩陣的轉置，I3表示三階單位矩陣；
| $ \begin{array}{l} \sum {({{p'}_k},q) = \frac{1}{N}} \sum\limits_{i = 1}^N {\left[ {\left( {{{({{p'}_k})}_i} - mean\left( {{{p'}_k}} \right)} \right)} \right.} \\ \quad\quad\quad\quad\quad\ \ \left. {{{({q_j} - mean\left( q \right))}^{\rm T}}} \right] \end{array} $ |
| $ {\rm{\Delta }} = {[\begin{array}{*{20}{c}} {{A_{23}}} & {{A_{31}}} & {{A_{12}}} \end{array}]^{\rm T}} $ |
| $ {A_{i,j}} = {\left[ {\sum {({{p'}_k},q)} - {{\left( {\sum {({{p'}_k},q)} } \right)}^{\rm T}}} \right]_{i,j}} $ |
式（6）中，mean（·）表示求均值。
設Q 最大特征值對應的特征向量為Qq ，Qq 用式（9）表示。
| $ {{\rm{Q}}_q} = \left[ {\begin{array}{*{20}{c}} {{q_0}} & {{q_1}} & {{q_2}} & {{q_3}} \end{array}} \right] $ |
求當前的旋轉矩陣R1k 和平移矩陣T1k，如式（10）和（11）所示。
| $ {{\rm{R}}_{1k}} = \left[ {\begin{array}{*{20}{c}} {{q_0} + {q_1} - {q_2}^2 - {q_3}^2} & {2({q_1}{q_2} - {q_0}{q_3})} & {2({q_1}{q_3} + {q_0}{q_2})}\\ {2({q_1}{q_2} + {q_0}{q_3})} & {{q_0}^2 + {q_1}^2 - {q_1}^2 - {q_3}^2} & {2({q_2}{q_3} - {q_0}{q_1})}\\ {2({q_1}{q_3} - {q_0}{q_2})} & {2({q_2}{q_3} + {q_0}{q_1})} & {{q_0}^2 + {q_3}^2 - {q_1}^2 - {q_2}^2} \end{array}} \right] $ |
| $ {T}_{1k} = center\left( q \right) - {R}_{1k}center\left( {{{p'}_k}} \right) $ |
④如式（12）所示，對p’k 所有頂點進行旋轉平移，實現與q 的配準。
| $ {p'_{\left( {k + 1} \right)}} = {R}_{1k}({p'_k}) + ({T}_{1k}) $ |
⑤重復上述迭代過程，直到前后兩次迭代的兩模型對應點間距離平方和均方根誤差的差值?e小于預設的閾值。
4 模型融合
本文采用德洛內區域生長（Delaunay-based region growing，DBRG）算法對配準以后由激光掃描圖像重構的牙冠與由口腔CT圖像重構的牙根進行融合，從而建立完整的牙齒三維模型[17]。搜索配準后激光掃描圖像重構的牙冠模型頂點在CT圖像重構的牙齒模型中的 1 至15階最近點，將這些點從牙齒模型中刪除后得到的頂點即是CT圖像重構的牙根頂點，從這些頂點可重構牙根模型。
采用DBRG算法實現牙冠與牙根的融合過程如下：
①將牙冠模型的頂點與牙根模型的頂點德洛內三角形化得到三角面片集合，記作T，從而建立兩模型邊界頂點間的拓撲信息，實現模型的拼接。
②在集合T 中選擇一個初始三角面片tstart 。初始三角面片從與其相鄰三角面片之間的二面角較小的面片中選擇。定義三角面片集合S，將初始三角面片tstart 加入S 中；定義三角面片邊的集合B，將初始切片的三條邊加入集合B 中。
③定義集合Ai ，將B 中每一條邊bi 的相鄰三角面片加入Ai 中。定義集合Ai_saved ，求集合Ai 中的各三角面片與邊bi 所從屬的三角面片之間的二面角，將夾角大于 155° 的三角面片加入集合Ai_saved 中。
④對B 中的所有的邊進行以上操作，得到集合AI_SAVED ，即所有Ai_saved 的并集。在集合AI_SAVED 中求對應最大二面角的三角面片AI_SAVED_MAX ，并從集合AI_SAVED 中剔除AI_SAVED_MAX 。
⑤檢測三角面片AI_SAVED_MAX 的拓撲關系是否正常。如果正常，則進行步驟⑥；若不正常，則返回步驟④。
⑥更新集合S 和B。將三角面片AI_SAVED_MAX 加入集合S 中，同時將不在集合B 中的AI_SAVED_MAX 的邊加入集合B，并剔除已在集合B 中的AI_SAVED_MAX 的邊。接著進行步驟③并以此往下進行，直到集合AI_SAVED 為空。
在上述過程中，155° 的限制是為了避免區域生長過程中三角面片回溯生長；檢測三角面片的拓撲關系正常的判斷標準為：待檢測面片AI_SAVED_MAX 與S 中所有的三角面片之間不會產生空間上的壓疊。
通過此算法從T 中篩選出來的三角面片集合S 即為融合得到的三維牙齒網格模型，但是會出現空洞，因此我們在局部孔洞處再次使用以上算法進行修補。
5 實驗結果與分析
5.1 分割重構結果
由某患者激光掃描圖像分割重構的牙冠模型如圖 5 左圖所示，由該患者口腔 CT 圖像分割重構的牙齒模型如圖 5 右圖所示。
圖5
某患者激光掃描圖像與口腔 CT 圖像分割重構結果
Figure5.
Results of segmentation and reconstruction from laser-scanned images and dental CT images of one subject
5.2 模型配準結果
激光掃描圖像重構的牙冠模型與口腔 CT 圖像重構的牙冠模型配準結果如圖 6 所示。其中，紫色為激光掃描圖像重構的牙冠模型，藍色為口腔 CT 圖像重構的牙冠模型。視覺上，配準后兩模型較好地重疊在一起。
圖6
不同視角下的牙冠配準結果 藍色：口腔 CT 圖像分割重構的牙冠模型；紫色：激光圖像分割重構的牙冠模型
Figure6.
Crown registration results from two different views Blue: crown model based on images; Purple: crown model based on laser-scanned images
為量化配準結果，本文對配準后激光掃描圖像重構的牙冠模型上每一個頂點 在口腔 CT 圖像重構的牙冠模型上搜索離最近的頂點qj ，并計算兩者之間的距離，以此計算平均配準誤差和最大點對的配準誤差（單位：mm）[18-19]。平均配準誤差和最大點對的配準誤差求解公式為：
| $ AD = \frac{{\sum\limits_{i = 1}^N {||{{({{p'}_{final}})}_i} - {q_i}|{|^2}} }}{N} $ |
| $ MD = \max (||{({p'_{final}})_i} - {q_i}||) $ |
其中，N 為激光牙冠模型頂點的個數。
經過對 5 位成年患者的圖像進行測試，測試結果為：AD=（0.26±0.03）mm，MD=（0.54±0.25）mm。 其中，某一位患者的配準結果如圖 7 所示，配準結果中存在少數點對之間距離較大，其原因主要為：①本文所用 CT 圖像的空間分辨率較低（0.30 mm）；②口腔 CT 圖像中牙齒分割結果存在一定誤差。
圖7
不同視角下的配準誤差圖
Figure7.
Error map of registration result from two different views
5.3 模型融合結果
本文模型融合主要發生在兩模型的牙頸處，牙根數據的完整性不會影響整個實驗方法的驗證。為降低實驗復雜度又不失合理性，我們采用部分口腔 CT 牙根數據參與實驗，結果如圖 8 所示， 重構結果可以用于臨床正畸治療中輔助醫師制定正畸治療方案以及矯治器輔助設計。
圖8
不同視角下的融合結果
Figure8.
Fusion results from two different views
6 結論
完整牙齒三維模型為正畸醫師進行診療提供了必要的輔助信息。目前，牙齒三維模型主要通過口腔 CT 圖像分割與重構得到。一方面，口腔 CT 圖像重構的模型精度較低，不適用于無托槽隱形矯治器的輔助設計；另一方面，在正畸治療的不同周期對患者進行重復口腔 CT 掃描會給患者帶來輻射傷害。為此，從提高口腔 CT 圖像重構的完整三維牙齒模型的精度和減少不必要的 CT 掃描輻射的角度出發，本文提出一種基于口腔 CT 圖像與激光掃描圖像融合重構牙齒三維模型的方法。該方法首先采用區域增長算法和改進的基于標記的快速分水嶺算法從激光掃描圖像重構的網格模型中分割出獨立的牙冠模型；采用混合水平集算法和 MC 算法從口腔 CT 圖像中分割重構出完整的牙齒三維模型。然后采用 PCA 算法和 ICP 算法將兩模型的牙冠區域進行配準。最后，采用 DBRG 算法將激光掃描圖像重構的牙冠與口腔 CT 圖像重構的牙根融合，建立完整牙齒三維模型。實驗結果驗證了該方法能有效地利用激光掃描圖像中牙冠信息與口腔 CT 圖像中牙根信息融合來重建完整牙齒的三維模型。利用該方法重建的牙齒模型牙冠部位精度更高，可用于無托槽隱形矯治器的輔助設計。同時，該方法只需在正畸治療前對患者進行一 次口腔 CT 掃描，在治療過程中進行激光掃描即可得到患者各治療周期完整的牙列模型，從而通過減少 CT 掃描次數來降低給患者帶來的輻射傷害。
引言
錯頜畸形不僅影響患者頜面美觀，還影響患者頜面發育及口腔功能，給患者的身體、社交及心理行為帶來不可忽視的傷害。錯頜畸形主要通過臨床口腔正畸治療進行矯治[1-2]。臨床上，正畸醫師需要利用牙列的三維信息進行錯頜畸形診斷、治療方案規劃以及矯治器設計。傳統的物理石膏模型制作過程復雜、制作周期長而且不含牙根。同時，物理石膏模型存在不易保存、攜帶不便及反復使用易損壞等缺點。隨著醫學成像及計算機技術的發展，利用口腔醫學圖像重建患者牙列的數字化三維模型已逐步應用于臨床口腔正畸輔助治療。
口腔領域常用的三維成像技術包括激光掃描（物理牙模激光掃描或口內激光掃描）和口腔計算機斷層掃描（computed tomography，CT）。激光掃描圖像主要包含患者牙冠的三維信息，其具有十微米級的圖像分辨率，所重構的模型精度較高。利用激光掃描圖像重建的患者牙冠模型被廣泛應用于無托槽隱形矯治器輔助設計中[3-7]。正畸治療的目標是要矯治患者牙齒使牙冠與牙根均排列整齊。在正畸治療前及治療過程中，醫師需要患者完整牙齒的三維信息配合診斷以及治療方案規劃。利用激光掃描圖像建立的牙冠三維模型在臨床應用中具有一定的局限性。
目前，完整牙齒的三維模型主要通過口腔 CT 圖像分割重構得到[8-9]。利用口腔 CT 圖像重建的模型雖然包含了完整的牙齒三維信息，但其精度相對較低，不適用于無托槽隱形矯正器的輔助設計[10]。同時，為了建立各周期不同牙頜形態下的牙齒三維模型，在治療過程中需要多次反復地對患者進行口腔 CT 掃描，給患者帶來輻射傷害。
本文提出了一種基于口腔 CT 圖像與激光掃描圖像融合的牙齒三維模型重構方法。該方法利用口腔 CT 圖像重構的牙根與激光掃描圖像重構的牙冠經配準融合后，建立完整的牙齒三維模型。所建立的模型牙冠部位精度更高，可適用于無托槽隱形矯治器的輔助設計。同時，利用該方法，患者只需在正畸治療前進行一次口腔 CT 掃描，在各治療周期中進行激光掃描即可建立各周期不同牙頜形態下完整牙齒的三維模型，從而通過減少 CT 掃描次數來降低給患者帶來的輻射傷害。
1 材料和方法
1.1 測試數據
本研究已通過中國科學院深圳先進技術研究院倫理審查委員會的審查和批準，并獲得受試者簽署的知情同意書。本文所采用的口腔 CT 圖像和激光掃描圖像均在中山大學光華口腔醫院獲取。本文共使用 5 位錯頜畸形患者的圖像進行實驗測試。所有患者的牙冠無缺損且牙根完整；部分患者存在拔牙情況。采集圖像受試者包括 3 名男性成年患者和 2 名女性成年患者，平均年齡為（23±1.41）歲。所采用的激光掃描圖像通過口內掃描儀（Konica Minolta，Aichi，Japan，50 μm）對患者進行口內掃描得到。所采用的口腔 CT 圖像由 CT 掃描機（NewTom VG，Italy）在 120 kV 的掃描電壓，5 mA 的掃描電流和曝光時間 6 s 情況下，對患者進行口腔 CT 掃描得到。CT 圖像的層厚為 1 mm，像素分辨率為 0.30 mm。
1.2 測試平臺
本研究在 Windows 8.1 系統下 MATLAB R2014a 軟件中通過自編代碼完成。
1.3 方法
基于口腔 CT 圖像與激光掃描圖像融合重建三維牙齒模型的流程如圖 1 所示。首先從激光掃描圖像中分割重構出牙冠模型；從口腔 CT 圖像中分割重構出完整的牙齒模型。然后，將由激光掃描圖像重構的牙冠模型配準到由 CT 圖像重構的牙齒模型的牙冠區域。最后，在配準模型上將激光掃描圖像重構出的牙冠與 CT 圖像重構的牙根融合重建，得到最終的三維牙齒模型。
圖1
基于口腔 CT 圖像與激光掃描圖像融合的牙齒三維模型重構流程圖
Figure1.
Diagram of 3D tooth model reconstruction based on fusion of dental CT images and laser-scanned images
2 模型分割與重構
2.1 激光掃描圖像的分割重構
激光掃描圖像獲取的是點云數據，本文先采用移動立方體（marching cube，MC）算法將激光掃描圖像點云數據重構為三角網格模型，然后采用一種改進的基于交互標記的快速分水嶺網格分割方法進行分割，得到獨立的激光牙冠模型[11]。
傳統基于交互標記的快速分水嶺網格分割方法以網格模型中三角面片的彎曲程度來設計高度函數，并以高度函數確定分割邊界[12]。該方法分割效率較高但往往存在過分割現象。同時，由于相鄰牙齒邊界處三角面片彎曲程度變化緩慢，利用該方法也難以實現相鄰牙齒的分割。本文提出一種改進的網格分割方法，其分割流程如圖 2 所示。對傳統方法的改進主要包括兩方面：①在使用快速分水嶺算法之前，在每顆牙齒的咬合區域使用區域生長算法進行初始分割，避免傳統方法在咬合區域出現過分割現象[13]。②同時利用三角網格模型中三角網格面片的彎曲程度和三角面片的面積設計高度函數以實現相鄰牙齒的分割。
圖2
由激光掃描圖像分割重建獨立的牙冠模型流程圖
Figure2.
Diagram of separated crown segmentation and reconstruction from laser-scanned images
對于給定的三角面片f1，f1 與其一階鄰域面片f2間的高度函數H（f1，f2）定義如式（1）所示：
| $ H({f_1},{f_2}) = {w_1}\frac{{area({f_1}) \! + \! area({f_2})}}{{\max \{ area({f_i})\} }} \! + \! {w_2}C({f_1},{f_2}) $ |
其中：fi 表示f1 的所有一階鄰域面片，area（·）表示三角面片的面積，max{·}表示取最大值算子，C（f1，f2）表示三角面片的彎曲程度[12]。w1 和w2 分別表示三角面片面積和三角面片彎曲程度的權重系數，根據模型差異需要調整，本文選擇w1=0.5，w2=60。
與傳統方法相比，本文提出的改進方法對牙冠網格模型的分割結果更準確，部分對比結果如圖 3 所示。
圖3
牙冠網格模型分割結果對比
Figure3.
Result comparison using different methods for crown mesh model segmentation
2.2 口腔 CT 圖像的分割重構
口腔 CT 圖像的分割重構流程圖如圖 4 所示。用戶首先需要手動選擇一張牙頸處的切片作為初始切片，并在這張初始切片上為每顆牙齒標記初始種子點。接著，從這張切片開始分別沿著牙冠方向和牙根方向自動實現各切片中獨立牙齒輪廓的分割。輪廓自動分割過程中，前一張切片的分割結果將作為當前待分割切片中牙齒的初始輪廓，利用混合水平集方法完成牙齒輪廓分割[8]。分割完成后，對分割出的每顆牙齒區域設置不同的灰度值，然后利用 MC 算法設置相應的閾值實現各獨立牙齒三維模型重構。
圖4
口腔 CT 圖像分割重構流程圖
Figure4.
Diagram of segmentation and reconstruction of dental CT images
3 模型配準
激光掃描圖像重構的牙冠模型與口腔 CT 圖像重構的牙齒模型之間的配準分為 2 個階段：粗配準和精配準。粗配準采用主成分分析（principal component analysis，PCA）算法，以實現兩種模型的大致對齊；精配準采用迭代最近點（iterative closest point，ICP）算法，對大致對齊后的兩模型的位置進行精細調整[14-15]。在進行配準操作前，本文用一平面從口腔 CT 圖像重構的牙齒模型中截取牙冠部分，所截取的牙冠將與激光掃描圖像重構的牙冠模型進行配準。
利用 PCA 算法進行粗配準的流程為：①分別計算兩種模型所有頂點的協方差矩陣Covlaser 和CovCT ；②分別計算Covlaser 和CovCT 的三個相互正交的單位特征向量，并以此建立兩種模型的貼體坐標系，坐標原點在各模型所有頂點的質心位置。③通過坐標系變換使兩種模型的貼體坐標系大致對齊，即求取粗配準的旋轉矩陣R0 和平移矩陣T0。設Covlaser 和CovCT 的單位特征向量矩陣分別為Rlaser 和RCT ，旋轉矩陣R0 和平移矩陣T0 的求解公式為：
| $ {R}_0 = {R}_{CT}{\rm {inv}} ({R}_{laser}) $ |
| $ {T}_0 = center\left( p \right) - {{R}_0 }center\left( q \right) $ |
其中，inv（·）表示矩陣的逆操作算子，p 和q 分別表示兩種模型所有頂點坐標構成的矩陣（以各頂點坐標為列向量），center（·）表示求質心。
利用 ICP 算法進行精配準的過程，即搜索使兩種模型中對應點距離的平方和均值誤差最小時的旋轉矩陣R1 和平移矩陣T1。
設粗配準后激光掃描圖像重構的牙冠模型頂點集合為p’，ICP 算法實現精配準的流程如下：
①利用k 維樹結構化粗配準后兩種模型的頂點集合[16]；
②對于p’ 中的任一頂點p’i ，在 CT 圖像重構得到的牙冠模型頂點集合q 中搜索與p’i 距離最近的頂點qj ，并求所有的點對（p’i ，qj ）之間距離平方和均方根誤差e：
| $ e =\sqrt { \frac{{\sum\limits_{i = 1}^N {||{{p'}_i} - {q_j}|{|^2}} }}{{ N}}} $ |
其中，N 為p’ 中頂點個數，||·||表示取模算子。
③利用四元數法求兩模型配準矩陣R1k 和T1k[15]；
中間矩陣Q 的求解如下：
| $ {\rm{Q}} = \left[ {\begin{array}{*{20}{c}} {tr\left( {\sum {({{p'}_k},q)} } \right)} & {{\Delta ^{\rm {\rm T}}}}\\ \Delta & {\sum {({{p'}_k},q)} + {{\left( {\sum {({{p'}_k},q)} } \right)}^{\rm T}} - {\rm tr}\left( {\sum {({{p'}_k},q)} } \right){{\rm{I}}_3}} \end{array}} \right] $ |
其中，tr（·）表示矩陣的跡，（·）T表示求矩陣的轉置，I3表示三階單位矩陣；
| $ \begin{array}{l} \sum {({{p'}_k},q) = \frac{1}{N}} \sum\limits_{i = 1}^N {\left[ {\left( {{{({{p'}_k})}_i} - mean\left( {{{p'}_k}} \right)} \right)} \right.} \\ \quad\quad\quad\quad\quad\ \ \left. {{{({q_j} - mean\left( q \right))}^{\rm T}}} \right] \end{array} $ |
| $ {\rm{\Delta }} = {[\begin{array}{*{20}{c}} {{A_{23}}} & {{A_{31}}} & {{A_{12}}} \end{array}]^{\rm T}} $ |
| $ {A_{i,j}} = {\left[ {\sum {({{p'}_k},q)} - {{\left( {\sum {({{p'}_k},q)} } \right)}^{\rm T}}} \right]_{i,j}} $ |
式（6）中，mean（·）表示求均值。
設Q 最大特征值對應的特征向量為Qq ，Qq 用式（9）表示。
| $ {{\rm{Q}}_q} = \left[ {\begin{array}{*{20}{c}} {{q_0}} & {{q_1}} & {{q_2}} & {{q_3}} \end{array}} \right] $ |
求當前的旋轉矩陣R1k 和平移矩陣T1k，如式（10）和（11）所示。
| $ {{\rm{R}}_{1k}} = \left[ {\begin{array}{*{20}{c}} {{q_0} + {q_1} - {q_2}^2 - {q_3}^2} & {2({q_1}{q_2} - {q_0}{q_3})} & {2({q_1}{q_3} + {q_0}{q_2})}\\ {2({q_1}{q_2} + {q_0}{q_3})} & {{q_0}^2 + {q_1}^2 - {q_1}^2 - {q_3}^2} & {2({q_2}{q_3} - {q_0}{q_1})}\\ {2({q_1}{q_3} - {q_0}{q_2})} & {2({q_2}{q_3} + {q_0}{q_1})} & {{q_0}^2 + {q_3}^2 - {q_1}^2 - {q_2}^2} \end{array}} \right] $ |
| $ {T}_{1k} = center\left( q \right) - {R}_{1k}center\left( {{{p'}_k}} \right) $ |
④如式（12）所示，對p’k 所有頂點進行旋轉平移，實現與q 的配準。
| $ {p'_{\left( {k + 1} \right)}} = {R}_{1k}({p'_k}) + ({T}_{1k}) $ |
⑤重復上述迭代過程，直到前后兩次迭代的兩模型對應點間距離平方和均方根誤差的差值?e小于預設的閾值。
4 模型融合
本文采用德洛內區域生長（Delaunay-based region growing，DBRG）算法對配準以后由激光掃描圖像重構的牙冠與由口腔CT圖像重構的牙根進行融合，從而建立完整的牙齒三維模型[17]。搜索配準后激光掃描圖像重構的牙冠模型頂點在CT圖像重構的牙齒模型中的 1 至15階最近點，將這些點從牙齒模型中刪除后得到的頂點即是CT圖像重構的牙根頂點，從這些頂點可重構牙根模型。
采用DBRG算法實現牙冠與牙根的融合過程如下：
①將牙冠模型的頂點與牙根模型的頂點德洛內三角形化得到三角面片集合，記作T，從而建立兩模型邊界頂點間的拓撲信息，實現模型的拼接。
②在集合T 中選擇一個初始三角面片tstart 。初始三角面片從與其相鄰三角面片之間的二面角較小的面片中選擇。定義三角面片集合S，將初始三角面片tstart 加入S 中；定義三角面片邊的集合B，將初始切片的三條邊加入集合B 中。
③定義集合Ai ，將B 中每一條邊bi 的相鄰三角面片加入Ai 中。定義集合Ai_saved ，求集合Ai 中的各三角面片與邊bi 所從屬的三角面片之間的二面角，將夾角大于 155° 的三角面片加入集合Ai_saved 中。
④對B 中的所有的邊進行以上操作，得到集合AI_SAVED ，即所有Ai_saved 的并集。在集合AI_SAVED 中求對應最大二面角的三角面片AI_SAVED_MAX ，并從集合AI_SAVED 中剔除AI_SAVED_MAX 。
⑤檢測三角面片AI_SAVED_MAX 的拓撲關系是否正常。如果正常，則進行步驟⑥；若不正常，則返回步驟④。
⑥更新集合S 和B。將三角面片AI_SAVED_MAX 加入集合S 中，同時將不在集合B 中的AI_SAVED_MAX 的邊加入集合B，并剔除已在集合B 中的AI_SAVED_MAX 的邊。接著進行步驟③并以此往下進行，直到集合AI_SAVED 為空。
在上述過程中，155° 的限制是為了避免區域生長過程中三角面片回溯生長；檢測三角面片的拓撲關系正常的判斷標準為：待檢測面片AI_SAVED_MAX 與S 中所有的三角面片之間不會產生空間上的壓疊。
通過此算法從T 中篩選出來的三角面片集合S 即為融合得到的三維牙齒網格模型，但是會出現空洞，因此我們在局部孔洞處再次使用以上算法進行修補。
5 實驗結果與分析
5.1 分割重構結果
由某患者激光掃描圖像分割重構的牙冠模型如圖 5 左圖所示，由該患者口腔 CT 圖像分割重構的牙齒模型如圖 5 右圖所示。
圖5
某患者激光掃描圖像與口腔 CT 圖像分割重構結果
Figure5.
Results of segmentation and reconstruction from laser-scanned images and dental CT images of one subject
5.2 模型配準結果
激光掃描圖像重構的牙冠模型與口腔 CT 圖像重構的牙冠模型配準結果如圖 6 所示。其中，紫色為激光掃描圖像重構的牙冠模型，藍色為口腔 CT 圖像重構的牙冠模型。視覺上，配準后兩模型較好地重疊在一起。
圖6
不同視角下的牙冠配準結果 藍色：口腔 CT 圖像分割重構的牙冠模型；紫色：激光圖像分割重構的牙冠模型
Figure6.
Crown registration results from two different views Blue: crown model based on images; Purple: crown model based on laser-scanned images
為量化配準結果，本文對配準后激光掃描圖像重構的牙冠模型上每一個頂點 在口腔 CT 圖像重構的牙冠模型上搜索離最近的頂點qj ，并計算兩者之間的距離，以此計算平均配準誤差和最大點對的配準誤差（單位：mm）[18-19]。平均配準誤差和最大點對的配準誤差求解公式為：
| $ AD = \frac{{\sum\limits_{i = 1}^N {||{{({{p'}_{final}})}_i} - {q_i}|{|^2}} }}{N} $ |
| $ MD = \max (||{({p'_{final}})_i} - {q_i}||) $ |
其中，N 為激光牙冠模型頂點的個數。
經過對 5 位成年患者的圖像進行測試，測試結果為：AD=（0.26±0.03）mm，MD=（0.54±0.25）mm。 其中，某一位患者的配準結果如圖 7 所示，配準結果中存在少數點對之間距離較大，其原因主要為：①本文所用 CT 圖像的空間分辨率較低（0.30 mm）；②口腔 CT 圖像中牙齒分割結果存在一定誤差。
圖7
不同視角下的配準誤差圖
Figure7.
Error map of registration result from two different views
5.3 模型融合結果
本文模型融合主要發生在兩模型的牙頸處，牙根數據的完整性不會影響整個實驗方法的驗證。為降低實驗復雜度又不失合理性，我們采用部分口腔 CT 牙根數據參與實驗，結果如圖 8 所示， 重構結果可以用于臨床正畸治療中輔助醫師制定正畸治療方案以及矯治器輔助設計。
圖8
不同視角下的融合結果
Figure8.
Fusion results from two different views
6 結論
完整牙齒三維模型為正畸醫師進行診療提供了必要的輔助信息。目前，牙齒三維模型主要通過口腔 CT 圖像分割與重構得到。一方面，口腔 CT 圖像重構的模型精度較低，不適用于無托槽隱形矯治器的輔助設計；另一方面，在正畸治療的不同周期對患者進行重復口腔 CT 掃描會給患者帶來輻射傷害。為此，從提高口腔 CT 圖像重構的完整三維牙齒模型的精度和減少不必要的 CT 掃描輻射的角度出發，本文提出一種基于口腔 CT 圖像與激光掃描圖像融合重構牙齒三維模型的方法。該方法首先采用區域增長算法和改進的基于標記的快速分水嶺算法從激光掃描圖像重構的網格模型中分割出獨立的牙冠模型；采用混合水平集算法和 MC 算法從口腔 CT 圖像中分割重構出完整的牙齒三維模型。然后采用 PCA 算法和 ICP 算法將兩模型的牙冠區域進行配準。最后，采用 DBRG 算法將激光掃描圖像重構的牙冠與口腔 CT 圖像重構的牙根融合，建立完整牙齒三維模型。實驗結果驗證了該方法能有效地利用激光掃描圖像中牙冠信息與口腔 CT 圖像中牙根信息融合來重建完整牙齒的三維模型。利用該方法重建的牙齒模型牙冠部位精度更高，可用于無托槽隱形矯治器的輔助設計。同時，該方法只需在正畸治療前對患者進行一 次口腔 CT 掃描，在治療過程中進行激光掃描即可得到患者各治療周期完整的牙列模型，從而通過減少 CT 掃描次數來降低給患者帶來的輻射傷害。
