基于ITK醫學圖像分割的應用

宋 靜， 何啟泠， 徐 慶， 王海江

(成都信息工程大學電子工程學院，四川成都610225)

0 引言

醫學影像是現代醫學診斷和輔助治療最重要的技術手段之一，使醫生能夠對人體解剖結構及功能的認識更加直觀與全面，影像診斷已經廣泛應用于醫療診斷、治療、術前計劃和術后監測的各個環節。隨著互聯網醫療的快速發展，醫學影像診斷相關的應用越來越朝著輕量級的方向發展，使醫學圖像可以在多種應用中展示與瀏覽而不是依賴于專用的圖像瀏覽軟件。通過各種醫學圖像處理方法幫助醫生提高診斷的準確性的同時，也在思考如何更加方便醫生進行診斷，如更方便與快捷地遠程影像診斷。圖像分割是圖像處理中的關鍵性技術，對圖像進行進一步分析和處理的前一過程基本都依賴于圖像分割。圖像分割是指根據圖像的不同特性(圖像的灰度、顏色、紋理等)分割成不同的區域，這一處理有助于提取目標圖像，如某一組織或是病變位置。

在互聯網醫療時代，大量的醫學圖像通過傳輸完成診斷，交流探討。在這一過程中對于圖像傳輸的質量和效率的控制直接影響互聯網與醫療的結合度以及用戶體驗。大量的醫學圖像對診斷有用的區域與圖像的實際大小總是存在一定的差值，就會有大量的無效數據加大存儲量與降低傳輸效率。另一方面在三維重建的過程中，圖像中的非有效成像區域對于三維可視效果也有極大的影響，改善這兩點正是研究的意義所在。ITK主要提供了醫學影像圖像分割和配準方面的功能，但是在可視化方面還要借助VTK(visualization toolkit)觀察結果以及進行交互顯示。

圖像分割在組織提取或區域提取方面已有很多研究，文中主要從基于圖像分割的實際應用方面出發，討論如何降低數據的傳輸量，提高圖像傳輸速率，這一點對于開發PACS系統和研究遠程醫療具有一定的意義。此外大量的醫學影像存在設備干擾影像，在單幅影像中，可能這些干擾影像還不影響醫生進行診斷，但是在對這些醫學影像進行三維重建或者一些特定處理時，將嚴重影響成像效果。文中也討論圖像分割在這方面的應用。

1 ITK簡介

ITK是一款醫學圖像處理軟件包，是一個開源，跨平臺的影像分析擴展軟件工具。ITK采用先進的多模態數據分割配準算法，是一個用于處理醫學圖像的開源軟件庫，其中有豐富的圖像分割與配準的算法程序，并且支持多種開發平臺和開發語言，方便開發者進行圖像處理方面的研究與應用。

ITK作為一個開源，跨平臺的圖像分析框架，里面大量的前沿算法實現可以廣泛用于圖像配準和分割［1］。ITK的主要特征有:(1)ITK提供了通用的方式表示圖像(任意維度)和面片(非結構化的meshes)，提供了大量用于進行分割和配準的算法，主要面向醫學應用。(2)ITK只提供最小化的文件接口而不提供可視化的用戶接口，留由其他庫提供。(3)大量使用泛型編程技術。(4)內存模型使用智能指針維護對象的引用計數，使用對象工廠實例化對象。(5)支持多線程并行處理。(6)使用命令/觀察者模式進行事件處理。(7)基于數據流的架構進行組織，數據被表示成數據對象，數據對象由處理對象處理，數據對象和處理對象連在一起形成Pipeline。(8)ITK主要采用管道模塊結構設計，在對于圖像處理后需要將圖像文件輸出，需要借助VTK顯示圖像。對于此類需求可通過ITK，VTK聯合編程實現。

許多圖像處理方法可以在MATLAB中實現，但ITK比MATLAB而言優點主要體現在幾個方面:(1)ITK提供了大量的完備的與可靈活調用的圖像處理算法，如圖像分割，圖像配準與濾波等。(2)有效地將ITK的圖像讀取與處理功能與VTK在圖像可視化方面的突出能力結合起來可以實現強大的圖像處理功能。(3)ITK作為一個開源的平臺，會隨著時間不斷的豐富并且應用會有更好的延伸性。文中的主要研究內容是基于圖像分割的實際應用，重點在應用性的拓展上而不是圖像處理算法的實現上，因此ITK平臺相對于MATLAB而言是更佳的選擇。

2 醫學圖像分割算法應用

用于診斷的各種掃描圖像中，包含大量的信息，如人體器官，病變組織及某些物體。為分析它們的形狀、形態和功能，或者對它們進行可視化處理，都需要先將其從所在的情景中提取出來。這一提取過程即為醫學圖像分割。醫學圖像分割算法包括閾值分割、區域增長、基于分水嶺的分割、FastMarching算法、LevelSet(水平集)等多種分割方法［2］。文中使用的方法主要是區域生長與閾值分割。由于篇幅有限，主要介紹使用到的兩種方法以及其在應用過程中所作的調整。為提取目標圖像區域，采用先對原圖進行二值化處理，再對圖像進行區域生長，得到最大的連通區域，除去原圖像中的設備干擾成像。

閾值分割就是確定一個或多個閾值，根據圖像中像素的灰度值與確定閾值的大小關系進行一個映射，從而達到分割效果［3］。設原圖像f(x，y)，分割后的圖像為g(x，y)，新得到的圖像為二值圖像。原圖和分割后的圖像與閾值關系為

在這一應用中，根據二值圖像可以大體確定圖像所在的區域，研究常見圖像的提取情況，根據二值圖像提取目標圖像如圖1～圖6所示(圖像均為提取數據在網頁顯示效果截圖)。

圖1 原圖

圖2 二值圖像

圖3 根據二值圖像提取圖

圖4 原圖

圖5 二值圖像

圖6 根據二值圖像提取

根據結果發現在提取無干擾成像的圖像時，僅根據二值圖像來提取就可以達到預想的目的，但是當圖像中有干擾成像時提取的情況還不夠理想，根據這一情形想到再提取圖像前還需要對二值圖像作進一步的處理，達到提取的目標數據在保證提取到完整圖像的同時要盡可能的減少數據量。通過分析大量醫學圖像可見醫學影像圖像的目標總是在圖像的中心而不是邊緣或角落。把目標當作一個完整的連通區域，干擾圖像就可被屏蔽，通過實驗發現ITK中的區域生長中的連通閾值提取的方法正好可以解決這一情況。

區域生長的基本思想是提取有相似性質的像素從而構成區域，對每個區域都要選取一個種子點，以和種子像素相似或相同的性質為準則在種子點周圍開始生長，合并滿足條件的點，然后將這些點再作為種子開始生長，直到沒有滿足條件的像素點為止停止生長。在這一過程中有3個關鍵點:分別是種子點的選取;生長準則的確定;區域生長的停止條件［4-5］。

種子點的選擇:應用中，由于要求在應用過程中是后臺自動處理生成，這一點決定的種子點是不能通過交互來得到的，因此文中選擇總是從圖像的中心開始生長。這一方法也會產生一個新的問題，即如果圖像的中心在二值化的過程中是一個不可生長的像素如何處理，針對這一情況，做了兩個幫助生長的處理。首先由根據圖像二值化得到的二值圖像進行統計，可得到有效像素點數與整個圖像中像素點數的比例，根據這一比例在圖像中心構造一個有效區域，確保生長順利進行。這一方法是基于醫學影像的目標區域總是在圖像的中心。第二個幫助生長的處理是因為在生長過程中發現有一些希望得到的小區域會在生長過程中丟失，于是就采取先對二值圖像作膨脹處理，確保在生長過程中不會丟失掉一些邊緣信息。表1為中心構造區域大小與圖像有效像素與所有像素的比值關系。

表1 有效區域比例與中心構造區域面積關系表

在實驗過程中，根據以上介紹的方法可以準確提取出目標圖像，并且在無干擾成像的情況下依然可以準確提取出目標圖像數據，如圖7～圖9所示。

圖7 原圖

圖8 連通區域圖

圖9 根據連通區域圖像提取圖

3 基于ITK中圖像封裝類的實現過程

用ITK進行醫學圖像的處理會比較方便，因為大量的圖像處理算法已經在ITK中封裝好，在編寫程序的過程中只需調用相應的庫文件就可加以應用。目的主要是為了討論如何提取出有效的數據信息，從而可以在網頁繪制過程中更快，在傳輸過程中更快。圖10為基于ITK實現這一圖像處理過程的流程圖。

圖10 圖像處理過程流程圖

3．1 圖像二值化處理

用二值化處理過程的目的主要是為后續的處理過程，如膨脹處理，目標連通區域生成［6］。選擇Otsu方法來生成二值圖像，作為一種常用的灰度圖像二值化的方法，用一個基于直方圖統計得到的灰度閾值來對圖像進行分割，從而把圖像分為目標與背景兩部分，可以幫助確定目標區域在圖像中的所在位置。在ITK中已經把這一方法封裝，直接調用方法即可，其閾值是基于統計自動得到，從而可自動對圖像進行分割。

用ITK中的itk::OtsuThresholdImageFilter類對輸入圖像進行二值化處理。其中把SetOutsideValue()設置為0，SetInsideValue()設置為255，分別代表處理后閾值范圍外/內的像素值;SetNumber of Histogram Bins()設置為256表示用于計算分割閾值的直方圖的灰度級數。

const unsigned int Dimension=2;//定義數據的維數

typedef signed short InputPixelType;//定義讀取對象的類型

typedef itk::Image＜InputPixelType，Dimension ＞InputImageType;//定義輸入圖像的類型

typedef itk::ImageFileReader＜InputImageType＞ReaderType;

ReaderType::Pointer reader=ReaderType::New();//建立一個讀取對象

reader-＞SetFileName(InFile);//讀取圖像文件

typedef itk::OtsuThresholdImageFilter＜InputImageType，InputImageType ＞FilterType;

FilterType::Pointer filter=FilterType::New();//定義FilterType指針

filter-＞SetInput(reader-＞GetOutput());//獲得讀取的CT圖像

const InputPixelType outsideValue=255;//設置分割的低閾值為255

const InputPixelType insideValue=0;//設置分割的高閾值為0

filter-＞SetOutsideValue(outsideValue);

filter-＞SetInsideValue(insideValue);

filter-＞SetNumberOfHistogramBins(256)

filter-＞Update();

3．2 二值化圖像膨脹處理

圖像膨脹就是用固定結構元素在二值圖像上移動，掃描圖像中的每一個元素，即用結構元素中的像素與圖像像素做“與”，如果相與結果全部為0，則該像素置0，反之置1［7-8］。在這里加入膨脹處理的原因是由二值圖像直接得到目標連通區域時可能會丟失一些邊緣，通過膨脹處理就可以使目標增大，這樣在生成連通區域圖時就盡可能的包含圖像本身的信息，盡量避免在由二值化到連通區域生成的過程中錯誤的丟失掉原本屬于目標區域的圖像。用ITK中的 itk:BinaryDilateImageFilter類對上一過程輸出的二值圖像進行膨脹。這一過程主要是為的保存目標圖像的一些邊緣信息。將結構元素的半徑設置為5，調用方法為structuringElement.SetRadius()。通過 GrayscaleDilateImageFilter類定義了一個膨脹處理對象grayscaleDilate，因此膨脹處理后的輸出為grayscaleDilate-＞GetOutput()。

原圖不經過膨脹處理和經過膨脹處理提取的圖像效果如圖11～13所示。

圖11 原圖

圖12 未經過膨脹處理的提取圖

圖13 經過膨脹處理的提取圖

可看出沒經過膨脹處理得到的結果可能會丟失掉一些邊緣信息，所以在處理過程中加入了膨脹處理這一過程，提高圖像提取時的完整性。

3．3 得到目標圖像連通區域

ITK中的ConnectedThreshold分割算法可以得到目標圖像的連通區域［9-11］。圖像中存在的噪聲有可能會降低濾波的效果因此先對圖像進行一個平滑處理，接下來選取種子點如上面介紹選取圖像的中心點作為種子點。還要通過SetLower()，SetUpper()兩個方法來設置圖像分割的高閾值和低閾值，這里設置時要特別注意取值范圍要包含二值圖像的有效值的灰度。SetReplaceValue()設置新的有效值。SetSeed()加載種子點索引。其中size代表原圖像的行列數。

typedef float InternalPixelType;//定義用于區域生長的對象類型

typedef itk::Image＜ InternalPixelType，Dimension ＞InternalImageType;//定義用于區域生長的圖像類型

typedef itk::ConnectedThresholdImageFilter＜InternalImageType，InternalImageType ＞ ConnectedFilterType;

ConnectedFilterType::PointerconnectedThreshold=ConnectedFilterType::New();

//定義ConnectedFilterType類型指針

connectedThreshold-＞SetInput(caster-＞GetOutput());

//將膨脹處理后輸出轉換數據類型后作為提取連通區域的輸入即caster-＞GetOutput()

connectedThreshold-＞SetLower(125);

connectedThreshold-＞SetUpper(500);

connectedThreshold-＞SetReplaceValue(255.0);

InternalImageType::IndexType index;//定義種子點

index［0］=size［0］ /2;//種子點設置為圖像的中心區域

index［1］=size［1］ /2;

connectedThreshold-＞SetSeed(index);

connectedThreshold-＞Update();

3．4 根據連通區域提取目標區域圖像數據

通過以上幾步處理就可得到目標區域的二值圖。掃描這二值圖像可以得到圖像起止點矩形區域的坐標(x1，y1，x2，y2)，根據這一矩形區域坐標索引從原圖中提取數據，提取出來的數據可直接在頁面繪制，這一數據為圖像的原始數據所以不會影響后續的圖像操作如調窗、反顯、偽彩等，后續如果還要進行圖像分割處理也不會有任何影響。只是通過這一操作節約前端處理的操作時間和傳輸時間。其中原始像素數據可用迭代法提取，速度會快于遍歷提取，主要代碼如下:

選取不同醫院設備提取的CT，MR(dicom圖像)進行實驗，處理前后數據量見表2。

表2 圖像數據提取前后大小對比

這一處理數據量的變化不是固定不變的，與目標區域和原圖的比例緊密相關，從總體看對數據量的減少能起到很大作用。這一處理還可以用在三維模型生成前的圖像處理，主要是去除在成像過程中設備成像的干擾。這里對ITK生成三維模型［12-13］就不贅述，實驗效果如圖14、15所示(由VTK顯示的面繪效果)。

圖14 未處理過的圖像模型

圖15 處理過后的圖像模型

4 結束語

ITK強大的圖像處理功能可以方便開發者快捷的編制程序，提高開發效率。圖像壓縮可以解決圖像傳輸效率問題，但是在醫學圖像的處理過程中對圖像的無損性要求相對更高，以上介紹的提取數據的方法對于醫學影像的目標區域而言始終保存原始數據，因為是從原始圖像中提取區域數據，對于單個圖像而言可能有更簡化的方法提取，但是對于要處理大量醫學圖像的軟件應用層次而言文中介紹的方法具有一定的適用性。無需交互特性使得應用范圍可以得到擴展。

［1］ 張芳，李強，王陽萍.ITK及其在醫學圖像分割中的應用［J］.微計算機信息，2008，24(18):304-306.

［2］ 徐彩云.圖像分割算法的研究綜述［J］.電腦知識與技術，2014，(11).

［3］ 王磊.多維Otsu方法在圖像分割中的研究［D］.濟南:山東師范大學，2009.

［4］ 李敏.數字人腦切片圖像自動分割算法的初步研究［D］.重慶:重慶大學，2010.

［5］ 鄧金蓮，李成貴.結合區域生長和 SUSAN的圖像分割方法［J］.電光與控制，2012，19(8):38-40.

［6］ 常江.車輛圖像局部識別［J］.信息網絡安全，2013，(11):94-96.

［7］ 黃穎，楊光瓊.結合小波系數的normalized cut分割算法［J］.計算機應用，2011，31(1):182-183.

［8］ 喬誠，李連東，王秀麗.ITK及其在分水嶺醫學圖像分割中的應用［J］.太原科技，2008，(2).

［9］ 牟春潔，張國英.基于區域邊界生長的圖像分割方法［J］.北京石油化工學院學報，2009，17(4):8-12.

［10］ 呂曉琪，任曉穎，賈東征.基于 ITK，VTK和MFC的DICOM圖像讀寫及顯示［J］.中國組織工程研究，2011，15(13):2416-2420.

［11］ 伍云智.基于VTK和ITK算法庫的研究與應用［D］.南昌:南昌大學，2010.

［12］ 袁杲，楊玲，朱小波，等.利用ITK和VTK集成實現三維醫學圖像的分割［J］.計算機應用與軟件，2009，26(2):248-250.

［13］ 陳保遠.基于ITK和VTK醫學圖像配準方法的研究與應用［D］.南昌:南昌大學，2012.