劉暢 周競(jìng)宇

摘要：在3D打印領(lǐng)域中，針對(duì)找到受損骨骼待修補(bǔ)部分的問(wèn)題上，提出一種利用模板進(jìn)行配準(zhǔn)找出需要修補(bǔ)部分的方法。這種模板可以是相似度高的模型，也可以是病人的鏡像骨。在完成方法過(guò)程中，利用到了最小盒算法和ICP配準(zhǔn)，前者是對(duì)模型的粗配準(zhǔn)，后者是對(duì)模型的精配準(zhǔn);在最后進(jìn)行剔除錯(cuò)誤部分和平滑，生成打印模型。實(shí)驗(yàn)表明，該方法能完整精確打印出待修補(bǔ)處三維模型，與缺損模型平滑連接。

關(guān)鍵詞：3D打印;骨骼修復(fù);三維重建;剛性配準(zhǔn);點(diǎn)云

中圖分類號(hào)：TP301.6 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2019）31-0231-03

1概述

隨著3D打印技術(shù)的飛速發(fā)展，醫(yī)學(xué)領(lǐng)域和計(jì)算機(jī)領(lǐng)域的快速結(jié)合，也帶來(lái)了更多優(yōu)秀人才去加入研究骨骼缺損修復(fù)這一領(lǐng)域。多方面領(lǐng)域人才的結(jié)合更快的推動(dòng)了骨骼修復(fù)的發(fā)展，生物醫(yī)療研究人員致力于生產(chǎn)優(yōu)質(zhì)的仿生材料去貼合人體構(gòu)造，而在計(jì)算機(jī)領(lǐng)域則是構(gòu)造更完美的模型去進(jìn)行打印。因此，如何設(shè)計(jì)更適應(yīng)骨骼缺陷的三維模型成了如今研究的重點(diǎn)。

在構(gòu)建契合破損骨骼的三維模型時(shí)，關(guān)鍵在于我們要得到精確的形狀，并使其能夠與正常骨平滑連接。在本文的研究中，我們應(yīng)需要找到合適的待匹配模板，而這種模板是自于模型庫(kù)或病人的對(duì)側(cè)骨。所有的后續(xù)研究都是在這樣的基礎(chǔ)上進(jìn)行的。

2總體流程

該修復(fù)方法首先重建受損骨骼模型嘲由已知的CT序列得到;接著對(duì)模型預(yù)處理，選取合適的模板;再通過(guò)對(duì)應(yīng)的點(diǎn)云模型求取變換矩陣組，最后進(jìn)行模型配準(zhǔn)和計(jì)算缺損部分，具體流程如圖1所示。重建過(guò)程在win10系統(tǒng)環(huán)境下完成，采用

2.1骨骼模型重建

根據(jù)原始CT圖像的橫斷面、冠狀面及矢狀面逐層分割，能更好地提高圖像的完整性和準(zhǔn)確性。CT數(shù)值反應(yīng)的是組織的密度，我們可以參照正常人體組織CT值表，通過(guò)調(diào)窗去除大部分非骨骼組織，實(shí)現(xiàn)粗分割;然后再利用分割算法[3]進(jìn)行精分割，最后對(duì)分割結(jié)果平滑處理。

2.1.1骨骼粗分割

通過(guò)查表可知，骨的平掃CT值在150～1000，因此我們對(duì)原始數(shù)據(jù)進(jìn)行范圍規(guī)定化，在此過(guò)程中，若CT值小于150，則灰度值置為0，;若CT值高于1000，則灰度值置為255;若CT值大于150而小于1000，則灰度值按照（1）得出，在（1）中g(shù)（i，j，k）為灰度值，f（i，j，k）為CT值。

g（i，j，k）=（f（i，j，k）-150/（1000-150）*255 （1）

2.1.2骨骼精分割

由于目標(biāo)具有相同的灰度特征，我們使用Otsu全局閾值分割方法方法是根據(jù)圖像上灰度值的分布將圖像分成背景和前景兩部分看待，前景就是我們要按照閾值分割出來(lái)的部分。背景和前景的分界值就是我們要求出的閾值。遍歷不同的閾值，計(jì)算不同閾值下對(duì)應(yīng)的背景和前景之間的類間方差，當(dāng)類間方差取得極大值時(shí)，此時(shí)對(duì)應(yīng)的閾值就是OTSU算法所求的閾值。

2.1.3平滑去噪

在進(jìn)行骨骼精分割以后，圖像中也存在少許噪點(diǎn)。此外，骨骼邊緣比較粗糙，會(huì)直接影響后續(xù)點(diǎn)云模型的精度。所以我們還需要進(jìn)行平滑處理。

步驟1形態(tài)學(xué)開運(yùn)算去噪：可選擇小半徑的圓形結(jié)構(gòu)元，對(duì)精分割結(jié)果進(jìn)行開運(yùn)算?？稍诒３止趋啦蛔兊那闆r下去除小噪點(diǎn)，并斷開不同骨骼間的偽連接。

步驟2高斯濾波平滑邊緣：選擇小半徑高斯濾波器對(duì)去噪結(jié)果進(jìn)行平滑，然后以閾值T=0.5分割平滑結(jié)果。

2.2求取缺損部分模型

求取缺損模型首先應(yīng)該是在適配模型庫(kù)中找到與缺損模型最相近的完整模型，由于我們需要得到最契合的缺損模型，因此本文提出的方法是在獲取模型點(diǎn)云的基礎(chǔ)上進(jìn)行配準(zhǔn)的。利用MarchingCubes算法獲取體模表面三角片，并計(jì)算三角片重心作為點(diǎn)云集。接著進(jìn)行點(diǎn)云降采樣，這樣做的目的是減少點(diǎn)云數(shù)量，提高配準(zhǔn)速率。

2.2.1點(diǎn)云粗配準(zhǔn)

在得到匹配模板和缺損模型的點(diǎn)云后，利用最小包圍盒算法進(jìn)行點(diǎn)云的粗配準(zhǔn)。最小包圍盒也叫外接最小矩陣，是一種求解離散點(diǎn)集最優(yōu)包圍空間的算法，基本思想是用體積稍大且特性簡(jiǎn)單的幾何體（稱為包圍盒）來(lái)近似地代替復(fù)雜的幾何對(duì)象。

首先，獲取已知點(diǎn)云的三個(gè)主方向，獲取質(zhì)心，計(jì)算協(xié)方差，獲得協(xié)方差矩陣，求取協(xié)方差矩陣的特征值和特長(zhǎng)向量，特征向量即為主方向。接著，利用已獲得的主方向和質(zhì)心，將輸入點(diǎn)云轉(zhuǎn)換至原點(diǎn)，且主方向與坐標(biāo)系方向重回，建立變換到原點(diǎn)的點(diǎn)云的包圍盒。最后，給輸入點(diǎn)云設(shè)置主方向和包圍盒，通過(guò)輸人點(diǎn)云到原點(diǎn)點(diǎn)云變換的逆變換實(shí)現(xiàn)。

通俗來(lái)說(shuō)，最小包圍盒算法是求解點(diǎn)集XOY平面內(nèi)投影點(diǎn)的最小面積包圍盒，即把空間化問(wèn)題轉(zhuǎn)化為平面問(wèn)題來(lái)求解。

2.2.2點(diǎn)云精配準(zhǔn)

由蘭州交通大學(xué)楊軍等人提出的改進(jìn)的ICP算法是本文用于點(diǎn)云精配準(zhǔn)的主要方法，在該方法中主要在初始值的敏感性以及整體與部分模型的配準(zhǔn)精度上做出了改進(jìn)。選擇改進(jìn)的ICP算法適用于精配準(zhǔn)是由于在初始配準(zhǔn)后還未達(dá)到想要的效果，因此選用KD-Tree搜索大幅度地提高了搜索效率，制定自適應(yīng)最有閾值則是在經(jīng)典的ICP算法上提高了配準(zhǔn)精度。

假如我們手上的數(shù)據(jù)有原始點(diǎn)云數(shù)據(jù)P和經(jīng)過(guò)平移加旋轉(zhuǎn)后的點(diǎn)云Q，那么ICP算法就是求出最優(yōu)的從P到Q的坐標(biāo)變換矩陣。從（2）中我們可以看到我們從P到Q的過(guò)程，其中旋轉(zhuǎn)為R.平移為T。

2.3計(jì)算缺損部分

為解決無(wú)法通過(guò)點(diǎn)云直接得到缺損部分的問(wèn)題，我們以點(diǎn)云模型為基準(zhǔn)建立初始采樣網(wǎng)格，后利用變換矩陣對(duì)網(wǎng)格進(jìn)行逆運(yùn)算，再對(duì)原始體數(shù)據(jù)進(jìn)行采樣，即可得到配準(zhǔn)后的體模型。

把適配模型網(wǎng)格稱為模板，缺損模型網(wǎng)格稱為實(shí)例，那么現(xiàn)在用模板減去實(shí)例就可以得到粗結(jié)果。當(dāng)然我們應(yīng)該在粗結(jié)果上進(jìn)行進(jìn)行錯(cuò)誤部分的提出，其中包括刪除扁平、細(xì)長(zhǎng)、顆粒狀部分，最后我們需要利用三維高斯核和球形中值濾波器對(duì)剩余部分進(jìn)行平滑處理。

3實(shí)驗(yàn)

在本實(shí)驗(yàn)中選用羊的脛骨作為實(shí)驗(yàn)對(duì)象，數(shù)據(jù)集一共83張516KB的CT圖像，CT數(shù)據(jù)如圖2所示。讀取整個(gè)序列，并根據(jù)dicom參數(shù)（ResaleSlope、RescaleIntercept、WindowCenter、WindowWidth、ShceLocationl進(jìn)行調(diào)窗與調(diào)序，獲得數(shù)值為0～255的灰度圖像;接著利用閾值分割提取所有骨骼，利用連通分量篩選目標(biāo)骨骼;再根據(jù)dicom參數(shù)（PixelSpacing、SliceTh-ickness），重建序列體數(shù)據(jù)，過(guò)程示例如圖3所示。

在選取到適匹配模板后，那就可以生成缺損模型和匹配到的完整模型的點(diǎn)云，如圖4所示，經(jīng)過(guò)點(diǎn)云的粗配準(zhǔn)即最小包圍盒算法以及精配準(zhǔn)ICP算法后的點(diǎn)云如圖5所示，在經(jīng)過(guò)這一步操作后我們可以得到變換矩陣，這一矩陣的得到是為了下一步做鋪墊。

在得到缺損模型前，需要用采樣網(wǎng)格對(duì)點(diǎn)云進(jìn)行包圍，在圖6中三維坐標(biāo)下豎直采樣網(wǎng)格代表缺損模型，傾斜采樣網(wǎng)格代表匹配到的完整模型，在已知變換矩陣的條件下對(duì)采樣網(wǎng)格進(jìn)行逆運(yùn)算，從而得到配準(zhǔn)后模型。

那么在這樣的情況下，就可以進(jìn)行兩模型的直接相減得到粗略的缺損模型，如圖7所示，在經(jīng)過(guò)剔除錯(cuò)誤部分和平滑后，就可以生成打印模型。MarchingCubes面繪制算法獲取三角片，并以光順優(yōu)化獲得更好的三角片，最后以sTL的格式導(dǎo)出，結(jié)果的缺損模型如圖8所示。