• Procedure DXFREAD
• Procedure SOLAO
  • Derivation of the procedure SOLAO
• Procedure PHASE0
• Procedure FOURIER
  • Derivation of the procedure FOURIER
• Procedure KAF
  • Derivation of the procedure KAF
• Procedure COREL

Doc. RNDr. Stanislav Barton, CSc.
Mendel university of Agriculture and Forestry in Brno
Agronomical faculty
Dept. of Automobile transport and principles of technics
Zemedelska 1; 61300 Brno
++ 05 4513 2127
barton@mendelu.cz

Identification and Recognition of Leaves and General 2-D Objects

Part 1:  - Transformation of DXF file into a parametric function

Theory:   The outline of an arbitrary complicated 2D object can be approximated by a polygon consisting of a finite number of line segments. A scanned image can be converted into line segments by the program Adobe Streamline or similar, and the resulting outline can be saved as DXF file = Data Exchange Format. The [x, y]  coordinates of polygon vertices can be found in this file. The simplest approximation of such a polygon can be done by a couple of parametrically dependent functions .The main question is the proper selection of such a  parameter an ordinary number of the verticies can be an example of the parameter. But for other manipulations this parameter will be very important. If the parameter is a length from the initial point on the perimeter, measured along the perimeter normalized to 2 , Fourier series can be used. As is shown in the section Fourier  and results from the orthogonality of goniometric functions, it is possible to compute Fourier coefficients independently, not only for x  and y  functions, but these coefficients are independent. The objects outline is converted into a couple of smooth continuous functions defined on the interval <0, 2 > . These functions can be described by the lists of coefficients and can be used to compute the similarity of two different objects, described by the coefficients of the linear correlation, for example. These coefficients will be computed in the first part.

Reading of the procedures

Procedure DXFREAD

The procedure reads the input DXF file and looks for x  and y  coordinates of the polygon describing the perimeter of the given object. Coordinates are read step by step as they are ordered in the input file. Coordinates are saved in the variables X  and Y  of type list.
Procedure's input is the name of the file, without extension, containing DXF data.  The output is two rlists of vertices, one for x coordinates and one for y.

Procedure SOLAO

Procedure orients the object in the coordinate system. The line of orientation  is chosen to minimize the moment of inertia (with respect to the axis of rotation in the xy  plane). So the parameters k  and q  minimize  , where k  and q  are parameters describing the general line , ds  = element of the curve (perimeter) and r  = distance of the element of the curve ds  from the line. The procedure finds two such lines, mutually perpendicular. One leads to a minimum of J(k,q)  the second one to a maximum. The whole object is shifted so that the intersection of these lines will be identical with the origin of the coordinate system [0,0] , and rotated so that line with the minimum moment of inertia will coincide with the x  axis. This leads to identical orientation of all objects. Because the object (polygon) is created from line segments, it is possible to use the following law to simplify computation
 => ,
where N  = count of found vertices, thus N-1  line segments, p = parameter describing line segment
The Maple functions makeproc  and optimize  were used to derive this procedure.
Procedures inputs are lists X  and Y  of coordinates of the vertices.  The procedure changes the values of X  and Y  inplace.

Derivation of the procedure SOLAO

>

restart;

>	with(linalg):  open library Linear Algebra

>	with(codegen,optimize,makeproc,cost);  Open library for code generation, automatic optimization and computer resource computation

>	Prx:=x=Kx*p+Qx; Pry:=y=Ky*p+Qy;  Parametric description of the line segment creating part of the perimeter. K = [Kx, Ky] is a vector corresponding to the line segment and  p = parameter , 0 <= p <= 1.

>	solx:=solve({subs(x=X[i],p=0,Prx),subs(x=X[i+1],p=1,Prx)},{Kx,Qx}); soly:=solve({subs(y=Y[i],p=0,Pry),subs(y=Y[i+1],p=1,Pry)},{Ky,Qy}); Computation of the components of the vector  K

>	S2 := (q-y+k*x)^2/(1+k^2);  Squared distance of the point [x,y]  from the line

>	S2s:=subs(Prx,Pry,S2);    [ x,y ]  moves along Prx, Pry  = one of line segments describing object, S2s  = the squared distance from the point on the line segment to the general line. Point position is given by the parameter p .

>	S2i:=int(S2s*L[i],p=0..1); Computation of the moment of inertia for one line segment

>	S2ii:=simplify(subs(solx,soly,S2i)); Kx and  Ky are described by  [x,y] coordinates of the initial and final point of the line segment. Substitution.

>	S2S1:=map(u->simplify(sum(u,i=1..N-1)),expand(S2ii));  Summation of individual integrals.

>	sigma:=[selectfun(S2S1,sum)[]];  Separation of individual independent summations

>	n:=nops(sigma);

>	Su:={seq(sigma[j]=Sigma||j,j=1..n)}; Substitution of summations leads to a simplified overview.

>	Bs:=map(u->rhs(u)=lhs(u),Su); Back substitution

>	S2S2:=simplify(subs(Su,S2S1)); Substituted integral

>	Sk:=numer(diff(S2S2,k)); Sq:=numer(diff(S2S2,q)); Normal equations. Their solution returns values of k  and q  minimizing S2S2  = moment of inertia of the object

>	sol:=allvalues(solve({Sk,Sq},{k,q})): Computed variables of k  and q . Output is not presented, because it is too long.

>	SOLA:=vector([subs(sol[1],[k,q])[],subs(sol[2],[k,q])[]]): Solution is converted into a vector.

>	Cost1:=cost(SOLA); Consumption of computer memory for the original solution = SOLA

>	SOLAp:=makeproc(SOLA,[seq(Sigma||j,j=1..n)]): Conversion of the solution into a procedure

>	SOLAO:=optimize(SOLAp): Procedure optimization

>	Cost2:=cost(SOLAo); Consumption of computer memory for the optimized solution = SOLAo

>	Cost1-Cost2; savings

>

Sense of the procedure

>	read `DXFREAD.pro`:

>	DXFREAD(`Oak1.dxf`); Input of the coordinates of the vertices

"Converting DXF file into lists of coordinates of vertices"

"102 vertices found"

>	P1:=plot(zip((x,y)->[x,y],X,Y),color=blue): Before centering and orientation

>	SOLAO(X,Y); application of the procedure

"Looking for center and axis of symmetry"

"Shift =[-3.846728402, -2.137677566]"

"Rotation = 2.201621101*degrees"

>	P2:=plot(zip((x,y)->[x,y],X,Y),color=red,thickness=2): Centered and orientated plot

>	plots[display]({P1,P2},scaling=constrained);

Procedure PHASE0

Parametric functions are used to describe a given object, therefore it is advantageous to set 0  to the value of the parameter so that the corresponding point will have always the same positional = phase angle. This procedure is searching for the point situated on the x  axis, and its x  coordinate is maximal. The procedure later reorders the lists of coordinates in that way, so this point will be the first one, but the shape of the object will be preserved.
Inputs and outputs of the procedure are lists of the coordinates of the vertices..

Sense of the procedure

>	P1:=plot(zip((x,y)->[x,y],X,Y)):

>	A1:=plots[arrow]([X[1],Y[1]],color=blue):  Displays original initial point

>	read `PHASE0.pro`:

>	PHASE0(X,Y);

"Setting 0 to initial phase"

"Zero point = [3.537382891, 0]"

>	A2:=plots[arrow]([X[1],Y[1]],color=red):  New initial point lays on the x  axis

>	plots[display]({P1,A1,A2},scaling=constrained);

Procedure FOURIER

Procedure Fourier  computes coefficients of the Fourier series substituting piecewise defined functions x(p)  and y(p) , 0 <= p <= L ,( L  = distance of the point [x(p),y(p)] , from the initial point, [x(0),y(0)]  measured along the perimeter). Fourier series will create smooth continuous functions defined on the interval <0,2 > .
The procedures entries are the lists of the coordinates of the vertices, the third parameter indicates the number of the members of the series, outputs are the lists of the Fourier coefficients

Derivation of the procedure FOURIER

>

restart;

>	q:=(f(p)-sum(c[i]*cos(p[i])+s[i]*sin(p[i]),i=0..3))^2; Squared residuals

>	q2:=expand(q):

>	q3:=eval(subs(seq(p[i]=i*p,i=0..3),q2)); The first four members will be used, for the sake of brevity.

>	q4:=map(u->simplify(int(u,p=0..2*Pi)),q3); Integral of the squared residuals

>	Var:={c[0],seq({c[i],s[i]}[],i=1..3)}; The set of the unknown Fourier coefficients

>	Eq:=map(u->diff(q4,u)=0,Var); The normal equations.

>	sol:=solve(Eq,Var); Coefficients formulas, can be simply generalized

>	Sol:={C0=1/2*int(f(p),p=0..2*Pi)/Pi,       Cn=int(f(p)*cos(n*p),p=0..2*Pi)/Pi,       Sn=int(f(p)*sin(n*p),p=0..2*Pi)/Pi}; Generalized solution

>	sol1:=subs('p=0..2*Pi'='p=X[i]..X[i+1]',f(p)=K[i]*p+Q[i],Sol); Generalization for the piecewise defined function one interval of the validity, all partial functions are line segments

>	sol2:=value(sol1); sol1  can be integrated

>	assign(sol2);

>	C0:=simplify(sum(C0,i=1..N-1)); and summarized for all line segments

>	Cn:=simplify(sum(Cn,i=1..N-1));

>	Sn:=simplify(sum(Sn,i=1..N-1));

>	kq:=solve({Y[i]=K[i]*X[i]+Q[i],Y[i+1]=K[i]*X[i+1]+Q[i]},{K[i],Q[i]});  and  can be computed from the coordinates of the initial and final point of the line segments.

Sense of the procedure

>	P1:=plot(zip((x,y)->[x,y],X,Y),color=blue): The line plot

>	FOURIER(X,Y,40);  Coefficients generation

"Computing 40 Fouriere's coefficients"

X0 = -.1328943774e-7

Cx[1] = 2.900225949, Sx[1] = .1419823020

Cx[40] = .4754955872e-3, Sx[40] = .3202506275e-2

Y0 = .3724225667e-7

Cy[1] = .6723654211e-1, Sy[1] = -1.715950597

Cy[40] = .1312131189e-2, Sy[40] = .1486127291e-2

>	Fx:=sum(Cx[i]*cos(i*x)+Sx[i]*sin(i*x),i=1..40):  Fourier series describing x(p)

>	Fy:=sum(Cy[i]*cos(i*x)+Sy[i]*sin(i*x),i=1..40):  Fourier series describing y(p)

>	P2:=plot([Fx,Fy,x=0..2*Pi],color=red,numpoints=600):  parametric plot

>	plots[display]({P1,P2},scaling=constrained,title="Graphical comaprison");

>	L:=[0.]: Computation of the paremeters coressponding to the vertices

>	for j from 2 to N do;   L:=[L[],L[j-1]+sqrt((X[j-1]-X[j])^2+(Y[j-1]-Y[j])^2)]; od:

>	L:=map(u->evalf(2*Pi*u/L[N]),L):

>	plot([zip((x,y)->[x,y],L,X),zip((x,y)->[x,y],L,Y),Fx,Fy],x=0..2*Pi,color=[blue,red,navy,brown],title="Oak (like a sine and cosine) functions",legend=["x lin. segm.","x(p) Fouriere","y lin. segm.","x(p) Fouriere"]);

>	R:=zip((x,y)->sqrt(x^2+y^2),X,Y):  Distances of the vertices

>	F:=map(u->evalf(subs(x=u,sqrt(Fx^2+Fy^2))),L):  Distances of the point corresponding with vertices

>	RE:=zip((x,y)->(1-y/x)*100,R,F):  The relative errors

>	Var:=sqrt(sum(RE[i]^2,i=1..N)/N):  The variance of the relative errors

>	plot([zip((x,y)->[x,y],L,RE),[[0,Var],[2*Pi,Var]],[[0,-Var],[2*Pi,-Var]]],color=[red,blue,blue],title="The relative error of the distance from the origin and its variance",labels=["Perimeter's parameter","%"]);

Increasing the number of the Fourier coefficients decreases the error of the approximation.

Processing

`"Converting DXF file into lists of coordinates of vertices"`

`"574 vertices found"`

`"Looking for center and axis of symmetry"`

`"Shift =[-12.21275220, -9.483011334]"`

`"Rotation = 3.418483572*degrees"`

`"Setting 0 to initial phase"`

`"Zero point = [9.226803319, 0]"`

`"Computing 50 Fouriere's coefficients"`

`X0 = -.2212253708e-6`

`Cx[1] = 8.672002591, Sx[1] = -1.017347828`

`Cx[50] = -.1191587246e-1, Sx[50] = -.1287842212e-1`

`Y0 = .1917817064e-6`

`Cy[1] = -.8152692550, Sy[1] = -6.792215980`

`Cy[50] = .2360530589e-1, Sy[50] = -.1289477692e-1`

Part 2:  - Computation of the transformation parameters, (mutual rotation, shift, rescaling and phase shift) of two different objects using Least Square Method and computing of the Coefficient of the Linear Correlation,.

Theory:  The Least Square Method is used to determine the parameters of the transformation. These parameters are:

      K  = rescaling of the first object

     Dx, Dy  = shift of the first object in x  and y

     a   = rotation of the second object
     f  = phase shift of the second object,

which will minimize .

With respect to ortogonality of individual terms of Fourier series creating Fx, Fy, Gx  and Gy,  it is possible to use a lot of simplifications leading to an exact analytical value of the integral mentioned above. The analytical solution of   a, Dx, Dy  a K  can be determined from this value. The value of f  can be found by its direct minimization.

If the exact values of transformation are known, it is possible to create new transformed functions FX, FY, GX  and GY , where:
      ,

      .

Similarity of the transformed objects will be described by the coefficient of the linear correlation computed below.

      ,      ,     

      ,      ,

      ,      ,

      ,

Reading of the procedures

Procedure KAF

Procedure KAF computes parameters of the transformation of two objects, resulting in a maximal similarity of both objects. The computation will be finished when the step size of the variable f , (see main page) is lower than a desired correctness = the third parameter.. The fourth, optional parameter is the number of lines connecting the corresponding points on perimeters of both transformed objects. If value of this parameter is 0 , both transformed objects are plotted without lines. If this parameter is missing, the  procedure will not create graphical output. The procedure creates new parametric dependent functions describing both transformed objects, saved as new global variables x1, y1, x2  and y2 .

Entries of the procedure are filenames containing coefficients of the Fourier equations describing the first and second object. The third parameter is the desired correctness step of f . The four optional parameter controls graphical output.

Derivation of the procedure KAF

>

restart;

Simplified functions describing x and y  coordinates of the first (Fx, Fy)  and second (Gx, Gy)  objects

>	Fx:=CFx[i]*cos(i*x)+CFx[j]*cos(j*x)+SFx[i]*sin(i*x)+SFx[j]*sin(j*x);

>	Fy:=CFy[i]*cos(i*x)+CFy[j]*cos(j*x)+SFy[i]*sin(i*x)+SFy[j]*sin(j*x);

>	Gx:=CGx[i]*cos(i*x)+CGx[j]*cos(j*x)+SGx[i]*sin(i*x)+SGx[j]*sin(j*x);

>	Gy:=CGy[i]*cos(i*x)+CGy[j]*cos(j*x)+SGy[i]*sin(i*x)+SGy[j]*sin(j*x);

The first object will be rotated through angle  a and shifted by Dx and  Dy.

>	FX:=Dx+expand(Fx*cos(a)+Fy*sin(a)); FY:=Dy+expand(-Fx*sin(a)+Fy*cos(a));

The phase of the second object will be shifted by f  and object will be rescaled by a factor of K  

>	GX:=expand(subs(x=x+f,K*Gx));

>	GY:=expand(subs(x=x+f,K*Gy));

Squared distance of the corresponding points of the first and second objects = squared residual. Output is too long, so it is not presented.

>	R2:=expand((FX-GX)^2+(FY-GY)^2):

>	assume(i,integer); assume(j,integer); Arguments of Fourier series contains always integer multiple of 2 .

>	R21:=map(u->simplify(Int(u,x=0..2*Pi)),R2): Integral of residuals. Values of Dx, Dy, K, a  and f , minimizing its value are to be found. Output is too long, it contains 262  integrals.

>	sigma:=selectfun(R21,Int); But only 15  of them are individual

>	Su:=[seq(sigma[k]=Sigma[k],k=1..nops(sigma))]; Substitution for the integrals

>	Bs:=map(u->rhs(u)=value(lhs(u)),Su); Back substitution, only 5 are not zero

>	R22:=subs(Su,R21): Integrals are substituted, still too long

>	R23:=subs(Bs,R22);  Values of the integrals are substituted. Simple and short.

>	W1:=selectfun(R23,`^`); Squared terms are extracted

>	W11:=select(has,W1,[i,j]);  Indexed variables are extracted

>	W12:=map(u->factor(u),collect(select(has,R23,W11),K));  , These terms are removed, and converted into a proper summation.

>	R24:=R23-expand(W12)+Pi*sum(C1x[k]^2+C1y[k]^2+S1x[k]^2+S1y[k]^2,k=1..N)             +Pi*K^2*sum(C2x[k]^2+C2y[k]^2+S2x[k]^2+S2y[k]^2,k=1..N);  This summation will replace old terms in R23 .

>	W1:=select(has,R24,CFx);  Simplification of terms containing  CFx

>	W12:=factor(select(has,W1,CGx));

>	W13:=-2*cos(a)*Pi*K*sum(C1x[k]*C2x[k]*cos(k*f),k=1..N);

>	W2:=W1-expand(W12)+W13;

>	W12:=factor(select(has,W2,CGy));

>	W13:=2*sin(a)*Pi*K*sum(C1x[k]*C2y[k]*cos(k*f),k=1..N);

>	W3:=W2-expand(W12)+W13;

>	W12:=factor(select(has,W3,SGy));

>	W13:=2*sin(a)*Pi*K*sum(C1x[k]*S2y[k]*sin(k*f),k=1..N);

>	W4:=W3-expand(W12)+W13;

>	W12:=factor(select(has,W4,SGx));

>	W13:=-2*cos(a)*Pi*K*sum(C1x[k]*S2x[k]*sin(k*f),k=1..N);

>	W5:=W4-expand(W12)+W13;

>	R25:=R24-W1+W5;

>	W1:=select(has,R25,SFx): Simplification of terms containing SFx, ..... Because it is the same method as before, outputs will not be presented.

>	W12:=factor(select(has,W1,CGy)):

>	W13:=-2*sin(a)*Pi*K*sum(S1x[k]*C2y[k]*sin(k*f),k=1..N):

>	W2:=W1-expand(W12)+W13:

>	W12:=factor(select(has,W2,SGy)):

>	W13:=2*sin(a)*Pi*K*sum(S1x[k]*S2y[k]*cos(k*f),k=1..N):

>	W3:=W2-expand(W12)+W13:

>	W12:=factor(select(has,W3,CGx)):

>	W13:=2*cos(a)*Pi*K*sum(S1x[k]*C2x[k]*sin(k*f),k=1..N):

>	W4:=W3-expand(W12)+W13:

>	W12:=factor(select(has,W4,SGx)):

>	W13:=-2*cos(a)*Pi*K*sum(S1x[k]*S2x[k]*cos(k*f),k=1..N):

>	W5:=W4-expand(W12)+W13:

>	R26:=R25-W1+W5:

>	W1:=select(has,R26,CFy):  Atd. s CFy.

>	W12:=factor(select(has,W1,CGx)):

>	W13:=-2*sin(a)*Pi*K*sum(C1y[k]*C2x[k]*cos(k*f),k=1..N):

>	W2:=W1-expand(W12)+W13:

>	W12:=factor(select(has,W2,CGy)):

>	W13:=-2*cos(a)*Pi*K*sum(C1y[k]*C2y[k]*cos(k*f),k=1..N):

>	W3:=W2-expand(W12)+W13:

>	W12:=factor(select(has,W3,SGx)):

>	W13:=-2*sin(a)*Pi*K*sum(C1y[k]*S2x[k]*sin(k*f),k=1..N):

>	W4:=W3-expand(W12)+W13:

>	W12:=factor(select(has,W4,SGy)):

>	W13:=-2*cos(a)*Pi*K*sum(C1y[k]*C2y[k]*sin(k*f),k=1..N):

>	W5:=W4-expand(W12)+W13:

>	R27:=R26-W1+W5:

>	W1:=select(has,R27,SFy):  Atd. s SFy.

>	W12:=factor(select(has,W1,CGy)):

>	W13:=2*cos(a)*Pi*K*sum(S1y[k]*C2y[k]*sin(k*f),k=1..N):

>	W2:=W1-expand(W12)+W13:

>	W12:=factor(select(has,W2,CGx)):

>	W13:=2*sin(a)*Pi*K*sum(S1y[k]*C2x[k]*sin(k*f),k=1..N):

>	W3:=W2-expand(W12)+W13:

>	W12:=factor(select(has,W3,SGx)):

>	W13:=-2*sin(a)*Pi*K*sum(S1y[k]*S2x[k]*cos(k*f),k=1..N):

>	W4:=W3-expand(W12)+W13:

>	W12:=factor(select(has,W4,SGy)):

>	W13:=-2*cos(a)*Pi*K*sum(S1y[k]*S2y[k]*cos(k*f),k=1..N):

>	W5:=W4-expand(W12)+W13:

>	R28:=R27-W1+W5:

>	sigma:=selectfun(R28,sum);  Collection of summations of the final term R28

>	Su:=[seq(sigma[l]=Sigma[l],l=1..nops(sigma))];  Substitution

>	Bs:=map(u->rhs(u)=lhs(u),Su);  and back substitution

>	R1:=subs(Su,R28);  Integral of the squared residuals is quite simple now.

>	R1x:=diff(R1,Dx); R1y:=diff(R1,Dy); Normal equations for Dx  and Dy .

>	Dx:=0: Dy:=0:  No shift. Notice that absolute terms X0  and Y0  were not used to create Fx, Fy, Gx  and Gy .

>	R1K:=collect(simplify(diff(R1,K)/2/Pi),[sin,cos]);  Normal equation for K

>	Ksol:=K=collect(solve(R1K,K),[sin,cos]);  Solution for K

>	R1a:=collect(simplify(diff(R1,a)/2/Pi),[sin,cos]);  Normal equation for a

>	Asol:=a=solve(R1a,a);  Solution for a

>	R1c:=collect(simplify(R1/Pi),[sin,cos,K]);  It is necessary to find the value of f  minimizing this term. It is simpler to solve its normal equation, because all summations contains f .

Procedure COREL

The procedure computes coefficients of linear correlation of two functions, R1(p)  and R2(p) . These functions are created from two pairs of input functions:  and , see main file for details. The Simpsons rule is used for numeric integration. As the coefficient of the linear correlation approaches 1, the similarity  increases. Zero value indicates no similarity of objects.

Entries of the procedure are four functions x1(p), y1(p), x2(p)  and y2(p)  describing parameter dependence of x  and y  coordinate of the first and second object. The fifth parameter controls the relative accuracy of the numeric integration.

>

Processing

"f = .2,  df = .1,  MIN = 4.74937257675374613"

"f = .16784667968750,  df = -.1220703125e-4,  MIN = 4.69629115004216838"

"_____________________________________________"

"Rotation of Maple1.sav = 8.0004480430310438605*degrees"

"Phase shift of Maple2.sav = 9.6176057627489623335*degrees"

"Scale of Maple2.sav = 2.7941995986574304523"

P1 = 7.912570760

V1 = 2.837610148

P2 = 7.785024297

V2 = 2.514006750

C12 = 1.944717629

Conclusion:

The presented algorithm can be used as an instrument to determine rate of similarity of simple two-dimensional objects. Its advantage is simplicity and possibility of conversion into source files of other programming languages. These procedures were designed for minimization of computer memory consumption, and so Maple special functions were not used, because these functions cannot be used out of Maple.

As examples we used a couple of leaf types: Maple, Oak, and Pothos. Linear coefficients describing similarity of individual leaves are shown in the following table.

It is evident that leaves of the same variety are the most similar each to other. Similarity depends on leaf segmentation, of course

The algorithms were derived in the program Maple 7, running on the computer Autocont, PIV, 512 MB RAM

 5. 2. 2002