AU2020103131A4

AU2020103131A4 - Leaf surface reconstruction and physically based deformation simulation based on the point cloud data

Info

Publication number: AU2020103131A4
Application number: AU2020103131A
Authority: AU
Inventors: Lin CAO; Zhidong HU; Ling Jiang; Yiduo Li; Lianfeng Xue; Ting YUN
Original assignee: Nanjing Forestry University
Current assignee: Nanjing Forestry University
Priority date: 2020-10-30
Filing date: 2020-10-30
Publication date: 2021-01-07
Anticipated expiration: 2028-10-30

Abstract

The invention discloses the leaf surface reconstruction and physically based deformation simulation based on the point cloud data. First, obtaining the accurate leaf boundary by polynomial fitting method from the scanned leaf point clouds, and calculating the main leaf vein of the leaf surface through computer graphics method; second, adopting the generalized bicubic tensor product based Bezier surface to fit the leaf point cloud data, so as to remove the noise points caused by wind shaking the scanned targets; finally, according to the mechanical model of solid mechanics, defining the different mechanical material properties of leaf vein and mesophyll, and building the nonlinear finite element deformation equation to calculate the leaf morphology after applying different load forces, and then the real deformation of leaves in natural environment can be simulated. Utilizing computer technique to analyze and acquire leaf properties automatically, the invention can accurately descript the natural morphological characteristics of leaf surface under different stress and physiological conditions. 1/6 L2.3 12 L2. oL2.6 L7 L2.9 Li Figure 1 The calculated scanning lines according to the edge points PI, ,3 33 1P4,1 2. a- P 24 P, , 44 L , L2,3 L26 P ~L2,9 Li Figure 2 The implementation diagram of leaf surface reconstruction using the generalized tensor product based on Bezier algorithm

Description

1/6

L2.3

12 L2. oL2.6 L7

L2.9

Li Figure 1 The calculated scanning lines according to the edge points

PI, ,3 33 1P4,1 2. a- P 24 P, , 44 L ,

L2,3

L26

P ~L2,9

Li Figure 2 The implementation diagram of leaf surface reconstruction using the generalized tensor product based on Bezier algorithm

Leaf surface reconstruction and physically based deformation simulation based

on the point cloud data

TECHNICAL FIELD

[01] The invention relates to reconstruction and physically based deformation simulation method for real leaf of broad-leaved trees based on laser scanned data collected by terrestrial laser scanner (TLS).

BACKGROUND

[02] Broad-leaved tree is any tree within the diverse botanical group of angiosperms, with flat, broad leaves and reticulate veins. Its leaves are evergreen or deciduous falling from the branches in autumn and winter. Moreover, the leaf shape varies with the species of tree, and the leaves may become twisted and distorted or the new growth will curl and become deform. With great economic value, some broad leaved trees are important timber species, some are precious wood, landscape plants, fruits, etc., and others are used as street trees or garden greening trees. The invention explores the leaf surface reconstruction and deformation of two tree species of Micheliafigo (Lour.) Spreng. and Cerasus x yedoensis 'Somei-yoshino'. The data is acquired by terrestrial laser scanner.

[03] The terrestrial laser scanner will not cause any damage to the measured targets and can accurately record the three-dimensional (3D) information of the targets in the form of point cloud data. Moreover, with the advancement of electronics and optics, the 3D laser scanner has become a powerful and incomparable measurement tool in metrology. Therefore, many foreign forestry researchers have carried out in depth research and discussion on the application of terrestrial 3D laser scanning technology in forestry. But so far, it has not been found to use discrete point cloud data for accurate leaf blade modelling and deformation simulation, which is mainly because the tree leaves have various physiological states resulting in different degree of leaf curling, and the external environment (e.g., wind blowing) has a continuous impact on the organs of trees, while it is difficult to obtain point cloud data without influence of external environmental disturbances. Therefore, the following problems need to be solved when using computer technique to automatically simulate the leaf surface construction and deformation.

[04] 1) When the terrestrial laser scanning is adopted for data acquisition, the integrity and reliability of the results are affected by the mutually occluded-vegetative elements and external environment influence such as wind causing the leaf blade vibration.

[05] 2) The luxuriant trees have numerous leaves with various shapes, inclination and azimuth angles. How to identify and segment the point cloud data belonging to different leaves from the whole scanned points of the target tree is a difficult problem to be solved.

[06] 3) The scanning data is easy to has deviation and noise due to the disturbance from external natural environment and the laser scanning coverage field decreased by occlusion. How to remove the deviation and depict the true deformed leaf surface from point clouds is a problem to be considered.

[07] 4) Laser scanning is commonly applied to obtain the discrete point cloud data for restoration of the target phonotypical traits. While the leaves are demonstrated as 3D curved surfaces, it is difficult to design a reasonable fitting algorithm with the purpose of converting the scanned points to real leaf surfaces. Aforementioned limitations are the obstacles to the use of computer theory in forest research, so automatically obtaining accurate forestry indicators from discrete laser point clouds is a critical issue to be solved.

SUMMARY

[08] Purpose of the invention:

[09] In order to overcome the shortcomings existing in the prior studies, the invention provides a method of leaf surface reconstruction and physically based deformation simulation based on the point cloud data, which aims at realizing an accurate and feasible leaf phonotypical trait acquisition and analysis platform with the help of laser scanner technology. Employing the computer graphics technology, the invention exhibits the feasibility to accurately describe the current state of leaf surface under different forest structure conditions and acquire forestry parameters automatically.

[010] Technical scheme:

[011] To achieve the above purpose, the technical scheme adopted by the invention is as follows:

[012] The leaf surface reconstruction and physically based deformation simulation based on the point cloud data comprises the following steps. First, collecting the 3D leaf point clouds of different broad-leaved tree species by high-precision laser scanner; second, obtaining the accurate leaf boundary by polynomial fitting from the point clouds of leaf edges, and calculating the main leaf vein of the leaf surface through computer graphics method; third, adopting the generalized bicubic tensor product based Bezier surface to reconstruct the leaf surface from the scanned points, so as to remove the interferences caused by wind shaking and eliminate the noise point clouds; finally, according to the mechanical model of solid mechanics, defining the different mechanical properties of leaf vein and mesophyll materials, and construct the nonlinear finite element deformation equation to calculate the leaf morphology after applying different external forces, so as to simulate the real deformation of leaves in natural environment.

[013] The specific implementation of the above method includes the following steps:

[014] (1) Data pre-processing

[015] (1.1) Scanning the whole tree and extracting two endpoints of leaf point

cloud, p,=(xY,,Z,) T , p, =(xy,,z,), wherein, p, is the point of the petiole tip

and p, is the apex of the leaf blade, i.e., the top of the midrib. Then calculating the line

L, between p, and p, to determine points on the main vein, i.e., midrib,

LIP: Pp+tx(p,-p,), wherein, t is the independent variable of the line. From the line

equation Li , it can be calculated that the normal vector k2 =-1/(p,-p,)

perpendicular to vector k =(p, - p) and aligned with leaf surface. Dividing the line

L, of the main vein into n equal parts and obtaining the n+1 points coupled with

normal vector k2 to make up n+1 perpendicular lines to form the leaf width, wherein

L 2 :p=pi+txk2 , i=1,2,3...,n+1. Paired with the region of each leaf surface

covered by the laser scanned points, the length of each line L2 should not exceed the

leaf surface region and hence the left end point p2 1 and right end point P 2 , of each

line L2 can be determined, so as to obtain the edge points on both sides of each leaf

surface.

[016] (1.2) Taking left and right end points-P2,1 , and p2 of each scanning line

L2 in step (1.1), and recording them respectively as p2 , : (xy1 ,z,) and P 2,i :

(xr,y,z,),wherein i =1,2,3...,n. Then, a screen strategy was performed to inspect

correctness of all the generated lines L2 . If the length of a line L2 is shorter than its

adjacent lines L 2,,- 1 and L2 ,, which means the left and right endpoints of line L2 is

closer to the leaf midrib than the endpoints of it adjacent lines. For this case, we considered that there is lack of integrity regarding the obtained scanned points on the line L2 due to occlusion or deviation caused by external environment. Hence, the

corresponding generated line L2 should be excluded.

[017] (1.3) Based on step (1.2), after removing all the false lines with biased results, renumbering all the remaining lines sequentially as L2,j j=1,2,3...,m and

m<! n .

[018] At the same time, the left and right points are marked as p2 ,: (x , y , z )

andP2 nj: (X ,yj,z). Then, the following results can be obtained as:

[019] the set of left endpoints of every line L 2 ,j is

[020] pe =se xfyz (int soZf in yz , =1,2,3. ,rmi.n

[021] the set of right endpoints of every line is

[022] p2, =x ,x, y, z ),...(x,yi, z ,j= 1,2,3... m. (2)

[023] The sets of the left and right endpoints of every line L constitute the

initial scanning contour of each leaf surface. For the left or right part contour of a leaf blade, the combination of polynomial curve fitting and intersection of different generation planes are designed to depict the accuracy edges of both sides of a real leaf surface.

[024] The specific steps are as follows:

[025] (1.3.1) For the left edge points p =(xL,yL,zL) of the leaf surface, using

the polynomial curve fitting method and taking the value y of the edge points as the

input variables to calculate the fitting coefficients and obtain the corresponding fitted

value for X and z' . The formulas are as follows:

n in in y )

[026] Y-1 >1 (3) n in in y |

[027] where YL x' yf) v 0+v y)+v,2 y +...v /_,(y and

L7:z = + y...v . Then, calculating the partial

derivatives of v and V in equation 3 to retrieve the values of coefficients

vYr,O',' ,2'",- v,...vv, for the polynomial curve fitting.

[028] (1.3.2) The smooth and unbiased left edge points of each leaf surface are

obtained after polynomial curve fitting and denoted asP'g,=-X-'L ZfL} . An

alternative explanation of our approach in (1.3.1) is using polynomial curve fitting method to obtain the two curves L, and L, on the X-Y and Y-Z planes, respectively.

The two curves L and L, are taken as the two generatrices and move along the

given paths (directrixes) of Z and X axes, respectively to form the two curved surfaces. The intersection of the two curved surfaces determines the leaf edge on the left side.

[029] The same operation is performed for the endpoints p', =(xR,yR,zR) on

the right side of each line L to depict the leaf edge on the right side.

[030] (2) Leaf surface fitting and reconstruction

[031] Based on traditional tensor product Bezier surface, the invention proposes generalized bicubic tensor product based on Bezier surface fitting method. The specific formula is as follows:

[032] Q(U, V)=I [py,j +(U -, )Di, + +(V - yJ) Ej, ]Bi,,(u)Byi,(V) i=0 j=0

(4)

[033] Wherein, I and yj represent the uniform intervals according to parameters

m and n , respectively. YO < Y, <...< , and YO < y5 <...< Y, satisfy the constrains

that 1 -z,=,- _ and y -y=,-y-. Setting m=n=3 to get x O=yO

, 1=yT =1/3, 2=y2= 2/3, 73 = Y=1. U,V {0 U 1;0: V!1}are sample values

and their values can be assigned manually, and u =(U- Y7)/(T,, -T)

, V=(V-yO)/(y, -YO). D and Ej are target partial vectors of (,,yj) in the x

direction and y -direction, respectively.

!a ~/(h 1 +h,)Fh, ((p _ - p h_1)+h_((p - p /h if i #0 and i m

[034] D = a p -po )/ho) if i=0

apj((p.,j - pr1,j,)/h,_)ifi=m

!i (fi +fj)( , 1 /-p,_&)/ +f p( -p P/j)f f j#0andj#n

[035] E, 1 = p 1 -J)/if) if j=0 -((p

S p((P, -P , if j=n

[036] where h=xj 1 -x,(i=0,1,...,m-1) and f=yj -yj(j=0,1,...,n-1), and

[037] a , E [0,1] are shape-related parameters and here a and p are

assigned to 0.8.

B, (u)= (1-u)' Bo, (v)= (1 -v)'

B10 (u)=3u(1-u)2 B,()31-2

[038] B (u) = and B B(,)(v) =3v(-v) B 2, (u)=3u 2(1-u) B ,,(v)= 3v 2(1-v)

(B 3 (u) (B3,(v)= V3

[039] The fitting efficacy using generalized tensor product based Bezier surface is better than the tensor product Bezier surface, which makes the fitting surface closer to the original surface of the studied leaf.

[040] Obtaining the Bezier surface through generalized tensor product and selecting 16 control points in an arrangement of 4*4 on the leaf surface, where every 4 points are selected from the 4 adjacent scanning lines in turn recorded as

Po,o, P0,1, PO,2 , Po,3

[041] p 1= 1,'' '''1,2I P1,3 . Bringing these 16 points into formula 4, so as P 2 ,0 , P 2 ,1I P 2 , 2 , P2,3

_P3 ,o, P 3 ,1 , P3,2, P 3 ,3]

to obtain the fitted surface. At last, repeating this process until all scanning points on all the scanning lines are processed to achieve the generation of the fitted leaf surface.

[042] (3) Deformation of the fitted leaf surface under stress

[043] (3.1) Using 0 to represent the state of the leaf before deformation, and

pe Q is a certain point on the leaf, wherein, Q c R'andp(p,,p,p, ). Using F to

represent the state of the leaf after deformation F cR 3 and p's F is a certain point on

the deformed leaf p'(p' Ip, p' , corresponding top. The displacement u is used to

describe the transformation from Q to F, that is u: Q->F , p' = p + u(p), wherein, if

p, p2 E Q , p, p2 e F , the spatial displacement d between pI, p 2 before deformation is

d p 2 - p ; in the same way, after deformation, the spatial displacement d between

p',pis d': d'= p' - p' . From the above formula, we can get that:

[044] d'=p2 +u(p 2 )-p 1 -u(p)=d+ u(p,+d)-u(p,)=(I+Vu)d (5)

[045] Wherein: I is identity matrix, I+Vu is deformation gradient, denoted as F =I+Vu . The decomposition of vector u into the three mutually orthogonal au au au ax ay az components (u,v, w) and Vu -

. ax ay az Ow Ow Ow ax ay az

[046] The deformation can be represented by the change between d, d', that is

[047] |d 2 -d2=d' Td'-d T d= d(FTF-I)d

[048] According to theory of elasticity, Green strain EE R 3 can be obtained from the above formula as follows:

1 1 E = -(F T F -I)= -(V+Vu T +Vu T Vu) (6)

[049] 2 2

[050] Besides, from the above formula, it can be noted that E is a symmetric matrix.

[051] For the leaf model which is composed of tetrahedron as the basic element, linear interpolation function is used to represent the displacement vector of each node, that is u=(u,v, w). For a tetrahedron composed of four nodes, each node has three degrees of freedom, along the x,y, z directions of the node coordinates. Let the four

vertices of tetrahedron be represe noted respectively asas,t, ) (t , t, ti",),

P,(t., tkyItk,, ) and P(tt , , and each vertex changes in the x , y and z

directions. The displacement mode expressed by the four vertices of tetrahedron is called element displacement mode or displacement function, which is generally expressed by a polynomial:

U5Fu =C1 ±+C1 2 x+C1 3 y+C1z4

[052] v=C2 1 +C22 x+C23 y+C 34z (7) (w=C 3 1 +C3 2 x+C33 y+C3 4 z

[053] From the above formula, the displacement field functions of four vertices are obtained as follows:

1 S,= ,(a,+bix+cy+dz) 6V'

Si = I(a 1 + bjx+ cjy+ djz)

[054] 6'(8) Sk = ,a+ bkx + cky+ dkz) 6V k 1 S, = V , (a,+b~x+cy + dz) 6 (

[055] Si is the interpolation function or shape function of the tetrahedron, Vis

1 t tj t,

11 t1 t. t the volume of tetrahedral element,V'- - ;t a", b, c, d,, -d, are 61 t k, t tkz

I t t t

respectively as follows:

tJ, ty j,z I ,y tI, tjx t,

[056] a = tk, tky tkz , b, 1 tky tkz C= tkx tkz

tiX tI, ti,z 1 t' tl t 1 t,

t, t., 1 tj, j,y I

d=- tk, tky 1

tix tiy 1

[057] As for a, bj, c j, dj, a k -.-and d, , they can be obtained by alternating

subscripts i, j, k, according to the right hand rule, for example,

t, tk,, t, tI ti, t tkx tkx ti,

[058] a = ti, yt, ti,z a2 =t, ti, t a,=t t t ti, i'y ti, tj,, j'y tj, k, k,x k,,

[059] From the above formula, we can calculate the displacement u:

ti

[060] u [IS, IS] ISk IS,] (9)

t j

[061] Wherein, ,i,f, I and i1 are 3x 1 displacement vectors of each vertex of

tetrahedron, I is identity matrix and S is 3x 3 displacement matrix for each vertex. Hence, the transformation for the current displacement u is retrieved.

[062] (3.2) Leaf deformation

[063] (3.2.1) Model representation

[064] The deformation model is represented as (V,G,P) , wherein

V={Vh=1,2,...,n} is tetrahedral set, which constitutes the leaf surface.

ER 3 |1i P=tjE 4} is the coordinate set of each vertex of a tetrahedron and

G={(I ,P-)|I-, ce V,i j} is a set of edge of tetrahedron. As for each vertexi,

N, ={P(,I )e G} is used to represent the set of the adjacent vertices of 1. At the

same time, for the convenience of calculation, it is assumed that the mass of all tetrahedrons in the deformation model is concentrated in their four vertices, so that the

mass mi of the point Ican be calculated as follows:

[065] mi= Z 1 -p V Vje N,

[066] Wherein, p, a constant, is tetrahedral density; V is the volume of the j

th tetrahedron with a vertex ofJ .

[067] (3.2.2) Definition of potential energy in deformation model

[068] Calculating the elastic potential energy of the deformation model based on following formula:

[069] T(P)= T(P)+T (P) (10)

[070] Wherein, 7(P) represents the potential energy that makes the deformed

model return to equilibrium or initial state, namely recovery potential energy; T,(P)

represents the energy required to maintain the deformation of the model, namely volume potential energy.

[071] (3.2.3) Definition of recovery potential energy

[072] Calculating the recovery potential energy based on following formula:

[073] T(p)= A L(P)(-R,(P);|2

[074] Wherein, L is differential operator of vertex It and r, is differential

coordinates of vertex I at equilibrium or static state. R(P) is rotation matrix of 13

and A is the Young's modulus of the deformation model. For the leaf solid model, Laplace operator and Laplace coordinate are used as differential operator and differential coordinate, and the calculation method of Laplace operator is as follows:

[075] 5,=A (v )=w w,f(P )-f(?)) (12) Vjc Ni

[076] Wherein, w= l, cot(Ok) is the weight of the side length, lkis the side

length of G(P,,P), opposite side of G(,POkis the dihedral angle of opposite side

1 G(P,,1, w = is the normalized weight of the centre vertex P-, V., is the (V.,)i

volume of the dual Voronoi domain of the central vertex P, as shown by the gray part

of the tetrahedron in Figure 7. f(P)=Z fJ ,(P) is a scalar field function, F() is a

piecewise linear basis function, $(I()= ,nP is the basis function at the central 3V vertex P, V' is the volume of tetrahedron, s is the area of the relative bottom triangle

PPP and hp is the normal direction of the bottom triangle pointing to the outside of

the tetrahedron. Schematic representation is shown in Figure 7.

[077] Formula (12) can be expressed as a matrix-form.

[078] S(S,6 , S tt) =ML ,P,) (13)

[079] wherein, 1(t1,,t ,t )andQo ,d,)arevectorrepresentationsofthe

three components of the initial tetrahedral vertex coordinates and Laplacian 1 coordinates, respectively, M is the diagonal matrix of the volume weightM 1-

L, is a symmetric sparse linear matrix of grid weights

- Ywi ,f i-=j kcNi

[080] (L,) = w , if Pe N,(e )

0, otherwise

[081] Then the Laplacian matrix can be expressed as L = M-L,.

[082] According to the properties of Laplacian coordinates, the Laplacian coordinates of volume meshes are invariant to the translation and are variant to the scaling and rotation.

[083] In order to calculate the recovery potential energy, the rotation matrix R,(P)

of each vertex p. is also calculated as follows:

[084] before =m.(P°before -before)

[085] Fafter= m(Pfter fer)

[086] Wherein, 'pbfo and 1afer represent the coordinate of the point I before

and after deformation, respectively. Then the linear transformation matrix A. of vertex

can be obtained by finding the minimum of the following equation

[087] |fter -A fore |2 jE N(i)

[088] When the minimum value of the above formula is obtained, we can get the following results:

[089] A = Af'A

[090] Wherein, Ar (Fafer (Fbefore ) and" ( before( or) )-1.

[091] Since Ar is a symmetric matrix without any rotation deformation, the

rotation matrix should be obtained by polar decomposition of the matrixAfr. At the

same time, setting the translation matrix of the model as T, so that R (TP)= TR,(P), that is,

[092] LTP-R,(TP)d =T(LiP-Rd,)

[093] It can be seen from the above that ,(TP)= f, (P) , meaning V,(P) is translation invariant.

[094] (3.3) Model deformation

[095] The deformation simulation of the model is controlled by the following Euler-Laplace motion equation:

[096] MP+ DJ+ Op= fe, (14)

[097] Wherein, M is the mass matrix of the system, D is the damping matrix of the system, T(P) is the potential energy of the system, f, is the external force on the

DV(P) system, is the internal force of the system. At the same time, the Jacobian Op 8 2 T (P) matrix of the internal force of the system is , namely, the stiffness matrix of the sPe P system K .For the damping matrix Dof the system, Rayleigh damping is adopted, so

[098] D= aM+§8K

[099] Wherein, KE R(3nx 3n) is the stiffness matrix of the system, M is the

mass matrix of the system and a and # are constant. In the solution, setting# =0 and

using implicit Euler method or implicit Newmark method to solve the above equation.

[0100] In general, dividing the volume model into many tetrahedral models, and applying external force to the model, so as to obtain the deformed leaf surface by solving the equation. According to the specific steps of our method, the deformation of any leaf surface under mechanical loading can be implemented, and the results can be validated by real trails with the selected leaf surface.

[0101] Beneficial effects:

[0102] The invention provides a method for leaf surface reconstruction and physically based deformation simulation based on the point cloud data. It can build an accurate and feasible 3D leaf surface modelling and analysis platform with the help of laser scanner technology. Besides, with the latest method of graphics and image science, the invention is benefit for describing various morphological structures of broad leaves under different stress and physiological conditions.

BRIEF DESCRIPTION OF THE FIGURES

[0103] Figure 1 is the calculated scanning lines according to the edge points.

[0104] Figure 2 is the implementation diagram of leaf surface reconstruction using the generalized tensor product based on Bezier algorithm.

[0105] Figure 3 the edge depiction and main vein generation for the leaf surface, carried out from step (3a) to step (3f).

[0106] Figure 4 is a process diagram of leaf surface reconstruction, guided by many fitted lines, carried out from step (4a) to step (4b).

[0107] Figure 5 is a set of diagrams showing polygonization and deformation of leaf surface under different external force, carried out from step (5a) to step (5f).

[0108] Figure 6 is another set of polygonization and deformation diagram of leaf surface under external stress, carried out from step (6a) to step (6f).

[0109] Figure 7 is the schematic representation showing the calculation of the Laplace operator.

DESCRIPTION OF THE INVENTION

[0110] The present invention is further described in combination with the attached figures.

[0111] A real leaf surface reconstruction and physically based deformation simulation for broad-leaved trees based on laser point cloud data comprises the following steps. Firstly, collecting the 3D leaf point clouds of different broad-leaved tree species by high-precision laser scanner; second, obtaining the accurate leaf boundary by polynomial fitting from the point clouds of leaf edges, and calculating the main leaf vein of the leaf surface through computer graphics method; thirdly, adopting the generalized bicubic tensor product based on Bezier surface to reconstruct the leaf surface from the scanned points, so as to remove the interferences caused by wind shaking and eliminate the noise point clouds; finally, according to the mechanical model of solid mechanics, defining the different mechanical properties of leaf vein and mesophyll materials, and construct the nonlinear finite element deformation equation to calculate the leaf morphology after applying different external forces, so as to simulate the real deformation of leaves in natural environment.

[0112] The specific implementation of the above method includes the following steps:

[0113] (1) Data pre-processing

[0114] (1.1) Scanning the whole tree and extracting two endpoints of leaf point

cloud, p,=xyz')e , p xy,z,)T, wherein, p, is the point of the petiole tip

Li between p, and p, to determine points on the main vein, i.e., midrib,

L : p=p,+tx(p, - p), wherein, t is the independent variable of the line. From the line equation L, , it can be calculated that the normal vector 2=-1/(P,-Pe) perpendicular to vector k =(p, - p) and aligned with leaf surface. Dividing the line

L 2 :p=pi+txk 2 , i=1,2,3...,n+1. Paired with the region of each leaf surface

covered by the laser scanned points, the length of each line L should not exceed the

leaf surface region and hence the left end point p 2 , and right end point P 2 ,, of each

surface.

[0115] (1.2) Taking left and right end points- p2 1, and p2 ,, of each scanning line

L2 in step (1.1), and recording them respectively as p2, : (x,,y,,z,) and p :

(x,y,,z,.)), wherein i =1,2,3...,n. Then, a screen strategy was performed to inspect correctness of all the generated lines L 2 . If the length of a line L2 is shorter than its

adjacent lines L2,,-, and L2 ,, which means the left and right endpoints of line L2 is

corresponding generated line L2 should be excluded.

[0116] (1.3) Based on step (1.2), after removing all the false lines with biased results, renumbering all the remaining lines sequentially as L 2 j, j =1,2,3...,m and

m < n . At the same time, the left and right points are marked as p2 ,k: (x,,y,z, )and

P2 ,,: (Xj,y,zj). Then, the following results can be obtained as:

[0117] the set of left endpoints of every line L, is

[0118] pe =xfgYtZ)iY Z ine jt=f1,y23...

[0119] the set of right endpoints of every line is

[0120] p2 =x y X2z),(x ,zf),(xiyz , j=1,2,3...,m. (2)

[0121] The sets of the left and right endpoints of every line L2 , constitute the

[0122] The specific steps are as follows:

[0123] (1.3.1) For the left edge points p =(xL,yL, zL) of the leaf surface, using

value for X and z' . The formulas are as follows:

m in in o

[0124] Y-1 (3) ninE8 qin yL)-z

[0125] wherein L 1 :x'L=p (yjj v -+v y +...+vnl _y and

L7:z =v y ... v_ . Then, calculating the partial

derivatives of v, and V in equation 3 to retrieve the values of coefficients

vr,O',v2'--v,rni , ... v v_,,_1,vOv for the polynomial curve fitting.

[0126] (1.3.2) The smooth and unbiased left edge points of each leaf surface are

obtained after polynomial curve fitting and denoted as P',,-={x',yL L . An

The two curves L and L, are taken as the two generatrices and move along the

given paths (directrixes) of Z and X axes, respectively to form the two curved surfaces, The intersection of the two curved surfaces determines the leaf edge on the left side.

[0127] The same operation is performed for the endpoints p, =(x'y,<z) on

the right side of each line L to depict the leaf edge on the right side.

[0128] (2) Leaf surface fitting and reconstruction

[0129] Based on traditional tensor product Bezier surface, the invention proposes generalized bicubic tensor product based on Bezier surface fitting method. The specific formula is as follows:

[0130] Q(U, V)=I [py,j +(U -Y,)Di, + +(V - Y I)Ej, ]Bj,,(u)B ,n(v) i=0 j=0

(4)

[0131] Wherein, Yj and y represent the uniform intervals according to parameters

m and n , respectively. YO < Y, <...< , and Yo < y <...< Y, satisfy the constrains

that -z,= -,_, and y -y7,=y7-y . Setting m=n=3 to get xo =yO

1=y=1/3,T2=y 2 =2/3, 3 =y 3 =1. U,V{0 UO V1;0 V 1}are sample values

and their value can be assigned manually, and u = (U - Y,) /(x, - Y,)

, v=(V-yo)/(j, -O). D and Ej are target partial vectors of (i,Y.) in the x

direction and y -direction, respectively.

!a j/h 1 +h,)Fh, ((p _ - p /h1)+h_ 1(( -p /h if i #0 and i m

[0132] D ao, ((PI-pO )/ho) if i=0

aLj ((p -p,j _)/h,)ifi=m

!i, (fi+f )&(, p 1 -p,_)/ + f(( p, -p P/f )# if j 0 andj#n

[0133] E,1 = p 1 -J)/if) if j=0 -((p

p p((P -p )/f_ if j= n

[0134] wherein hi =x -x (i= 0,1,...,m -1) and f=y - yi (j= 0,1,..., n -1), and

[0135] a,l,Ge [0,1] are shape-related parameters and here a y and #,j are

assigned to 0.8.

Bo, (u)= (1-u)' B, (v)= (1 -v)'

[0136] Bu= B101(u)=3u(1-u)2 '" and B B ()3(-) B(,)(v) =3v(-v) B 2,(u)=3u 2 (1-u) B, (v)=3v 2(1-v)

(B 3,(u) (B3, (v)= V3

[0137] The fitting efficacy using generalized tensor product based on Bezier surface is better than the tensor product Bezier surface, which makes the fitting surface closer to the original surface of the studied leaf.

[0138] Obtaining the Bezier surface through generalized tensor product and selecting 16 control points in an arrangement of 4*4 on the leaf surface, where every 4 points are selected from the 4 adjacent scanning lines in turn recorded as

P0,0 I P0,1 I P0,21 PO9,3

p 1 ,0 ,p 1 ,1 ,p 1 ,2 ,p 1 ,3

[0139] p1 = . Bringing these 16 points into formula 4, so as P2 ,0, P 2 ,1, P 2 ,2 , P 2 ,3

_P3 ,0 , P 3 ,1 , P 3 ,2 , P3 ,3 ] to obtain the fitted surface. At last, repeating this process until all scanning points on all the scanning lines are processed to achieve the generation of the fitted leaf surface.

[0140] (3) Deformation of the fitted leaf surface under stress

[0141] (3.1) Using Q to represent the state of the leaf before deformation, and p e Q

is a certain point on the leaf, wherein, Q c R3 and p:(p,, p,, pz). Using F to represent

the state of the leaf after deformation F cR 3 and p' F is a certain point on the

deformed leaf p'p',p',pf) corresponding to p . The displacement u is used to

describe the transformation from Q to F, thatis u:->, p' = p+ u(p), wherein, if

p, p2 EQ ,p' e F , the spatial displacement d between pI, p 2 before deformation is

p',pis d': d'= p' - p' . From the above formula, we can get that:

[0142] d p=p 2 +u(p2)-p 1 -u(p) =d+u(p,+d)-u(p,)=(I+Vu)d (5)

[0143] Wherein: I is identity matrix, I+Vu is deformation gradient, denoted as F=I+Vu . The decomposition of vector u into the three mutually orthogonal

au au au ax ay az

components (u,v, w) and Vu= ax ay az aw aw aw ax ay az

[0144] The deformation can be represented by the change between d, d', that is

[0145] |d1 2 -Id2=d'Td'-d Td=dT (FTF-I)d

[0146] According to theory of elasticity, Green strain EE R 3 can be obtained from the above formula as follows:

1 1 E=-(FT F-I)=-(V+Vu T +Vu T Vu) (6)

[0147] 2 2

[0148] Besides, from the above formula, it can be noted that E is a symmetric matrix.

[0149] For the leaf model which is composed of tetrahedron as the basic element, linear interpolation function is used to represent the displacement vector of each node, that is u=(u,v, w). For a tetrahedron composed of four nodes, each node has three

degrees of freedom, along the x,y, z directions of the node coordinates. Let the four

vertices of tetrahedron be represe noted respectively asas,t, ) (t , t, ti",),

and n(, (,t, , ) t ,tt and each vertex changes in the x , y and z

[0] u =C1 ±+C1 2 x+C1 3 y+C14 z

[0150] v=C21 +C 22 x+C23 y+C 34z (7) (w=C 3 1 +C3 2 x+C33 y+C3 4 z

[0151] From the above formula, the displacement field functions of four vertices are obtained as follows:

1 S, = , (a + bix+ ciy+ diz) 6V'

Si = I(a,+bjx+cy+dz)

[0152] 6Va (8) 1 Sk = ---- a, + bkx + cky + dkz) 6V k 1 S, = , (a,+b~x+cy + dz) 6V (

[0153] Si is the interpolation function or shape function of the tetrahedron, Vis

1 t t t. 11 tj tj t, the volume of tetrahedral element, V'=- ; aj, bi, ci, d, - d, are 61 tkx tky tkz

1 t t t

respectively as follows

x, y, z1 1 yj zi xj 1 z xj y 1

[0154] a, = xk Yk Zk ,b,=-1 Yk Zk, Cik = x 1 zk, di= -Xk Yk 1 x, Y, Z, 1 y, Z, X, 1 z x, y, 1

[0155] As for aj,b, c, d j, -and d,, they can be obtained by alternating subscripts

i, j,k,according to the right hand rule, for example,

Xk Yk Zk X1 y, z, Xi Y, z,

[0156] a1 = x, y, z, a = x, y, zi a,= zj xj yj X, Yi Z, X1 Y1 z1 Xk Yk Zk

[0157] From the above formula, we can calculate the displacement u: t.

[0158] u= I S S S ] -' (9) tk

ti

[0159] Wherein,,T, I and , are 3x1 displacement vectors of each vertex of

[0160] (3.2) Leaf deformation

[0161] (3.2.1) Model representation

[0162] The deformation model is represented as (V,G,P) , wherein, V=(1,2,...,n) is tetrahedral set, which constitutes the leaf surface.

P={pie R 3 1i4 is the coordinate set of each vertex of a tetrahedron and

G={(1,P-)|I,=e V,i #j} is a set of edge of tetrahedron. As for each vertex i,

{P I(I-,P)e G} is used to represent the set of the adjacent vertices of 1 . At the same time, for the convenience of calculation, it is assumed that the mass of all tetrahedrons in the deformation model is concentrated in their four vertices, so that the

mass mi of the point P can be calculated as follows:

11

[0163] mi = -pVf 4 Vje N

[0164] Wherein, p1 , a constant, is tetrahedral density; V' is the volume of the j

th tetrahedron with a vertex of 1' .

[0165] (3.2.2) Definition of potential energy in deformation model

[0166] Calculating the elastic potential energy of the deformation model based on following formula:

[0167] T(P)= T(P)+ T (P) (10)

[0168] Wherein, V,(P) represents the potential energy that makes the deformed

model return to equilibrium or initial state, namely recovery potential energy; V,(P)

represents the energy required to maintain the deformation of the model, namely volume potential energy;

[0169] (3.2.3) Definition of recovery potential energy

[0170] Calculating the recovery potential energy based on following formula:

[0171] T(p)= 2 | I Li(P)- R(P)r ||2 ic~v

[0172] Wherein, Li is differential operator of vertex 1. and r, is differential

coordinates of vertex P at equilibrium or static state. R,(P) is rotation matrix of I

[0173] , = A,(Jw)=w w(f(,)- f(1p)) (12) VjE Ni

[0174] Wherein, wj= - Z cot(O,) is the weight of the side length, lk is the side

lengthof G(P,I), opposite side of G(,P),Okis the dihedral angle of opposite side

G(P,,1P, w, =( is the normalized weight of the centre vertex '-, Vda, is the

volume of the dual Voronoi domain of the central vertex 13, as shown by the gray part

of the tetrahedron in Figure 7. f(P)=Z fl,(P) is a scalar field function, F() is a

-s n piecewise linear basis function,,(13)= ( is the basis function at the central 3V' vertex P, V' is the volume of tetrahedron, s is the area of the relative bottom triangle

P-PP and f is the normal direction of the bottom triangle pointing to the outside of

the tetrahedron. Schematic representation is shown in Figure 7.

[0175] Formula (12) can be expressed as a matrix-form.

[0176] ,8 5( ,,3 )= I(ti, t , t,) = M- (P, P,) (13)

[0177] Wherein, -(t,,,t,,t) and are vector representations of

the three components of the initial tetrahedral vertex coordinates and Laplacian 1 coordinates, respectively, M is the diagonal matrix of the volume weightM = w. L, is a symmetric sparse linear matrix of grid weights

- wi, I i-j ke Ni

[0178] (L, ), = i w, if Pg- EN(i) 0, otherwise

[0179] Then the Laplacian matrix can be expressed as L = M-1L.

[0180] According to the properties of Laplacian coordinates, the Laplacian coordinates of volume meshes are invariant to the translation and are variant to the scaling and rotation.

[0181] In order to calculate the recovery potential energy, the rotation matrixR1 (P)

of each vertex ] is also calculated as follows:

[0182] F 'et"r"=m (Pef"t""-peore

[0183] F ai = my(Pair -pIta)

[0184] Wherein, .hef"rand .""" represent the coordinate of the point 1 before

and after deformation, respectively. Then the linear transformation matrix A of vertex

can be obtained by finding the minimum of the following equation

[0185] || 1 ' 4 .I'r|2 jc N(i)

[0186] When the minimum value of the above formula is obtained, we can get the following results:

[0187] AA ,"rr

[0188] Wherein, A'=( '"(F; or) and Arr ( before (Fbefre) )

[0189] Since Ar is a symmetric matrix without any rotation deformation, the

rotation matrix should be obtained by polar decomposition of the matrix Af'. At the

same time, setting the translation matrix of the model as T, so that R (TP)= TR.(P),

that is,

[0190] LITP - R,(TP)d, = T(LP - R,d )

[0191] It can be seen from the above thatV,(TP)=V,(P) , meaning V,(P) is

translation invariant.

[0192] (3.3) Model deformation

[0193] The deformation simulation of the model is controlled by the following Euler-Laplace motion equation:

[0194] MP+ DP+ P= f, (14)

[0195] Wherein, M is the mass matrix of the system, D is the damping matrix of

the system, V(P) is the potential energy of the system, f, is the external force on the

DV(P) system, is the internal force of the system. At the same time, the Jacobian ap 2 matrix of the internal force of the system is - P), namely, the stiffness matrix of the aPcP system K. For the damping matrix D of the system, Rayleigh damping is adopted, so

D = aM +/K.

[0196] Wherein, Ke R(3nx 3n) is the stiffness matrix of the system, M is the

mass matrix of the system and a and p are constant. In the solution, setting§3=0 and

[0197] In general, dividing the volume model into many tetrahedral models, and applying external force to the model, so as to obtain the deformed leaf surface by solving the equation. According to the specific steps of our method, the deformation of any leaf surface under mechanical loading can be implemented, and the results can be validated by real trails with the selected leaf surface.

[0198] Table 1. The comparison regarding the calculated area of the reconstructed leaf surface versus the manual measurement value from the studied leaf in real word

The studied Leaf area measured by LI- Area of the reconstructed Deviations leaf 3000C (cm2) leaf surface (cm 2

) 1 54.3 51.7 4.79% 2 48.1 50.2 4.18% 3 43.7 45.6 4.17%

[0199] Table 2. The number of the scanned points for each leaf and the time consuming for each step of leaf blade reconstruction from point clouds

Number of Number of Time-consuming (s) The studied scanned scanned points Edge Leaf eaf surface leaf points after Be e depiction surface reconstruction surface fitting dpcin fitting 1 1507 1078 0.257 1.631 0.207 2 1326 928 0.225 1.245 0.182 3 1073 789 0.197 1.101 0.176

[0200] Although the invention has been described with reference to specific examples, it will be appreciated by those skilled in the art that the invention may be embodied in many other forms, in keeping with the broad principles and the spirit of the invention described herein.

[0201] The present invention and the described embodiments specifically include the best method known to the applicant of performing the invention. The present invention and the described preferred embodiments specifically include at least one feature that is industrially applicable

Claims

THE CLAIMS DEFINING THE INVENTION ARE AS FOLLOWS:

1. A real leaf surface reconstruction and physically based deformation simulation of broad-leaved trees based on laser point cloud data, characterized by comprising the following steps. Firstly, collecting the 3D leaf point clouds of different broad-leaved tree species by high-precision laser scanner; second, obtaining the accurate leaf boundary by polynomial fitting from the point clouds of leaf edges, and calculating the main leaf vein of the leaf surface through computer graphics method; third, adopting the generalized bicubic tensor product based on Bezier surface to reconstruct the leaf surface from the scanned points, so as to remove the interferences caused by wind shaking and eliminate the noise point clouds; finally, according to the mechanical model of solid mechanics, defining the different mechanical properties of leaf vein and mesophyll materials, and construct the nonlinear finite element deformation equation to calculate the leaf morphology after applying different external forces, so as to simulate the real deformation of leaves in natural environment.

2. The real leaf surface reconstruction and physically based deformation simulation of broad-leaved trees based on laser point cloud data according to claim 1, characterized by including following specific steps:

(1) Data pre-processing

(1.1)Scanning the whole tree and extracting two endpoints of leaf point cloud,

T Pe =(Xe Ye, z) , P, =(x, y, z)T,wherein, p, is the point of the petiole tip and p, is

the apex of the leaf blade, i.e., the top of the midrib. Then calculating the line L,

between p, and p, to determine points on the main vein, i.e., midrib,

LI: p-p,+tx(p, -pe), wherein, t is the independent variable of the line. From the line

equation Li , it can be calculated that the normal vector k2 =-1/(p,-Pe)

perpendicular to vector k, =(p, - pe) and aligned with leaf surface. Dividing the line

L 2 :p=pL+txk 2 , i=1,2,3...,n+1. Paired with the region of each leaf surface covered by the laser scanned points, the length of each line L should not exceed the leaf surface region and hence the left end point p2 , and right end point P2 ,, of each line L2 can be determined, so as to obtain the edge points on both sides of each leaf surface.

(1.

2) Taking left and right end points- p2 1 and p2,, of each scanning line L2

in step (1.1), and recording them respectively as p2 : (x 1 ,y 1,z 1 )and p :

(Xi,YiZ), wherein i=1,2,3...,n. Then, a screen strategy was performed to inspect

adjacent lines L2 , and L, , which means the left and right endpoints of line L2 is

corresponding generated line L2 should be excluded.

(1.

3) Based on step (1.2), after removing all the false lines with biased results, renumbering all the remaining lines sequentially as L 2 , j= 1,2,3..., m and m < n

. At the same time, the left and right points are marked asp 2 : (x,,y,,zb)and

p2 , : (x, y,, z, . Then, the following results can be obtained as:

the set of left endpoints of every line L is

p = , ),(x2,yz j=1,2,3...,m. (1)

the set of right endpoints of every line is

p, =, x2y z)(x iz2),...,(xM MJXRYR )(2, R,y R,z , j ,m. (2) =1,2,3 ...

The sets of the left and right endpoints of every line L constitute the initial

scanning contour of each leaf surface. For the left or right part contour of a leaf blade, the combination of polynomial curve fitting and intersection of different generation planes are designed to depict the accuracy edges of both sides of a real leaf surface.

The specific steps are as follows:

(1.3.1) For the left edge points p', =(x , y , z ) of the leaf surface, using the

polynomial curve fitting method and taking the value y of the edge points as the input

variables to calculate the fitting coefficients and obtain the corresponding fitted value

for x' andzL. The formulas are as follows:

mn (y-xn (3) m n'<py

where L ( Y ) Y~, +v",(Y +VX2 (Y m ) 2 "'Y,,n- ()nn and

L :z vy .. v y . Then, calculating the partial

derivatives of v and Vz in equation 3 to retrieve the values of coefficients

vx,I,,v, Iv, 2 ,...v, -i'v,,v, 1,vz 2,2,...vzrn-i for the polynomial curve fitting.

(1.3.2) The smooth and unbiased left edge points of each leaf surface are obtained after polynomial curve fitting and denoted as I' ' L }. An

The two curves L,, and L,, are taken as the two generatrices and move along the

The same operation is performed for the endpoints p, =(x,) on the

right side of each line L to depict the leaf edge on the right side.

(2) Leaf surface fitting and reconstruction

Based on traditional tensor product Bezier surface, the invention proposes generalized bicubic tensor product based on Bezier surface fitting method. The specific formula is as follows:

m n Q(U, V)= $$[pi,jy+(U -,)Di,jy+(V -y,)Ei, ]Bi,,(u)B ,n(v) (4) i=0 j=0

Wherein, Yj and yJ represent the uniform intervals according to parameters m

and n respectively. YO < Y, <...< ,, andyO < y <...< Y, satisfy the constrains that

,-(_= -- _,and y, -yj=y-yj . Setting m=n=3 to get To=yo=0,

71 = y, 1/3, T2 =i 2 = 2/3, 3=y=1. U,V {0< U 1; 0: V! 1} are sample values

and their value can be assigned manually, and u = (U - Y,) /(,, - Y,)

v= (V-yo)/(y, - jO). D, and E are target partial vectors of (2,y.) in the x

direction and y -direction, respectively.

!aj/(h +hi) )h((pi,- p)j h,_)+h p ((p 1+1- p1)/h)] if i # 0 and i # m

D, = ao (( pi-po ho) if i=0

aLj ((prn -pj )/h,_j)ifi=m

! fi V + fj )[)F'(p ,-p 1 f_)+f ((Plj-p)f f j0andjn

E0 = p -p fO) if j=0

§fj (pi~ n- )/fn-_ if j= n

wherein h , 1 -zT(i=0,1,...,m-1) and f =y -yj (j=0,1 n -1), and

a1 ,§j e [0,1] are shape-related parameters and here a and #,j are assigned to 0.8.

Bo,(u)=(1-u)3 B (v)=(1-v)' 2 B1, (u= 3 u(1-u)2 (V) BndB 1, (v)= 3v(1-v) (U)( 3 B 2 (v 3V 2 ( V)

[B3,, (U) =u B3, ,(V)= v'

The fitting efficacy using generalized tensor product based on Bezier surface is better than the tensor product Bezier surface, which makes the fitting surface closer to the original surface of the studied leaf.

Obtaining the Bezier surface through generalized tensor product and selecting 16 control points in an arrangement of 4*4 on the leaf surface, where every 4 points are selected from the 4 adjacent scanning lines in turn recorded as

0O,0 I P0,1 I P0,21 PO9,3

P1,0O3P1,1 'P1,2 p 1 ,3 p, = IO . P I P1,2Bringing these 16 points into formula 4, so as to P2 ,0, P 2 ,1, P 2 ,2 , P 2 ,3

P3 ,0 P3 ,1 I P 3 ,2 P 3 ,3

obtain the fitted surface. At last, repeating this process until all scanning points on all the scanning lines are processed to achieve the generation of the fitted leaf surface.

(3) Deformation of the fitted leaf surface under stress

(3.1) Using Q to represent the state of the leaf before deformation, and p E-?

is a certain point on the leaf, wherein, Q c R3 and p(p,, p,, p). Using F to represent

the state of the leaf after deformation F R3 and p'e F is a certain point on the

deformed leaf p'(p', p', p) , corresponding to p . The displacement u is used to

describe the transformation from Q to F, that is u:Q ->F, p' = p+ u(p), wherein, if

pp2 E Q , p , p2 F , the spatial displacement d between p, p 2 before deformation is

d =p 2 - p ; in the same way, after deformation, the spatial displacement d between

p i,p ' isd' =p' - p' . From the above formula, we can get that:

d' = p 2 +u(p2)-p1 -u(p)=d+u(p,+d)-u(p)=(I+Vu)d (5)

Wherein: I is identity matrix, I+Vu is deformation gradient, denoted as F=I+Vu . The decomposition of vector u into the three mutually orthogonal au au au ax ay az

components (u,v, w) and Vu a z .y ax ay az aw aw aw ax ay az

The deformation can be represented by the change between d, d', that is

|d'|2 -d|2 -d T d = d(FT F-I)d = d' T 'd'

According to theory of elasticity, Green strain EE R"' can be obtained from the above formula as follows:

E= (FTF- I)= (V+VuT+VuTVu) (6) 2 2

Besides, from the above formula, it can be noted that E is a symmetric matrix.

For the leaf model which is composed of tetrahedron as the basic element, linear interpolation function is used to represent the displacement vector of each node, that is u= (u,v, w). For a tetrahedron composed of four nodes, each node has three degrees of freedom, along the x,y, z directions of the node coordinates. Let the four

vertices of tetrahedron be represented respectively as 1 ",) , ,) (titit (t , t ,,)

IP(tX,tIk,tk,) and P e t d(t t ,), and each vertex changes in the x , y and z

!u =C +C U 1I x+1yC4 12 X±C13 Y±C4z v= C2 1 +C2 2x+C2 3y+C3 4z (7)

w= C3 1 +C 3 2 x+C33 y+C3 4z

From the above formula, the displacement field functions of four vertices are obtained as follows:

1 S.= ,(a,+b~x+c~y+dz) 6V' 1 Sj= 6V -(a+b x+cy+djz) (8) 1 Sk= (ak+bkx+cky+dkz) 6V' 1 S,= ,(a,+blx+c,y+dz) 6V'

Si is the interpolation function or shape function of the tetrahedron, Vis the

1 t1j t' t, 11 t t. t volume of tetrahedral element, V'= - 1 ; a,,b,,c,, d 1 ,---d are 61 tkx tky tkz

1 t t t

respectively as follows:

tjX t t , t t, t 1 t~ tx t 1,y a =tkx t t , b=-1t t c =t 1 t ,di=-t t 1 t1' t' tl~ 1 tI tl tl, 1 tl tl tI 1

As for aj,bj,cj, dj, a k ... and d,, they can be obtained by alternating subscripts

i, j,k,according to the right hand rule, for example,

k,x kyk ,x ,y ,,z i=x iy t i,z

aj = ,t,x t'y tl,z ak =~ i yt' i,z a, = j,x j'x tj'x t.'I,x ti.l,y t, I,z t.j' j,x t.j,y U' J,zkx tk~x tk~x kx tk~x kx

From the above formula, we can calculate the displacement u:

t .

u =ivJ=[JISiIS IISk ] (9) tk t

Wherein, ii,i, and i are 3x1 displacement vectors of each vertex of

tetrahedron, I is identity matrix and S is 3x 3 displacement matrix for each vertex.

Hence, the transformation for the current displacement u is retrieved.

(3.2) Leaf deformation

(3.2.1) Model representation

The deformation model is represented as (V,G,P) , wherein

V={JVh=1,2,-..,n} is tetrahedral set, which constitutes the leaf surface.

P={IER 3 1i 4} is the coordinate set of each vertex of a tetrahedron and

G={(1I,)),IP e V,i# j} is a set of edge of tetrahedron. As for each vertexi,

N, ={P I(IP,P)e G} is used to represent the set of the adjacent vertices ofJ). At the

mass mi of the point I can be calculated as follows:

1 4 VjcN,

Wherein, p1 , a constant, is tetrahedral density; Vf is the volume of the j th

tetrahedron with a vertex of 1.

(3.2.2) Definition of potential energy in deformation model

Calculating the elastic potential energy of the deformation model based on following formula:

T(P)= T(P) +T, (P) (10)

Wherein, V,(P) represents the potential energy that makes the deformed

model return to equilibrium or initial state, namely recovery potential energy; V,(P) represents the energy required to maintain the deformation of the model, namely volume potential energy;

(3.2.3) Definition of recovery potential energy

Calculating the recovery potential energy based on following formula:

T(p)= ||L[(P)-R (P)|2 (11)

Wherein, Li is differential operator of vertex 1' and r, is differential

coordinates of vertex 1] at equilibrium or static state. R(P) is rotation matrix of 13

,5, = A,(wJ)=w w1,(f(P)-f(17)) (12) VjE Ni

in Wherein, w.= l1cot(O,) is the weight of the side length, lk is the side

length of G(P, PI), opposite side of G( , Pj), 0, is the dihedral angle of opposite side

1 G(P, ,),1 w, = is the normalized weight of the centre vertex P, V., is the (Vdual) i

of the tetrahedron in Figure 7. f(P)= fD(P) is a scalar field function, (-) is a

-s n piecewise linear basis function, (1_)= ," is the basis function at the central 3V' vertex P , V' is the volume of tetrahedron, s is the area of the relative bottom triangle

PPPj and hp is the normal direction of the bottom triangle pointing to the outside of

the tetrahedron. Schematic representation is shown in Figure 7.

Formula (12) can be expressed as a matrix-form.

-5(5(o ,y, ,) )= P(ti, ,ti, )=M-I,(17,PPi) (13)

Wherein, ,ta and(,, , are vector representations of the

three components of the initial tetrahedral vertex coordinates and Laplacian 1 coordinates, respectively, M is the diagonal matrix of the volume weightM,

L, is a symmetric sparse linear matrix of grid weights

-Y wi, i-j ke Ni

(L,),j= wj, if- iGNJ,

) 0, otherwise

Then the Laplacian matrix can be expressed as L = M-1L.

According to the properties of Laplacian coordinates, the Laplacian coordinates of volume meshes are invariant to the translation and are variant to the scaling and rotation.

In order to calculate the recovery potential energy, the rotation matrixR,(P)

of each vertex ]is also calculated as follows:

F heftreM (pem"re Pelore)

= airM (palter- paler

Wherein, pbere and ar represent the coordinate of the point P before and

after deformation, respectively. Then the linear transformation matrix A of vertex

can be obtained by finding the minimum of the following equation

|F"" - A-F, c 2

jcN(i)

When the minimum value of the above formula is obtained, we can get the following results:

A, A|"rAr"

Wherein, A|" = (Fl"ter (1ber") ) and A,"= ( 1bfor (F b°")y-I je Ni e

Since Arr is a symmetric matrix without any rotation deformation, the rotation

matrix should be obtained by polar decomposition of the matrix Ap/. At the same time,

setting the translation matrix of the model as T, so that R,(TP)= TR,(P), that is,

LTP - RI(TP)d, = T(LP - Rd,)

It can be seen from the above that V,(TP)= V(P) , meaning V,(P) is

translation invariant

(3.3) Model deformation

The deformation simulation of the model is controlled by the following Euler Laplace motion equation:

MP+DP+ T(P) (14)

Wherein, M is the mass matrix of the system, D is the damping matrix of the

system, T(P) is the potential energy of the system, f,, is the external force on the

system, is the internal force of the system. At the same time, the Jacobian C9P

a 2 T (P) matrix of the internal force of the system is , namely, the stiffness matrix of the a8aP system K. For the damping matrix D of the system, Rayleigh damping is adopted, so

D = aM +8K

Wherein, KE R(3nx 3n) is the stiffness matrix of the system, M is the

mass matrix of the system and a and # are constant. In the solution, setting§3=0 and

In general, dividing the volume model into many tetrahedral models, and applying external force to the model, so as to obtain the deformed leaf surface by solving the equation. According to the specific steps of our method, the deformation of any leaf surface under mechanical loading can be implemented, and the results can be validated by real trails with the selected leaf surface.