基于隨機森林算法的小鼠micro- CT影像中骨骼關節特征點定位

屠睿博 陳中華 王洪凱

(大連理工大學生物醫學工程系，遼寧 大連 116024)

基于隨機森林算法的小鼠micro- CT影像中骨骼關節特征點定位

屠睿博 陳中華 王洪凱#*

(大連理工大學生物醫學工程系，遼寧 大連 116024)

隨著小動物成像技術的發展，技術人員每天需要處理的小動物影像數量急劇增長，這使得自動化的小動物圖像分析方法成為迫切的需求。在小鼠圖像分析方面，小鼠靈活多變的身體姿態給自動化的圖像分析帶來困難。基于隨機森林算法實現小鼠micro- CT圖像中骨骼關節點的自動定位，為解決小鼠影像中身體姿態的自動識別打下基礎。該算法主要分3步：先通過分類隨機森林算法得到小鼠骨骼關節點的粗定位，再通過回歸隨機森林算法進一步減小定位誤差，最后通過圖匹配的方法在備選點中挑選正確位置上的關節點。對49例不同身體姿態的小鼠全身三維micro- CT圖像進行測試，全身關節點定位的成功率為98.27%，定位誤差的中值為0.68 mm。同時驗證聯合使用分類與回歸隨機森林的必要性，并探究訓練數據的數量對不同骨關節的識別效果的影響。研究為小鼠micro- CT影像中身體姿態的識別提供一種新方法，為后續的自動化圖像配準、圖像分割以及自動化圖像測量提供重要的定位信息。

小動物影像分析； 骨關節點定位； 隨機森林； 模式識別； 顯微CT

引言

在癌癥研究和新藥試制的臨床階段，小動物成像起著重要的作用。近些年來，隨著小動物成像研究的增多和小動物成像設備的發展[1- 4]，小動物影像的數量劇烈增長。與此同時，研究人員處理分析小動物影像的壓力也隨之增大，不得不每天從事繁冗、單一的基本操作。因此，對小動物影像自動分析系統的需求愈加迫切。在小鼠全身影像的自動分析算法中，數字解剖圖譜的全身配準是最基礎的一步。通過配準，圖譜中包含的豐富的器官解剖結構信息可被映射到個體圖像中，從而為后續的精確器官分割以及各種影像參數的測量打下基礎。

目前，阻礙小鼠圖譜與個體影像精確配準的最大障礙之一，就是圖譜與個體之間的姿態差異。由于小鼠的骨骼關節靈活，不同個體之間姿態差異可以很大。為解決這一問題，研究人員嘗試了多種方案，主要可分為3類：非線性變形配準、活動關節配準、統計形狀配準。Kesner等通過人工定義PET影像中的特征點實現圖譜與PET影像的配準是一種非線性變形配準的方法[5],但是非線性變形在姿態差異較大情況下容易導致器官和骨骼的扭曲變形。Baiker等通過定義小鼠活動骨骼圖譜的方式，在一定程度上克服了小鼠姿態變化對配準的影響[6- 8]，但是仍然無法自動配準較大的姿態差異。筆者前期研究通過統計訓練數據形狀進行建模，實現對身體姿勢和器官形態同時有差異圖像的配準[9- 10]，但是這種統計形狀模型只適合配準訓練樣本集所包含的身體姿態。

為了解決大幅度自由姿態變形的小鼠圖譜配準問題，骨關節點的定位至關重要。近年來，隨機森林算法在醫學影像和人體表面掃描影像中實現快速準確的關節定位。Shotton使用分類隨機森林算法和Mean- Shift聚類方法在人類體表掃描深度圖中預測關節點位置，由隨機森林算法得到全身關節點的類別判斷[11]，Mean- shift算法實現最后關節點的篩選；由于體表掃描影像的獲取較為便利，而且他們還通過人體數字模型模擬方式產生巨大的訓練數據，所以可以提供充足的訓練數據。這樣龐大的數據量，對于臨床CT圖像的收集是不可能的。不同于體表掃描的深度圖，CT圖像數據的干擾因素也很多，比如設備的差異帶來像素分布范圍和圖像分辨率的不同。因此在CT圖像數據上使用隨機森林算法，后期對分類結果的進一步處理是必要的。Donner等針對CT圖像，在隨機森林算法得到的定位基礎上，進一步根據相鄰關節點之間的相對位置關系來最終確定每個關節點的最優位置[12]，在人體CT影像的骨骼關節定位方面取得了良好的效果。

在小動物成像領域，微型斷層成像(micro- computed tomography, micro- CT)技術具有高分辨率、方便監測、易于建立多模態影像、提供精確的解剖信息等優點[13], 因此獲得了廣泛的應用。在小鼠的micro- CT影像中，骨組織的對比度好，包含較為清晰的全身關節圖像特征。因此本文著重研究小鼠micro- CT影像中骨關節的自動識別問題。雖然在人體三維CT影像領域，骨關節的定位已經取得較好的效果，但是在小鼠的骨關節定位領域還未見相關的研究。本研究針對小鼠micro- CT影像，使用隨機森林的分類算法首先得到對醫學圖像關節點的初步預判，隨后使用回歸隨機森林對已有的分類結果進行優化，得到高精準度的關節點備選點，最后通過基于訓練數據的結構圖進行關節備選點的全局優化，在已有的關節識別結果中挑選最優的備選點。以下章節將先簡要介紹隨機森林算法的原理，然后重點介紹本研究方法，并給出實驗結果和總結。

1 方法

1.1 隨機森林算法簡介

隨機森林是由多顆決策樹并聯而成的集成分類系統，每一棵決策樹都是一個能力較弱的模式識別分類器(簡稱“弱分類器”)，而多棵樹的集成則構成了一個強分類器[14]。

在20世紀80年代，早期決策樹的研究工作主要在“分類回歸樹”(CART)[15]中體現，研究者則主要集中于使用已有數據構造最優的決策樹。在20世紀90年代，學者們發現使用集成學習器(比如一般的弱習器)可以得到更高的準確率[16]。集成方法結合決策樹便產生了決策森林，首先應用在手寫數字的識別上[17]。20世紀90年代末，Breiman提出“隨機森林”的概念[18]，并且實現更廣泛的使用。

為了解隨機森林的工作原理，首先需要了解森林中每一棵決策樹的原理。如圖1所示，待識別數據的n維特征向量v從樹狀結構最頂層的根節點輸入，并逐層向下走，每到達一個節點就根據該節點所保存的判別式進行一次判別，決定繼續走向左下還是右下節點，直到到達最底層某個葉子節點(見圖中方塊節點)，由最終所達到的葉子節點給出最終輸出結果。圖中虛線表示了特征向量從根節點到葉子節點所走過的一個可能的路徑。

在樹中每一個節點，判別過程如下：隨機抽取v中的d個特征構成一個d維(d

圖1 決策樹原理Fig.1 The decision tree

根據應用目的的不同，隨機森林可以分為分類森林(classification forest)和回歸森林(regression forest)[19]。分類森林的作用是解決傳統的模式分類問題，當輸入特征到達決策樹最底層的葉子節點后，以葉子節點中所保存的訓練樣本中最大概率類別作為分類的結果。回歸森林則是為了解決函數擬合的問題，根據輸入的特征向量v來預測一個連續的輸出值y=f(v)，整個森林就相當于是對函數f(v)的擬合；在回歸森里中，當輸入特征到達決策樹最底層的葉子節點后，會以葉子節點中所保存的訓練樣本的條件概率函數來預測該葉子的y值。舉例說明，在本研究的應用中，分類森林用來識別圖像中的像素屬于哪個關節點，而回歸森林用來預測像素到關節點中心的偏移量。

分類森林與回歸森林在決策樹每個節點的判別原理上是基本相似的，只是葉子節點的輸出原理不同。關于分類森林和回歸森林的詳細介紹推薦參考文獻[20]。對于這兩類隨機森林的訓練和測試方法，本研究均采用文獻[20]所發布的sherwood C++函數庫，因此在后面章節中將不做累述，本研究重點介紹針對小鼠micro- CT圖像的算法設計。

1.2 算法總體流程

所處理的小鼠micro- CT影像都是將二維斷層影像在三維空間中堆疊重構成的三維圖像。處理的過程是逐個掃描三維圖像的體素，對每個體素提取100維的特征向量作為分類森林和回歸森林的輸入。分類森林被用來識別每個像素位于哪個關節點附近，從而在每個關節點附近得到一群備選點組成的點云，并得到每個點屬于每個關節的概率；進一步，回歸森林會計算備選點云中每個點到關節中心的偏移量，更加精確地定位關節的中心位置；針對分類森林和回歸森林可能產生錯誤識別的問題，用圖匹配的方法、利用相鄰關節之間的相對位置關系做約束，篩選掉可能產生的錯誤識別點，得到最終的定位結果。

本方法中的分類森林、回歸森林和匹配圖模型均通過有監督的訓練得到。在訓練階段，每幅訓練圖像中的關節點位置由人類專家手工標定。如圖2所示，算法分為測試階段和訓練階段。訓練數據是標記關節的小鼠三維CT圖像，用于訓練分類、回歸隨機森林和匹配圖。測試階段以小鼠的膝關節識別為例說明。先通過分類隨機森林得到初步測試結果，再通過回歸隨機森林縮小識別結果范圍，最后通過匹配圖篩選正確的關節識別結果。

1.3 分類森林

分類隨機森林算法是關節點定位的第一步，對小鼠全身micro- CT圖像進行處理，預測每一個體素屬于不同類別的概率，得到初步的類別概率預測結果。

1.3.1 分類森林的數據預處理

分類森林的訓練數據是n幅不同身體姿態的小鼠全身CT圖像，其中骨關節點的位置坐標已經被人為標記。為了滿足訓練正樣本的數量，需要擴充每幅圖像中關節標記點的個數。由于關節點的標記是人工標記的，并且關節是一定范圍內存在的，所以擴充關節點標記的個數是有必要的。另外，為了準確識別關節并且增加算法的魯棒性，也需要擴充每個關節的標記點。擴充正樣本的方式是用已標記關節點為中心、邊長為5個體素的正方體范圍內的體素作為關節點標記點，這樣每幅圖像有20個關節，每個關節有125個標記點，由原來的每幅圖像20個正樣本擴充到2 500個正樣本。

圖2 算法的流程Fig.2 Algorithm Workflow

除了訓練需要的正樣本，算法需要提供負樣本，即背景類。背景樣本的采樣方式有很多，這里只介紹其中的兩種。第一種方式是背景采樣，在去除標定的關節點后，對圖像剩余的空間進行均勻采樣，或者隨機采樣一定數量的負樣本。第二種方式是在標記的關節點周圍采樣負樣本點。注意，對于負樣本采樣的數量，應該遵循負樣本類別的數量與正樣本數量近似相等的原則。因為當負樣本數量遠大于正樣本數量時，隨機森林的訓練將是一種不平衡訓練，影響對各類別的預測。同時，若采用第二種負樣本采樣方式，負樣本的采樣空間應該距離正樣本的標記空間一定的距離，以免出現對有歧義的空間進行標記，導致隨機森林學習失敗、分類能力減弱。

為了有效控制負樣本的存在范圍，采用簡單的閾值分割得到骨骼，并對分割結果用半徑為5 mm的球體算子做形態學膨脹，以膨脹后的結果為骨骼遮罩，從中提取負樣本進行訓練。在后期算法測試的過程中，也對測試圖像生成同樣的遮罩并從中尋找骨骼關節點，實驗證明這種方法比全局采樣得到負樣本的方法識別結果敏感性要好。同時，骨骼遮罩的方式還有利于縮小關節點搜索范圍，提升算法速度，并防止小鼠其他部位影響識別效果。

1.3.2 分類森林的訓練

特征設計是隨機森林的基礎，特征的好壞直接決定后續識別能力的強弱。體素灰度值作為特征，這種簡單的特征不僅可以提升算法的速度，而且由隨機森林在計算機視覺中應用的經驗得知，以體素或者像素值作為特征符合隨機森林的學習習慣，可以得到很好的學習效率和效果。本研究采用隨機森林圖像分類算法中常用的雙體素差值特征，如圖3所示，對于一個體素的位置u，定義其特征為d維特征向量，有

(1)

(2)

式中，δi1和δi2為隨機產生的三維偏移向量，函數代表像素灰度值。

圖3 雙體素差值特征Fig.3 Double voxel difference feature

這種定義實際上是將體素u的第i維特征定義為其周圍隨機抽取的兩個體素的灰度差值。這種差值可以摒除因為圖像整體明暗偏差造成的影響。隨機抽取的過程遵循高斯分布，其均值為0，方差為7個體素(實驗嘗試確定)，共抽取100次，從而形成100維的特征向量v(u)。這里引入隨機性，因為隨機森林要求每一棵樹要盡可能不同，這樣可以避免學習得到的森林只能識別某一特定模式的圖像，從而具備較強的魯棒性，可以識別多種姿態小鼠CT圖像的骨關節點。

最后，關于分類森林訓練的最優化方法，采用訓練數據劃分后的信息增量作為判決條件，隨機地計算若干次劃分結果，并選取可以最大信息增益的參數作為本次訓練的參數[20]。

1.3.3 分類森林的測試

由于分類森林的訓練是基于骨骼遮罩，并提取了100維的特征向量，因此在測試過程中要使用相同的閾值提取骨骼遮罩，并遵照訓練過程為每個節點確定的偏移量來提取100維像素特征。

1.4 回歸森林

1.4.1 數據預處理

回歸森林使用的原始訓練數據同分類隨機森林相同。由于回歸隨機森林的輸出是連續的向量，所以在訓練階段需要給算法提供每個訓練樣本點指向關節中心點的偏移向量。所有訓練數據是采用以關節點為中心，21為邊長的正方體范圍內的點，并以每個點指向關節點的三維向量為訓練的輸出。因此每幅圖像20個關節，每個關節有9 261個訓練數據，這要比之前的分類隨機森林的訓練數據要多，因此回歸隨機森林的訓練耗時更長。與分類隨機森林不同，回歸隨機森林不需要提供負樣本的數量，因為不需要提供指向背景類的向量。

1.4.2 回歸隨機森林的訓練

回歸隨機森林訓練目的是根據圖像的局部特征學習出某一位置指向關節點的向量。所采用的輸入特征與分類隨機森林相同，定義同式(1)、(2)，同樣是通過隨機抽取形成100維的特征向量v(u)。在訓練方法上，回歸隨機森林也是采用信息增量作為判決條件，它與分類隨機森林的不同點主要體現在兩個方面：其一，由于分類森林的輸出是不連續的類別變量，而回歸森林的輸出是連續的實數變量，因此在訓練每個節點的時候，分類森林采用離散的信息熵，而回歸森林采用連續的信息熵；其二，在葉子節點的訓練方式上，分類森林選取葉子節點中出現次數最多的訓練樣本類別作為輸出類別，而回歸森林則首先獲取葉子節點中每個訓練樣本指向關節中心的位移向量，然后對這些向量擬合一個三維的條件高斯分布，以像素特征v(u)為條件，計算該條件下的高斯分布均值作為輸出的位移向量。受篇幅所限，推薦讀者參考文獻[18]進一步了解算法細節。

本研究針對20個不同的關節點訓練20個不同的森林，預測某一點到某一個關節的距離只需輸入到相應的森林中進行預判。

1.4.3 回歸隨機森林的測試

在回歸隨機森林的測試階段，需要分別預測每一個體素測試距離20種不同類別關節點的偏移向量。測試時，需要先統計分類森林結果中不同類別識別概率中的最大概率，只針對那些識別概率大于最大概率一半的點進行測試，從而摒除不必要的干擾點，減小計算量。

得到初步的偏移向量后，移動測試點以遞歸的方式繼續進行回歸測試。當遞歸的次數過少時由于偏移向量指示范圍有限而導致誤差較大，而當遞歸次數過多時，將會導致測試點的過度移動反而加大定位的誤差。因此，本研究選用遞歸3次進行測試。得到的測試結果是每一個類別的備選點經過3次移動后的新位置。

1.5 使用圖匹配方法篩選關節備選點

回歸隨機森林得到的結果是某一個關節點的多個備選點，需要進一步從這些備選點中挑選出一個最優點作為最終的定位結果，本研究使用基于隨機場的圖匹配方法來實現選取最優點。

1.5.1 匹配圖的訓練

若想在多個位置中挑選出一個關節點，除了考慮該點的識別概率，還應該考慮與它相鄰的關節點的相對位置關系，因此需要根據訓練數據學習出對某一關節點識別最有影響力的幾個相鄰關節點。

假設一共有N個關節點，每一個關節點有N-1個關節點與其相連，首先需要根據訓練數據找到對每個關節點影響較大的幾個關節點。倘若兩個關節點在多幅訓練數據中分布的相對位置關系是穩定的，則認為二者之間有較強的相互影響，其相對位置關系對于關節備選點的篩選十分有利。

關節點之間相對位置關系的穩定性可以通過熵的概念來描述，相對位置關系穩定的兩個關節點之間的熵最小。在N幅訓練數據中，第i幅圖像的兩個關節點s、t(s是待訓練關節點，t是與s相鄰關節點)，二者的三維坐標分別是Isi和Iti，t相對于s的偏移向量是ΔIs,ti=Iti-Isi。由于訓練數據的數目是有限的，為了以有限訓練樣本模擬無窮多樣本中關節點之間的相對位置分布，假設每個訓練樣本都可擴充為無窮多個相似的樣本，以三維向量ΔIs,ti為均值構造一個三維高斯分布，gs,ti=N(ΔIs,ti,σg)來代表s與t之間的相對位置分布，其中分布半徑σg是一個經驗參數(其取值對最終結果影響很小，本研究中設為2個像素)。對于N幅訓練圖像，會得到N個gs,ti高斯分布，對這N個分布的求均值并計算均值熵，即可作為對關節s與t之間相對位置穩定性的衡量，有

(3)

計算s同其他所有關節點之間的熵并排序，選出最小的前T個作為對s影響最大的關節點(本研究選取T=2)，用于在后面的圖匹配過程中幫助挑選關節s的最優備選點。

1.5.2 圖匹配方法的測試

圖匹配的作用是為每個關節點在多個備選點中挑選出最優的點作為最終結果，其挑選過程不僅利用了相鄰關節點之間的相對位置關系作為先驗知識，還利用了之前分類隨機森林的每個備選點屬于該關節的概率。圖匹配算法的代價函數定義為

(4)

這個典型的隨機場問題可以通過經典的圖割方法[21]來求解，如果優化問題的規模足夠小，也可以通過遍歷所有組合的方式來求解。這里解釋一下一元項和二元項的定義。一元項的計算是關節識別概率的信息量的計算，有

(5)

(6)

(7)

(8)

在使用圖割方法求解該隨機場時，需要采用多標簽的圖割方法[22]，這種方法采用迭代機制，通過多次雙標簽的圖割運算來近似求得最優的多標簽分割結果。在本研究中，骨骼關節l對應隨機場中的節點，關節對e代表隨機場中的邊，式(4)中的一元項對應圖割中的t- link，二元項對應圖割中的n- link。由于每個關節具有多個備選點，因此所有關節的備選點可以組成一個備選點集∏；假設∏中共有K個點，則編號1～K對應了圖割問題中的K個標簽，從而將備選點選取問題轉化成了多標簽的圖割問題。需要注意的是，如果某個備選點k是屬于關節i的備選點，則須將k與其他關節之間的t- link值設為無窮大，以避免圖割算法將k賦值給其他關節。

1.6 算法評價方法

本研究所采用的訓練數據是49幅不同身體姿態的小鼠的全身micro- CT圖像，采集自美國加州大學洛杉磯分校Crump分子影像中心。圖像由256×256×496個體素構成，每個體素的尺寸為0.204 mm×0.203 mm×0.202 mm。筆者訓練的分類森林包含50棵樹，每棵樹深度為15層；回歸森林包含20棵樹，深度為10層；以上參數通過實驗嘗試得到。被識別的骨骼關節點包括20種，分別是鼻尖、脊柱與頭骨連接處、胸椎上端，胸椎下端、胸骨柄上端、胸骨柄下端、骨盆后延、骨盆前端結合處、左肩、左肘、左手、右肩、右肘、右手、左胯、左膝、左踝、右胯、右膝、右踝。

測試采用leave- one- out的交叉驗證方式進行，即每測試一幅圖像，就以另外48幅作為訓練數據。對于被測試圖像，人類專家的手工標記相當于衡量算法精度的金標準；對于訓練圖像，手工標記則是訓練數據的一部分。

2 結果

圖4展示了多幅不同姿態的下小鼠骨骼關節定位結果。結果圖以骨骼閾值分割結果的半透明三維面繪制的方式展現，突出小鼠的骨骼結構，算法給出的關節點以圓點的形式標出。部分圖像中的CT床由于像素值較高，同骨骼一起被繪制出，但是我們的算法并沒受到CT床的干擾。可以看出算法對不同姿態的小鼠均實現了較好的關節定位，圖中第一行中間兩個小鼠圖像中的CT床沒有被去除，第二行左側兩個圖像沒有小鼠的足部，但是這些因素都不影響小鼠關節的識別。

圖4 多種姿態小鼠的關節識別結果。小鼠的骨骼以閾值分割的半透明三維面繪制方式展現，識別出的骨骼關節圓點形式標記圖像上Fig.4 Bone joints localization results of various body postures. The skeletons are displayed as surface rendering of the threshold bones, the detected joints are marked as dots.

2.1 算法準確性分析

這里采用均值誤差和中位數誤差評估定位關節點的準確率。

圖5展示算法對每個關節識別的均值誤差和中位數誤差，由于踝關節的訓練數據和測試數據不足沒有將踝關節加入考慮。定位誤差定義為算法定位的位置與人類專家標記位置之間的歐氏距離。誤差計算是基于49幅圖像的測試結果統計得到，可以看到大部分關節的中值誤差在0.8 mm以下，均值誤差在1.0 mm以下，所有關節點的中值誤差為0.63 mm，均值誤差為0.97 mm。(中值誤差為3.16個體素，均值誤差為4.80個體素)，說明本方法對于小鼠micro- CT圖像骨關節有較好識別和定位能力。

圖5 關節定位誤差Fig.5 Bone joints localization errors

圖5中誤差較大的兩個關節是左膝和右膝，主要原因并不是算法對這兩個關節本身的識別能力弱，而是由于部分測試圖像的踝關節缺失，會造成算法把一部分膝關節備選點識別為踝關節，使膝關節的識別精度受到影響。

2.2 算法識別成功率分析

這里首先定義識別失敗，當某一關節點識別到其他關節的位置時記作識別失敗，當識別的關節點距離正確位置與其他錯誤位置的距離相等或者大于到錯誤關節點的距離也記作識別失敗。當識別的關節點在正確位置，以及相對于錯誤位置的距離遠大于正確關節點的距離，定義為識別成功。

算法對絕大多數關節均能正確識別，唯一的錯誤識別來自少數缺少踝關節圖像：當小鼠的腳沒有被掃描到圖像中，膝蓋的識別會受到影響，因為膝蓋和踝關節相關性比較高，在通過匹配圖篩選關節點時，有踝關節可以幫助篩選膝關節的備選點。

如圖6所示，橫坐標為定義的關節各類別，識別成功率為成功識別的關節個數與同類關節全部參與識別的關節個數的比值。可以發現，由于部分數據的踝關節缺失的干擾，使膝關節的識別率也受到影響，但是單個關節的最低識別率仍達92%。綜合考慮所有關節，一共49×20=980個關節，有17個錯誤識別的關節(包含識別錯誤的踝關節)，占所有關節的1.73%，即全身關節點定位的成功率為98.27%。

圖6 關節識別成功率Fig.6 Bone joints recognition rate

2.3 回歸隨機森林的必要性探究

本研究方法中，回歸森林用于進一步改善關節識別的精度。但是，回歸森林的使用會延長算法執行時間，因此有必要驗證回歸森林是否真的十分必要。這里將單獨使用分類森林結合匹配圖進行關節的定位，并同分類森林+回歸森林+匹配圖進行關節定位進行比較，探究回歸森林發揮的作用。

圖7 關于是否使用回歸森林的誤差對比。(a)“分類森林+回歸森林+匹配圖”的誤差；(b)“分類森林+匹配圖”的誤差Fig.7 Identification error comparison between with and without the regression forest. (a) The errors of “classification forest +regression forest+graph matching”; (b) The errors of “classification forest + graph matching”.

圖7(a)、(b)所示分別是“分類森林+回歸森林+匹配圖”和“分類森林+匹配圖”的結果箱線圖，縱坐標表示定位誤差，橫坐標為骨關節類別。可以看到，圖7(a)中各關節識別誤差的四分位間距超過5個體素的只有一個關節點(上胸骨柄)，其他識別誤差結果穩定在上下四分位相差為2～3個體素，證明本研究所用方法識別小鼠骨關節穩定性高。而圖7(b)中多數關節的中值誤差是圖7(a)的2倍左右，且上下四分位的間距也更大，說明如果沒有回歸森林，大部分關節的定位精度和穩定性都會下降。

2.4 訓練數據量對算法精度的影響

這里研究訓練數據數量和訓練數據小鼠姿態對小鼠骨關節定位的影響。如圖8所示，49幅測試圖像中的身體姿態主要呈現為3類：最左側的姿態雙腿張開，中間的腿部收縮，右側的軀干扭曲，這三種姿勢定義為1、2、3類，其他各類非主流姿勢定義為4類。

圖8 不同姿態的小鼠。(a)正常放置姿態；(b)脊柱與后肢蜷縮姿態；(c)脊柱側向彎曲姿態。Fig.8 Mice with different postures. (a)Normal posture; (b)Spine and lower limbs curled up; (c)Spine lateral-bending.

作為測試，首先對以上4類姿態每類選出6幅圖像，共24幅圖像。然后，從第1類姿態的6幅中選出一幅作為測試圖像，其他23幅作為訓練圖像。一共進行了5次實驗，第1次只用了6幅屬于第2類姿態的訓練圖像，第2次增加了6幅第3類姿態的圖像，第3次再增加6幅第4類姿態的圖像。第4次實驗，從已有的18幅訓練圖像中去掉6幅(每類姿態2幅)，增加5幅1類姿態的訓練圖像。第5次實驗，使用全部23幅測試圖像。

圖9展示了這5次實驗的結果，圖的橫軸為實驗的序號，縱軸為識別的誤差。由于小鼠大多數姿態中身體軀干部分變化很少，而四肢的姿態變化很大，所以將軀干與四肢分開進行誤差分析。通過觀察箱線圖的四分位間距，可以發現，位于四肢關節的識別穩定性收斂比軀干關節穩定性收斂快，這是由于軀干脊柱上各個脊椎的相似性大，其脊柱周圍的灰度值分布復雜，算法更難精確定位脊柱上的關節點，而四肢關節特征明顯辨識度強，識別的穩定性收斂更快。

圖9 訓練數據量與識別誤差箱線圖。(a)全身關節誤差分析；(b)四肢關節誤差分析；(c)軀干關節誤差分析Fig.9 Training data number and localization error.(a)Whole- body joints error analysis; (b)Limbs joints error analysis; (c)Trunk joints error analysis

對比5次實驗，前3次實驗的訓練數據數目每次增加6幅圖像，但是都不包含與測試圖像相似的小鼠姿態。得出結論是：訓練數據越多，識別的誤差越小。第4次與第3次的訓練樣本數量相當，但是第4次包含了與測試小鼠姿態相同的圖像，對比實驗3、4的結果可以得出結論：當訓練數據包含測試小鼠的姿態，識別的誤差將減小并且更加穩定。所以訓練數據應該盡可能包含多種小鼠姿態。第5次實驗包含全部23幅訓練圖像，而誤差卻沒有比第4次明顯改善，對于軀干關節的穩定性甚至比第4次略差，說明訓練樣本的數目已經基本達到了飽和。

3 討論

從實驗結果可以看出，本研究方法對各個關節的定位誤差中位數可以控制在0.63 mm，相當于3.16個體素，關節識別的成功率達到92%以上。定位精度主要受回歸森林影響，而成功率主要由分類森林決定。影響定位精度和成功率的主要因素都是部分圖像中踝關節的缺失。因此，今后的研究有必要重點解決踝關節缺失的問題。

通過對比“分類森里+回歸森林+圖匹配”和“分類森林+圖匹配”，可以明顯看出后者在精度和穩定性方面的優勢。但是，兩種方法中位數誤差之間的差距在兩個體素以內，并且脊柱和頭骨連接處、骨盆前和右手聯合使用回歸森林的效果略差于單獨使用分類森林。所以在嚴格要求識別精度的情況下，有必要使用回歸森林，但是考慮到定位時間和相對低精度的要求，可以考慮只使用分類森林結合匹配圖的方式進行定位。

通過5次不斷增加訓練樣本的實驗，可以看到訓練樣本數目對于精度和穩定性的重要作用，訓練樣本越多，精度和穩定性越好。為了使算法的精度和穩定性達到飽和區，本方法所需的訓練樣本數最小應該在18以上。還可看到，在訓練集中包含于測試樣本相同姿態的訓練數據對于提高精度和穩定性也很重要。因此，為了達到對不同姿態的良好定位精度，應盡量在訓練集中包含各種身體姿態。在未來的研究工作中，需要收集更多不同姿態的影像，將其加入訓練集，從而進一步提高算法對不同姿態的普適性。

本研究僅涉及到了micro- CT這一種小動物成像模式，為了進一步拓寬算法的適用范圍，今后還將考慮對核磁、PET、SPECT等影像模式開發類似的關節點定位算法。由于其他影像模式對于骨組織的對比度分辨率不如CT好，后續研究的難度可能會更大。

雖然本研究的小鼠關節點定位問題對于小動物影像分析具有較重要意義，但是國內外尚無同類研究。正如本文引言部分所介紹，現有的各種小鼠micro- CT圖像分析方法多側重于全身圖譜的配準，但是其配準精度又都受到小鼠姿態變化的限制，足以說明本研究方法對于現有研究有促進意義。對比人類臨床CT圖像中的關節定位精度，Donner等[4]采用隨機森林算法得到的中值誤差為2.23個體素，與本研究相當。

4 結論

本研究提出了一種較為準確穩定的算法用于定位小鼠micro- CT圖像中的骨骼關節點。本方法使用分類和回歸隨機森林實現了對CT圖像關節點備選點的選擇，然后應用圖匹配算法根據關節相對位置關系來篩選備選關節點。基于49幅包含不同身體姿態的測試圖像，本研究方法的中值誤差為0.63 mm，均值誤差為0.97 mm(中值誤差為3.16個體素，均值誤差為4.80個體素)。回歸森林的使用，對于提高定位精度有重要作用，在訓練集中包含各種不同的身體姿態對于提高精度十分關鍵。本研究方法對于腳踝缺失的CT數據的識別效果不夠理想，但是對完整的CT數據包括小鼠全身的圖像識別成功率很好。本研究為小鼠micro- CT影像的姿態識別和自動分析提供了新的方法途徑。

(致謝：衷心感謝美國加州大學洛杉磯分校Crump分子影像中心的Arion Chatziiannou教授為本研究提供小動物CT影像數據。)

[1] Zhang Guanglei, Liu Fei, Pu Huangsheng, et al. A direct method with structural priors for imaging pharmacokinetic parameters in dynamic fluorescence molecular tomography [J]. IEEE Transactions on Biomedical Engineering, 2014, 61(3): 986- 990.

[2] Wang Luyao, Zhu Jun, Liang Xiao, et al.Performance evaluation of the Trans- PET?BioCaliburn?LH system: A large FOV small- animal PET system [J]. Physics in Medicine and Biology, 2015, 60(1): 137- 50.

[3] Wang Huina, Liu Chengbo, Gong Xiaojing, et al. In vivo photoacoustic molecular imaging of breast carcinoma with folate receptor- targeted indocyanine green nanoprobes [J]. Nanoscale, 2014, 6(23): 14270- 14279.

[4] 楊昆, 李真, 楊永鑫. 小動物正電子發射斷層成像儀探測器發展 [J]. 中國生物醫學工程學報, 2014, 33(2): 218- 226.

[5] Kesner AL, Dahlbom M, Huang SC, et al. Semiautomated analysis of small- animal PET data [J]. Journal of Nuclear Medicine, 2006, 47(7): 1181- 1186.

[6] Baiker M, Vastenhouw B, Branderhorst W, et al. Atlas- driven scan planning for high- resolution micro- SPECT data acquisition based on multi- view photographs: a pilot study [C]// Michael IM, Kenneth HW, eds. SPIE Medical Imaging 2009: Visualization, Image- Guided Procedures, and Modeling. Lake Buena Vista: SPIE, 2009: 72611L1- 72611L8.

[7] Baiker M, Staring M, L?wik CWGM, et al. Automated Registration of whole- body follow- up micro CT data of mice [C]//Medical Image Computing and Computer- Assisted Intervention- MICCAI 2011. Berlin: Springer Berlin Heidelberg, 2011: 516- 523.

[8] Baiker M, Milles J, Dijkstra J, et al. Atlas- based whole- body segmentation of mice from low- contrast Micro- CT data [J]. Med Image Anal, 2010, 14(6): 723- 737.

[9] Wang H, Stout DB, Chatziioannou AF. A method of 2D/3D registration of a statistical mouse atlas with a planar X- ray projection and an optical photo [J]. Med Image Anal, 2013, 17(4): 401- 416.

[10] Wang H, Stout DB, Chatziioannou AF. Estimation of mouse organ locations through registration of a statistical mouse atlas with micro- CT images [J]. IEEE Trans Med Imaging, 2012, 31(1): 88- 102.

[11] Shotton J, Sharp T, Kipman AA, et al. Real- time human pose recognition in parts from single depth images [J]. Communications of the ACM, 2013, 56(1): 116- 124.

[12] Donner R, Menze BH, Bischof H, et al. Global localization of 3D anatomical structures by pre- filtered Hough forests and discrete optimization [J]. Med Image Anal, 2013, 17(8): 1304- 1314.

[13] 林修煅. Micro- CT成像系統及其應用研究 [D].西安：西安電子科技大學,2012.

[14] 方匡南，吳見彬，朱建平，等. 隨機森林方法研究綜述 [J].統計與信息論壇, 2011, 26(3):32- 38.

[15] Van Ryzin J, Breiman L, Friedman JH, et al. Classification and Regression Trees [J]. Journal of the American Statistical Association, 1986,33(1):128- 128.

[16] Schapire RE. The strength of weak learnability[J].Machine Learning, 1990, 5(2): 197-227.

[17] Amit Y, Geman D. Shape quantization and recognition with randomized trees [J]. Neural Computation, 1997, 9(7): 1545- 1588.

[18] Breiman L. Random forests [J]. Machine Learning, 2001, 45(1): 5- 32.

[19] 李欣海. 隨機森林模型在分類與回歸分析中的應用 [J]. 應用昆蟲學報, 2013, 50(4): 1190- 1197.

[20] Criminisi A, Jamie Shotton. Decision Forests for Computer Vision and Medical Image Analysis[M]. Berlin: Springer Science & Business Media, 2013:340- 366.

[21] Pushmeet K, Philip HST. Efficiently solving dynamic markov random fields using graph cuts [C]// IEEE International Conference on Computer Vision (ICCV 05). Beijing: IEEE, 2005: 922- 929..

[22] Boykov Y, Veksler O, Zabih R. Fast approximate energy minimization via graph cuts[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2001, 23(11):1222- 1239.

Bone Joints Localization in Mouse Micro- CT Images Using Random Forests Algorithm

Tu Ruibo Chen Zhonghua Wang Hongkai#*

(Department of Biomedical Engineering, Dalian University of Technology, Dalian 116024, Liaoning, China)

Along with the rapid development of imaging techniques for small animals, more and more images obtained from small animals need to be analyzed per day, therefore automated image analysis method has become an urgent demand. For mice images, the significant inter- subject posture variations become a major difficulty for automated image analysis. In this paper, an automatic bone joint localization method was developed for mouse micro- CT images, so as to help with posture identification of mouse body. The proposed method was composed of three steps: (1) classification random forests for rough joint localization, (2) aggregating the results of classification through regression forest, and (3) picking up landmarks in the right position by the mapping graph. The method achieved automatic bone joint localization for 49 test images of different body postures. The median localization error of the whole body CT images was 0.68 mm. The success rate of localization was 98.27%. We also demonstrated the necessity of combining classification and regression random forest and discussed the influence on localization with different number of training data. With this new method for mouse micro- CT posture identification was expected to provide helpful information for the subsequent image registration, segmentation and measurements.

small animal image analysis；bone joint localization；random forest；pattern recognition；micro- CT

10.3969/j.issn.0258- 8021. 2017. 03.001

2016-08-19， 錄用日期:2016-10-21

國家自然科學基金(61571076)；國家自然科學基金青年基金(81401475)；遼寧省自然科學基金(2015020040)

R318

A

0258- 8021(2017) 03- 0257- 10

*通信作者(Corresponding author)，E- mail: wang.hongkai@dlut.edu.cn