苗 騰 ， 趙春江， 陸聲鏈， 郭新宇

（1．中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院，北京 100083；2．國家農(nóng)業(yè)信息化工程技術(shù)研究中心，北京 100097）

植物器官的真實(shí)感建模與繪制方法是計算機(jī)圖形學(xué)領(lǐng)域的研究熱點(diǎn)，其中對于葉片的繪制更是重中之重。植物葉片的形態(tài)具有多樣性、不規(guī)則性的特點(diǎn)，如何逼真地模擬其結(jié)構(gòu)特征和形態(tài)細(xì)節(jié)是一項巨大的挑戰(zhàn)。植物葉片三維建模的應(yīng)用需求主要來自兩方面：一是作為影視動畫以及游戲的自然場景元素，二是農(nóng)學(xué)研究中作為虛擬農(nóng)業(yè)試驗(yàn)的虛擬實(shí)驗(yàn)對象。在動畫和游戲制作領(lǐng)域，對于形態(tài)結(jié)構(gòu)細(xì)節(jié)的精細(xì)程度要求不高，更多的是考慮模型的視覺效果以及繪制過程的實(shí)時性。而對于農(nóng)學(xué)研究來說，仿真結(jié)果需要真實(shí)地反映器官的生理特征和品種特性，所以生成的模型要達(dá)到較高的精細(xì)度。

L系統(tǒng)作為通用的植物拓?fù)浣Y(jié)構(gòu)描述語言，在植物模擬中應(yīng)用較廣，很多研究也將其應(yīng)用到葉片的建模中[1-2]，但是對于具有隨機(jī)邊緣的結(jié)構(gòu)，生成的結(jié)果單調(diào)呆板。一些學(xué)者利用邊緣輪廓構(gòu)建葉片的模型[3-4]，這類方法交互性強(qiáng)，可以設(shè)計任意外形的葉片，但是對葉片的自然彎曲等細(xì)節(jié)的表達(dá)存在缺陷。文獻(xiàn)[5]提出了利用真實(shí)三維數(shù)據(jù)重構(gòu)葉片的方法，可以得到自然界中葉片的真實(shí)形態(tài)，但是操作過程十分復(fù)雜。除了對葉片幾何形態(tài)的描述研究之外，一些學(xué)者還從葉片紋理[6-7]以及光照[8-9]等方面對葉片進(jìn)行模擬，在像素空間得到了很好的真實(shí)感。一些學(xué)者也注意到，植物葉片具有豐富的形態(tài)結(jié)構(gòu)，利用上面的方法難以進(jìn)行如卷曲、萎蔫等特定形態(tài)的植物葉片建模，并提出了專門的方法[10]。但這些方法的操作過程都比較復(fù)雜，交互設(shè)計不夠自然。

本文提出了一種新的植物葉片三維形態(tài)交互式設(shè)計方法，利用拉普拉斯變形技術(shù)對葉片曲面進(jìn)行交互式編輯，從而達(dá)到設(shè)計各種形態(tài)的三維植物葉片造型的目的。實(shí)驗(yàn)證明，本方法可以快速構(gòu)建植物葉片的三維模型，同時對真實(shí)葉片的細(xì)節(jié)特征保持較好，避免了失真。

1 葉片曲面輪廓特征點(diǎn)的獲取

本文以黃瓜葉片為例，說明特征點(diǎn)獲取過程。

圖1 利用Fastscan三維掃描儀探筆進(jìn)行數(shù)據(jù)采集

利用三維掃描儀Fastscan獲取葉片的輪廓中軸點(diǎn)數(shù)據(jù)，如圖1所示。用Fastscan的探筆點(diǎn)取需要獲得的輪廓點(diǎn)以及葉脈主脈中軸點(diǎn)（本文采用黃瓜葉片的主葉脈作為中軸線），只獲取少量具有明顯特征的點(diǎn)即可，這樣將大大減少采集數(shù)據(jù)時的工作量。在點(diǎn)取過程中，需要分別記錄下左、右輪廓以及中軸上選取點(diǎn)的數(shù)量，以便系統(tǒng)可以對數(shù)據(jù)進(jìn)行分類。如有需要，可以在Fastscan配套的軟件Polhemus里查看、編輯點(diǎn)信息。

除了利用掃描儀獲得輪廓中軸點(diǎn)數(shù)據(jù)外，還可以基于葉脈形態(tài)特征計算這些點(diǎn)信息。葉脈形態(tài)特征是通過葉脈實(shí)測數(shù)據(jù)的統(tǒng)計學(xué)規(guī)律總結(jié)出來的，為此進(jìn)行了一系列的數(shù)據(jù)采集并作統(tǒng)計學(xué)分析。選生長均勻一致的 15株進(jìn)行取樣，在每株的同一葉位上選取葉片進(jìn)行數(shù)據(jù)采集，這樣選取的葉片生理年齡基本相同。將葉片盡量展平直至可以近似為平面，測量特征較明顯的一級葉脈數(shù)據(jù)（長度以及角度）。如圖 2，測量 OA、OB、OC、OE、AF、CG、DI、BH的長度，以及∠COA、∠EOC、∠DOE、∠BOD、∠OAF、∠OCG、∠ODI、∠OBH的角度值。為方便分析，對葉脈進(jìn)行編號，OE為1號葉脈，OD +DI為2號，OB+BH為3號，OC+CG為4號，OA+AF為5號。

圖2 測量葉脈數(shù)據(jù)

得到數(shù)據(jù)之后，對其進(jìn)行處理，首先處理長度數(shù)據(jù)。設(shè)OE為1，對其他葉脈長度數(shù)據(jù)做歸一化，得到各個葉脈長度之間的比例系數(shù)，之后計算比例系數(shù)的相關(guān)性，發(fā)現(xiàn)數(shù)據(jù)間的相關(guān)性較高，所以將得到的比例系數(shù)近似作為葉脈長度之間的比例因子，這樣在系統(tǒng)中只需要輸入葉脈OE的長度，就可以近似得出其他的葉脈長。相關(guān)性數(shù)據(jù)如表1所示。之后處理角度數(shù)據(jù)，對角度值按區(qū)間進(jìn)行劃分，選取出現(xiàn)概率最大的角度范圍作為各個角度的取值范圍。

表1 各編號葉脈長度數(shù)據(jù)相關(guān)性

計算出葉脈長度以及角度值之后，可以根據(jù)二維平面幾何簡單地計算出葉片關(guān)鍵特征點(diǎn)E、G、F、H、I的位置，之后再根據(jù)各個點(diǎn)之間的邊緣形態(tài)特征計算出輪廓點(diǎn)，如EG之間用直線擬合，F(xiàn)O之間用二次曲線擬合。為了避免呆板，在擬合中會加入適當(dāng)?shù)碾S機(jī)處理。而中軸上的點(diǎn)只需要在OE線段上選取即可。

本方法支持手動添加并調(diào)整輪廓中軸點(diǎn)的功能。用戶可以在操作界面上添加點(diǎn)信息，也可以隨意拖動點(diǎn)，這樣可以較靈活的修改曲面信息。圖3為交互編輯數(shù)據(jù)點(diǎn)的操作過程。

圖3 交互編輯數(shù)據(jù)點(diǎn)

2 葉片曲面網(wǎng)格生成

以輪廓點(diǎn)以及中軸點(diǎn)為控制點(diǎn)，利用三次樣條插值曲線分別擬合出左、右輪廓線以及中軸線，之后按照一定的間隔在輪廓以及中軸線中獲得一系列的采樣點(diǎn)；將左右輪廓采樣點(diǎn)分別與中軸采樣點(diǎn)順次相連形成多邊形，然后根據(jù)DeWall算法[11]實(shí)現(xiàn)三維空間中不規(guī)則多邊形的三角化。為了使曲面更加光滑，可以對曲面進(jìn)行細(xì)分，本文采用自適應(yīng)細(xì)分方法[12]進(jìn)行細(xì)分。

獲得三角網(wǎng)格之后，計算每個點(diǎn)的紋理坐標(biāo)。本文簡單的將葉片曲面投影到一個二維平面，然后與葉片紋理做映射。該二維平面通過主成分分析方法計算，計算步驟如下：

1）計算所有頂點(diǎn)的重心坐標(biāo)；

2）構(gòu)建離散度矩陣，并通過離散度矩陣計算協(xié)方差矩陣；

3）對協(xié)方差矩陣進(jìn)行特征分解；

4）其中最大特征值對應(yīng)的特征向量為第1主成分，第2大特征值對應(yīng)的特征向量為第2主成分；

5）第1主成分向量和第2主成分向量組成的平面確定為該二維平面。

真實(shí)世界中，黃瓜葉片的正反面是不同的，為了模擬這種情況，本文渲染時繪制兩遍曲面，第1次繪制時開啟背面剔除，并映射葉片正面紋理；第2次繪制時開啟正面剔除，映射葉片背面紋理。這樣會大大增加葉片的真實(shí)感。

3 葉片曲面的變形計算

本文采用拉普拉斯變形方法對葉片曲面進(jìn)行變形。拉普拉斯方法速度較快，能夠達(dá)到實(shí)時交互，而且其交互性比較簡單，用戶可以在三角網(wǎng)格上選取一系列操作點(diǎn)，然后拖動任意一個操作點(diǎn)，之后整個曲面中非操作點(diǎn)將產(chǎn)生變形，而操作點(diǎn)盡量保持不變。同時，拉普拉斯方法在保持網(wǎng)格局部細(xì)節(jié)方面有天然優(yōu)勢。

3.1 網(wǎng)格處理

由于本文方法得到的網(wǎng)格會出現(xiàn)有兩個邊都在葉片邊緣的三角形，首先要對這種三角片進(jìn)行處理。處理流程如下：

1）建立每個頂點(diǎn)的信息索引，包括頂點(diǎn)與相鄰點(diǎn)的索引關(guān)系，以及與所在面的索引關(guān)系；

2）遍歷頂點(diǎn)信息，如果頂點(diǎn)只在一個三角形中，則該三角形為待處理的三角形，該頂點(diǎn)為孤立點(diǎn)；

3）分割待處理的三角形；

4）更新頂點(diǎn)信息索引。

分割待處理的三角形的方法如圖4所示。其中圖4(a)為要處理的三角形在網(wǎng)格中的拓?fù)涫疽鈭D，虛線為葉片的輪廓邊。對圖4(a)情況進(jìn)行處理，如圖4(b)所示，取BC中點(diǎn)E，然后連接AE、DE。

圖4 網(wǎng)格處理

3.2 拉普拉斯坐標(biāo)

設(shè)v1, v2,…vn為網(wǎng)格上的頂點(diǎn)位置，i'表示與vi相鄰的頂點(diǎn)索引集，則頂點(diǎn)vi的拉普拉斯坐標(biāo)用如下式計算

其中，ωij表示vi與vj所在邊相對于vi的權(quán)重，本文取ωij= c otαij+ c ot βij。這樣表示的拉普拉斯坐標(biāo)可以看作為vi點(diǎn)的平均法向曲率的近似值。任意點(diǎn)的拉普拉斯坐標(biāo)可以看作是其相鄰點(diǎn)集合的權(quán)重平均，所以它在一定程度上表達(dá)了該點(diǎn)的局部幾何信息。拉普拉斯坐標(biāo)的計算公式表達(dá)成矩陣形式如下

其中， L ( x, y, z)為 n×3階拉普拉斯坐標(biāo)矩陣，V (x, y, z)為n×3階頂點(diǎn)坐標(biāo)矩陣，L為n×n階拉普拉斯系數(shù)矩陣，L中的對角線元素為非對角線元素為-ωij(i，j=1,2…n)，當(dāng) vi、vj不相鄰時， -ωij為 0。L 矩陣的秩為n-1，且在變形中該矩陣保持不變。

3.3 拉普拉斯變形系統(tǒng)

構(gòu)造一個誤差函數(shù)，用來表達(dá)變形前后幾何信息的差異，誤差函數(shù)可以用下式表示

ui為操作點(diǎn)的坐標(biāo)，Vhandle 表示操作點(diǎn)的集合，ki為調(diào)控參數(shù)。我們的目的是找到最合適的點(diǎn)，使得上式最小。解這一問題最后可以歸結(jié)為解超定稀疏線性方程組的問題

V(x, y, z)為待求的未知頂點(diǎn)坐標(biāo)矩陣，H為操作點(diǎn)系數(shù)矩陣，該矩陣元素hii=1，當(dāng)且僅當(dāng)vi點(diǎn)為操作點(diǎn)，其余元素都為0。h(x, y, z)為操作點(diǎn)變形之后的坐標(biāo)與ki的乘積矩陣。每次操作控制點(diǎn)進(jìn)行變形時，都更新b，之后解線性方程組得到變形之后的點(diǎn)坐標(biāo)。本文利用稀疏矩陣的Cholesky分解法解該方程組。

4 實(shí)驗(yàn)結(jié)果與討論

上述系統(tǒng)在VC2005和OpenGL圖形引擎編程環(huán)境下實(shí)現(xiàn)，并在配置為1.96G HZ的CPU ,DDR 2G內(nèi)存以及ATI RHD2400顯卡的PC機(jī)上進(jìn)行實(shí)驗(yàn)。圖5為一個黃瓜葉片曲面網(wǎng)格形成過程：圖5(a)為初始輸入的輪廓以及中軸點(diǎn)（由掃描儀獲得），圖 5(b)為由樣條插值生成的輪廓中軸曲線，圖5(c)為生成的初始網(wǎng)格，圖5(d)在圖5(c)的基礎(chǔ)上進(jìn)行了一次細(xì)分，圖 5(e)消除了有兩條邊在輪廓上的三角形，紅色線框指出了經(jīng)過處理的網(wǎng)格。圖6為葉片的變形過程：圖6(a)為未變形時的葉片曲面，圖6(b)為用鼠標(biāo)選取操作點(diǎn)的過程，紅色的點(diǎn)為操作點(diǎn)，圖6(c)、(d)、(e)為經(jīng)過變形之后的葉片曲面，其中由藍(lán)色圓框圈中的點(diǎn)為在變形過程中拖動的操作點(diǎn)。該葉片曲面網(wǎng)格包含 2600個點(diǎn)，每次變形計算用時 219毫秒，能夠滿足實(shí)時性交互速度。

利用本文方法可以較容易地生成葉片在自然情況下發(fā)生的變形現(xiàn)象。圖7為利用本文方法對黃瓜葉片萎蔫的模擬：其中圖7(a)為現(xiàn)實(shí)中黃瓜葉片萎蔫的形態(tài)，圖7(b)為生成的黃瓜葉片曲面，圖7(c)和圖7(d)為對圖7(b)進(jìn)行變形操作后得到的葉片萎蔫模擬效果。圖8為對西瓜葉片自然狀態(tài)下的模擬：圖8(b)為西瓜葉片曲面，圖8(c)和圖8(d)為變形之后得到的西瓜葉片在自然狀態(tài)下的形態(tài)模擬。

葉面積是農(nóng)業(yè)生產(chǎn)上的重要指標(biāo)，利用葉面積可以計算作物冠層光分布，從而估算植物的光合作用強(qiáng)度[13]。本文對系統(tǒng)生成的葉片面積與實(shí)際的葉片面積進(jìn)行了比較。分析利用實(shí)測數(shù)據(jù)作為輪廓中軸輸入生成的葉面積誤差，其中輸入的左右邊緣點(diǎn)各為22個，中軸點(diǎn)為5個。實(shí)際葉片的面積利用三維掃描儀得到的結(jié)果來計算，而對于系統(tǒng)生成曲面的面積通過累加所有三角面片面積來獲得，且每個三角形的面積用海倫公式計算。計算結(jié)果如表2所示，我們?nèi)〕稣`差最大的5組數(shù)據(jù)，可以看出誤差都限制在8%之內(nèi)。對于根據(jù)葉脈形態(tài)方法生成的葉片面積，也用上述方法計算誤差。取生長期的葉片作為測試，得到的誤差在14%左右。

用實(shí)測數(shù)據(jù)作為系統(tǒng)輸入得出的葉片在葉面積上與實(shí)際相差較小，可以為光分布計算打下良好基礎(chǔ)，但是需要通過對實(shí)際葉片進(jìn)行掃描，前期工作量相對第二種方法要大；對于利用葉脈形態(tài)作為系統(tǒng)輸入的方法，其中的參數(shù)是根據(jù)一定量的數(shù)據(jù)進(jìn)行統(tǒng)計學(xué)分析得到的，可以反映各個時期葉片的總體形態(tài)特征，但是會有一定誤差，誤差基本在14%以內(nèi)。但是利用該方法可以只通過輸入主脈長度就可以迅速構(gòu)建出葉片曲面，而不需要對每個葉片進(jìn)行三維數(shù)字化，如果需要迅速構(gòu)建葉片形態(tài)，利用這種方法最方便。

圖5 葉片曲面生成過程

圖6 葉片變形過程

圖7 黃瓜葉片萎蔫模擬

圖8 自然狀態(tài)的西瓜葉片模擬

表2 生成葉片曲面面積與真實(shí)葉片面積的比較結(jié)果

5 結(jié) 論

本文介紹了一種新的三維植物葉片設(shè)計方法，通過拉普拉斯網(wǎng)格變形技術(shù)對葉片的曲面網(wǎng)格進(jìn)行交互式編輯，利用本方法生成的葉片不僅在形態(tài)上具有較強(qiáng)的真實(shí)感，同時能夠很好地保持葉片的面積特征。需要說明的是，本文的方法不一定要以葉片輪廓及中軸點(diǎn)數(shù)據(jù)作為輸入生成網(wǎng)格曲面，也可以直接在已有的網(wǎng)格曲面上進(jìn)行變形操作，能夠應(yīng)用于所有的植物葉片，因此具有較好的普適性。

本文的系統(tǒng)在渲染流程中只應(yīng)用了簡單的紋理以及基本的光照，并沒有考慮散射光以及葉片材質(zhì)特性對渲染結(jié)果的影響。在今后的工作中，將加入對該部分的模擬，以求得到更真實(shí)的效果。

[1]Rozenberg G ,Lindenmayer A. Developmental systems with locally catenative formulas [J]. Acta Informatica,1973, 2(3): 214-248.

[2]Rodkaew Y, Lursinsap C, Fujimoto, et al. Modeling leaf shapes using L-systems and genetic algorithms [C]//Proceedings of NICOGRAPH2002, Tokyo, Japan,2002: 73-78.

[3]Mündermann L, MacMurchy P, Pivovarov J, et al.Modeling lobed leaves[C]//Proceedings of Computer Graphics International, Tokyo, Japan, 2003: 60-65.

[4]Lu S L, Guo X Y, Zhao C J, et al. Modeling of organs and structures generating for crops [C]// International Symposium on Crop Modeling and Decision Support,Nanjing, China, April, 2008: 669-673.

[5]Loch B, Belward J, Hanan J. Application of surface fitting techniques for the representation of leaf surfaces [C]//Zerger A, Argent R M(eds). International Congress on Modelling and Simulation, Melbourne: MODSIM Press, 2005: 1272-1278.

[6]Rodkaew Y, Siripant S, Lursinsap C, et al. An algorithm for generating vein images for realistic modeling of a leaf [C]//Proceedings of Computational Mathematics and Modeling (CMM2002), Bangkok,Thailand, 2002: 223-228 .

[7]Runions A, Fuhrer M, Lane B, et al. Modeling and visualization of leaf venation patterns [J]. ACM Transactions on Graphics, 2005, 24(3): 702-711.

[8]Wang L F, Wang W L, Dorsey J, et al. Real-time rendering of plant leaves [J]. ACM Transactions on Graphics, 2005, 24(3): 712-719.

[9]Baranoski G V G, Rokne J G. Efficiently simulating scattering of light by leaves [J]. The Visual Computer,2001, 17(8): 491-505.

[10]Sung M H, Simpson B, Gladimir V G B. Interactive venation-based leaf shape modeling [J]. Computer Animation and Virtual Worlds, 2005, 16(3-4):415-427.

[11]Cignoni P, Montani C, Scopigno R, et al. A fast divide and conquer delaunay triangulation algorithm in ed [J]. CAD, 1998, 30(5): 333-341.

[12]陸聲鏈, 郭新宇, 李長鋒. 三維植物葉片精確建模和繪制技術(shù)研究[J]. 中國圖象圖形學(xué)報, 2009,14(4): 731-737.

[13]溫維亮, 孟 軍, 郭新宇, 等. 基于輻射照度的作物冠層光分布計算系統(tǒng)設(shè)計[J]. 農(nóng)業(yè)機(jī)械學(xué)報,2009, 40: 190-193.