本發(fā)明涉及古生物化石技術領域，更為具體的是涉及一種古生物化石的三維重建方法。

背景技術：

化石是古生物學研究的基本材料，其主要信息都攜帶在化石的形狀之中，化石標本是一類特殊的物體，是生物在死亡以后經過搬運或否，然后埋葬，在經石化過程而成。在埋藏和石化過程中，軟組織腐爛降解，硬組織礦化成巖石，由于生物的結構很復雜，直到成為化石以后有的形態(tài)暴露在外面而一目了然，而有些結構形態(tài)被包藏在骨骼內部而無法直接看到。如果是數量較多的標本，可通過實體切片進行觀察，但當標本數量很少，屬于珍稀標本時，只能進行無創(chuàng)解剖，即通過透視方法獲取化石標本內部的圖像，并根據這些圖像提供的三維信息來研究化石內部的構造形態(tài)。

目前，獲取標本內部形態(tài)信息的途徑主要通過CT掃描成像的方法，而CT掃描成像由一系列二維X光透視掃描圖像組成，可以通過計算機程序進行三維圖像的重建，在進行影像分割提取工作時，需要將每一張二維切片影像進行觀察，鑒別不同的器官、組織或結構成分的界線。

但化石保存狀況受諸多因素影響，由于化石標本在保存中經過了埋藏或者搬運，在軟組織腐爛講解以后，有些空穴會被圍巖或礦物質填充，使得標本內部的密度分布發(fā)生與原結構不同的現象。有時化石標本內部或有變形，或因充填物的密度與骨骼接近而無法區(qū)分骨骼與充填物的界線，使得標本的影像色度分布復雜，難以使用現有的自動選擇工具進行圖形分割，只能使用手工分割，不僅使得提取的三維圖形進度下降，并且分割的工作量很大。

技術實現要素：

為了解決上述問題，本發(fā)明提供了一種計算過程簡單并能快速對化石內部復雜結構進行重建，且能準確區(qū)分骨骼和充填物界線的古生物化石的三維重建方法。

本發(fā)明采用的技術方案是：一種古生物化石的三維重建方法，包括以下步驟：

S1，向古生物化石投射線狀激光進行CT掃描；

從至少兩個角度連續(xù)采集射線狀激光照射的古生物化石的圖像信息，得二維CT掃描切片影像；

S2，將二維CT掃描切片影像轉換為BMP格式；

S3，對BMP格式的二維CT掃描切片影像進行預處理；

所述預處理包括圖像平滑、圖像增強、圖像插值、圖像分割過程；

S4，對經預處理后的二維CT掃描切片影像進行壓縮處理得壓縮CT掃描切片數據；

S5，對壓縮CT掃描切片數據進行圖像處理；

S6，將經圖像處理后的數據進行非線性優(yōu)化與幾何擬合得擬合數據；

S7，對擬合數據進行數據可視化計算，得3D可視化模型；

S8，對3D可視化模型進行動態(tài)模擬計算。

作為上述方案的進一步設置，所述圖像平滑采用中值濾波法。

作為上述方案的進一步設置，所述圖像增強采用拉普拉斯銳化法。

作為上述方案的進一步設置，所述圖像插值采用線性加權平均法。

作為上述方案的進一步設置，所述圖像分割包括圖像二值化、邊緣檢測兩個步驟。

作為上述方案的進一步設置，所述邊緣檢測采用Canny邊緣算子。

作為上述方案的進一步設置，所述圖像處理包括2D和3D圖像處理。

作為上述方案的進一步設置，所述可視化模型進行動態(tài)模擬計算前先進行圖像修復。

本發(fā)明計算過程簡單、運算速度高，可對含有充填物的生物化石影像的邊界進行識別與自動化分割。

附圖說明

圖1為本發(fā)明的結構示意圖。

具體實施方式

下面結合附圖及實施例對本發(fā)明做進一步描述。

如圖1所示，一種古生物化石的三維重建方法，包括以下步驟：

S1，向古生物化石投射線狀激光進行CT掃描；

從至少兩個角度連續(xù)采集射線狀激光照射的古生物化石的圖像信息，得二維CT掃描切片影像；

如果被測物體較大，可采取分塊多測站的形式對古生物化石進行掃描，然后進行量測或對應對點云進行拼接，為保證拼接的準確性，應使相鄰測站掃描的影像具有一定的重疊度。

S2，將二維CT掃描切片影像轉換為BMP格式；

S3，對BMP格式的二維CT掃描切片影像進行預處理；

所述預處理包括圖像平滑、圖像增強、圖像插值、圖像分割過程；

由于激光掃描所得到的數據是全離散的矢量距離點，也即是點云，在掃描的過程之中，掃描速度、設備精度、被測古生物化石的表面情況和操作者的熟練程度等都會對測量數據造成影響，使得得到的數據可能帶有很多離散點和小振幅噪聲，會影響重建后模型的質量。

因此，用計算機對圖像做后處理，對獲取的圖像進行增強信噪比的工作，即濾除圖形的噪聲和干擾，突出感興趣對象區(qū)域或邊緣，便于邊緣提取和圖像分割。

1、圖像平滑

對于噪聲干擾的圖像，可用圖像平滑方法來濾除噪聲。在頻率域中，低頻相應于空間的大物體，高頻相應于細節(jié)和邊緣，與噪聲相混淆。一個較好的平滑方法或濾波方法應該是既能消掉這些噪聲又不使圖像的邊緣輪廓和線條變模糊。

中值濾波在一定條件下可以克服線性濾波器如最小均方濾波、均值濾波等帶來的圖像細節(jié)模糊，而且對抑止圖像中的脈沖干擾和椒鹽噪聲特別有效。因為這些噪聲在圖像中往往以孤立點的形式出現，與之對應的像素又比較少，所以采用中值濾波能有效地去除這些噪聲，達到圖像平滑的目的。

中值濾波法一般采用一個含有奇數個點的滑動窗口，對該滑動窗口內的諸像素灰度排列，用其中值代替窗口中心像素原來的灰度值，或者說用局部中值代替局部平均值。具體來說，就是假設有一個N維奇數的離散序列a₁，a₂，…，a_n，其中值為m，那么用它來代替窗口中心所對應像素的灰度。例如，取一維序列的元素為N＝5個，數值分別為80,90,200,110,1200計算的中值m＝110。將m＝110放在這5個數中間，對序列重新排序后得80,90,110,120,200。顯然，如果200是一個噪聲點，通過中值濾波后，圖像質量就會有很大的改進。

一般情況下，序列為N的一維窗口的中值為：

median(x₁，x₂，...，x_N)＝m (I)

如果N為偶數，則取兩個中間數的平均值作為中值。

2、圖像增強

圖像增強的目的是為了改善圖像的視覺效果，或者是為了更便于人或計算機的分析和處理。為了使一幅圖像的邊緣更加鮮明，有時也需要對圖像進行尖銳化增強處理，因而圖像尖銳化技術也常常用做圖像增強的工具。在圖像增強的過程中突出了一部分信息，同時可能會壓制另一部分信息。

本發(fā)明的圖像增強采用拉普拉斯銳化法，原理如下：

對于給定的二維離散圖像∫(x，y)，可以計算出其一階差分為：

二階差分為：

從而根據拉普拉斯算子：

計算出：

對于因擴散引起的圖像模糊，可以用式(3.12)進行銳化：

其中，k是與擴散有關的函數。如果令k＝1，則有：

g(x，y)＝5∫(x，y)-∫(x-1，y)-∫(x+1，y)-∫(x，y+1)-∫(x，y-1) (7)

3、圖像插值

當斷層圖像間的距離比斷層圖像內像素間的距離大得多時，就需要用圖像插值方法在原來的斷層圖像之間再插值生成一些中間斷層圖像。因為許多體數據三維重建和三維顯示方法要求體數據是各向同性的。因而斷層圖像就要插值成各向同性，即經插值后的斷層圖像序列中斷層間距等于斷層圖像內像素間距。然而，圖像插值是一個具有很大任意性的問題。為了使圖像插值成為一個確定的、可解的問題，通常引入下面三個約束條件：

(1)插值圖像要與原始斷層圖像相似；

(2)插值圖像與兩個原始斷層圖像的相似度應該分別和它與這兩個斷層圖像的距離成反比關系；

(3)插值圖像序列應該呈現出從一幅原始斷層圖像到另一幅原始斷層圖像的漸變過程。

最簡單的插值方法是對上下兩個相鄰的斷層圖像進行加權平均，產生一組插值圖像。當斷層間距與斷層圖像內像素間距相差不是很大時，采用這種插值方法是可行的。

但是，當斷層間距與像素間距相差很大時，這種方式生成的插值圖像模糊不清。其原因是，兩個相鄰斷層圖像中處于同一位置的像素不一定對應同一物質，他們的加權平均沒有什么意義。對斷層間距較大的插值，一種較好的方法是基于匹配的圖像插值，這種方法能夠更好的構造中間斷層圖像。關于斷層圖像的匹配是模式識別和計算機視覺領域所研究的主要問題之一。圖像匹配的任務就是要在兩幅圖像之間建立對應關系，也就是尋找圖像間的變換。傳統(tǒng)的匹配方法一般包括特征選取、特征匹配和整幅圖像之間匹配映射關系的確定三大步驟，匹配的可靠性和準確性很大程度取決于是否能夠找到魯棒的方法去完成特征提取和特征匹配。

本發(fā)明采用線性加權平均法來進行圖像差值，線性加權平均的圖像插值方法描述如下：

設S_k(i，j)、S_k+1(i，j)分別是第k層和第k+1層切片圖像。他們之間的插值圖像可表示為：

S_λ(i，j)＝(1-λ)·S_k(i，j)+λ·S_k+1(i，j) (8)

其中λ＝d₁/(d₁+d₂)，d₁、d₂分別是插值圖像到第k、k+1層圖像的距離。顯然，當λ＝0時，S_λ(i，j)＝S_k(i，j)，當λ＝1時，S_λ(i，j)＝S_k+1(i，j)。給出一組{λ_i|λ_i∈(0，1)，i＝1，2，...，n}，就相應的得n個插值圖像。為得到等間距的插值圖像，參數序列{λ_i|i＝1，2，...，n}應該取

4、圖像分割

對CT圖像進行基于邊緣檢測的圖像分割處理，包括兩個步驟：圖像二值化和邊緣檢測。

(1)圖像二值化

在進行圖像處理時，由于普通計算機的內存和計算能力都是有限的，從而為了有效降低需要處理的數據量，同時清晰分辨圖像邊界輪廓，可采用圖像二值化的方法對圖像樣本進行預處理。二值圖像有以下幾方面特點：

①計算二值圖像特性的算法非常簡單，容易理解和實現，并且計算速度很快。

②二值圖像所需的內存少，對計算設備要求低。

③許多二值圖像處理技術也可以用于灰度圖像的處理。

二值化指將灰度圖像轉換為二值圖像，二值圖像指圖像中的每個象素只取兩個離散的值之一，即非此即他。用式(9)表示：

上式中，∫(x，y)表示一幅數字圖像，x、y表示該圖像中某象素的坐標值，而0和1表示該象素的象素取值。

由灰度級直方圖(Grey Level Histogram)確定整體閾值

設規(guī)范化灰度值g范圍為0≤g≤1，g＝0為最黑，g＝1為最白。M為灰度級數目，p(g_k)為第k級灰度的概率。n_k是在圖像中出現的灰度級為k的次數，n為圖像中像素的總數。則有：

p(g_k)＝n_k/n 0≤g≤1，k＝1，2，...，M (10)

通常稱以p(g_k)為縱坐標，以g_k為橫坐標的圖形為灰度直方圖。閾值應該取在兩個峰值的波谷處，波谷越深陡二值化效果越好。

根據灰度直方圖選取閾值T將圖像分成C₁和C₂兩個類，其中C₁類像素灰度值比T值高，C₂類像素灰度值比T值低，計算類間方差σ²_B和類內方差σ²_W：

類間方差：

類內方差：

其中，w₁和w₂分別是類C₁和類C₂的像素數，u₁和u₂分別是類C₁和類C₂中像素的灰度平均值，u是所有像素的灰度平均值，σ²₁和σ²₂分別是類C₁和類C₂的灰度方差。

使分離η(T)為最大值的閾值T即為最佳閾值，即：

(2)邊緣提取

將CT圖像二值化后，根據圖像本身的情況可以對其進行邊緣檢測。其目的是把圖像中人們感興趣的部分分離出來，突出想要的目標，減少信息量。

圖像的邊緣是圖像的最基本特征，所謂邊緣是指其周圍像素灰度有階躍變化或屋頂變化的那些像素的集合，本發(fā)明采用Canny邊緣算子來計算：

Canny檢測主要檢測階躍性邊緣，Canny算子的原理如下，

式(14)定義了G(X)：

分別對式(14)求x偏導，求y偏導得：

利用兩個微分模板估計梯度矢量

其幅度邊緣為強度梯度方向角為

對像素邊緣強度和梯度方向角應用非最大值抑制技術，具體如下：

①計算梯度方向上與3×3邊界相交的兩個虛擬像素Q₁和Q₂。(i，j)是3×3模塊的中心像素，Q₁和Q₂是(i，j)梯度方向上與3×3邊界相交的虛擬像素。

②用所Q₁在邊框上的兩個與Q₁相鄰像素的梯度強度線性內插Q₁點的邊緣強度用Q₂所在邊框上的兩個與Q₂相鄰像素的梯度強度線性內插Q₂點的邊緣強度

③如果中心像素(i，j)的邊緣強度同時大于和則保留此中心像素為候選邊緣，否則極為非邊緣像素；

應用Hysteris閾值判斷是否是邊緣像素，具體如下：

①對每個候選邊緣，如果則標記為邊緣像素；

②對剩下的每個候選邊緣，如果在其3×3鄰域，至少已有一個鄰像素是邊緣，那么標記為邊緣；

③迭代進行第二步，知道沒有新的邊緣像素被標記出位置，則剩下的像素被標記為非邊緣像素；經過以上步驟，圖像就變成二值化的邊緣圖像。

S4，對經預處理后的二維CT掃描切片影像進行壓縮處理得壓縮CT掃描切片數據；

由于測量精度的提高，形成的數據可能會很大，而且其中會包括大量的冗余數據，必須對數據進行壓縮處理，因此，采取先確定整個數據在整個三維空間的最外區(qū)域，然后，將整個區(qū)域細分，在每一個小區(qū)域內求小區(qū)域內的點云，然后根據點云曲率的特點壓縮點云，主要步驟如下：

①確定原始點云數據最外區(qū)域；

②區(qū)域細化；

③細化區(qū)域內點云曲率的計算；

④計算小區(qū)域內每兩個的點云曲率的差值；

⑤確定曲率差值的閾值，保留小閾值的點，刪除大于閾值的亮點中曲率較小的點；

⑥壓縮結果不滿足要求時，循環(huán)步驟4和5，第二次壓縮至數據壓縮比在70％～80％。

S5，對壓縮CT掃描切片數據進行圖像處理；

由于投影數據形成的圖像為2D圖像，不能直觀的顯示古生物化石的效果，因此對濾波反投影數據進行2D和3D圖像處理，使之可以形成3D圖像

S6，將經圖像處理后的數據進行非線性優(yōu)化與幾何擬合得擬合數據；

S7，對擬合數據進行數據可視化計算，得3D可視化模型；

由于發(fā)掘的古生物化石在形成過程中可能出現破碎、扭曲和不完整的情況，因此需要對可視化模型進行圖像修復，圖像修復采用Criminisi算法，其算法如下：

(1)讀入可視化模型的圖像與圖像掩碼；

(2)計算待修復區(qū)域的邊界，初始化置信度和數據項；

(3)計算邊界上像素點的置信度，數據項，以及鄰域相關信息，根據優(yōu)先權計算公式計算各點的優(yōu)先權；

(4)以優(yōu)先權最高的像素點為中心，形成初始待修復模板，再根據模板自適應大小方法，修改模板大小，若模板超過一定閾值大小，記錄其邊界；

(5)按照新的模板匹配方式，計算每個已知模板與待修復塊的直方圖相交距離和兩者間的SSD值，以搜尋最佳匹配模板；

(6)將最佳匹配模板的信息拷貝到待修復模塊中的未知部分，更新置信度和優(yōu)先權等信息；

(7)若邊界不為空，轉到第三步繼續(xù)執(zhí)行；否則，轉到下一步；

(8)根據第三步記錄的較大模塊的邊界，對其采用Telea算法再進行一次快速修復，以消除邊界效應，然后算法結束。

S8，對3D可視化模型進行動態(tài)模擬計算，使形成的古生物化石可視化模型可以以動態(tài)形式呈現。

上面結合附圖對本發(fā)明的具體實施方式作了詳細說明，但是本發(fā)明并不限于上述實施方式，在本領域普通技術人員所具備的知識范圍內，還可以在不脫離本發(fā)明宗旨的前提下作出各種變化。