基于隨機森林算法的小鼠micro- CT影像中骨骼關(guān)節(jié)特征點定位

屠睿博 陳中華 王洪凱

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

基于隨機森林算法的小鼠micro- CT影像中骨骼關(guān)節(jié)特征點定位

屠睿博 陳中華 王洪凱#*

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

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

小動物影像分析； 骨關(guān)節(jié)點定位； 隨機森林； 模式識別； 顯微CT

引言

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

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

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

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

1 方法

1.1 隨機森林算法簡介

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

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

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

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

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

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

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

1.2 算法總體流程

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

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

1.3 分類森林

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

1.3.1 分類森林的數(shù)據(jù)預處理

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

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

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

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

1.3.2 分類森林的訓練

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

(1)

(2)

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

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

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

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

1.3.3 分類森林的測試

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

1.4 回歸森林

1.4.1 數(shù)據(jù)預處理

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

1.4.2 回歸隨機森林的訓練

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

本研究針對20個不同的關(guān)節(jié)點訓練20個不同的森林，預測某一點到某一個關(guān)節(jié)的距離只需輸入到相應(yīng)的森林中進行預判。

1.4.3 回歸隨機森林的測試

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

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

1.5 使用圖匹配方法篩選關(guān)節(jié)備選點

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

1.5.1 匹配圖的訓練

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

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

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

(3)

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

1.5.2 圖匹配方法的測試

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

(4)

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

(5)

(6)

(7)

(8)

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

1.6 算法評價方法

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

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

2 結(jié)果

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

圖4 多種姿態(tài)小鼠的關(guān)節(jié)識別結(jié)果。小鼠的骨骼以閾值分割的半透明三維面繪制方式展現(xiàn)，識別出的骨骼關(guān)節(jié)圓點形式標記圖像上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 算法準確性分析

這里采用均值誤差和中位數(shù)誤差評估定位關(guān)節(jié)點的準確率。

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

圖5 關(guān)節(jié)定位誤差Fig.5 Bone joints localization errors

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

2.2 算法識別成功率分析

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

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

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

圖6 關(guān)節(jié)識別成功率Fig.6 Bone joints recognition rate

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

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

圖7 關(guān)于是否使用回歸森林的誤差對比。(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)所示分別是“分類森林+回歸森林+匹配圖”和“分類森林+匹配圖”的結(jié)果箱線圖，縱坐標表示定位誤差，橫坐標為骨關(guān)節(jié)類別?？梢钥吹剑瑘D7(a)中各關(guān)節(jié)識別誤差的四分位間距超過5個體素的只有一個關(guān)節(jié)點(上胸骨柄)，其他識別誤差結(jié)果穩(wěn)定在上下四分位相差為2～3個體素，證明本研究所用方法識別小鼠骨關(guān)節(jié)穩(wěn)定性高。而圖7(b)中多數(shù)關(guān)節(jié)的中值誤差是圖7(a)的2倍左右，且上下四分位的間距也更大，說明如果沒有回歸森林，大部分關(guān)節(jié)的定位精度和穩(wěn)定性都會下降。

2.4 訓練數(shù)據(jù)量對算法精度的影響

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

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

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

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

圖9 訓練數(shù)據(jù)量與識別誤差箱線圖。(a)全身關(guān)節(jié)誤差分析；(b)四肢關(guān)節(jié)誤差分析；(c)軀干關(guān)節(jié)誤差分析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次實驗的訓練數(shù)據(jù)數(shù)目每次增加6幅圖像，但是都不包含與測試圖像相似的小鼠姿態(tài)。得出結(jié)論是：訓練數(shù)據(jù)越多，識別的誤差越小。第4次與第3次的訓練樣本數(shù)量相當，但是第4次包含了與測試小鼠姿態(tài)相同的圖像，對比實驗3、4的結(jié)果可以得出結(jié)論：當訓練數(shù)據(jù)包含測試小鼠的姿態(tài)，識別的誤差將減小并且更加穩(wěn)定。所以訓練數(shù)據(jù)應(yīng)該盡可能包含多種小鼠姿態(tài)。第5次實驗包含全部23幅訓練圖像，而誤差卻沒有比第4次明顯改善，對于軀干關(guān)節(jié)的穩(wěn)定性甚至比第4次略差，說明訓練樣本的數(shù)目已經(jīng)基本達到了飽和。

3 討論

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

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

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

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

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

4 結(jié)論

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

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

[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] 楊昆, 李真, 楊永鑫. 小動物正電子發(fā)射斷層成像儀探測器發(fā)展 [J]. 中國生物醫(yī)學工程學報, 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成像系統(tǒng)及其應(yīng)用研究 [D].西安：西安電子科技大學,2012.

[14] 方匡南，吳見彬，朱建平，等. 隨機森林方法研究綜述 [J].統(tǒng)計與信息論壇, 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] 李欣海. 隨機森林模型在分類與回歸分析中的應(yīng)用 [J]. 應(yīng)用昆蟲學報, 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