基于直方圖熵的體數據分類算法研究

摘 要：體數據的分類用于確定體素的可見性，在三維體繪制中起著重要的作用。提出一種基于熵的體數據分類算法。首先根據累計直方圖將體數據的直方圖進行分段，然后根據熵判別式在每個分段中計算一個閾值作為阻光度傳遞函數的分段點，再根據阻光度傳遞函數計算出體素的阻光度值，完成體數據的分類。以工業CT體數據為對象進行實驗，其結果表明，所提出的算法較好地實現了體數據的分類，體繪制結果清晰，且能夠實現試件的模擬拆卸。

關鍵詞：體數據分類； 熵； 累計直方圖； 體繪制

中圖分類號：TP391 文獻標志碼：A 文章編號：1001-3695(2008)08-2410-03

Classification of volume data based on entropy of histogram

LI Jin， ZHOU Lu-lu， YU Hong， LIANG Hong

( College of Automation， Harbin Engineering University， Harbin 150001， China)

Abstract:The classification step used to assign the appropriate opacity to each voxel is very important in the volume rende-ring. This paper proposed a new classification algorithm for volume data， which was based on the entropy. Firstly， partitioned the histogram of the volume data into different subsections through calculating the accumulated histogram of volume data according to the number of object classes. Secondly， in each subsection computed a threshold based on the entropy. Finally， assigned the opacity of voxel by a transfer function， which was split into subsection by these thresholds. Experimental results from industrial CT images were presented and it indicated that the classification of volume data was obtained with the better accuracy by the proposed algorithm. The industrial components are shown explicitly in the volume rendering results and successfully disassembled in the computer simulated.

Key words：classification of volume data; entropy; accumulated histogram; volume rendering

0 引言

體繪制^[1，2]主要處理定義在三維或更高維網格上的標量和向量數據，是科學可視化最活躍的研究子領域之一^[3]。與傳統的計算機圖形學相比，三維體繪制的優點在于不僅能描述物體的表面，而且能揭示物體內部復雜的結構，使人們看到通常情況下看不到的物體內部結構。三維體繪制通常包含兩個主要的步驟，即體數據的分類和體數據的繪制。體數據的分類主要是給體數據中的每個體素分配一個阻光度，以確定體素的可見性；體數據的繪制是根據視線以及體素的可見性將三維體數據投影到二維圖像平面，顯示感興趣的特征。因而，體數據的分類對三維體繪制的結果有很大的影響，只有對三維體數據進行準確的分類，把體數據中蘊涵的各種各樣的物質區分開，才能經過后續處理，有選擇、有目的對體數據進行顯示，得出合理的圖像。

多年來人們一直對體數據的分類問題進行研究，提出了一些用于體繪制的體數據分類算法，如使用特征距離圖譜的分類算法^[4]、基于小波域隱式馬爾可夫模型的分類算法^[5]、基于模糊分析的分類算法^[6]等，但體數據的分類一直是一個沒有徹底解決好的問題^[7]。本文提出了一種基于熵的體數據分類算法，較好地實現了工業CT體數據的分類。

1 體數據的構成及阻光度函數模型

1.1 體數據的構成

本文對規則體數據進行分類，三維體數據來自一系列連續的二維CT圖像，這些圖像也被稱為切片。切片按照圖1中的規則排列起來構成三維體數據。假設三維體數據由n個切片構成，每個切片的尺寸為k×l，(x，y)j 是切片j上的一個像素， 其灰度用pj(x，y)表示，則體素(x，y，z)的灰度值v(x，y，z)滿足

v(x，y，i)=pi(x，y)； x，y，i∈I； 0＜i≤n； 0＜x≤k； 0＜y≤l(1)

1.2 阻光度傳遞函數模型

體數據的分類是通過阻光度傳遞函數實現的，為了靈活地顯示物體中任意感興趣的部分，阻光度傳遞函數通常為分段函數。設v(x，y，z)、α(x，y，z)分別為體素(x，y，z)的灰度值和阻光度值。假設在體數據中有n種待顯示的物質，第m(m=1，2，…，n)種待顯示的物質的灰度范圍是vm~vm+1(vm＜vm+1)，與vm對應的阻光度值為αm。本文采用的阻光度傳遞函數能突出表面強度和物體實體感。其模型如式(2)所示^[1]。

α(x，y，z)=|v(x，y，z)|αm[v(x，y，z)-vm）/（vm+1-vm)]+

αm+1[vm+1-v(x，y，z)/（vm+1-vm)]

vm≤v(x，y，z)≤vm+1

0 其他(2)

其中：v(x，y，z)為梯度算子，可以用下式計算：

v(x，y，z)=[v(x+1，y，z)-v(x-1，y，z)]/2

[v(x，y+1，z)-v(x，y-1，z)]/2

[v(x，y，z+1)-v(x，y，z)]/2(3)

2 分類算法描述

從式(2)可以看出，在同一灰度區間，體素的阻光度值與該體素的灰度與梯度有關，同時，受阻光度傳遞函數分段點的影響。本文提出的算法首先由二維CT圖像建立體數據，假定體數據中包含n種物質，則需要計算出n-1個分段點將阻光度傳遞函數分為n段；然后計算體數據的直方圖，并根據其累計直方圖將體數據的直方圖劃分成n-1段；最后根據體數據直方圖的熵在每一個分段中計算出一個閾值作為阻光度傳遞函數的分段點，利用式(2)對每個體素進行阻光度賦值，完成體數據的分類。其流程如圖2所示。

3 體數據直方圖的分段

3.1 體數據直方圖的定義

對于按照1.1節中的方式構成的體數據，本文定義其直方圖h1(i)為式（4）的形式：

h1(i)=∑x+1r=x-1∑y+1s=y-1∑z+1t=z-1g(i)/k×l×n; i=0，1，…，255（4）

其中：∑x+1r=x-1∑y+1s=y-1∑z+1t=z-1g(i)表示灰度為i的體素個數；g(i)表示灰度為i的體素的標記。

g(i)=1v(r，s，t)=i； i∈[0，255]

0其他 (5)

體數據的直方圖通常存在許多毛刺，在利用直方圖進行體數據相關特征（如直方圖的熵、矩等）計算時容易影響其準確性，所以在很多情況下要將直方圖進行平滑，提高直方圖特征計算的精度。本文所采用的計算方法如式（6）所示^[8]。

h(i)=∑Wl=1[h1(i-1)×k(l，σ)+h1(i+l)×k(l，σ)]（6）

其中:h1(i)是原始直方圖，可由式（4）求出；h(i)是平滑后的直方圖；W是平滑區域的單邊寬度；k(l，σ)是高斯函數為

k(l，σ)=exp(-(l×σ)2/2)（7）

3.2 基于累計直方圖的體數據直方圖分段方法

對于存在n種物質的體數據，其直方圖被n-2個點j1，j2，…，jn-2劃分成n-1個部分。圖3給出了一個例子，體數據包含五種物質，即n＝5。其中（a）為原始體數據直方圖，存在一些毛刺，不夠平滑；（b）給出了按照式（6）的方法濾波后的平滑體數據直方圖，并標記了直方圖的三個分段點j1、j2、j3。

本文提出了一種基于累計直方圖的體數據直方圖劃分方法。首先找出體數據直方圖的非零端點fmin和fmax，然后將灰度區間(fmin ， fmax)平均劃分成n-1個部分，每一個部分可以表示成區間的形式

(fp， fp+1)； p=1，…，n-1； fp＜fp+1（8）

端點fp可以由式（9）計算：

fp=fminp=1

p×(fmax+fmin)/np=2，…，n-2

fmaxp=n-1（9）

再用式（10）計算每個區間的累計直方圖

hp=∑fp+1i=fph(i)；p=1，2，…，n-1（10）

若hp＜1/(m+2)(p=0，1，…，m-1)，則該分段需要被擴展，并重新計算hp，直到滿足式（11）并保證分段區間相互不重疊。

hp≥1/(m+2)； p=1，2，…，n-1（11）

圖4給出了直方圖分段方法的流程和分段擴展的方法。

4 阻光度傳遞函數分段點的計算

熵是熱力學中的概念，在信息論中用熵可以作為信息量的度量，基于熵的理論已應用于很多領域，如語音識別^[9]、醫療診斷^[10，11]等。本文將熵的概念應用于體數據的分類，其可以看做是體數據統計特性的一種表現形式，是體數據信息量的度量。由式(4)中體數據直方圖的定義可知，體數據的直方圖可以看做是體素灰度的概率分布，根據熵的定義，體數據包含信息量的大小可以由下式計算：

H=-∑fmaxi=fminh(i) ln h(i)（12）

其中：fmin、fmax分別為體數據中體素的最小灰度值和最大灰度值；h(i)為體數據的直方圖。

KSW熵方法是由J. N. Kapur等人^[12]提出的，最早應用于二維圖像的單閾值分割，是基于兩個分布假定原理的方法。本文利用KSW熵方法在體數據直方圖的每個分段中自動找到合適的閾值，作為阻光度傳遞函數的分段點，可以實現體數據的二類及多類分類。

設體數據的直方圖為h(i)，某分段區間(fp，fp+1)內直方圖的最大灰度值和最小灰度值用fpmax、fpmin表示，滿足

fp≤fpmin＜fpmax≤fp+1（13）

則h(i)的灰度值為i(fpmin＜i＜fpmax)的體素出現的概率，假設分段區間(fp，fp+1)中的閾值為tp，tp將體數據直方圖的該分段區間分成兩部分，分別用A和B表示，設A用來描述灰度值為i(fmin≤i≤tp)的體素的分布，B用來描述灰度值為i(tp＜i≤fmax)的灰度分布，則它們的概率分布為

A:h(fpmin)/Pt，h(fpmin+1)/Pt，…，h(t)/Pt

B:h(t+1)/(Mt-Pt)，h(t+2)/

(Mt-Pt)，…，h(fpmax)/(Mt-Pt)（14）

其中:Pt=∑tpi=f pminh(i)（15）

Mt=∑f pmaxi=f pminh(i)（16）

定義與這兩個概率分布相關的熵為

HA(tp)=ln Pt+Ht/Pt（17）

HB(t)=ln(Mt-Pt)+(Hmax-Ht)/(Mt-Pt)（18）

式中，Hmax=∑f pmaxi=f pminh(i)ln h(i)（19）

Ht=-∑tpi=f pminh(i)ln h(i)（20）

Pt、Mt同式（15）（16）。

該體數據分段的總熵H(t)為HA(t)與HB(t)之和，即

H(tp)=ln Pt(1-Pt)+Ht/Pt+(Hmax-Ht)/(1-Pt)（21）

求使得H(tp)最大的tp值即為分段區間(fp，fp+1)的閾值。計算出每個分段區間的閾值，由tp(p=1，2，…，n)構成阻光度傳遞函數的分段點。

5 實驗結果和結論

本文給出了兩組工業CT體數據的分類實驗結果。一組體數據來自于帶有未焊透缺陷的鋼板，對鋼板進行CT掃描，得到15幅二維CT圖像，每張圖像的尺寸為512×512；另一組體數據來自于電鉆，對電鉆進行CT掃描，得到134幅二維CT圖像，每張圖像的尺寸為512×512。由這些圖像按照1.1節的方式構成體數據，按照本文方法對體數據進行分類后，用shear-warp算法進行三維體繪制，得到實物的三維構型。

鋼板體數據的分類實驗結果如圖5 所示。其中：(a)為鋼板的實物圖像，該鋼板帶有焊縫未焊透缺陷；(b)為鋼板經掃描后得到的二維CT圖像;(c)為鋼板體數據的近似阻光度曲線；(d) 為按(c)中阻光度值進行體繪制得到的鋼板不同側面的三維構型。電鉆體數據的分類實驗結果如圖6 所示。其中：(a)為電鉆的實物圖像；(b)為電鉆經掃描后得到的二維CT圖像；(c)為電鉆體數據的近似阻光度曲線；(d)(e)分別為按(c)中阻光度值進行體繪制得到的電鉆整體不同側面的整體的三維構型和電鉆內部鋼制結構不同側面的三維構型。

從圖5(b)可以看出，鋼板的二維CT圖像帶有嚴重的偽影，因此，其體數據中也包含偽影。利用本文的算法求得阻光度傳遞函數的分段點為126，可以較準確地實現該體數據的分類，其近似阻光度值如(c)所示。由于在同一種物質中，體素的梯度相近，在同一種物質中阻光度近似為線性變化。從(d)可以看出，利用本文的算法對體數據進行分類，其體繪制結果不受偽影的影響，且清晰地顯示出鋼板的構型及焊縫的未焊透缺陷。

從圖(b)可以得出，電鉆體數據中包含三種物質，利用本文算法求得阻光度傳遞函數的分段點為14和109，實現了該體數據的多種物質的劃分，將體數據中表示不同結構的體素區分開，可以通過改變阻光度來顯示感興趣的結構。(c)給出了顯示電鉆整體結構和內部鋼制結構的阻光度值，在同一種物質中阻光度近似為線性變化。從(d)(e)可以看出，利用本文的算法對體數據進行分類，其體繪制結果較為清晰，且能夠單獨顯示電鉆內部的鋼制結構，較好地實現了電鉆的模擬拆卸。

參考文獻：

[1]LEVOY M.Volume rendering:display of surfaces from volume data[J].IEEE Computer Graphics and Applications， 1988，8(5):29-37.

[2]DREBIN R A， CARPENTER L， HANRAHAN P. Volume rendering[J]. Computer Graphics， 1988， 22: 65-74.

[3]唐澤圣. 三維數據場可視化 [M]. 北京:清華大學出版社， 1999.

[4]SHIN B S. An efficient classification and rendering method using tagged distance maps[J]. Visual Computer， 2004， 20(8-9): 540-553.

[5]LIN Xue-yan， LIU Zheng-guang. Classification and rendering of brain MR imaging based on wavelet-domain hidden Markov mode[C]//Proc of International Conference on Imaging: Technology and Applications for 21st Century. Beijing: [s.n.]， 2005:342-343.

[6]PARVEEN， RUNA， TODD-POKROPEK， et al. Classification of MR brain tissues using fuzzy estimation[J]. Nuclear Science Sympo-sium Conference Record， 2006， 4: 2613-2619.

[7]管偉光. 體視化技術及其應用[M].北京:電子工業出版社，1998:2-4.

[8]楊靜宇，曹雨龍. 計算機圖像處理及常用算法手冊[K].南京:南京大學出版社，1997:173-209.

[9]ARONOWITZ H， BURSHTEIN D. Efficient speaker recognition using approximated cross entropy[J]. IEEE Trans on Audio， Speech and Language Processing， 2007， 15(7): 2033-2043.

[10]ABASOLO D， HORNERO R， ESCUDERO J， et al. Approximate entropy and mutual information analysis of the electroencephalogram in Alzheimer’s disease patients[C]//Proc of the 3rd International Conference on Advances in Medical，Signal and Information Processing. 2006:1-4.

[11]NATWONG B， SOORAKSA P， PINTAVRIOOJ C，et al. Wavelet entropy analysis of the high resolution ECG[C]//Proc of the 1st IEEE Conference on Industrial Electronics and Applications. 2006:1-4.

[12]KAPUR J N， SAHOO P K， WONG A K C.A new method for gray-level picture thresholding using the entropy of the histogram[J]. Computer Vision， Graphics， and Image Processing，1985，29(3): 273-285.

注：本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文