基于直方圖熵的體數(shù)據(jù)分類算法研究

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

關(guān)鍵詞：體數(shù)據(jù)分類； 熵； 累計直方圖； 體繪制

中圖分類號：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]主要處理定義在三維或更高維網(wǎng)格上的標量和向量數(shù)據(jù)，是科學可視化最活躍的研究子領(lǐng)域之一^[3]。與傳統(tǒng)的計算機圖形學相比，三維體繪制的優(yōu)點在于不僅能描述物體的表面，而且能揭示物體內(nèi)部復雜的結(jié)構(gòu)，使人們看到通常情況下看不到的物體內(nèi)部結(jié)構(gòu)。三維體繪制通常包含兩個主要的步驟，即體數(shù)據(jù)的分類和體數(shù)據(jù)的繪制。體數(shù)據(jù)的分類主要是給體數(shù)據(jù)中的每個體素分配一個阻光度，以確定體素的可見性；體數(shù)據(jù)的繪制是根據(jù)視線以及體素的可見性將三維體數(shù)據(jù)投影到二維圖像平面，顯示感興趣的特征。因而，體數(shù)據(jù)的分類對三維體繪制的結(jié)果有很大的影響，只有對三維體數(shù)據(jù)進行準確的分類，把體數(shù)據(jù)中蘊涵的各種各樣的物質(zhì)區(qū)分開，才能經(jīng)過后續(xù)處理，有選擇、有目的對體數(shù)據(jù)進行顯示，得出合理的圖像。

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

1 體數(shù)據(jù)的構(gòu)成及阻光度函數(shù)模型

1.1 體數(shù)據(jù)的構(gòu)成

本文對規(guī)則體數(shù)據(jù)進行分類，三維體數(shù)據(jù)來自一系列連續(xù)的二維CT圖像，這些圖像也被稱為切片。切片按照圖1中的規(guī)則排列起來構(gòu)成三維體數(shù)據(jù)。假設(shè)三維體數(shù)據(jù)由n個切片構(gòu)成，每個切片的尺寸為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 阻光度傳遞函數(shù)模型

體數(shù)據(jù)的分類是通過阻光度傳遞函數(shù)實現(xiàn)的，為了靈活地顯示物體中任意感興趣的部分，阻光度傳遞函數(shù)通常為分段函數(shù)。設(shè)v(x，y，z)、α(x，y，z)分別為體素(x，y，z)的灰度值和阻光度值。假設(shè)在體數(shù)據(jù)中有n種待顯示的物質(zhì)，第m(m=1，2，…，n)種待顯示的物質(zhì)的灰度范圍是vm~vm+1(vm＜vm+1)，與vm對應(yīng)的阻光度值為αm。本文采用的阻光度傳遞函數(shù)能突出表面強度和物體實體感。其模型如式(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)可以看出，在同一灰度區(qū)間，體素的阻光度值與該體素的灰度與梯度有關(guān)，同時，受阻光度傳遞函數(shù)分段點的影響。本文提出的算法首先由二維CT圖像建立體數(shù)據(jù)，假定體數(shù)據(jù)中包含n種物質(zhì)，則需要計算出n-1個分段點將阻光度傳遞函數(shù)分為n段；然后計算體數(shù)據(jù)的直方圖，并根據(jù)其累計直方圖將體數(shù)據(jù)的直方圖劃分成n-1段；最后根據(jù)體數(shù)據(jù)直方圖的熵在每一個分段中計算出一個閾值作為阻光度傳遞函數(shù)的分段點，利用式(2)對每個體素進行阻光度賦值，完成體數(shù)據(jù)的分類。其流程如圖2所示。

3 體數(shù)據(jù)直方圖的分段

3.1 體數(shù)據(jù)直方圖的定義

對于按照1.1節(jié)中的方式構(gòu)成的體數(shù)據(jù)，本文定義其直方圖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的體素個數(shù)；g(i)表示灰度為i的體素的標記。

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

0其他 (5)

體數(shù)據(jù)的直方圖通常存在許多毛刺，在利用直方圖進行體數(shù)據(jù)相關(guān)特征（如直方圖的熵、矩等）計算時容易影響其準確性，所以在很多情況下要將直方圖進行平滑，提高直方圖特征計算的精度。本文所采用的計算方法如式（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是平滑區(qū)域的單邊寬度；k(l，σ)是高斯函數(shù)為

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

3.2 基于累計直方圖的體數(shù)據(jù)直方圖分段方法

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

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

(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）計算每個區(qū)間的累計直方圖

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）并保證分段區(qū)間相互不重疊。

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

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

4 阻光度傳遞函數(shù)分段點的計算

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

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

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

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

設(shè)體數(shù)據(jù)的直方圖為h(i)，某分段區(qū)間(fp，fp+1)內(nèi)直方圖的最大灰度值和最小灰度值用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)的體素出現(xiàn)的概率，假設(shè)分段區(qū)間(fp，fp+1)中的閾值為tp，tp將體數(shù)據(jù)直方圖的該分段區(qū)間分成兩部分，分別用A和B表示，設(shè)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）

定義與這兩個概率分布相關(guān)的熵為

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）。

該體數(shù)據(jù)分段的總熵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值即為分段區(qū)間(fp，fp+1)的閾值。計算出每個分段區(qū)間的閾值，由tp(p=1，2，…，n)構(gòu)成阻光度傳遞函數(shù)的分段點。

5 實驗結(jié)果和結(jié)論

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

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

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

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

參考文獻：

[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]唐澤圣. 三維數(shù)據(jù)場可視化 [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]管偉光. 體視化技術(shù)及其應(yīng)用[M].北京:電子工業(yè)出版社，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.

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