一種簡單可靠的昆蟲撲翼運動圖像序列自動分析方法

魏 榛,楊基明,C.P.Ellington

(1.中國科學技術大學,安徽合肥 230027;2.劍橋大學,英國 CB2 3EJ)

0 引 言

昆蟲飛行的運動學描述一直是撲翼飛行空氣動力學研究所需的重要信息。高速攝影則是獲得昆蟲撲翼運動細節的絕佳工具。如今,昆蟲觀測中高速相機拍攝速度已達到1000fps以上,一次可存儲上千幀昆蟲運動的數字圖像。此時采用人工逐幀提取圖像中昆蟲形態學特征的方法[1]來分析數據顯然不可取,開發一種能夠自動分析海量昆蟲撲翼運動圖像的方法十分必要。

然而,昆蟲運動圖像的復雜性和經典運動圖像識別技術的局限性,使得實現這樣的自動分析非常不容易。首先,昆蟲圖像輪廓十分復雜：真實昆蟲翼面不僅具有揮拍角、迎角和偏斜角的變化還具有扭轉和彎曲變形[2-3],而且昆蟲身體還有其自身的運動姿態。其次,撲翼運動的某些時刻翼面圖像和昆蟲身體圖像間會發生部分遮擋。再次,昆蟲翼面本身的透明性使得翼面圖像和背景十分接近。再則,每幀圖像都需要同時提取多個形態學特征：左右翼面各自的翼尖、翼根點,前、后緣輪廓,身體的頭部和腹部端點等。最后,經典的運動圖像識別方法直接處理昆蟲圖像問題大都有局限：如基于背景提取的模板匹配方法[4-5]雖然可以適應遮擋、透明和噪聲等問題,但較難適應昆蟲圖像輪廓的柔性大變形;而光流法[7]可以克服柔性大變形影響,卻無法較好處理遮擋、透明和噪聲等問題。

為了構建一種切實可用的,無需人工干預,不針對特定昆蟲,運算量也較小的昆蟲撲翼圖像自動分析方法,經過對昆蟲撲翼運動圖像的內在特點的分析,獨辟蹊徑地提出了一種完全不同于經典運動圖像識別方法的新方法,并嘗試了使用程序合成的虛擬昆蟲運動圖像和真實昆蟲的運動圖像測試該方法。

1 基本假設和圖像處理流程

昆蟲圖像的形態學特征提取是獲得昆蟲飛行的運動學描述十分基本而關鍵的步驟：單相機觀測系統得到昆蟲形態學特征后即可使用對稱假設[1,8]等求出運動學描述;多相機系統得到每個相機視角的昆蟲形態學特征后,即可通過對應特征的空間匹配[5-6,9-10]求出運動學描述。因此,這里不討論具體的觀測系統,而是解決單個視角昆蟲圖像上的自動昆蟲形態學特征提取這一基本而關鍵的難題。為了實現所提出的處理方法,對于選取的昆蟲圖像需要滿足一些簡單的假設。

首先,圖像序列中昆蟲的身體和翼面間的遮擋不能太嚴重。嚴重遮擋的時候,不便于同時提取兩翼面的形態學特征。其次,相鄰兩幀圖像間的翼面位置變化要大于昆蟲身體的位置變化。再次,圖像中翼面前緣輪廓線與翼面后緣輪廓線只在翼尖點和翼根點相交,即不考慮翼面嚴重扭轉的情形。最后,假設昆蟲身體部分只發生剛性運動,且昆蟲頭部端點到翼根點的距離要小于昆蟲腹部端點到翼根點的距離。

選取合適的相機角度和拍攝速度,上述假設對于果蠅,黃蜂,鷹蛾等昆蟲都較容易滿足。得到圖像序列之后,首先需要獲得每幀圖像中昆蟲部分的二值圖,然后再根據得到的二值圖計算昆蟲的形態學特征。與經典運動圖像處理方法截然不同：這里采用了相鄰三幀圖像的巧妙運算來實現昆蟲翼面和身體的二值圖提取,借助了昆蟲內在特征和上述假設,實現從二值圖逐一自動分析所需的多個昆蟲形態學特征。

2 昆蟲翼面和身體二值圖的提取

為得到昆蟲翼面和身體的二值圖,利用了相鄰三幀圖像的差別。此方法既克服了經典運動圖像識別中直接二值圖提取時無法應對昆蟲翼面圖像與背景反差太小的問題,又克服了背景靜止圖像對二值圖提取時的干擾。其獨到過程如下圖1所示。

圖1 翼面二值圖提取原理Fig.1 Principle of wing binary image extraction

翼面二值圖提取的操作原理如圖 1所示：圖1(a)表示從i-1到i+1時刻連續的三幀昆蟲翼面輪廓關系;圖 1(b)中,把第i-1與i幀的灰度圖像相減并二值化得到的圖像標號記作i-(i-1),其中白色部分為相鄰兩幀間翼面因重疊被減去的部分。然后使用圖像形態學算法可以填充這個白色區域部分,得到的圖像標號記作i-1/2;標號i+1/2的圖像提取的方法完全類似。最后,對標號i-1/2和i+1/2的翼面二值圖進行“與”操作,即得到i時刻昆蟲翼面部分的二值圖。

此方法巧妙地使用了兩次相鄰幀的相減圖像的二值化操作來提取中間幀的二值圖。因為是從相減圖像中恢復原圖像信息,所以不需要使用背景圖像的灰度估計等,也就避免了翼面與背景反差太小時經典運動圖像處理方法的直接二值化困難。同時,經典的圖像處理方法只使用一次相鄰幀的相減圖像的二值化操作來估計運動目標的大體位置,這實質上并沒有充分地利用幀間相減圖像的二值圖所包含的全部信息。使用兩次相鄰幀的相減圖像的二值化操作,則充分利用了相減圖像的二值圖信息,可以直接還原出中間幀的二值圖,因而與圖像中翼面如何變形無關,這也避免了經典圖像處理方法中模板匹配時不適合柔性大變形的難處。當然,圖1是在理想情況下的操作過程,實際使用中還必須考慮真實圖像二值化的方法和圖像形態學處理的技巧。

相鄰兩幀圖像相減時,所得灰度圖像反映了相減圖像間的灰度變化;對于昆蟲運動圖像而言,相鄰兩幀間的最大差別來自于翼面運動所帶來的灰度變化。由于相鄰兩幀圖像中翼面變化并不大,翼面灰度與背景灰度差別也不大,所以相減所得的圖像的灰度動態范圍較小,因此先要增強相減圖像上翼面運動區域所代表的灰度特征。設相減圖像中一共有k階灰度,第i階灰度值為gi,具gi灰度值的像素個數為Ni,增強之后的加權灰度值Gi由下式計算

公式(1)所構造的灰度值權重函數W是一個與像素個數比的對數成反比的函數,即相減圖像上具有某個灰度值的像素個數占圖像總像素個數的比例越高,就越不可能是翼面運動有關的部分所在的像素,而相減圖像上那些灰度值較大且像素個數較少的像素則最有可能反映翼面的運動區域,在公式(1)中也取得了較大的灰度權重而被圖像增強。顯然這與相減圖像的物理直觀一致,比經典圖像處理方法中的簡單二值化方法更能強化圖像中所包含的運動和輪廓信息。

增強處理之后,相減圖像能夠充分地反映翼面運動區域的灰度特征,再取閾值進行相減圖像的二值化,就可得到較完整的相減圖像的二值圖,即圖1中的標號為i-(i-1)的圖像或標號(i+1)-i的圖像。

以標號i-(i-1)的圖像為例,分析得到標號i-1/2的圖像的形態學處理技巧。真實圖像由于翅脈的灰度反差與翼面部分十分不一致或相減圖像的二值化閾值選取偏差,所得到的標號i-(i-1)的圖像一般會出現連通區域破碎。這時圖像形態學處理進行適當修補,處理過程如圖2所示。

圖2 中間過程二值圖的形態學操作Fig.2 Morphological operation of the binary image in process

首先通過連通域分割得到二值圖的分塊結果如圖2按面積由大到小編號為1,2,3…;然后,借助輪廓操作對每個分塊區域求取與編號1區域的外切線連通域;最后填充生成的連通域,再“并”操作即得到標號i-1/2的圖像。經典的運動圖像處理方法中,對于標號i-(i-1)的圖像一般不會采取如圖2的深入處理,因為經典的運動圖像處理方法并非基于相鄰三幀圖像的兩次圖像相減的,而在獨辟蹊徑的巧妙方法中,圖2的圖像形態學操作使得最終標號i-1/2的圖像與標號i+1/2的圖像的“與”操作成為可能。

經過“與”操作之后,第i幀圖像中的昆蟲翼面二值圖已經成功提取出來。不難發現,本方法巧妙地回避了正面解決昆蟲圖像輪廓柔性大變形和昆蟲圖像質量不佳等問題,獨辟蹊徑把輪廓提取問題簡化成兩次相鄰幀的相減圖像的二值化操作為核心的簡單處理過程,因而能夠處理經典運動圖像處理方法所不適合的難題。

由于昆蟲身體不透明,因而與圖像背景反差鮮明;所以得到第i幀的昆蟲翼面二值圖之后,直接在第i幀原圖像中扣除翼面二值圖部分,再對原圖像中翼面二值圖所在的局部區域進行二值化即可得到昆蟲身體部分的二值圖。第i幀的昆蟲翼面二值圖已經獲得,所以不僅昆蟲所在區域已經確定,而且扣掉翼面原圖部分之后,昆蟲身體的二值圖提取就不再有來自翼面灰度的干擾,所以難度大大降低。

3 昆蟲形態學特征提取

得到各幀的昆蟲翼面和身體二值圖后,即可提取昆蟲形態學特征。經典圖像處理方法中,基于背景提取的模板匹配方法如果已經完成匹配,那么就可以從用于匹配的參照模板圖像上直接獲得昆蟲的形態學特征。然而模板匹配方法首先面對此難題就很難實現圖像匹配;就算得到匹配圖像,也是基于特定昆蟲的匹配,無法適應非特定昆蟲的識別要求。而另一類經典運動圖像處理方法如光流法,且不計較難以處理具有遮擋、透明和噪聲等問題的昆蟲圖像,即使得到輪廓結果,仍然也需要設計方法來進行昆蟲形態學特征的分析,并無捷徑可走,因此與本方法相比并無優勢。根據所述方法已經分離出第i幀的翼面二值圖和昆蟲身體二值圖,因此,昆蟲形態學特征的提取成為可能。昆蟲形態學特征提取過程如圖3所示。

圖3 昆蟲形態學特征提取步驟Fig.3 Step of insect morphology feature extraction

由于使用了常見昆蟲的內在特征作為形態學特征提取的依據,因此對于滿足這些特征的昆蟲,這種形態學提取方法都是可靠的。也正因為此,所以所提供的形態學特征提取方法是不針對特定昆蟲的,也是無需人工干預的,這大大優于經典的匹配的方法。

4 測試結果

采用了由程序合成的虛擬昆蟲運動圖像和真實的昆蟲運動圖像來完成測試上述的方法。程序合成的虛擬昆蟲運動圖像用于驗證結果的正確性,真實的昆蟲運動圖像用于驗證該方法的實用性。

程序合成的虛擬昆蟲圖像使用了MATLAB生成,基于理想相機模型的投影關系,模擬了俯視相機角度的懸停飛行的鷹蛾多個揮拍周期的運動。虛擬圖像如圖4所示。

由于圖像序列為程序合成,所以預先知道虛擬昆蟲的形態學特征,因此能夠對比本方法自動處理得到的形態學特征。形態學特征點的比較結果以X,Y方向的誤差像素大小來衡量,形態學特征輪廓的比較結果以輪廓的平均誤差像素來衡量。此處合成了600×800分辨率,90幀,共3個撲翼周期的灰度圖像用于測試,最后的誤差結果如圖5所示。

圖4 程序合成的虛擬昆蟲運動圖像截圖Fig.4 Motion video clip of the virtual insect generated by program

圖5 虛擬昆蟲處理結果的誤差分布Fig.5 Error distribution of virtual insect result

圖5只列出左翼的誤差結果,右翼結果類似。虛擬圖像中昆蟲身體長度約為174個像素,揮拍幅度150°左右,上圖可見,虛擬昆蟲頭部和腹部端點識別結果平均誤差在0.6%的昆蟲身體像素左右;翼尖識別結果除去少數幾幀因為翼尖位置太靠近身體部分而在X方向上識別誤差較大外,平均誤差在4.2%的昆蟲身體像素左右;翼根點識別使用了擬合方法,因此誤差呈現較規律的分布,平均誤差約在2.9%的昆蟲身體像素左右;翼面前緣和后緣輪廓平均誤差都控制在較小范圍,大約在0.4%的昆蟲身體像素左右。綜上可知,在虛擬昆蟲圖像上本方法能夠給出較好的識別結果。因為虛擬昆蟲圖像不能很好地模擬真實高速攝影圖像效果中的背景噪聲,圖像輪廓模糊等情況,因此需要真實昆蟲圖像作為進一步的測試。

測試所使用的真實昆蟲運動圖像是在進行偏航角調整的機動飛行過程中的長尾管蚜蠅(Eristalis tenax)俯視圖像,一共150幀,600×800分辨率,拍攝速度為4700f/s。圖6為真實昆蟲運動圖像部分幀的截圖。

使用本處理方法進行自動分析,即可得到昆蟲的形態學特征。處理程序運行于MAT LAB V7.5,一共處理得到148幀圖像中的運動特征,平均每幀畫面運算耗時8.45s,圖7列出了圖6中的圖像幀的處理結果。

從圖6的圖像可見,昆蟲左右翼的運動非常的不對稱;除前緣之外的昆蟲翼面灰度都十分接近背景灰度;幀編號fn=100的圖中右翼面較為靠近昆蟲身體。而圖7所得提取結果在這些情況下都能較為準確地判別昆蟲的特征點和翼面輪廓。由此可知,獨辟蹊徑所采用的二值圖和形態學特征提取方法對于真實昆蟲圖像能夠自動分析給出較好的結果,這已經克服了經典運動圖像處理方法在昆蟲運動圖像處理上的難度,具有實用價值。由各幀所得的形態學特征數據進一步可以繪制出圖像坐標系下昆蟲身體頭尾連線角度和兩翼面與身體夾角的變化曲線,如圖8所示。

圖6 真實昆蟲運動圖像截圖(fn,fn為幀編號)Fig.6 Motion video clip of real insect

圖7 真實昆蟲圖像的形態學特征提取結果Fig.7 Morphology feature extraction result of real insect

圖8 真實昆蟲形態學特征角度變化Fig.8 Morphological feature angle curve of real insect

圖8中的 α,θL,θR并非嚴格意義下的昆蟲身體空間姿態角和翼面揮拍角,但卻能夠充分反映圖像上昆蟲姿態的變化。圖8可知,一方面在接近30ms的時間內昆蟲翼面完成了5個周期的揮拍運動,且左右翼面極不對稱的運動使得昆蟲身體獲得了約50°的偏航角調整;另一方面圖8曲線是直接使用圖像提取的形態學特征進行繪圖的結果,未經過后期濾波等處理,可知根據形態學特征來計算的角度變化較為平滑,能夠充分地反映圖像中的昆蟲的運動過程,再次證明了所提供的方法能夠自動提取較為準確可靠的形態學特征。

5 結 論

程序合成的虛擬昆蟲運動圖像和真實昆蟲的運動圖像的測試結果表明：該方法簡單可靠,無需人工干預,不針對特定昆蟲,運算量也較小,并能夠提供較好的識別結果。這給昆蟲圖像的形態學特征自動識別難題提供了一個切實可行的解決方法,極大地節省了昆蟲撲翼運動圖像的分析時間。

該方法獨辟蹊徑,采用了相鄰三幀圖像的巧妙運算來實現昆蟲翼面和身體的二值圖提取,借助了昆蟲內在特征和簡單假設,實現從二值圖逐一自動分析所需的多個昆蟲形態學特征,完全不同于模板匹配或光流法等經典的運動圖像分析方法。事實證明該方法能夠較好地適應于昆蟲這種復雜變形體運動,并能夠在翼面透明,昆蟲身體有遮擋,高速攝影本底噪聲大的環境中,順利地自動提取二值圖并分析形態學特征。采用本文所提供的方法可以更加高效便捷地分析昆蟲飛行的運動學特征,這對于昆蟲飛行的空氣動力學研究將很有幫助。

[1] ELLINGTON C P.The aerodynamics of hovering insect flight[J].III.Kinematics Philosophical Transactions of the Royal Society of London Series B-Biological Sciences,1984,(305)：41-46.

[2] CHENG P,HU J,ZHANG G F,et al.Deformation measurements of dragonfly's wings in free flight by using windowed fourier transform[J].Optics and Lasers in Engineering,2008,(46)：157-161.

[3] WALKER S M,WALKER S M,THOMAS L R,et al.Deformable wing kinematics in free-flying hoverflies[J].Journal of the Royal Society Interface,2010,(7)：131-142.

[4] HENIKOFF J,SHAPIRO L G.Representative pat-terns for model-based matching[J].Pattern Recognition,1993,(26)：1087-1098.

[5] FONTAINE E I,ZABALA F,DICKINSON M H,et al.Wing and body motion during flight initiation in Drosophila revealed by automated visual tracking[J].JEB,2009,(212)：1307-1323.

[6] 趙創新,徐進良,張永立,等.基于單攝像機的昆蟲自由飛行參量三維重構[J].光學學報,2006,(26)：61-66.

[7] BARRON J L,FLEE D J,BEAUCHEMIN S S,et al.Performance of optical-flow techniques[J].International Journal of Computer Vision,1994,(12)：43-77.

[8] WILLMOTT A P,ELLINGTON C P.Measuring the angle of attack of beating insect wings：Robust three dimensional reconstruction from two dimensional images[J].JEB,1997,(200)：2693-2704.

[9] WALLACE I D,LAWSON N J,HARVEY A R,et al.High-speed photogrammetry system for measuring the kinematics of insect wings[J].Applied Optics,2006,(45)：4165-4173.

[10]FRY S N,DICKINSON M H.Kinematics and aerodynamics of free flight maneuvers in Drosophila[J].A-merican Zoologist,2001,(41)：1447-1447.