部分頜骨缺損、囊腫和需要種植牙的患者，在口腔診斷過程中醫生為了觀察其完整的牙列信息，需要患者拍攝全景X光片，或者需要醫生手動繪制牙弓線用以生成全景圖。針對患者拍攝全景X光片將增加額外負擔以及醫生手動分割牙弓線耗時過長的問題，本文基于錐形束計算機斷層掃描（CBCT）提出了一種全自動全景圖重建方法，使用V型網絡（VNet）將牙齒與背景預分割，生成對應的二值圖，再利用貝塞爾曲線定義最佳牙弓曲線，生成口腔全景圖。此外，本研究還解決了將牙頜骨誤識別為牙弓、生成的牙弓線覆蓋牙弓區域不全面和魯棒性低等問題，可為今后牙科診斷智能化以及提高醫生工作效率奠定理論基礎。
引用本文： 侯昌鵬, 朱赴東, 章高華, 呂震, 劉云峰, 朱偉東. 基于預分割和貝塞爾函數的口腔全景圖重建方法. 生物醫學工程學雜志, 2023, 40(5): 894-902. doi: 10.7507/1001-5515.202302036 復制
0 引言
隨著居民生活水平和維護口腔健康意識的不斷提高，牙科領域的手術量也在逐年上漲。研究報道，世界上近90%的人患有不同程度的牙齒疾病，不同年齡階層的人都存在牙科治療的需求[1]。目前，在牙科治療的臨床實踐中，通常需要不同類型的醫學成像技術來輔助診斷。例如X射線成像、三維（three dimensional，3D）口腔內掃描圖像和錐形束計算機斷層掃描圖像（cone beam computerized tomography，CBCT）等技術均可用于輔助診斷、手術和制定治療計劃。在所有成像技術中，CBCT是提供完整牙齒和牙槽骨特征信息的唯一成像方式，常用于頜面外科腫瘤、囊腫、牙體牙髓病、顳下頜關節病和骨髓炎等的診斷，還可輔助口腔種植牙手術的術前指導以及阻生齒術前的診斷及指導[2-6]。此外，口腔全景圖則是展示患者牙齒陣列信息和周邊組織情況的重要視圖。口腔全景圖能提供包含牙齒和頜骨的全局視角，更有利于展示牙齒的細節特征和整體走向，幫助醫生更加準確地判斷牙齒的具體健康狀況，通常用于評估患者的智齒、齲齒發展狀態和規劃植牙手術方案。對比而言，原始的CBCT圖像，無法顯示全局的牙齒排列情況，因此從CBCT圖像中高效精準地生成口腔全景圖在數字牙科領域研究中至關重要。
目前，全景圖的獲取方式有兩種：一是拍攝三維的電子計算機斷層掃描（computerized tomography，CT）或CBCT圖像后，醫生手動勾畫牙弓線來生成口腔全景圖。二是拍攝口腔全景片，這是口腔外的一種X射線成像膠片[7]。在實際牙科診斷過程中，患有頜骨缺損、囊腫等情況的患者大部分都需要拍攝CBCT圖像，但再次拍攝口腔全景片，會給患者造成額外的輻射負擔。針對上述問題，本文基于CBCT圖像提出了一種先分割牙齒組織后定義牙弓線的口腔全景圖生成方法。
在醫學圖像計算領域，許多學者針對牙齒的自動分割，開展了大量研究。但由于牙齒排列緊湊以及牙齒與頜骨組織密度相近的特點，如何從醫學影像中精準地分割牙齒的邊界與體積，一直都是困擾學者的難題。傳統的口腔圖像分割有基于區域的、基于閾值的、基于聚類的、基于邊緣跟蹤的，以及基于分水嶺算法的等眾多方法[8-12]。但是，由于傳統方法只在單一場景下表現良好，因此無法直接應用于復雜場景[13-15]。基于卷積神經網絡（convolutional neural networks，CNN）的深度學習方法，以任務為導向在大規模數據中進行學習，具備提取代表性和預測性特征的強大能力，在各個領域展現了良好的應用前景。
近年來，深度學習在計算機視覺和醫學圖像計算方面取得了巨大成功，涌現了一系列可用于牙齒和骨骼結構分割的神經網絡[16]。田素坤等[17]提出了一種基于多層三維CNN的多層級檢測方法，用于分割和識別牙齒類型，在其自采集的牙科數據集上分割精度保持在85.8%以上，單牙分割的準確性為89.8%。Koch等[18]證明了U型網絡（U-Net）可利用數據對稱性、網絡集成和測試時間增強來提高分割性能，在公開的森林（Silva）數據集上使用戴斯相似系數（Dice similarity coefficient，DSC）評估，達到了93.4%的精度[19]。Cui等[20]利用生成性對抗網絡為牙齒分割開發綜合語義信息，在公開的肺部結節行號數據庫（lung nodule database，LNDb）數據集上的交并比（intersection over union，IOU）為90.4%。
不少學者在自動化生成口腔全景圖方面同樣做出了卓越貢獻。Bae等[21]通過對牙齒分割和咬合平面的確定，從CBCT圖像的咬合平面中生成牙齒與偏移平面之間的交點。通過三次B樣條曲線近似擬合交點曲線，最終確定牙弓形態生成全景圖[22]。Yun等[23]提出了一種高對比度全景圖重建方法，首先從包含牙齒的軸向切片生成的軸向最大強度投影（maximum intensity projection，MIP）圖像中檢測牙弓厚度，減少全景圖像重建中的非感興趣組織[24]；然后，在圖像合成階段，提出了一種新的合成算法，減少非感興趣組織對圖像對比度的影響；最后，將圖像增強算法應用于合成圖像，以提高最終全景圖像的細節對比度。Amorim等[25]提出了一種使用CT圖像自動化重建二維全景牙科圖的方法，該方法使用貝塞爾曲線和優化算法，最終可得到患者下頜區域的縱向視圖，從而避免拍攝額外的全景X光片。Amorim等[25]將該方法應用于5名患者，雖然其中一些患者存在牙齒缺失的情況，但試驗最終仍然獲得了對比度良好的平滑全景圖像。
綜上所述，基于醫學圖像的自動化重建口腔全景圖還處于實驗探索階段，存在自動化程度低、魯棒性差等缺點。本文基于CBCT圖像提出了一種先分割牙齒組織，后定義牙弓曲線的全景圖重建方法。將從以下內容展開研究：① 基于V型網絡（VNet）模型，針對醫學分割中小體積目標漏檢的問題，選用DSC來構成損失函數，有效解決上下采樣過程造成的特征丟失。② 采用預分割后的二值圖像堆，定位牙弓線繪制圖層，并與分水嶺算法處理的圖像進行融合，獲得包含牙弓與頜骨的二值圖。③ 優化提取牙弓區域骨架算法，解決骨架曲線分叉導致的牙弓線提前終止問題。④ 基于形態學算子、分割和濾波方法的使用，改善了過去口腔全景圖牙齒重疊、牙弓區域覆蓋不全、圖像模糊等問題。
1 方案整體設計
口腔全景圖重建方法，主要分為三個階段：
（1）牙齒分割，對CBCT圖像使用VNet模型進行體積分割，以大量殘差連接代替普通的卷積堆疊，提取感興趣的牙齒區域[26]。
（2）牙弓骨架曲線提取，從分割后的二值圖像堆中，定位牙弓區域的繪制圖層，再提取出能表達牙弓走勢的骨架曲線。
（3）定義牙弓線，隨機選擇控制點，進行貝塞爾曲線迭代，定義最佳牙弓曲線，生成口腔全景圖。
1.1 牙齒預分割
首先將輸入的三維CBCT圖像裁剪為256 × 256 × 256的大小，再對裁剪后的原始CBCT圖像進行窗寬和窗位截斷，將亨氏單位（hounsfield units，Hu）在0～2 500內的值，進行均值為0，方差為1的歸一化處理，對其余像素點進行截斷處理。
本文選用VNet模型作為主干網絡分割牙齒與背景區域，網絡結構如圖1所示，使用不同大小的卷積核進行卷積，并使用累加像素特征方式對圖像進行編解碼。針對牙齒小目標容易漏檢的問題，網絡使用跨層特征重構方法，在圖1中通過水平連接表示，在訓練過程中捕捉可能在編碼路徑中丟失的細粒度細節，以此提高了預測的質量。分割網絡的輸出大小與原始圖像相同，通過柔性最大（softmax）函數進行二分類，輸出圖像上每個像素點歸屬牙齒的概率。
圖1
VNet網絡結構圖
Figure1.
VNet network structure diagram
口腔CBCT圖像中，感興趣的牙齒組織通常只占較小的區域，這很容易導致在訓練中陷入局部極小值，從而產生一個預測結果偏向背景的網絡。并且，如果患者體內存在金屬植入體，會導致圖像產生嚴重的金屬偽影，令牙齒組織更加難以檢測。因此，在過去的牙齒分割方法中，存在金屬偽影的CBCT圖像常出現漏檢牙齒或僅檢測到半顆牙等情況。
針對分割任務中經常出現的上述問題，本文基于DSC來構成戴斯損失函數（Dice loss），在圖像正負樣本比例失衡的場景下Dice loss有著不錯的解決性能。該損失函數中像素的損失值不僅和預測值有關，還和周圍的鄰近像素點相關，并且對訓練結果的影響不會隨著目標圖像的大小產生變化。將其作為損失函數，可以在訓練過程中更注重對目標區域的檢測。Dice loss的計算公式如式(1)所示：
 |
式中，LDice代表Dice loss，?i代表預測分割體積，yi代表真實輪廓標簽，N代表參與訓練的像素點數，i代表不同像素點的編號。
但只訓練Dice loss容易不穩定，尤其是小目標的情況下，而且極端場景下會導致梯度飽和現象。本文結合了Dice loss和二元交叉熵損失(cross entropy loss)構成最終的損失函數，其損失如式(2)～式(3)所示：
 |
 |
式中，Lcross代表二元交叉熵損失函數，L代表最終的損失函數。
1.2 牙弓區域的骨架曲線提取
將最大灰度投影圖作為牙弓線繪制層，計算過程如圖2所示。將分割結果轉為二值圖像堆作為輸入，構建三維矩陣。選擇三維矩陣中含有牙齒體素最多的切片，其所在位置具有牙齒尖端的切片，并且橫斷面中含有頜骨區域的有牙區域面積最大，將該圖層作為牙弓線繪制圖像層。
圖2
最大灰度投影圖計算過程
Figure2.
The calculation process of the maximum grayscale projection map
取出原始圖像堆中位于牙弓線繪制圖層的切片，經過分水嶺算法處理，得到了包含頜骨特征的二值圖。從預分割后的二值圖像堆中挑出位于牙弓線繪制圖層的切片，與分水嶺算法分割后的二值圖進行特征融合，此處的特征融合使用二進制中的“或”操作，將兩幅二值圖上值為1的像素進行保留，得到包含牙齒和頜骨特征的二值圖，整體流程如圖3所示。
圖3
提取全景圖流程圖
Figure3.
Flow chart of panorama extraction
部分患者存在缺牙現象，會導致圖像存在非連通區域。為了打通非連通區域，本文使用形態擴張算子對牙齒和頜骨進行膨脹，膨脹算子設為50，膨脹后使用高斯濾波器進行平滑處理，濾波器的標準差為2，高斯核大小為3×3。得到牙齒與牙頜骨的最大連通區，膨脹后的圖像可打通連續缺牙四顆的非連通區域。經過上述初步處理得到了平滑的聯通區，使用中軸算法進行骨架曲線的提取，生成一條穿越所有牙齒的骨架曲線，再通過貝塞爾曲線定義牙弓線，如圖3中紅色曲線所示，即可生成口腔全景圖[27]。
通過上述方法得到的牙弓骨架曲線，并不是一條光滑的曲線，并且部分圖像生成的骨架曲線還存在分叉等情況，如圖4所示的紅圈部分。
圖4
分叉骨架曲線
Figure4.
Bifurcated skeleton curve
為了消除分叉，考慮到二值圖是由像素組成的連續矩陣，且值只有0和1，本文采用深度優先搜索算法（depth first search，DFS），遍歷骨架，提取長度最大的連續曲線，具體步驟如下：
（1）從左上角的起始點出發，每次朝著八個方向（上、下、左、右、右上、右下、左上、左下）前進，并使用額外的空間存儲走過的像素點，經過的點不再重復，以此實現剪枝。
（2）到達新的點位后判斷像素值是否為0，為0，則需返回上一個像素點，朝另一個方向前進（像素為0則表述當前點不包含在骨架曲線內）；為1，則判斷當前經過的路程是否為全局最長，通過判定就替換全局最長路徑。
（3）遍歷所有像素值為1的點位后，輸出更新后的骨架曲線。
1.3 貝塞爾曲線定義牙弓線
消除骨架曲線的分叉后，曲線還是由不規則線段構成，而牙弓線的建立必須符合以下要求：
（1）牙弓線必須穿越牙弓區域的中心，能如實地反映患者牙弓的形狀。
（2）牙弓線必須具備連續的性質，能準確地計算牙弓線上任意點的坐標，這樣才能生成不失真的口腔全景圖。
（3）口腔全景圖是牙齒陣列沿垂直于上、下頜骨所在曲面方向投影生成的二維展開圖，所以牙弓線必須是能夠計算任意位置法向量的光滑曲線。
（4）由于CBCT圖像的二維切片之間是離散化的，在生成口腔全景圖時，需要在牙弓線上進行等弧長取點采樣的操作，所以定義牙弓線的曲線必須滿足等弧長劃分的要求。
為實現上述條件，必須采用能夠描述牙弓走勢的數學曲線模型進行擬合。考慮到貝塞爾曲線的特點，本文決定采用貝塞爾曲線來定義牙弓線，貝塞爾曲線的特點如下：
（1）貝塞爾曲線，是開始于起點結束于終點的曲線，其余中間的點都是作為曲線形狀和曲率的控制點。
（2）貝塞爾曲線的所有控制點都在曲線上。
（3）一條貝塞爾曲線可切分成兩段或多段子曲線，每一段子曲線仍然屬于貝塞爾曲線。
（4）貝塞爾曲線中某一個控制點的位置發生改變，整條曲線的形狀都會發生變化。由于牙弓線是受到整體牙弓特征的影響，所以貝塞爾曲線的這條性質對于定義牙弓線尤為重要。
下頜角的左右兩端，可以作為貝塞爾曲線的起點和終點，位于牙齒上的點位作為約束貝塞爾曲線的控制點，貝塞爾曲線的定義如式(4)所示：
 |
式中，Pi代表貝塞爾曲線的控制點；t代表曲線參數，范圍從0～1；曲線越接近1，得到的點就越接近曲線的末端；在n次貝塞爾曲線中，曲線都是連續的；B(t)代表貝塞爾曲線的函數方程；i代表曲線控制點的編號。
為了提高生成全景圖的質量，生成的貝塞爾曲線應該盡可能地貼近骨架曲線。那么就需要優化初始的貝塞爾曲線，去定義一條盡可能接近骨架的曲線，以最小化骨架和貝塞爾曲線之間的距離。本文采取的方式是通過對控制點集合不斷迭代的方式去最小化貝塞爾曲線和骨架之間的距離，如式(5)所示：
 |
式中，f（Pi）代表骨架曲線與貝塞爾曲線距離的計算函數，Γa是骨架曲線長度坐標，Φ表示曲線骨架和貝塞爾曲線之間的距離。min表示最小化，dΓ表示骨架曲線的微分算子。
由于圖像表示采用離散表示法，因此不能應用于連續函數。因此，優化函數必須適用于離散分析，離散的公式如式(6)所示：
 |
式中，bi是貝塞爾曲線的第i個控制點，ei是骨架的第i點，這是離bi最近的點，m是計算離散點的個數，想要更全局的優化曲線，m的取值要盡可能的大。為了獲得優化的貝塞爾曲線，本文使用順序最小二乘規劃（sequential least squares programming，SLSQP）算法，選擇n個隨機控制點不斷迭代，得到優化后的貝塞爾曲線[28]。
需要特別注意的是，第一個和最后一個控制點是固定的，因為貝塞爾曲線經過處于牙頜骨位置的起始點和終點。例如，如果貝塞爾曲線有5個控制點，則有3個點是可調節的變量。
2 實驗與數據分析
2.1 數據與評價指標
本研究試驗數據在浙江大學醫學院附屬口腔醫院，使用CBCT影像設備進行采集，共采集150名成年人的口腔CBCT數據。將其隨機劃分成6：2：2的訓練集、測試集（用于驗證牙齒預分割精度）以及驗證集（用于評估口腔全景圖重建方法）。患者均為經過醫生篩選，并符合試驗條件的受試者。試驗研究已通過浙江大學醫學院附屬口腔醫院倫理委員會批準［（2020年度）倫審研第（37）號］，所有受試者均由其本人在簽署知情同意書的情況下，自愿參加本試驗。
本文主要采用DSC、敏感度（sensitivity，SEN）、IOU和對稱位置的平均表面距離（average symmetric surface distance，ASSD）四種指標進行算法評價，各評價指標定義如式（7）～式（10）所示：
 |
 |
 |
 |
式中，真陽性（true positive，TP）表示模型正確分割出牙齒標簽的樣本數；真陰性（true negative，TN）表示模型正確分割出背景區域的樣本數；假陽性（false positive，FP）表示模型將背景誤分為牙齒對象的樣本數；假陰性（false negative，FN）表示模型將牙齒區域誤認為背景的樣本數。式（10）中S(A)代表模型A中的表面體素集合，S(B)代表模型B中的表面體素集合，SA表示S(A)內任意體素，SB表示S(B)內任意體素，d(SA, S(B))表示集合S(A)內的任意體素到集合S(B)的最短距離，d(SB, S(A))表示集合S(B)內的任意體素到集合S(A)的最短距離。
2.2 實驗結果與分析
模型訓練環境為計算機操作系統Windows 10.0（Microsoft Inc, 美國），使用深度學習網絡框架Pytorch 1.8.0（Meta Platforms Inc, 美國）。主要硬件包括中央處理器Intel i9 11900K（Intel Corporation, 美國）以及顯示加速卡GTX3090（NVIDIA Corporation, 美國），使用小批量隨機梯度下降法來訓練模型，初始學習率為1 × 10?3，衰減率為1 × 10?4，隨著迭代次數增多，不斷地降低學習速率，動量設置為0.9，訓練周期為10 000個迭代周期。在30組測試集上的分割評價結果如表1所示。
表1
分割效果評價表
Table1.
Segmentation effect evaluation table
分割結果如圖5所示，本文對牙齒的預分割為語義分割，精度較高，能完整地將牙齒組織和背景區分開，滿足后續牙弓線生成的使用條件。
圖5
牙齒分割結果
Figure5.
Tooth segmentation results
為了評價貝塞爾曲線定義牙弓線的效果，對含30份CBCT患者圖像的驗證集，進行牙弓線點集合的成本函數評估，如表2所示列出了成本函數評估結果。
表2
成本函數評估結果
Table2.
Cost function evaluation results
對于采用3個控制點定義的貝塞爾曲線，成本函數最大不超過16.32 mm，平均（4.21 ± 2.54）mm。對于采用7個控制點定義的貝塞爾曲線，成本函數總是小于5.32 mm，平均（1.67 ± 0.87）mm。在這種情況下，貝塞爾曲線有時并沒有覆蓋到所有的牙齒。但是貝塞爾曲線在7個控制點下相比于3個控制點更接近藍色的骨架曲線，如圖6所示。對于11個和15個控制點，成本函數的最大偏差為1.95 mm，平均（0.94 ± 0.46）mm和1.53 mm，平均（0.73 ± 0.31）mm。在這兩種情況下，貝塞爾曲線能覆蓋所有牙齒且非常接近骨架曲線。這些結果表明，11個和15個控制點，比3個和7個控制點更好。
圖6
牙弓線對比實驗圖
Figure6.
Comparison experiment diagram of dental arch lines
本研究選擇來自驗證集中具有代表性的五位患者的CBCT圖像，其中兩位患者下頜部分經過修復手術，下頜骨部分缺失；一位患者則經歷過種植牙手術，存在種植體，而金屬種植體會在CBCT圖像上產生金屬偽影；最后兩位患者經歷過齲齒修復手術。對于每位患者，均使用3、7、11、15個控制點進行對比測試。
如圖6所示，從3個控制點的貝塞爾曲線，可以觀察到橙色的貝塞爾曲線在一些重要的牙齒區域并沒有與牙齒重合，也就是漏掉了部分牙齒，未覆蓋的牙齒區域在圖6中由紅色虛線框標出。
貝塞爾曲線重建的全景圖對比如圖7所示。由3個控制點獲得的全景圖像缺少牙齒特征信息，五位患者門牙區域的信息均不明確，在圖7中由紅色虛線框標出。控制點增加到7個點，所有的牙齒都出現在全景圖中，然而還是存在圖像模糊的問題，牙根部位的細節難以觀察。實驗結果表明，對于11個和15個控制點，生成的全景圖像非常相似，且11個和15個控制點的全景圖像比3個和7個控制點的全景圖像更加柔和，顯示的細節也更多。
圖7
口腔全景對比圖
Figure7.
Comparison of oral panorama images
對于部分下頜骨缺失的患者，本方法仍然能夠找到牙齒尖端切片所在的位置，并定位剩余頜骨區域，將所有牙齒展示在口腔全景圖中。對于存在種植體的患者以及經歷齲齒修復的患者，本文方法均能較好地克服金屬植入物造成的偽影影響，繪制符合要求的牙弓線。
本文相比于Amorim等[25]提出的牙弓線提取算法更能適用復雜圖像的情況，Amorim等[25]在提取牙弓線繪制層時，采用的方式是以Hu值1 500為界限，對圖像進行二值化，然后將有牙面積最大的切片作為繪制層。但是簡單地根據Hu值作為判別牙齒組織的條件存在較大缺點，例如會將牙槽骨部分也判斷為有牙區域，在選擇繪制層的時候有可能選擇存在牙槽骨的切片層，如圖8所示。本文使用VNet模型預分割了牙齒與背景，準確地獲取有牙面積最大的圖層，解決了該問題。
圖8
牙弓線繪制圖層對比
Figure8.
Comparison of drawing layers of dental arch lines
在骨架曲線提取上，本文的算法也更具優勢。Amorim等[25]使用的骨架提取算法無法解決局部分叉的情況，也難以更好地表達牙弓的形狀特征，而且在進行點采樣時，采樣曲線只從起始點到分叉的末尾即結束，導致生成的牙弓線未經過所有牙齒，生成的全景圖不完善。如圖9所示，文獻[25]方法獲得的牙弓線并未覆蓋所有的牙弓，導致生成的全景圖僅覆蓋了1顆右側邊緣的牙齒。
圖9
算法效果對比圖
Figure9.
Algorithm effect comparison chart
3 結論
針對口腔領域手動勾畫牙弓線，分割全景圖耗時過長的問題，本文提出了一種全自動重建全景圖的方法，采用預分割牙齒的方法定位牙弓線繪制圖層，解決了準確定位牙弓線繪制層、缺牙情況下無法擬合牙弓以及骨架曲線分叉的問題。實驗結果表明，該方法非常有效且魯棒性很高，對于存在頜骨缺失、金屬植入物、缺牙等情況均能生成符合要求的口腔全景圖，可以作為牙科領域的重要輔助工具。且本文方法的計算成本低，對于硬件的要求不高，可為每位患者重建一次全景圖。本文研究結果表明，11個和15個控制點的貝塞爾曲線可以作為生成全景牙齒圖像的選用策略，有望為口腔科醫生提供相關醫學診斷參考。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：侯昌鵬負責論文撰寫，構思設計實驗；朱赴東負責課題的構思與設計；章高華負責課題背景調研，參與實驗；呂震負責指導實驗與儀器設計；劉云峰負責文章主持和計劃安排，指導實驗；朱偉東負責仿真與數據分析指導及論文審閱修訂。
倫理聲明：本研究通過了浙江大學醫學院附屬口腔醫院醫學倫理委員會的審批[批文編號：(2020年度)倫審研第(37)號]。
0 引言
隨著居民生活水平和維護口腔健康意識的不斷提高，牙科領域的手術量也在逐年上漲。研究報道，世界上近90%的人患有不同程度的牙齒疾病，不同年齡階層的人都存在牙科治療的需求[1]。目前，在牙科治療的臨床實踐中，通常需要不同類型的醫學成像技術來輔助診斷。例如X射線成像、三維（three dimensional，3D）口腔內掃描圖像和錐形束計算機斷層掃描圖像（cone beam computerized tomography，CBCT）等技術均可用于輔助診斷、手術和制定治療計劃。在所有成像技術中，CBCT是提供完整牙齒和牙槽骨特征信息的唯一成像方式，常用于頜面外科腫瘤、囊腫、牙體牙髓病、顳下頜關節病和骨髓炎等的診斷，還可輔助口腔種植牙手術的術前指導以及阻生齒術前的診斷及指導[2-6]。此外，口腔全景圖則是展示患者牙齒陣列信息和周邊組織情況的重要視圖。口腔全景圖能提供包含牙齒和頜骨的全局視角，更有利于展示牙齒的細節特征和整體走向，幫助醫生更加準確地判斷牙齒的具體健康狀況，通常用于評估患者的智齒、齲齒發展狀態和規劃植牙手術方案。對比而言，原始的CBCT圖像，無法顯示全局的牙齒排列情況，因此從CBCT圖像中高效精準地生成口腔全景圖在數字牙科領域研究中至關重要。
目前，全景圖的獲取方式有兩種：一是拍攝三維的電子計算機斷層掃描（computerized tomography，CT）或CBCT圖像后，醫生手動勾畫牙弓線來生成口腔全景圖。二是拍攝口腔全景片，這是口腔外的一種X射線成像膠片[7]。在實際牙科診斷過程中，患有頜骨缺損、囊腫等情況的患者大部分都需要拍攝CBCT圖像，但再次拍攝口腔全景片，會給患者造成額外的輻射負擔。針對上述問題，本文基于CBCT圖像提出了一種先分割牙齒組織后定義牙弓線的口腔全景圖生成方法。
在醫學圖像計算領域，許多學者針對牙齒的自動分割，開展了大量研究。但由于牙齒排列緊湊以及牙齒與頜骨組織密度相近的特點，如何從醫學影像中精準地分割牙齒的邊界與體積，一直都是困擾學者的難題。傳統的口腔圖像分割有基于區域的、基于閾值的、基于聚類的、基于邊緣跟蹤的，以及基于分水嶺算法的等眾多方法[8-12]。但是，由于傳統方法只在單一場景下表現良好，因此無法直接應用于復雜場景[13-15]。基于卷積神經網絡（convolutional neural networks，CNN）的深度學習方法，以任務為導向在大規模數據中進行學習，具備提取代表性和預測性特征的強大能力，在各個領域展現了良好的應用前景。
近年來，深度學習在計算機視覺和醫學圖像計算方面取得了巨大成功，涌現了一系列可用于牙齒和骨骼結構分割的神經網絡[16]。田素坤等[17]提出了一種基于多層三維CNN的多層級檢測方法，用于分割和識別牙齒類型，在其自采集的牙科數據集上分割精度保持在85.8%以上，單牙分割的準確性為89.8%。Koch等[18]證明了U型網絡（U-Net）可利用數據對稱性、網絡集成和測試時間增強來提高分割性能，在公開的森林（Silva）數據集上使用戴斯相似系數（Dice similarity coefficient，DSC）評估，達到了93.4%的精度[19]。Cui等[20]利用生成性對抗網絡為牙齒分割開發綜合語義信息，在公開的肺部結節行號數據庫（lung nodule database，LNDb）數據集上的交并比（intersection over union，IOU）為90.4%。
不少學者在自動化生成口腔全景圖方面同樣做出了卓越貢獻。Bae等[21]通過對牙齒分割和咬合平面的確定，從CBCT圖像的咬合平面中生成牙齒與偏移平面之間的交點。通過三次B樣條曲線近似擬合交點曲線，最終確定牙弓形態生成全景圖[22]。Yun等[23]提出了一種高對比度全景圖重建方法，首先從包含牙齒的軸向切片生成的軸向最大強度投影（maximum intensity projection，MIP）圖像中檢測牙弓厚度，減少全景圖像重建中的非感興趣組織[24]；然后，在圖像合成階段，提出了一種新的合成算法，減少非感興趣組織對圖像對比度的影響；最后，將圖像增強算法應用于合成圖像，以提高最終全景圖像的細節對比度。Amorim等[25]提出了一種使用CT圖像自動化重建二維全景牙科圖的方法，該方法使用貝塞爾曲線和優化算法，最終可得到患者下頜區域的縱向視圖，從而避免拍攝額外的全景X光片。Amorim等[25]將該方法應用于5名患者，雖然其中一些患者存在牙齒缺失的情況，但試驗最終仍然獲得了對比度良好的平滑全景圖像。
綜上所述，基于醫學圖像的自動化重建口腔全景圖還處于實驗探索階段，存在自動化程度低、魯棒性差等缺點。本文基于CBCT圖像提出了一種先分割牙齒組織，后定義牙弓曲線的全景圖重建方法。將從以下內容展開研究：① 基于V型網絡（VNet）模型，針對醫學分割中小體積目標漏檢的問題，選用DSC來構成損失函數，有效解決上下采樣過程造成的特征丟失。② 采用預分割后的二值圖像堆，定位牙弓線繪制圖層，并與分水嶺算法處理的圖像進行融合，獲得包含牙弓與頜骨的二值圖。③ 優化提取牙弓區域骨架算法，解決骨架曲線分叉導致的牙弓線提前終止問題。④ 基于形態學算子、分割和濾波方法的使用，改善了過去口腔全景圖牙齒重疊、牙弓區域覆蓋不全、圖像模糊等問題。
1 方案整體設計
口腔全景圖重建方法，主要分為三個階段：
（1）牙齒分割，對CBCT圖像使用VNet模型進行體積分割，以大量殘差連接代替普通的卷積堆疊，提取感興趣的牙齒區域[26]。
（2）牙弓骨架曲線提取，從分割后的二值圖像堆中，定位牙弓區域的繪制圖層，再提取出能表達牙弓走勢的骨架曲線。
（3）定義牙弓線，隨機選擇控制點，進行貝塞爾曲線迭代，定義最佳牙弓曲線，生成口腔全景圖。
1.1 牙齒預分割
首先將輸入的三維CBCT圖像裁剪為256 × 256 × 256的大小，再對裁剪后的原始CBCT圖像進行窗寬和窗位截斷，將亨氏單位（hounsfield units，Hu）在0～2 500內的值，進行均值為0，方差為1的歸一化處理，對其余像素點進行截斷處理。
本文選用VNet模型作為主干網絡分割牙齒與背景區域，網絡結構如圖1所示，使用不同大小的卷積核進行卷積，并使用累加像素特征方式對圖像進行編解碼。針對牙齒小目標容易漏檢的問題，網絡使用跨層特征重構方法，在圖1中通過水平連接表示，在訓練過程中捕捉可能在編碼路徑中丟失的細粒度細節，以此提高了預測的質量。分割網絡的輸出大小與原始圖像相同，通過柔性最大（softmax）函數進行二分類，輸出圖像上每個像素點歸屬牙齒的概率。
圖1
VNet網絡結構圖
Figure1.
VNet network structure diagram
口腔CBCT圖像中，感興趣的牙齒組織通常只占較小的區域，這很容易導致在訓練中陷入局部極小值，從而產生一個預測結果偏向背景的網絡。并且，如果患者體內存在金屬植入體，會導致圖像產生嚴重的金屬偽影，令牙齒組織更加難以檢測。因此，在過去的牙齒分割方法中，存在金屬偽影的CBCT圖像常出現漏檢牙齒或僅檢測到半顆牙等情況。
針對分割任務中經常出現的上述問題，本文基于DSC來構成戴斯損失函數（Dice loss），在圖像正負樣本比例失衡的場景下Dice loss有著不錯的解決性能。該損失函數中像素的損失值不僅和預測值有關，還和周圍的鄰近像素點相關，并且對訓練結果的影響不會隨著目標圖像的大小產生變化。將其作為損失函數，可以在訓練過程中更注重對目標區域的檢測。Dice loss的計算公式如式(1)所示：
|
式中，LDice代表Dice loss，?i代表預測分割體積，yi代表真實輪廓標簽，N代表參與訓練的像素點數，i代表不同像素點的編號。
但只訓練Dice loss容易不穩定，尤其是小目標的情況下，而且極端場景下會導致梯度飽和現象。本文結合了Dice loss和二元交叉熵損失(cross entropy loss)構成最終的損失函數，其損失如式(2)～式(3)所示：
|
|
式中，Lcross代表二元交叉熵損失函數，L代表最終的損失函數。
1.2 牙弓區域的骨架曲線提取
將最大灰度投影圖作為牙弓線繪制層，計算過程如圖2所示。將分割結果轉為二值圖像堆作為輸入，構建三維矩陣。選擇三維矩陣中含有牙齒體素最多的切片，其所在位置具有牙齒尖端的切片，并且橫斷面中含有頜骨區域的有牙區域面積最大，將該圖層作為牙弓線繪制圖像層。
圖2
最大灰度投影圖計算過程
Figure2.
The calculation process of the maximum grayscale projection map
取出原始圖像堆中位于牙弓線繪制圖層的切片，經過分水嶺算法處理，得到了包含頜骨特征的二值圖。從預分割后的二值圖像堆中挑出位于牙弓線繪制圖層的切片，與分水嶺算法分割后的二值圖進行特征融合，此處的特征融合使用二進制中的“或”操作，將兩幅二值圖上值為1的像素進行保留，得到包含牙齒和頜骨特征的二值圖，整體流程如圖3所示。
圖3
提取全景圖流程圖
Figure3.
Flow chart of panorama extraction
部分患者存在缺牙現象，會導致圖像存在非連通區域。為了打通非連通區域，本文使用形態擴張算子對牙齒和頜骨進行膨脹，膨脹算子設為50，膨脹后使用高斯濾波器進行平滑處理，濾波器的標準差為2，高斯核大小為3×3。得到牙齒與牙頜骨的最大連通區，膨脹后的圖像可打通連續缺牙四顆的非連通區域。經過上述初步處理得到了平滑的聯通區，使用中軸算法進行骨架曲線的提取，生成一條穿越所有牙齒的骨架曲線，再通過貝塞爾曲線定義牙弓線，如圖3中紅色曲線所示，即可生成口腔全景圖[27]。
通過上述方法得到的牙弓骨架曲線，并不是一條光滑的曲線，并且部分圖像生成的骨架曲線還存在分叉等情況，如圖4所示的紅圈部分。
圖4
分叉骨架曲線
Figure4.
Bifurcated skeleton curve
為了消除分叉，考慮到二值圖是由像素組成的連續矩陣，且值只有0和1，本文采用深度優先搜索算法（depth first search，DFS），遍歷骨架，提取長度最大的連續曲線，具體步驟如下：
（1）從左上角的起始點出發，每次朝著八個方向（上、下、左、右、右上、右下、左上、左下）前進，并使用額外的空間存儲走過的像素點，經過的點不再重復，以此實現剪枝。
（2）到達新的點位后判斷像素值是否為0，為0，則需返回上一個像素點，朝另一個方向前進（像素為0則表述當前點不包含在骨架曲線內）；為1，則判斷當前經過的路程是否為全局最長，通過判定就替換全局最長路徑。
（3）遍歷所有像素值為1的點位后，輸出更新后的骨架曲線。
1.3 貝塞爾曲線定義牙弓線
消除骨架曲線的分叉后，曲線還是由不規則線段構成，而牙弓線的建立必須符合以下要求：
（1）牙弓線必須穿越牙弓區域的中心，能如實地反映患者牙弓的形狀。
（2）牙弓線必須具備連續的性質，能準確地計算牙弓線上任意點的坐標，這樣才能生成不失真的口腔全景圖。
（3）口腔全景圖是牙齒陣列沿垂直于上、下頜骨所在曲面方向投影生成的二維展開圖，所以牙弓線必須是能夠計算任意位置法向量的光滑曲線。
（4）由于CBCT圖像的二維切片之間是離散化的，在生成口腔全景圖時，需要在牙弓線上進行等弧長取點采樣的操作，所以定義牙弓線的曲線必須滿足等弧長劃分的要求。
為實現上述條件，必須采用能夠描述牙弓走勢的數學曲線模型進行擬合。考慮到貝塞爾曲線的特點，本文決定采用貝塞爾曲線來定義牙弓線，貝塞爾曲線的特點如下：
（1）貝塞爾曲線，是開始于起點結束于終點的曲線，其余中間的點都是作為曲線形狀和曲率的控制點。
（2）貝塞爾曲線的所有控制點都在曲線上。
（3）一條貝塞爾曲線可切分成兩段或多段子曲線，每一段子曲線仍然屬于貝塞爾曲線。
（4）貝塞爾曲線中某一個控制點的位置發生改變，整條曲線的形狀都會發生變化。由于牙弓線是受到整體牙弓特征的影響，所以貝塞爾曲線的這條性質對于定義牙弓線尤為重要。
下頜角的左右兩端，可以作為貝塞爾曲線的起點和終點，位于牙齒上的點位作為約束貝塞爾曲線的控制點，貝塞爾曲線的定義如式(4)所示：
|
式中，Pi代表貝塞爾曲線的控制點；t代表曲線參數，范圍從0～1；曲線越接近1，得到的點就越接近曲線的末端；在n次貝塞爾曲線中，曲線都是連續的；B(t)代表貝塞爾曲線的函數方程；i代表曲線控制點的編號。
為了提高生成全景圖的質量，生成的貝塞爾曲線應該盡可能地貼近骨架曲線。那么就需要優化初始的貝塞爾曲線，去定義一條盡可能接近骨架的曲線，以最小化骨架和貝塞爾曲線之間的距離。本文采取的方式是通過對控制點集合不斷迭代的方式去最小化貝塞爾曲線和骨架之間的距離，如式(5)所示：
|
式中，f（Pi）代表骨架曲線與貝塞爾曲線距離的計算函數，Γa是骨架曲線長度坐標，Φ表示曲線骨架和貝塞爾曲線之間的距離。min表示最小化，dΓ表示骨架曲線的微分算子。
由于圖像表示采用離散表示法，因此不能應用于連續函數。因此，優化函數必須適用于離散分析，離散的公式如式(6)所示：
|
式中，bi是貝塞爾曲線的第i個控制點，ei是骨架的第i點，這是離bi最近的點，m是計算離散點的個數，想要更全局的優化曲線，m的取值要盡可能的大。為了獲得優化的貝塞爾曲線，本文使用順序最小二乘規劃（sequential least squares programming，SLSQP）算法，選擇n個隨機控制點不斷迭代，得到優化后的貝塞爾曲線[28]。
需要特別注意的是，第一個和最后一個控制點是固定的，因為貝塞爾曲線經過處于牙頜骨位置的起始點和終點。例如，如果貝塞爾曲線有5個控制點，則有3個點是可調節的變量。
2 實驗與數據分析
2.1 數據與評價指標
本研究試驗數據在浙江大學醫學院附屬口腔醫院，使用CBCT影像設備進行采集，共采集150名成年人的口腔CBCT數據。將其隨機劃分成6：2：2的訓練集、測試集（用于驗證牙齒預分割精度）以及驗證集（用于評估口腔全景圖重建方法）。患者均為經過醫生篩選，并符合試驗條件的受試者。試驗研究已通過浙江大學醫學院附屬口腔醫院倫理委員會批準［（2020年度）倫審研第（37）號］，所有受試者均由其本人在簽署知情同意書的情況下，自愿參加本試驗。
本文主要采用DSC、敏感度（sensitivity，SEN）、IOU和對稱位置的平均表面距離（average symmetric surface distance，ASSD）四種指標進行算法評價，各評價指標定義如式（7）～式（10）所示：
|
|
|
|
式中，真陽性（true positive，TP）表示模型正確分割出牙齒標簽的樣本數；真陰性（true negative，TN）表示模型正確分割出背景區域的樣本數；假陽性（false positive，FP）表示模型將背景誤分為牙齒對象的樣本數；假陰性（false negative，FN）表示模型將牙齒區域誤認為背景的樣本數。式（10）中S(A)代表模型A中的表面體素集合，S(B)代表模型B中的表面體素集合，SA表示S(A)內任意體素，SB表示S(B)內任意體素，d(SA, S(B))表示集合S(A)內的任意體素到集合S(B)的最短距離，d(SB, S(A))表示集合S(B)內的任意體素到集合S(A)的最短距離。
2.2 實驗結果與分析
模型訓練環境為計算機操作系統Windows 10.0（Microsoft Inc, 美國），使用深度學習網絡框架Pytorch 1.8.0（Meta Platforms Inc, 美國）。主要硬件包括中央處理器Intel i9 11900K（Intel Corporation, 美國）以及顯示加速卡GTX3090（NVIDIA Corporation, 美國），使用小批量隨機梯度下降法來訓練模型，初始學習率為1 × 10?3，衰減率為1 × 10?4，隨著迭代次數增多，不斷地降低學習速率，動量設置為0.9，訓練周期為10 000個迭代周期。在30組測試集上的分割評價結果如表1所示。
表1
分割效果評價表
Table1.
Segmentation effect evaluation table
分割結果如圖5所示，本文對牙齒的預分割為語義分割，精度較高，能完整地將牙齒組織和背景區分開，滿足后續牙弓線生成的使用條件。
圖5
牙齒分割結果
Figure5.
Tooth segmentation results
為了評價貝塞爾曲線定義牙弓線的效果，對含30份CBCT患者圖像的驗證集，進行牙弓線點集合的成本函數評估，如表2所示列出了成本函數評估結果。
表2
成本函數評估結果
Table2.
Cost function evaluation results
對于采用3個控制點定義的貝塞爾曲線，成本函數最大不超過16.32 mm，平均（4.21 ± 2.54）mm。對于采用7個控制點定義的貝塞爾曲線，成本函數總是小于5.32 mm，平均（1.67 ± 0.87）mm。在這種情況下，貝塞爾曲線有時并沒有覆蓋到所有的牙齒。但是貝塞爾曲線在7個控制點下相比于3個控制點更接近藍色的骨架曲線，如圖6所示。對于11個和15個控制點，成本函數的最大偏差為1.95 mm，平均（0.94 ± 0.46）mm和1.53 mm，平均（0.73 ± 0.31）mm。在這兩種情況下，貝塞爾曲線能覆蓋所有牙齒且非常接近骨架曲線。這些結果表明，11個和15個控制點，比3個和7個控制點更好。
圖6
牙弓線對比實驗圖
Figure6.
Comparison experiment diagram of dental arch lines
本研究選擇來自驗證集中具有代表性的五位患者的CBCT圖像，其中兩位患者下頜部分經過修復手術，下頜骨部分缺失；一位患者則經歷過種植牙手術，存在種植體，而金屬種植體會在CBCT圖像上產生金屬偽影；最后兩位患者經歷過齲齒修復手術。對于每位患者，均使用3、7、11、15個控制點進行對比測試。
如圖6所示，從3個控制點的貝塞爾曲線，可以觀察到橙色的貝塞爾曲線在一些重要的牙齒區域并沒有與牙齒重合，也就是漏掉了部分牙齒，未覆蓋的牙齒區域在圖6中由紅色虛線框標出。
貝塞爾曲線重建的全景圖對比如圖7所示。由3個控制點獲得的全景圖像缺少牙齒特征信息，五位患者門牙區域的信息均不明確，在圖7中由紅色虛線框標出。控制點增加到7個點，所有的牙齒都出現在全景圖中，然而還是存在圖像模糊的問題，牙根部位的細節難以觀察。實驗結果表明，對于11個和15個控制點，生成的全景圖像非常相似，且11個和15個控制點的全景圖像比3個和7個控制點的全景圖像更加柔和，顯示的細節也更多。
圖7
口腔全景對比圖
Figure7.
Comparison of oral panorama images
對于部分下頜骨缺失的患者，本方法仍然能夠找到牙齒尖端切片所在的位置，并定位剩余頜骨區域，將所有牙齒展示在口腔全景圖中。對于存在種植體的患者以及經歷齲齒修復的患者，本文方法均能較好地克服金屬植入物造成的偽影影響，繪制符合要求的牙弓線。
本文相比于Amorim等[25]提出的牙弓線提取算法更能適用復雜圖像的情況，Amorim等[25]在提取牙弓線繪制層時，采用的方式是以Hu值1 500為界限，對圖像進行二值化，然后將有牙面積最大的切片作為繪制層。但是簡單地根據Hu值作為判別牙齒組織的條件存在較大缺點，例如會將牙槽骨部分也判斷為有牙區域，在選擇繪制層的時候有可能選擇存在牙槽骨的切片層，如圖8所示。本文使用VNet模型預分割了牙齒與背景，準確地獲取有牙面積最大的圖層，解決了該問題。
圖8
牙弓線繪制圖層對比
Figure8.
Comparison of drawing layers of dental arch lines
在骨架曲線提取上，本文的算法也更具優勢。Amorim等[25]使用的骨架提取算法無法解決局部分叉的情況，也難以更好地表達牙弓的形狀特征，而且在進行點采樣時，采樣曲線只從起始點到分叉的末尾即結束，導致生成的牙弓線未經過所有牙齒，生成的全景圖不完善。如圖9所示，文獻[25]方法獲得的牙弓線并未覆蓋所有的牙弓，導致生成的全景圖僅覆蓋了1顆右側邊緣的牙齒。
圖9
算法效果對比圖
Figure9.
Algorithm effect comparison chart
3 結論
針對口腔領域手動勾畫牙弓線，分割全景圖耗時過長的問題，本文提出了一種全自動重建全景圖的方法，采用預分割牙齒的方法定位牙弓線繪制圖層，解決了準確定位牙弓線繪制層、缺牙情況下無法擬合牙弓以及骨架曲線分叉的問題。實驗結果表明，該方法非常有效且魯棒性很高，對于存在頜骨缺失、金屬植入物、缺牙等情況均能生成符合要求的口腔全景圖，可以作為牙科領域的重要輔助工具。且本文方法的計算成本低，對于硬件的要求不高，可為每位患者重建一次全景圖。本文研究結果表明，11個和15個控制點的貝塞爾曲線可以作為生成全景牙齒圖像的選用策略，有望為口腔科醫生提供相關醫學診斷參考。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：侯昌鵬負責論文撰寫，構思設計實驗；朱赴東負責課題的構思與設計；章高華負責課題背景調研，參與實驗；呂震負責指導實驗與儀器設計；劉云峰負責文章主持和計劃安排，指導實驗；朱偉東負責仿真與數據分析指導及論文審閱修訂。
倫理聲明：本研究通過了浙江大學醫學院附屬口腔醫院醫學倫理委員會的審批[批文編號：(2020年度)倫審研第(37)號]。
