Reconstructing Surfaces from Fundamental Form Coefficients

by Tali Avigad
Department of Mathematics and Computer Science,
Bar-Ilan University, Ramat-Gan, Israel
Email : avigad@macs.biu.ac.il

Reconstruct_patch is a procedure for calculating a numerical approximation of a surface given its fundamental form coefficients(1) .

The procedure is based on the fundamental theorem of surfaces(2), which says that if certain compatibility conditions are satisfied, then it is possible to reconstruct a surface that is unique up to a Euclidean motion. The procedure first checks the validity of the compatibility conditions and a suitable message is printed. Afterwards the procedure calculates a numerical approximation by solving the Gauss-Weingarten(3) equations using the dsolve[dverk78] function.

The input:

reconstruct_patch( E , F , G , L , M , N , u , v , min_u , min_v , max_u , max_v , no_of_steps_u , no_of_steps_v)

E , F , G - functions of u and v which represent the first fundamental form coefficients.
L, M, N - functions of u and v which represent the second fundamental form coefficients.
min_u , min_v , max_u , max_v - the size of the grid.
no_of_steps_u - the number of steps in u-direction .
no_of_steps_v - the number of steps in v-direction .

The output:

A message regarding the validity of the compatibility conditions.
A plot of the approximated surface at the grid points.

Mathematical background

# ***********************************************************************

# PROC. reconstruct_patch approximate a patch using the sequence of column-rows

# ***********************************************************************

> reconstruct_patch := proc(E,F,G,L,M,N,u,v,min_u,min_v,max_u,max_v,no_of_steps_u,no_of_steps_v)

> local a11, a12, a22, b1_1, b2_1, b1_2, b2_2, g1_11, g2_11, g1_12, g2_12, g1_22, g2_22, X, sys_u, sys_v, init1, init2, init3, initu1, initu2, initu3, initv1, initv2, initv3, res_u_sys1, res_u_sys2, res_u_sys3, res_v_sys1, res_v_sys2, res_v_sys3, i, j, const_v, res1, res2, res3;

> global delta, Eu, Ev, Fu, Fv, Gu, Gv, E0, F0, G0, alpha11, alpha12, alpha22, beta1_1, beta2_1, beta1_2, beta2_2,gamma1_11, gamma2_11 , gamma1_12 , gamma2_12 , gamma1_22 , gamma2_22,gamma2_22u,gamma2_12v, gamma1_22u, gamma1_12v, step_size_u, step_size_v, find_init_val, check_compatibility_conditions, e, f, g;

# ***************************************************************

# variable initialization

# ***************************************************************

> e:=E;

> f:=F;

> g:=G;

> alpha11:=L;

> alpha12:=M;

> alpha22:=N;

> delta:= simplify(E*G-F^2);

> if (delta=0) then print('EG-FF_is_zero'); fi;

> Eu := diff(E,u);

> Ev := diff(E,v);

> Fu := diff(F,u);

> Fv := diff(F,v);

> Gu := diff(G,u);

> Gv := diff(G,v);

> beta1_1:=simplify((M*F-L*G)/delta);

> beta1_2:=simplify((N*F-M*G)/delta);

> beta2_1:=simplify((L*F-M*E)/delta);

> beta2_2:=simplify((M*F-N*E)/delta);

> gamma1_11:= simplify((G*Eu-2*F*Fu+F*Ev)/(2*delta));

> gamma1_12:= simplify((G*Ev-F*Gu)/(2*delta));

> gamma1_22:= simplify((2*G*Fv-G*Gu-F*Gv)/(2*delta));

> gamma2_11:= simplify((2*E*Fu-E*Ev-F*Eu)/(2*delta));

> gamma2_12:= simplify((E*Gu-F*Ev)/(2*delta));

> gamma2_22:= simplify((E*Gv-2*F*Fv+F*Gu)/(2*delta));

> gamma2_22u:= diff(gamma2_22,u);

> gamma2_12v:= diff(gamma2_12,v);

> gamma1_22u:= diff(gamma1_22,u);

> gamma1_12v:= diff(gamma1_12,v);

#=======================================================

# Proc check_compatibility_conditions()

#=======================================================

> check_compatibility_conditions := proc()

> local cond4, cond5, cond6, flag1, flag2, flag3, flag4, flag5, flag6, result;

> result:=1;

> interface(showassumed=0);

> assume(u,real); assume(v,real);

> cond4:=simplify((alpha11*gamma1_12+alpha12*(gamma2_12-gamma1_11)-alpha22*gamma2_11)-diff(alpha11,v)+diff(alpha12,u));

> cond5:=simplify((alpha11*gamma1_22+alpha12*(gamma2_22-gamma1_12)-alpha22*gamma2_12)-diff(alpha12,v)+diff(alpha22,u));

> cond6:=simplify((f*(gamma2_22u-gamma2_12v+gamma1_22*gamma2_11-gamma1_12*gamma2_12) +e*(gamma1_22u-gamma1_12v+gamma1_22*gamma1_11+gamma2_22*gamma1_12-gamma1_12*gamma1_12-gamma2_12*gamma1_22))-alpha11*alpha22+alpha12^2);

> flag4:=is(numer(cond4)=0);

> flag5:=is(numer(cond5)=0);

> flag6:=is(numer(cond6)=0);

> if (flag4 and flag5 and flag6)=false then print('The_compatibility_conditions_are_not_satisfied'); result:=0; elif (flag4=FAIL or flag5=FAIL or flag6=FAIL) then print('The_compatibility_conditions_may_not_be_satisfied'); fi;

> RETURN(result);

> end;

#============================================================================

# Proc find_init_val() This is an annoying procedure to put the output of dsolve in the right order!

#============================================================================

> find_init_val := proc(lst)

> local i,surf_val,Xu_val,Xv_val,N_val;

> for i from 1 to nops(lst) do

> if op(1,op(i,lst))='Nor_v(v)' then N_val:=op(2,op(i,lst)) fi;

> if op(1,op(i,lst))='U(v)' then Xu_val:=op(2,op(i,lst)) fi;

> if op(1,op(i,lst))='y(v)' then surf_val:=op(2,op(i,lst)); fi;

> if op(1,op(i,lst))='Nor_u(u)' then N_val:=op(2,op(i,lst)); fi;

> if op(1,op(i,lst))='diff(x(u),u)' then Xu_val:=op(2,op(i,lst)) fi;

> if op(1,op(i,lst))='x(u)' then surf_val:=op(2,op(i,lst)); fi;

> if op(1,op(i,lst))='diff(y(v),v)' then Xv_val:=op(2,op(i,lst)) fi;

> if op(1,op(i,lst))='V(u)' then Xv_val:=(op(2,op(i,lst))) fi;

> od;

> RETURN([surf_val,Xu_val,Xv_val,N_val]);

> end;

~~~~~~~~~~~~~~~~~~~~

# Main procedure

~~~~~~~~~~~~~~~~~~~~

> if (check_compatibility_conditions()=0) then RETURN() fi; u:='u'; v:='v';

> step_size_u:=abs(max_u-min_u)/(no_of_steps_u-1);

> step_size_v:=abs(max_v-min_v)/(no_of_steps_v-1);

>

> X:=[];

> a11:=simplify(subs(u=min_u,alpha11));

> a12:=simplify(subs(u=min_u,alpha12));

> a22:=simplify(subs(u=min_u,alpha22));

> b1_1:=simplify(subs(u=min_u, beta1_1));

> b1_2:=simplify(subs(u=min_u, beta1_2));

> b2_1:=simplify(subs(u=min_u, beta2_1));

> b2_2:=simplify(subs(u=min_u, beta2_2));

> g1_11:=simplify(subs(u=min_u, gamma1_11));

> g1_12:=simplify(subs(u=min_u, gamma1_12));

> g1_22:=simplify(subs(u=min_u, gamma1_22));

> g2_11:=simplify(subs(u=min_u, gamma2_11));

> g2_12:=simplify(subs(u=min_u, gamma2_12));

> g2_22:=simplify(subs(u=min_u, gamma2_22));

> E0:=(subs({u=min_u,v=min_v},E));

> F0:=(subs({u=min_u,v=min_v},F));

> G0:=(subs({u=min_u,v=min_v},G));

>

> sys_v:={(D@@2)(y)(v)=g1_22*U(v)+g2_22*D(y)(v)+a22*Nor_v(v), D(U)(v)=g1_12*U(v)+g2_12*D(y)(v)+a12*Nor_v(v), D(Nor_v)(v)=b1_2*U(v)+b2_2*D(y)(v)};

> initv1:={ D(y)(min_v)=F0/sqrt(E0), y(min_v)=0, Nor_v(min_v)=0 , U(min_v)= sqrt(E0) };

> initv2:={ D(y)(min_v)=sqrt((E0*G0-F0^2)/E0), y(min_v)=0, Nor_v(min_v)=0 , U(min_v)=0 };

> initv3:={ D(y)(min_v)=0, y(min_v)=0, Nor_v(min_v)=1 , U(min_v)= 0};

> res_v_sys1:=dsolve(sys_v union initv1 ,{y(v),Nor_v(v),U(v)},type=numeric, method=dverk78, output=procedurelist):

> res_v_sys2:=dsolve(sys_v union initv2 ,{y(v),Nor_v(v),U(v)},type=numeric, method=dverk78, output=procedurelist):

> res_v_sys3:=dsolve(sys_v union initv3 ,{y(v),Nor_v(v),U(v)},type=numeric, method=dverk78, output=procedurelist):

>

> for i from 1 to no_of_steps_v do

> const_v:=min_v+step_size_v*(i-1);

> a11:=simplify(subs(v=const_v,alpha11));

> a12:=simplify(subs(v=const_v,alpha12));

> a22:=simplify(subs(v=const_v,alpha22));

> b1_1:=simplify(subs(v=const_v, beta1_1));

> b1_2:=simplify(subs(v=const_v, beta1_2));

> b2_1:=simplify(subs(v=const_v, beta2_1));

> b2_2:=simplify(subs(v=const_v, beta2_2));

> g1_11:=simplify(subs(v=const_v, gamma1_11));

> g1_12:=simplify(subs(v=const_v, gamma1_12));

> g1_22:=simplify(subs(v=const_v, gamma1_22));

> g2_11:=simplify(subs(v=const_v, gamma2_11));

> g2_12:=simplify(subs(v=const_v, gamma2_12));

> g2_22:=simplify(subs(v=const_v, gamma2_22));

> sys_u:={(D@@2)(x)(u)=g1_11*D(x)(u)+g2_11*V(u)+a11*Nor_u(u), D(V)(u)=g1_12*D(x)(u)+g2_12*V(u)+a12*Nor_u(u), D(Nor_u)(u)=b1_1*D(x)(u)+b2_1*V(u)};

> init1:=find_init_val(res_v_sys1(const_v));

> init2:=find_init_val(res_v_sys2(const_v));

> init3:=find_init_val(res_v_sys3(const_v));

> initu1:={x(min_u)=init1[1], D(x)(min_u)=init1[2], V(min_u)=init1[3], Nor_u(min_u)= init1[4] };

> initu2:={x(min_u)=init2[1], D(x)(min_u)=init2[2], V(min_u)=init2[3], Nor_u(min_u)= init2[4] };

> initu3:={x(min_u)=init3[1], D(x)(min_u)=init3[2], V(min_u)=init3[3], Nor_u(min_u)= init3[4] };

> res_u_sys1:=dsolve(sys_u union initu1 ,{x(u),Nor_u(u),V(u)},type=numeric, method=dverk78, output=procedurelist):

> res_u_sys2:=dsolve(sys_u union initu2 ,{x(u),Nor_u(u),V(u)},type=numeric, method=dverk78, output=procedurelist):

> res_u_sys3:=dsolve(sys_u union initu3 ,{x(u),Nor_u(u),V(u)},type=numeric, method=dverk78, output=procedurelist):

> for j from 1 to no_of_steps_u do

> res1:=find_init_val(res_u_sys1(min_u+step_size_u*(j-1)));

> res2:=find_init_val(res_u_sys2(min_u+step_size_u*(j-1)));

> res3:=find_init_val(res_u_sys3(min_u+step_size_u*(j-1)));

> X:=[op(X),[res1[1],res2[1],res3[1]]];

> od;

> od;

> with(plots):

> print(polygonplot3d (X,style=point,axes=frame,title='The_approximated_patch'));

> RETURN(X);

> end:

>

>

# ----------------------------------------------------------------------------------

# EXAMPLE 1- The unit sphere :-

# X(u,v) = [ cos(u)*sin(v) , sin(u)*sin(v) , cos(v) ]

# E:=sin(v)^2, F:=0, G:=1, L:=sin(v)^2, M:=0, N:=1

# ----------------------------------------------------------------------------------

> min_u:=Pi/2: min_v:=Pi/2: max_u:=Pi/2+1: max_v:=Pi/2+1: u_steps:=21: v_steps:=11:

> patch1:=reconstruct_patch( sin(v)^2 , 0 , 1 , sin(v)^2 , 0 , 1 , u , v , min_u , min_v , max_u , max_v , u_steps , v_steps):

Warning, the name changecoords has been redefined

We compute the average error of the numerical procedure as a check.

The exact solution is

> X1:=[ -cos(u)*sin(v)+cos(min_u)*sin(min_v) , -cos(v)+cos(min_v) ,-sin(u)*sin(v)+sin(min_u)*sin(min_v)];

Compute values of the exact solution at the grid points

> h_u:=(max_u-min_u)/(u_steps-1):

> h_v:=(max_v-min_v)/(v_steps-1):

> x1:=[]:

> for i from 1 to v_steps do for j from 1 to u_steps do x1:=[op(x1),evalf(subs({v=min_v+(i-1)*h_v,u=min_u+h_u*(j-1)}, X1))]:od:od:

Compute average error (in each component) from the exact patch :

> average:=[0,0,0]:

> for i from 1 to u_steps*v_steps do

> average:= evalm(average + evalm(patch1[i]-x1[i])):

> od:

> evalf(evalm((average/(u_steps*v_steps))));

# -------------------------------------------------------------------------- #

# PATCH2 - A cylinder :- #

# X(u,v) = [ sin(u) , v , cos(u) ] #

# E:=1, F:=0, G:=1, L:=-1, M:=0, N:=0 #

# -------------------------------------------------------------------------- #

> min_u:=Pi/2: min_v:=Pi/2: max_u:=Pi/2+1: max_v:=Pi/2+1: u_steps:=21: v_steps:=11:

> patch2:=reconstruct_patch(1,0,1,-1,0,0,u,v, min_u , min_v , max_u , max_v , u_steps , v_steps):

Again compute average error. The exact solution is

> X2:=[-cos(u)+cos(min_u),v-min_v,sin(u)-sin(min_u)];

> h_u:=(max_u-min_u)/(u_steps-1): h_v:=(max_v-min_v)/(v_steps-1): x2:=[]:

> for i from 1 to v_steps do for j from 1 to u_steps do x2:=[op(x2),evalf(subs({v=min_v+(i-1)*h_v,u=min_u+h_u*(j-1)}, X2))]:od:od:

> average:=0:

> for i from 1 to u_steps*v_steps do average:= evalm(average +evalm(evalm(patch2[i]-x2[i]))) od:

> evalf(evalm(average/(u_steps*v_steps)));

# ------------------------------------------------------------------------------------------------------------------------ #

# PATCH3 - Torus :- #

# X(u,v) = [(4+sin(v))*cos(u) , (4+sin(v))*sin(u) , cos(v)] #

# E:=(4+sin(v))^2, F:=0, G:=1, L:=(4+sin(v))*sin(v), M:=0, N:=1 #

# ------------------------------------------------------------------------------------------------------------------------ #

> min_u:=Pi/2: min_v:=Pi/2: max_u:=Pi/2+1: max_v:=Pi/2+1: u_steps:=21: v_steps:=11:

> patch3:=reconstruct_patch((4+sin(v))^2,0,1,(4+sin(v))*sin(v),0,1,u,v, min_u , min_v , max_u , max_v , u_steps , v_steps):

# ----------------------------------------------------------------------------- #

# PATCH4- The unit sphere coefficients with a change in F :- #

# E:=sin(v)^2, F:=0.0001*cos(v), G:=1, L:=sin(v)^2, M:=0, N:=1 #

# ----------------------------------------------------------------------------- #

> min_u:=Pi/2: min_v:=Pi/2: max_u:=Pi/2+1: max_v:=Pi/2+1: u_steps:=21: v_steps:=21:

> patch4:=reconstruct_patch(sin(v)^2 , 0.0001*cos(v) , 1 , sin(v)^2 , 0 , 1 ,u,v, min_u , min_v , max_u , max_v , u_steps , v_steps):

# --------------------------------------------------------------------------- #

# PATCH5- The unit sphere coefficients with a change in E :- #

# E*:=sin(v)^2+sin(v)^4 , F:=0, G:=1, L:=sin(v)^2, M:=0, N:=1 #

# --------------------------------------------------------------------------- #

> min_u:=Pi/2: min_v:=Pi/2: max_u:=Pi/2+1: max_v:=Pi/2+1: u_steps:=21: v_steps:=21:

> patch5:=reconstruct_patch(sin(v)^2+sin(v)^4 , 0 , 1, sin(v)^2 , 0 , 1 ,u,v, min_u , min_v , max_u , max_v , u_steps , v_steps):