基于明膠實(shí)驗(yàn)的彈創(chuàng)空腔序列圖像處理技術(shù)*

沈 杰，王 濤，廖振強(qiáng)

(南京理工大學(xué)機(jī)械工程學(xué)院，南京 210094)

0 引言

空腔形成理論指出，高速飛行的殺傷元侵徹人體組織時(shí)，在組織中形成的瞬時(shí)空腔會(huì)經(jīng)多次膨脹與壓縮，使空腔周圍組織在極短的時(shí)間內(nèi)經(jīng)歷多次擠壓與牽拉[1]。在瞬時(shí)空腔膨脹波及的范圍內(nèi)，人體組織受到猛烈扯動(dòng)，發(fā)生明顯的位移，甚至?xí)毫押颓袛?，該區(qū)域范圍是組織器官受到嚴(yán)重創(chuàng)傷的區(qū)域[2]，因此研究瞬時(shí)空腔膨脹所能達(dá)到的最大空間范圍具有重要意義。為了檢測(cè)出彈創(chuàng)空腔最大膨脹波及范圍的邊緣輪廓，針對(duì)單幀圖像表示的瞬時(shí)空腔輪廓精確提取困難，且不能有效反映最大膨脹波及范圍等問題，文中提出序列瞬時(shí)空腔輪廓“或”運(yùn)算疊加法，算法重點(diǎn)在最終的空腔輪廓疊加結(jié)果，對(duì)單幅圖片處理的輪廓精度要求可以降低，而且在程序設(shè)計(jì)上更容易實(shí)現(xiàn)。

1 瞬時(shí)空腔最大膨脹波及范圍的定義

如圖1所示，設(shè)侵徹時(shí)刻從t1到t2時(shí)，彈創(chuàng)空腔輪廓由實(shí)線膨脹到虛線處，則在侵徹距離l1處，對(duì)應(yīng)t1和t2時(shí)刻的空腔直徑分別為Dl1t1和Dl1t2，在空腔膨脹周期內(nèi)，此距離處的空腔直徑不斷變化，故在此侵徹距離上必然存在一個(gè)空腔直徑的最大值Dl1max，同樣在侵徹距離為l2處也必然存在一個(gè)空腔直徑的最大值Dl2max，依次類推，即可得到一組不同侵徹距離處空腔直徑的最大值，最后根據(jù)得出的這一組數(shù)值可以繪制出一個(gè)空腔輪廓，該輪廓即為瞬時(shí)空腔最大膨脹波及范圍。

圖1 瞬時(shí)空腔最大膨脹波及范圍的定義

2 基于“或”運(yùn)算的瞬時(shí)空腔最大膨脹波及范圍輪廓提取

文獻(xiàn)[3]提出了一種通過拍攝瞬時(shí)空腔在兩個(gè)正交方向的投影圖像，將空腔截面進(jìn)行橢圓近似對(duì)空腔三維重建，為定量研究武器致傷效應(yīng)提供了依據(jù)。但是，由于瞬時(shí)空腔最大膨脹波及范圍是瞬時(shí)空腔在整個(gè)膨脹過程中所形成的，并不是某個(gè)時(shí)刻的瞬時(shí)空腔。傳統(tǒng)的圖像處理方法需對(duì)序列圖像中的每幀圖像進(jìn)行精確的輪廓邊緣檢測(cè)，比較各個(gè)侵徹距離處的輪廓大小，得到一組最大直徑值。因?yàn)橐Ａ羲休喞獢?shù)據(jù)作比較，對(duì)每幅圖片處理的輪廓精度要求較高。由于受實(shí)驗(yàn)環(huán)境的限制較大，得到的序列圖片大多不夠清晰，噪聲較多，對(duì)每幀空腔圖像進(jìn)行精確的輪廓邊緣檢測(cè)尚無很好的通用算法。文獻(xiàn)[4]探討了創(chuàng)傷彈道空腔X射線投影圖像的邊緣檢測(cè)技術(shù)，常規(guī)的邊緣檢測(cè)算子Sobel算子、Log算子、Canny算子不能有效的檢測(cè)出封閉邊緣。最后提出的基于簡(jiǎn)化M-S模型的改進(jìn)水平集方法，也不適用于文中的高速攝影圖像。

文中提出的算法先對(duì)單幀圖像進(jìn)行Canny算子邊緣檢測(cè)得到空腔輪廓，再對(duì)序列空腔輪廓圖進(jìn)行“或”運(yùn)算疊加，在疊加的過程中進(jìn)行適當(dāng)?shù)某胩幚恚コ嘤嗟膫芜吘?，盡可能的提高輪廓邊緣檢測(cè)的精度，最后通過濾波去噪和連接算法得到所需的輪廓，為空腔容積計(jì)算提供數(shù)據(jù)基礎(chǔ)。

3 單幀圖像瞬時(shí)空腔輪廓提取

3.1 Canny算子邊緣檢測(cè)

明膠經(jīng)脫色處理后灰度均勻，瞬時(shí)空腔膨脹時(shí)輪廓邊緣與明膠有較大的灰度差異。文中用信噪比大、檢測(cè)精度高的Canny[5]算子檢測(cè)邊緣。由于Canny算子首先對(duì)圖像進(jìn)行高斯平滑，對(duì)噪聲有一定的抑制能力。同時(shí)使用高閾值判斷邊緣點(diǎn)，用低閾值連接邊緣點(diǎn)，避免了大量虛假邊緣的產(chǎn)生，可以比較精確的檢測(cè)出邊緣。

圖2 高速攝像拍攝的原始圖像

圖3 Canny算子邊緣檢測(cè)

殺傷元在明膠中飛行時(shí)，所形成的空腔內(nèi)部也產(chǎn)生灰度變化，因此使用Canny算子檢測(cè)出來的邊緣除了空腔的輪廓邊緣還含有很多非感興趣邊緣，如在瞬時(shí)空腔內(nèi)部檢測(cè)到的邊緣，如圖3。而瞬時(shí)空腔的膨脹范圍只與其外輪廓有關(guān)。為了減少偽邊緣帶來的影響，文中在用Canny算子檢測(cè)后的邊緣圖像基礎(chǔ)上提出空腔外輪廓提取算法。

3.2 空腔外輪廓提取

將一幅經(jīng)Canny算子檢測(cè)后的邊緣圖像表示成大小為m×n的矩陣A:

1)產(chǎn)生一幅與邊緣圖像大小相等的臨時(shí)圖像Btemp，并將所有元素值置為0。

2)對(duì)第j列元素，從第1行向第m行檢測(cè)，記錄檢測(cè)到的第一個(gè)值為1的元素，記下該元素的位置(iup，jup)。從第m行向第1行檢測(cè)，記錄檢測(cè)到的第一個(gè)值為1的元素，記下該元素的位置(idown，jdown)。

3)在臨時(shí)矩陣Btemp中將位置為(iup，jup)，(idown，jdown)的元素置為1。

4)遍歷每列元素，重復(fù)2)、3)過程，最終得到的矩陣Btemp即為空腔輪廓外邊緣圖像。

圖4 單幀圖像輪廓提取

經(jīng)外輪廓提取后，得到的單幀瞬時(shí)空腔外輪廓如圖4所示。

4 空腔輪廓圖像“或”運(yùn)算

經(jīng)輪廓提取后，得到了單幀圖像的瞬時(shí)空腔輪廓圖。將處理得到的序列空腔輪廓圖像進(jìn)行“或”運(yùn)算疊加，多幀圖像的瞬時(shí)空腔輪廓由于存在重疊、相鄰，表現(xiàn)為輪廓變粗。

在殺傷元侵徹明膠時(shí)，由于明膠抖動(dòng)或其他外界條件的影響，有時(shí)檢測(cè)出來的單幀瞬時(shí)空腔輪廓并不精確，如圖5(c)空腔下邊緣輪廓并不能很好的檢測(cè)出來，同時(shí)不可避免的還存在一些偽邊緣。但是由于處理的不是單幀圖像，而是連續(xù)多幀圖像的“或”運(yùn)算疊加，少數(shù)不精確的空腔輪廓可以由多數(shù)相對(duì)精確的輪廓彌補(bǔ)，隨著序列圖像幀數(shù)的增多，空腔輪廓繼續(xù)變粗，表現(xiàn)得更突出。偽邊緣由于隨機(jī)性，在對(duì)瞬時(shí)空腔輪廓圖進(jìn)行“或”運(yùn)算時(shí)，偽邊緣重疊、相鄰的概率小，呈現(xiàn)為一些細(xì)小的邊緣，為了防止偽邊緣被錯(cuò)誤的加強(qiáng)，可設(shè)置一定的除噪周期，選取合適的結(jié)構(gòu)元素對(duì)“或”運(yùn)算疊加后的圖像進(jìn)行形態(tài)學(xué)開運(yùn)算，消除不重要的圖像細(xì)節(jié)，同時(shí)保證圖像的全局幾何形狀不失真。

圖6 連續(xù)30幀圖像疊加后進(jìn)行輪廓提取

5 瞬時(shí)空腔最大膨脹波及范圍輪廓曲線濾波去噪

“或”運(yùn)算疊加圖經(jīng)過開運(yùn)算后偽邊緣得到了抑制，再對(duì)其進(jìn)行空腔輪廓提取，由空腔外輪廓提取算法可知，所得的輪廓曲線即為瞬時(shí)空腔最大膨脹波及范圍輪廓。由于“開”操作以及Canny算子本身的局限，最終得到的空腔輪廓不連續(xù)且不可避免存在噪聲點(diǎn)如圖6(c)。提取得到的空腔輪廓由上、下兩條輪廓邊緣組成，每條輪廓上都存在相應(yīng)的噪聲點(diǎn)，噪聲點(diǎn)表現(xiàn)為偏離空腔輪廓，傳統(tǒng)的對(duì)曲線濾波方法有中值濾波、均值濾波、頻域?yàn)V波，針對(duì)各種濾波方法所呈現(xiàn)出來的問題，文中提出一種基于斜率濾點(diǎn)的方法，取輪廓曲線上含有噪聲點(diǎn)的一部分?jǐn)?shù)據(jù)進(jìn)行對(duì)比分析，結(jié)果表明基于斜率濾點(diǎn)的方法最有效。

5.1 中值濾波與均值濾波

通過一個(gè)采樣窗口從原始數(shù)據(jù)中提取奇數(shù)個(gè)數(shù)據(jù) y(i－ k)，y(i－ k+1)，…，y(i)，y(i+1)，…，y(i+k)，中值濾波是對(duì)這些數(shù)據(jù)按大小進(jìn)行排序，用排序后的中間值作為輸出值Y(i)，均值濾波是將提取數(shù)據(jù)的平均值作為輸出值Y(i)。將采樣窗口遍歷所有數(shù)據(jù)。為了使新的數(shù)據(jù)與原來的數(shù)據(jù)個(gè)數(shù)相同，先對(duì)原來數(shù)據(jù)兩端擴(kuò)展k個(gè)數(shù)據(jù)，其值為端部的數(shù)據(jù)。

圖7 中值濾波和均值濾波

中值濾波法可以去除一些噪聲點(diǎn)，特別是當(dāng)采樣窗口中的噪聲點(diǎn)個(gè)數(shù)小于采樣窗口長(zhǎng)度一半時(shí)，濾波效果顯著，但對(duì)于寬度大的噪聲曲線段，需加大窗口，但是，若盲目的加大窗口，對(duì)于無噪聲段曲線，由于窗口的加大，會(huì)使曲線失真。對(duì)于此類無法預(yù)知寬度的噪聲，中值濾波并不能完全適用。同樣均值濾波有一定的濾波能力，使曲線趨于平滑，但是同樣不能去除寬度大的噪聲曲線段。

5.2 頻域?yàn)V波

在信號(hào)處理中，噪聲能量一般集中于高頻，可設(shè)計(jì)低通濾波器，將高頻段的噪聲濾掉，使低頻分量順利通過而有效地阻止高頻分量，通過頻域處理去除噪聲點(diǎn)。

1)空腔輪廓點(diǎn)是一系列的離散點(diǎn)，對(duì)離散點(diǎn)進(jìn)行離散傅里葉變換，得到F(u):

2)建立低通濾波器H(u); (3)

3)G(u)=F(u)H(u); (4)

4)對(duì)G(u)進(jìn)行離散傅里葉反變換;

圖8 頻域?yàn)V波

圖8分別選用巴特沃斯低通濾波器和高斯低通濾波器濾波。由于噪聲并不只是一些階躍的點(diǎn)，有時(shí)也表現(xiàn)為偏移正常曲線的光滑曲線段，光用低通濾波的方法，并不能有效地去除噪聲。

5.3 基于斜率濾波

文中提出基于斜率的濾波方法，分別對(duì)上、下兩條輪廓曲線從曲線上某起始點(diǎn)向左右方向?yàn)V噪聲點(diǎn)，可達(dá)到預(yù)期效果。

具體實(shí)現(xiàn)算法步驟如下:

1)對(duì)每條輪廓線，從連接良好，不存在噪聲處的輪廓線上選取兩相鄰邊緣點(diǎn) Pi、Pj，j=i+1。

2)設(shè)置斜率閾值Δk。

3)計(jì)算Pi與左邊相鄰點(diǎn)Ph的斜率絕對(duì)值Ki，Pj與右邊相鄰點(diǎn)Pk的斜率絕對(duì)值Kj;

4)判斷Ki是否大于 Δk，若是，則 Ph為噪聲，去除。否則Ph為輪廓邊緣點(diǎn)，保留，同時(shí)i=i－1;判斷Kj是否大于 Δk，若是，則 Pk為噪聲，去除。否則 Pk為輪廓邊緣點(diǎn)，保留，同時(shí)j=j+1。

5)重復(fù)3)、4)過程，直到Pj為空腔輪廓最右邊的邊緣點(diǎn)，Pi為空腔輪廓最左邊的邊緣點(diǎn)。

6)通過插值的方法將原來噪聲點(diǎn)的位置補(bǔ)上插值點(diǎn)實(shí)現(xiàn)輪廓邊緣連接。

圖9 基于斜率濾波

該算法不但能有效的去除偏離曲線的單個(gè)噪聲點(diǎn)，而且對(duì)一些偏離曲線的大寬度噪聲曲線段也能較好的去除，同時(shí)對(duì)于準(zhǔn)確的邊緣點(diǎn)不做處理，最大限度的保留輪廓曲線的真實(shí)形狀。

6 實(shí)驗(yàn)結(jié)果與分析

對(duì)瞬時(shí)空腔劇烈膨脹成球形的序列圖像處理如圖10，由于拍攝條件好，空腔輪廓與明膠有較大的灰度差別，用σ=4的高斯濾波器先對(duì)圖像平滑，在前5幀疊加時(shí)，沒有帶入偽邊緣，疊加了10幀后，產(chǎn)生了細(xì)小的偽邊緣，設(shè)定開運(yùn)算周期為10，即每疊加10幀就對(duì)疊加后的邊緣圖進(jìn)行開運(yùn)算，隨著幀數(shù)的增加，空腔輪廓加強(qiáng)越明顯，偽邊緣也不斷地去除，對(duì)90幀空腔輪廓圖像或運(yùn)算結(jié)果如圖10(c)，最終提取的空腔輪廓如圖10(d)，該實(shí)驗(yàn)證明用文中提出的“或”運(yùn)算疊加法能有效獲得瞬時(shí)空腔的最大膨脹范圍輪廓。

圖10 瞬時(shí)空腔膨脹呈球形

圖11 瞬時(shí)空腔膨脹呈錐形

瞬時(shí)空腔膨脹呈錐形的序列圖像處理如圖6。該實(shí)驗(yàn)中各幀圖像包含的噪聲較多，因此在圖6(a)所示的30幀空腔輪廓“或”運(yùn)算疊加結(jié)果中可以清楚地看到若干偽邊緣存在，進(jìn)行周期除噪后還有一些噪聲點(diǎn)殘留，使得圖6(c)所示的最終提取的空腔輪廓并不完整，據(jù)此計(jì)算的空腔容積將產(chǎn)生較大的誤差。用文中提出的基于斜率的曲線濾波算法得到的最終空腔輪廓如圖11所示，對(duì)比可知，最終空腔輪廓得到了很好的改善。

7 結(jié)束語(yǔ)

文中針對(duì)殺傷元造成組織嚴(yán)重?fù)p傷區(qū)域的形成特點(diǎn)，提出了序列瞬時(shí)空腔輪廓“或”運(yùn)算疊加處理方法，使瞬時(shí)空腔最大膨脹波及范圍輪廓的提取重點(diǎn)落在最終疊加結(jié)果圖的處理上，避免了需要精確檢測(cè)單幀圖像輪廓，通過基于斜率的濾波方法可以得到最終較為精確的輪廓，為空腔容積計(jì)算提供了數(shù)據(jù)基礎(chǔ)。若對(duì)實(shí)驗(yàn)進(jìn)行正交拍攝，采用此方法也可對(duì)彈創(chuàng)瞬時(shí)空腔最大膨脹波及范圍進(jìn)行三維重建。文中提出的方法可為槍彈的致傷效果提供一定的評(píng)價(jià)參考，也為各種殺傷武器的研究設(shè)計(jì)提供指導(dǎo)。

[1]安波，蔣浩征，李楊.高速投射物對(duì)生物致傷時(shí)瞬時(shí)空腔的理論研究[J].北京理工大學(xué)學(xué)報(bào):英文版，2001，10(3):272－277.

[2]Fernando Spencer Netto，Dylan Pannell，Homer C Tien.Hollow-point ammunition and handguns:The potential for large temporary cavities[J]. Injury Extra，2008，39(2):50－52.

[3]林勇，王經(jīng)瑾，宋征，等.明膠彈創(chuàng)空腔閃光 X射線投影圖像三維重建[J].清華大學(xué)學(xué)報(bào):自然科學(xué)版，2002，42(12):1576－1578.

[4]賀成，王濤，廖振強(qiáng)，等.創(chuàng)傷彈道空腔圖像邊緣檢測(cè)技術(shù)研究[J].計(jì)算機(jī)工程與設(shè)計(jì)，2011，32(1):248 －250.

[5]Rafael C Gonzalez，Richand E Woods，Steven L Eddins.數(shù)字圖像處理(MATLAB版)[M].阮秋琦，譯.北京:電子工業(yè)出版社，2005.