[發(fā)明專利]面向激光點(diǎn)云數(shù)據(jù)的闊葉樹真實(shí)葉片建模與形變方法

有效

• 申請(qǐng)?zhí)枺?/td>: 201410436293.8
• 申請(qǐng)日: 2014-09-25
• 公開（公告）號(hào): CN105654543B
• 公開（公告）日: 2019-01-22
• 發(fā)明（設(shè)計(jì)）人: 薛聯(lián)鳳;云挺;嵇俊
• 申請(qǐng)（專利權(quán)）人: 南京林業(yè)大學(xué)
• 主分類號(hào): G06T17/00
• 分類號(hào): G06T17/00
• 代理公司: 北京匯捷知識(shí)產(chǎn)權(quán)代理事務(wù)所(普通合伙) 11531
• 代理人: 李宏偉

1.一種面向激光點(diǎn)云數(shù)據(jù)的闊葉樹真實(shí)葉片建模與形變方法，其特征在于：首先運(yùn)用高精度的激光掃描儀獲取不同闊葉樹種的三維葉面點(diǎn)云數(shù)據(jù)；其次從掃描獲取的葉子點(diǎn)云中采用多項(xiàng)式擬合的方法獲取精確的樹葉邊界，并結(jié)合計(jì)算機(jī)圖形學(xué)的方法計(jì)算葉面的主葉脈；接著采用雙三次廣義張量積Bezier曲面對(duì)葉面點(diǎn)云數(shù)據(jù)進(jìn)行擬合，進(jìn)而去除由于風(fēng)吹抖動(dòng)造成的掃描誤差并消除噪聲點(diǎn)云點(diǎn)；再次，根據(jù)固體力學(xué)受力模型，定義葉脈與葉肉不同的力學(xué)屬性，并構(gòu)造非線性有限元的受力形變方程，計(jì)算施加不同載荷力后的葉子形態(tài)，進(jìn)而模擬樹葉在自然環(huán)境中真實(shí)的形變；該方法包括如下步驟：

(1)數(shù)據(jù)獲取

(11)掃描獲取整株樹木，提取葉面點(diǎn)云的兩個(gè)端點(diǎn)，p_e＝(x_e，y_e，z_e)^T，p_s＝(x_s，y_s，z_s)^T，分別認(rèn)為p_e是葉尾點(diǎn)，p_s是葉脈的頂端，通過計(jì)算p_e和p_s之間的連線L₁來確定主葉脈上的點(diǎn)，L₁：p＝p_e+t×(p_s-p_e)，t為連線的自變量，同時(shí)確定與矢量相垂直的法向量將主葉脈L₁等間隔n等份后取線上面的n+1個(gè)等分點(diǎn)，這些等分點(diǎn)與法向量構(gòu)成了葉子寬度的n+1條掃描線L_2，i，其中i＝1，2，3...，n+1，定位求取L_2，i左右兩端的端點(diǎn)p_2，li和p_2，ri，即可以獲取葉面在不同部分的邊緣點(diǎn)；

(12)對(duì)步驟(11)的多條L_2，i掃描線，取每條掃描線的左右兩端的端點(diǎn)p_2，li和p_2，ri，分別記為p_2，li：(x_li，y_li，z_li)和p_2，ri：(x_ri，y_ri，z_ri)，其中i＝1，2，3...，n+1；除去兩側(cè)的兩條掃描線，若某條掃描線的左端點(diǎn)相對(duì)于相鄰兩條掃描線的左端點(diǎn)向主葉脈內(nèi)凹，且右端點(diǎn)相對(duì)于相鄰兩條掃描線的右端點(diǎn)也向主葉脈內(nèi)凹，則認(rèn)為該掃描線存在遮擋或由誤差造成，剔除該掃描線；

(13)基于步驟(12)，剔除完所有存在遮擋或由誤差造成的掃描線后，對(duì)所有掃描線按序重新編號(hào)為p_2，lj和p_2，rj，左右兩端的端點(diǎn)分別記為p_2，lj：(x_lj，y_lj，z_lj)和p_2，rj：(x_rj，y_rj，z_rj)；最終得到：

左端點(diǎn)的集合為p_2l＝{(x_l1，y_l1，z_l1)，(x_l2，y_l2，z_l2)，...，(x_lj，y_lj，z_lj)，...}；

右端點(diǎn)的集合為p_2r＝{(x_r1，y_r1，z_r1)，(x_r2，y_r2，z_r2)，...，(x_rj，y_rj，z_rj)，...}左端點(diǎn)的集合和右端點(diǎn)的集合構(gòu)成了葉面的初始掃描輪廓；對(duì)于一片葉子的左/右邊部分輪廓，采用不同投影面擬合與求交的方法定位真實(shí)葉子的邊緣；具體步驟如下：

(131)對(duì)于左/右邊部分輪廓上的邊緣點(diǎn)p_2，j＝(x，y，z)，將(x_i，y_i，z_i)i＝1，2，3...，n作為葉子幾何邊緣的參數(shù)，運(yùn)用多項(xiàng)式曲線擬合的方法，把邊緣點(diǎn)y值作為輸入?yún)?shù)，計(jì)算擬合系數(shù)來求取相應(yīng)的x′和z′，公式為：

整體算式表示為：

x≈x′＝v_x(y)＝v_x1yⁿ+v_x2y^n-1+v_x3y^n-2+...+v_xn′-1y+v_xn′

z≈z′＝v_z(y)＝v_z1yⁿ+v_z2y^n-1+v_z3y^n-2+...+v_zn′-1y+v_zn′

其中，v_x1v_x2v_x3...v_xn′-1v_xn′為計(jì)算得到的多項(xiàng)式擬合系數(shù)；

(132)得到新的葉面邊緣點(diǎn)為P′_edge＝{x′_l，y_l，z′_l；x′_r，y_r，z′_r}，從而得到平滑和無偏差的葉子邊緣點(diǎn)；

在這一部分中，也可以理解為運(yùn)用了多項(xiàng)式擬合分別把葉子邊緣點(diǎn)投影到X-Y平面和Y-Z平面上，接著對(duì)兩個(gè)投影面求交后，從而定位真實(shí)的擬合后的葉子邊緣；

(2)葉面擬合與葉面重建

在傳統(tǒng)的張量積Bezier曲面的基礎(chǔ)上，提出了雙三次廣義張量積Bezier曲面擬合方法；廣義張量積Bezier曲面公式為：

其中，和分別代表按照參數(shù)m、n等間隔分割：取m＝n＝3，得到U，V{0≤U≤1；0≤V≤1}是采樣值，采用精度根據(jù)m、n的大小設(shè)置，D_i，_j和E_i，j分別為處的x方向的目標(biāo)偏導(dǎo)矢和y方向的目標(biāo)偏導(dǎo)矢；

張量積得Bezier曲面，在單只葉面以4*4的16個(gè)控制點(diǎn)，從每條掃描線中依次選取4個(gè)點(diǎn)，分別為p_0，j，p_1，j，p_2，j，p_3，j，再選擇相鄰的四條線段，共取16個(gè)點(diǎn)構(gòu)成網(wǎng)格的16個(gè)點(diǎn)，把這16個(gè)點(diǎn)代入廣義張量積Bezier曲面公式p_i，j中，進(jìn)而采用廣義雙三次Bezier曲面對(duì)其擬合，最后，重復(fù)上述過程直到所有掃描線上的掃描點(diǎn)都遍歷一遍；

(3)葉面受力形變

(31)用Ω來表示葉面形變前的狀態(tài)，p∈Ω是葉面上的某個(gè)點(diǎn)p：(p_x，p_y，p_z)；用Γ表示樹葉形變后的狀態(tài)且p′∈Γ是形變后樹葉上的某個(gè)點(diǎn)p′(p′_x，p′_y，p′_z)，與p對(duì)應(yīng)；根據(jù)Ω，用位移u來描述Γ，即u：Ω→Γ，p′＝p+u(p)，其中：設(shè)p₁，p₂∈Ω，p′₁，p′₂∈Γ，形變前p₁，p₂之間的空間位移d表示為：d＝p₂-p₁，同理，形變后p′₁，p′₂之間的空間距離位移d′：d′＝p′₂-p′₁，由上式可得：

其中：I為單位矩陣，為形變梯度，用F來表示，即

通過d，d′之間的變化來表示其形變：

|d′|²-|d|²＝d′^Td′-d^Td＝d^T(F^TF-I)d

根據(jù)彈性力學(xué)知識(shí)，可由上式得出格林應(yīng)變E∈R^3×3為：

從上式可看出E是對(duì)稱矩陣，對(duì)于以四面體為基本單元組成的葉面模型而言，采用線性插值函數(shù)來表示每個(gè)節(jié)點(diǎn)的位移向量u＝(u，v，w)；對(duì)于一個(gè)由4個(gè)節(jié)點(diǎn)組成的四面體，每個(gè)節(jié)點(diǎn)有3個(gè)自由度，分別沿節(jié)點(diǎn)坐標(biāo)的x，y，z方向，設(shè)四面體四個(gè)頂點(diǎn)分別表示t_i(x_i，y_i，z_i)，t_j(x_j，y_j，z_j)，t_k(x_k，y_k，z_k)，t_l(x_l，y_l，z_l)，通過采用四面體的四個(gè)頂點(diǎn)來表現(xiàn)位移形式的模式稱為單元位移模式或位移函數(shù)，一般用一次多項(xiàng)式表示：

由上式，得到4個(gè)單元位移場(chǎng)函數(shù)為：

S_i為單元的插值函數(shù)或形函數(shù)，V為四面體單元的體積，而a_i，b_i，c_i，d_i，…d_l為：至于a_j，b_j，c_j，d_j，…和d_l各項(xiàng)按右手定則輪換下標(biāo)i，j，k，l即可求得，例如：

由上面式子可以得到：

其中：是四面體每個(gè)頂點(diǎn)的位移向量，I為單位矩陣，它們組成的列向量稱為四面體單元每個(gè)頂點(diǎn)的位移列陣，在三維空間中，Piola-Kirchhoff應(yīng)力σ被計(jì)算出，其表示為3×3的對(duì)稱矩陣，如下所示：

(32)樹葉變形

(321)模型表示

將形變模型表示為(V，G，P)，其中V＝(1，2，…，n)表示頂點(diǎn)集，G＝{(i，j)|i，j∈V，i≠j}是邊的集合，P＝{P_i∈R³|1≤i≤n}是模型中每個(gè)頂點(diǎn)的坐標(biāo)集合；對(duì)于每個(gè)頂點(diǎn)i，用N_i＝{j|(i，j)∈G}表示i的領(lǐng)域的集合；同時(shí)，為了方便計(jì)算，設(shè)定形變模型中所有四面體的質(zhì)量都被平均分配到四個(gè)頂點(diǎn)上，則m_i表示i點(diǎn)的質(zhì)量，有：

其中，ρ_j為四面體密度，為常量；表示與i相鄰的第j個(gè)四面體的體積；使用r向量表示頂點(diǎn)的初始位置，p向量表示頂點(diǎn)形變后的位置，則r_i和p_i分別表示頂點(diǎn)i的初始位置和形變后的位置；

(322)定義形變模型中的勢(shì)能

按下式計(jì)算形變模型的彈性勢(shì)能：

V(P)＝V_r(P)+V_v(P)

其中，V_r(P)表示使形變后的模型恢復(fù)到平衡或靜止?fàn)顟B(tài)的勢(shì)能，即恢復(fù)勢(shì)能；V_v(P)表示保持模型形變所需要的能量，即體積勢(shì)能；

(323)恢復(fù)勢(shì)能的定義

按下式計(jì)算恢復(fù)勢(shì)能：

其中，L_i表示i節(jié)點(diǎn)的微分算子，d_i表示i節(jié)點(diǎn)在平衡或靜止?fàn)顟B(tài)時(shí)的微分坐標(biāo)；

R_i(P)是i節(jié)點(diǎn)的旋轉(zhuǎn)矩陣，λ為形變模型的楊氏模量；對(duì)于樹葉這種固體模型，采用拉普拉斯算子和拉普拉斯坐標(biāo)作為微分算子和微分坐標(biāo)，計(jì)算方法如下：

其中，為邊長(zhǎng)的權(quán)值，l_k為邊e_ij的對(duì)邊e_pq的邊長(zhǎng)，θ_k為對(duì)邊e_pq的二面角，是中心頂點(diǎn)V_i的歸一化權(quán)值，V_dual是中心頂點(diǎn)V_i的對(duì)偶的Voronoi體域的體積，為標(biāo)量場(chǎng)函數(shù)，Φ_i(·)是分段線性基函數(shù)，為中心頂點(diǎn)V_i處的基函數(shù)，v為四面體體積，S_i為相對(duì)的底三角形的面積，為底三角形的指向四面體外的法線方向；

將表示成矩陣形式為：(δ^(x)，δ^(y)，δ^(z))＝L(P^(x)，P^(y)，P^(z))＝M^-1L_s(P^(x)，P^(y)，P^(z))，其中，(P^(x)，P^(y)，P^(z))和(δ^(x)，δ^(y)，δ^(z))分別表示初始網(wǎng)格頂點(diǎn)坐標(biāo)和拉普拉斯坐標(biāo)的三個(gè)分量的向量表示形式，M是體積權(quán)值的對(duì)角線矩陣L_s是網(wǎng)格權(quán)值的對(duì)稱稀疏線性矩陣：

則拉普拉斯算子矩陣可以表示為：L＝M^-1L_s；

根據(jù)網(wǎng)面拉普拉斯坐標(biāo)性質(zhì)可知，體網(wǎng)格的拉普拉斯坐標(biāo)也具有平移不變性，不具有旋轉(zhuǎn)不變性和縮放不變性；

為了計(jì)算恢復(fù)勢(shì)能，還要計(jì)算每個(gè)頂點(diǎn)i的旋轉(zhuǎn)矩陣R_i(P)，計(jì)算方式如下：

r_ij＝m_j(r_j-r_i)

p_ij＝m_j(p_j-p_i)

r_i和p_i分別表示節(jié)點(diǎn)i在變形前后的坐標(biāo)，則節(jié)點(diǎn)i的線性變換矩陣A_i通過求解下式的極小值獲得：

當(dāng)獲得上式的極小值后，可得：

其中，

由于為對(duì)稱矩陣，但此矩陣不包含任何的旋轉(zhuǎn)形變，所以旋轉(zhuǎn)矩陣應(yīng)通過矩陣進(jìn)行極分解獲得；同時(shí)設(shè)定模型的旋轉(zhuǎn)矩陣為T，則R_i(TP)＝TR_i(P)，即：L_iTP-R_i(TP)d_i＝T(L_iP-R_id_i)

由前述可知V_r(TP)＝V_r(P)，即V_r(P)是旋轉(zhuǎn)不變的；

(4)模型形變

模型形變模擬由下面的歐拉-拉普拉斯運(yùn)動(dòng)方程控制：

其中，m是系統(tǒng)的質(zhì)量矩陣，D是系統(tǒng)的阻尼矩陣，V(P)是系統(tǒng)的勢(shì)能，f_ext是系統(tǒng)所受的外力，是系統(tǒng)的內(nèi)力，同時(shí)系統(tǒng)的內(nèi)力的雅克比矩陣為即系統(tǒng)的剛度矩陣K，對(duì)于系統(tǒng)的阻尼矩陣D，采用Rayleigh阻尼：

D＝αM+βK

其中，K∈R^3n，3n是系統(tǒng)的剛度矩陣，M是系統(tǒng)的質(zhì)量矩陣，α，β是常數(shù)；求解過程中設(shè)定β＝0，采用隱式的歐拉方法或者采用隱式的Newmark方法求解上式；最后將體模型剖分成四面體模型，并對(duì)該模型施加外力，通過求解方程獲取受力變形后的葉片。