TWI571806B - Tooth body image joining method - Google Patents

Description

牙體影像的接合方法

本創作是關於一種牙體影像的接合方法，特別是指可根據複數牙體影像的重覆區域彼此疊合而接合影像的方法。

在建立口腔內牙體的影像時，牙科醫生係會持一探頭伸入患者的口腔中以進行牙體的拍攝。以下顎為例，若要取得整排牙齒的舌側面影像，則牙科醫生需要對下顎的牙齒進行拍攝一系列的影像，例如沿著患者左邊的大臼齒往右邊的大臼齒方向之軌跡拍攝。於拍攝完成後，係由與該探頭連線之電腦裝置接收該些影像，該電腦裝置可分辨出該些影像彼此間的關聯性，將該複數影像接合而得到完整的一下顎牙齒舌側面影像。

然而，由於牙科醫生進行拍照時，手持探頭的拍照角度不盡相同，因此當該複數影像彼此接合後，牙齒輪廓有可能產生不連續的狀況，例如牙縫之線條產生錯位，造成不精確的牙齒影像。

有鑒於習知接合影像作法會產生不連續的狀況，因此本創作的主要目的是提供一種牙體影像的接合方法，用以解決接合影像不連續的問題。

本創作牙體影像的接合方法包含有： 由一探頭取得一第一影像與一第二影像； 由一影像處理裝置接收並儲存該些影像，以計算該第一影像與第二影像之間的重覆區域以及該第二影像的一第一位移量，以利用該第一位移量移動該第二影像； 擷取該第二影像中的齒面特徵點以及於非齒面特徵區域的樣本點Bsi； 計算該第二影像齒面特徵點與樣本點Bsi在第一影像中的對應點Asi； 根據Asi與Bsi建立一旋轉矩陣與該第二影像的一第二位移量，以根據該旋轉矩陣與第二位移量移動該第二影像與第一影像接合。

根據本創作的方法，係根據第一位移量以初步移動該第二影像，再根據該旋轉矩陣與第二位移量以進行第一影像與第二影像彼此相對位置的微調，透過此微調手段可將第一影像與第二影像對齊，避免牙齒輪廓不連續的情況。

請參考圖1所示，執行本創作方法的系統包含有一探頭10與一影像處理裝置20，該影像處理裝置20可為電腦，係連線到該探頭10以進行資料傳輸。圖2為本創作方法的流程圖。

首先該探頭10係取得複數牙體的影像，各影像的三維資料係以空間座標參數所表示(步驟101)。該探頭10可為手持式探頭或固定式探頭，以手持式探頭為例，牙科醫師將探頭10伸入患者的口腔內以進行牙齒的拍攝作業，其中該探頭10內設有一姿態感測器11，該姿態感測器11可為陀螺儀，用以產生該影像當下的空間座標參數，在直角座標系統下，空間座標參數包含有x軸旋轉索引值(yaw)、y軸旋轉索引值(pitch)、z軸旋轉索引值(roll)、x座標(dx)、y座標(dy)與z座標(dz)，是以，每張拍攝的影像都有對應的空間座標參數。以下顎牙齒為例，該探頭10可沿著患者左側的大臼齒、小臼齒、犬齒至門牙方向的軌跡對不同牙齒的舌側面依序拍攝，進而得到複數具有連續關係的牙體影像。

以探頭10拍攝得到複數影像後，由該影像處理裝置20接收並儲存該些影像與對應的空間資料參數(步驟102)，於取得該些影像後，該影像處理裝置20負責將該些影像彼此接合。請參考圖3與圖4所示，本創作以兩個接續拍攝的第一影像31與第二影像32之接合動作為例說明，該兩影像31、32各自拍到相同的牙齒40，其中第一影像31與第二影像32為點雲資料(point cloud)。

該影像處理裝置20可在直角座標系統進行影像接合的運算。在進行接合之前，請參考圖5與圖6，首先該影像處理裝置20計算該第一影像31與第二影像32之間的重覆區域310、320以及第二影像32的一第一位移量(步驟103)。所述重覆區域310、320是指該兩影像31、32所拍到相同物體(即牙齒40)的分佈區域。在直角座標系統中，該第一影像31與第二影像32主要是排列在X軸上，因此若要接合該兩影像31、32，則是令該第一影像31朝第二影像32沿著X軸平移，或令第二影像32朝第一影像31沿著X軸平移，直到第一影像31與第二影像32的重覆區域310、320彼此重疊，本較佳實施例是令第二影像32朝第一影像31沿著X軸平移。

假設平移前該第二影像32包含一參考點，該參考點的座標為Bi，則平移後該參考點新座標為Bnewi，則Bnewi=Bi+Dx，Dx為第二影像32的第一位移量，其中i=1,2,3,…,N_B ，N_B 為第二影像32點雲資料的總點數。Dx是在[D_ini -d,D_ini +d]範圍內搜尋取得，[D_ini -d,D_ini +d]之範圍是由使用者操作該影像處理裝置20而設定，D_ini 為初始位移量，d為搜尋範圍。

關於第一位移量Dx之建立，本創作是根據三維物體(即第一影像31與第二影像32)的深度圖(depth image)，深度圖可以D(x,y)表示之，並計算深度圖的梯度大小(gradient magnitude)，本創作是將梯度沿著Y軸方向投影為例，進而分別計算該第一影像31與第二影像32於X座標的梯度累計值。是以，該兩影像31、32共產生兩個梯度累計值。其中請參考圖5，拍攝第一影像31時，因拍攝死角或唾液反光而導致第一影像31有缺損部分311，缺損部分311係不列入運算，其梯度累計設為0。

本創作深度圖的梯度大小表示為，透過濾波器計算如下，本創作是採用Sobel濾波器：

其中

假設在X座標的梯度累計值表示為：

則該第一影像31與第二影像32的梯度累計值則可分別表示為Acc_A (x)與Acc_B (x)，該第一影像31與第二影像32的梯度累計值相關性可表示如下： Î[D_ini -d,D_ini +d]

該第一位移量為。

請參考圖7所示，即是將該第二影像32根據第一位移量平移後，與第一影像31初步接合的示意圖，由圖5~圖7中可見，該第一影像31與第二影像32的重覆區域310、320已大致重疊，惟牙縫處41(由較密集的點所構成的區域代表牙縫)產生不連續的情形。

為了微調第一影像31與第二影像32的相對位置以達到對齊目的時，首先擷取該第二影像32中的齒面特徵點(步驟104)，接著擷取該第二影像32在非齒面特徵區域的樣本點(步驟105)，如此一來，因為是利用擷取出來的特徵點與樣本點進行後續計算，不是影像中的每一點資料進行計算，除了可以降低運算資料量，還可以提升較具代表性特徵點的比例，以提高運算速度並增進比對的精確度。

在第104步驟中，齒面特徵點為齒面的三維幾何特徵點，即具有明顯幾何變化的特徵點。本創作可利用曲率計算三維幾何變化。以S(p)(shape index)表示幾何特徵，考慮點雲資料中的任一點p，則

S(p)係介於0與1之間，其中k₁ (p)與k₂ (p)分別為該p點的兩個主曲率(principal curvature)，且k₁ (p)³k₂ (p)；S(p)為0時代表p的周邊形狀是球杯狀(spherical cup)，S(p)為1時代表p的周邊形狀是球冠狀(spherical cap)，至於中間值0.5則表示p點為馬鞍點(saddle point)，至於其他的值則介於這些形狀之間。由於牙齒上包含有凹陷處、凸起處與馬鞍處，故S(p)可分別表示牙齒的凹陷處、凸起處與馬鞍處，其中S(p)為1~0.9時代表p的周邊形狀是球冠狀，S(p)為0~0.1時代表p的周邊形狀是球杯狀，S(p)為0.45~0.55時代表p點為馬鞍點，以上S(p)之數值範圍僅是較佳實施例而已，不用以限制本創作。是以，對該第二影像32的點雲資料進行S(p)的計算，可取得該第二影像32特徵點。

於齒面特徵點取完後，進行第105步驟以在非齒面特徵區域擷取樣本點，本創作可係利用均勻取樣手段(uniform sampling)對該第二影像32的非齒面特徵區域擷取樣本點，該非齒面特徵區域係指該第二影像32中齒面特徵點以外的區域。本較佳實施例中，係定義三維座標系的Z軸為深度，因此沿著X軸與Y軸，並於每間隔距離△x與△y進行取樣，即可完成樣本點的擷取動作。

根據第104步驟與第105步驟分別取得齒面特徵點與樣本點後，接著計算第一影像31與第二影像32的對應關係。以第二影像32為例，取其特徵點與樣本點的聯集後，其聯集內包含有一點Bsi，i=1,2,3…N_BS ，N_BS 為第二影像32齒面特徵點數量與樣本點數量的總和，若第一影像31包含有一點Asi，i=1,2,3…N_AS ，N_AS =N_BS ，Asi為Bsi的對應點(corresponding point)，則(Asi, Bsi)形成兩影像的對應關係。因為探頭10在進行口內掃描拍照時，口腔內被遮蔽的地方無法被掃描到，例如從頰側或舌側拍攝時凹陷處易被遮蔽，導致第二影像32的一點Bsi無法與第一影像31對應而產生誤差，因此本創作以內插手段修補第一影像31與第二影像32因遮蔽所產生的空缺處，以解決資料缺失造成的對應困難。本創作之內插手段可以是線性內插(linear interpolation)、三次內插(cubic interpolation)或其他內插方式，所述內插計算為公知常識，再此不贅述。

計算第一影像31與第二影像32的對應關係時，以第二影像32中任一Bsi點為例，將該Bsi與第一影像31的所有點雲資料進行Ak-d tree (Approximated k-d tree, Ak-d tree)演算(可參考Michael Greenspan等人所著的”Approximate K-D tree Search for Efficient ICP”，如附件)，計算Bsi與第一影像31所有點雲資料的距離，進而決定該第一影像31中距離該Bsi最近的點作為對應的Asi，演算時係限制兩對應點Asi、Bsi間的距離在一門檻距離內，藉以忽略距離過遠的對應點，避免影響後續的估算，並節省運算量(步驟106)。前述實施例是先計算第二影像32的特徵點與樣本點後，再計算出與Bsi對應的Asi；於另一實施例中，亦可先計算第一影像31的特徵點與樣本點後，再計算出與Asi對應的Bsi。本創作是以Ak-d tree演算法計算最接近的Asi與Bsi，但不以此為限。

當Asi與Bsi之對應關係建立後，該影像處理裝置20接著進行剛體轉換(rigid transformation)以微調第一影像31與第二影像32的位置(步驟107)，所述剛體轉換包含旋轉與位移之動作。各第一影像31與第二影像32透過合適的旋轉與位移可彼此疊合，而旋轉與位移係透過三維點的互協方差矩陣(cross covariance matrix)M_AB 進行計算，互協方差矩陣如下：其中

為第一影像31齒面特徵點座標與樣本點座標的平均值，為第二影像32齒面特徵點座標與樣本點座標的平均值，前述互協方差矩陣建立後，產生如下的4×4矩陣：

其中tr(M_AB )代表矩陣M_AB 主對角線的和，E(M_AB )最大特徵值(eigenvalue)對應的特徵向量(eigenvector)為使得均方差(mean square error)為最小的旋轉。若E(M_AB )p=λp之關係成立，p與l分別為E(M_AB )的特徵向量與特徵值，可由特徵多項式(characteristic polynomial)的根計算特徵值，再代入E(M_AB )p=lp，解線性方程組求得特徵向量p。本創作以電腦執行數值分析方法，以Householder法將E(M_AB )轉換為三對角線矩陣(tridiagonal matrix)，再以QL演算法取得特徵值與特徵向量，可參考Numerical Recipes in C(ISBN:0 521 43108 5)書中所提到的函式tred2()與tqli()計算特徵向量與特徵值(第108~109頁、第113~115頁、第469~480頁)，以下說明Householder演算流程。

要找到特徵值與特徵向量，需先化簡矩陣，若原始矩陣是對稱矩陣，則可化簡為三對角線矩陣(tridiagonal matrix)，也就是除了主對角線與其上下相鄰的兩對角線之外，元素均為零的矩陣，再以疊代的方式求解。我們使用的化簡方法為Householder法，此法將一個n乘n的對稱矩陣化簡為tridiagonal matrix，首先考慮一Householder矩陣P，此矩陣的形式為P=1-2w×w^T ，w為一個長度為一的向量，也就是|w|² =1，由於P² =1-4w×w^T +4w× (w^T ×w)×w^T =1，所以P^-1 =P;此外，P為對稱矩陣，P^T =P，可得P^-1 =P^T ，所以P為正交矩陣(orthogonal matrix)。

把P改寫為其中

假設向量x為要化簡矩陣A的第一行(column)，令，為單位向量，則

其中為的第一個元素。這樣的結果顯示P矩陣把矩陣A的第一行(column)化簡為第一元素不為零且其餘元素為零的向量。如果略過矩陣A第一行的第一元素，而以其第二元素之後的部分形成向量x來建構Householder矩陣，若以此Householder矩陣建立一n´n矩陣：則

則也有類似的效果，此時第一行第二元素之後均為零，為扣除第一列與第一行所得的矩陣；若()，則可以讓第一列第二元素之後變成零，而得到

此時A₁ 的第一列與第一行的第二元素之後為零，為扣除第一列與第一行所得的矩陣，由於為正交矩陣，所以直接將上式的寫成。接著，重複上述步驟，也就是取矩陣A第二行第三元素之後的部分形成向量來建構Householder矩陣，若以此Householder矩陣建立一n´n矩陣：

則可得前兩行與前兩列除三個主對角線之外均為零的矩陣，重複上述步驟，可建立tridiagonal matrix。

而為方便運算，可重新整理流程，可化簡如下： , =

從A的最後一行開始建Householder矩陣，演算流程如下： 第m次疊代，m=1,2,3,…,n-2，進行下列步驟： 1.     建立, for i=n-m+1=n,n-1,…,3 其中，且前的正負號由的正負號決定，應設定成與的正負號一樣。此外，如果很小，則略過本次的疊代。 2.     取得向量，其中常數3.     建立向量，其中常數4.     有了向量q，可進行矩陣化簡

完成上述步驟之後，就可以將主對角三個對角線以外的元素變成零而得到tridiagonal matrix，而能夠對得到的tridiagonal matrix求取特徵向量，取得的特徵向量再加以反轉換，也就是乘上即可得矩陣A的特徵向量。以本創作而言，上述的矩陣A以E (M_AB )代入執行。

由於E(M_AB )為4×4矩陣，所以特徵向量長度為4，表示為 [q1 q2 q3 q4]，可轉換表示為旋轉矩陣R：

於取得旋轉矩陣R後，利用第一影像31與旋轉後的第二影像32計算一第二位移量T，。是以，該第二影像32的新座標為，並評估鄰近Asi與之誤差是否夠落入一由使用者設定的容許範圍，若誤差落在容許範圍，表示該第一影像31與第二影像32可完整對齊接合，如圖8所示，牙縫處42(由較密集的點所構成的區域代表牙縫)已經沒有錯位的情形，牙縫處42為平滑的線條，則該影像處理裝置20儲存第一影像31與第二影像32的最佳接合位置的計算結果。若鄰近Asi與之誤差落超出容許範圍，係回復執行第106步驟以再次計算旋轉矩陣R與第二位移量T，直到誤差落到容許範圍。或者，若重覆執行第106步驟的次數達到一門檻值，則以產生最小的誤差的該次第二位移量T為最後的第二位移量T。

若欲完成半顎牙齒影像，則可於第一影像31與第二影像32完成接合後，可根據前述步驟進行第二影像32與一連續的第三影像的進行接合，並依此類推，直到完成單顎所有牙齒的接合，得到一精確的半顎牙齒影像。

10‧‧‧探頭
11‧‧‧姿態感測器
20‧‧‧影像處理裝置
31‧‧‧第一影像
310‧‧‧重覆區域
311‧‧‧缺損部分
320‧‧‧重覆區域
32‧‧‧第二影像
40‧‧‧牙齒
41‧‧‧牙縫處
42‧‧‧牙縫處

圖1：本創作之系統方塊示意圖。 圖2：本創作較佳實施例之流程示意圖。 圖3：本創作中第一影像的示意圖。 圖4：本創作中第二影像的示意圖。 圖5：本創作中第一影像的重覆區域示意圖。 圖6：本創作中第二影像的重覆區域示意圖。 圖7：本創作中第一影像與第二影像根據第一位移量接合的示意圖。 圖8：本創作中第一影像與第二影像進一步由第二位移量接合的示意圖。

Claims (6)

1. 一種牙體影像的接合方法，包含：由一探頭接續拍攝以取得一第一影像與一第二影像，該第一影像與該第二影像為具有連續關係的牙體影像；由一影像處理裝置接收並儲存該些影像，以計算該第一影像與第二影像之間的重覆區域以及該第二影像的一第一位移量，以利用該第一位移量移動該第二影像，使該第一影像與該第二影像的重覆區域重疊，其中，該第一位移量係根據該第一影像與該第二影像的梯度大小與梯度累計值計算而得；擷取該第二影像中的齒面特徵點以及於非齒面特徵區域的樣本點Bsi；計算該第二影像的齒面特徵點與樣本點Bsi在該第一影像中的對應點Asi；根據Asi與Bsi建立一旋轉矩陣與該第二影像的一第二位移量，以根據該旋轉矩陣與第二位移量移動該第二影像與第一影像接合。
2. 如請求項1所述牙體影像的接合方法，取得該第一位移量的步驟包含：根據深度圖分別計算該第一影像與第二影像的梯度大小；根據該第一影像與第二影像的梯度大小計算該第一影像與第二影像的梯度累計值；根據該第一影像與第二影像的梯度累計值計算該第一位移量。
3. 如請求項2所述牙體影像的接合方法，在取得齒面特徵點的步驟中，齒面特徵點為齒面的三維幾何特徵點，根據幾何特徵S(p)： p為該第二影像中的任一點，S(p)係介於0與1之間，k1(p)與k2(p)分別為該p點的兩個主曲率(principal curvature)，且k1(p)k2(p)；齒面特徵點係指S(p)為1~0.9、0~0.1或0.45~0.55的點。
4. 如請求項3所述牙體影像的接合方法，在取得樣本點的步驟中，係利用均勻取樣手段(uniform sampling)對該第二影像的非齒面特徵區域擷取樣本點。
5. 如請求項1至4中任一項所述牙體影像的接合方法，產生該旋轉矩陣的步驟包含： 產生一矩陣M_AB，，其中 ，，為該第一影像特徵點座標與樣本點座標 的平均值，為該第二影像齒面特徵點座標與樣本點座標的平均值；產生4×4矩陣E(M_AB)， 對E(M_AB)透過Householder與QL演算法進行特徵值與特徵向量的計算以得出該旋轉矩陣。
6. 如請求項5所述牙體影像的接合方法，該第二位移量為 ，其中R為該旋轉矩陣，該第二影像的特徵點或樣本點與該第一影像接合後的位置表示為。