##This Maple sheet proofs the equivalence of (2) and (4) in the Eurographics 2013 paper # Global Selection of Stream Surfaces # # J. Martinez Esturo, M. Schulze, C. RÃ¶ssl, H. Theisel # University of Magdeburg restart; with(LinearAlgebra): #Jacobian matrix J := Matrix([[u_x,u_y,u_z], [v_x,v_y,v_z], [w_x,w_y,w_z]]); #Velocity field V := ; #Tangent of seed curve Sdot := ; #Second derivative vector of seed curve Sddot := ; #Without loss of generality we assume than V is in the direction of the x-axis, and that Sdot is in the x-y plane. #This gives a normal in the direction of the z-axis v := 0; w := 0; wsdot := 0; #We use this assumption to set the orientation of the normal, which simplifies computations later on. Any alternative of #assume(u>0,vsdot>0); #assume(u<0,vsdot<0); #assume(u<0,vsdot>0); #works equally well for the proof. assume(u>0,vsdot<0); #First and second order partials of stream surface Xs := Sdot; Xt := V; Xss := Sddot; Xst := Multiply(J,Sdot); Xtt := Multiply(J,V); #Classic differential geometric coefficients of first fundamental form E := Multiply(Transpose(Xs),Xs); F := Multiply(Transpose(Xs),Xt); G := Multiply(Transpose(Xt),Xt); #Unnormalized normal HNormal := CrossProduct(Xs,Xt); #Unit normal NNormal := HNormal / sqrt(HNormal[1]^2 + HNormal[2]^2 + HNormal[3]^2); #More classic differential geometric coefficients of second fundamental form L := simplify(Multiply(Transpose(NNormal),Xss)); M := simplify(Multiply(Transpose(NNormal),Xst)); N := simplify(Multiply(Transpose(NNormal),Xtt)); #Now we compute the principal directions P1,P2. K := Matrix([[dt^2,-ds*dt,ds^2],[E,F,G],[L,M,N]]); dt1 := eval(solve(Determinant(K)=0,dt)[1],ds=ds1); dt2 := eval(solve(Determinant(K)=0,dt)[2],ds=ds2); ds1 := 1; ds2 := 1; #Unnormalized principal directions HP1 := ds1*Xs + dt1*Xt; HP2 := ds2*Xs + dt2*Xt; #Principal directions P1 := HP1 / sqrt(HP1[1]^2 + HP1[2]^2 +HP1[3]^2); P2 := HP2 / sqrt(HP2[1]^2 + HP2[2]^2 +HP2[3]^2); #Compute the angle alpha #Normalized velocity Vn := V / Norm(V,2); #Squared cosine of angle between between V and P1: cosp1q := Multiply(Transpose(P1),Vn)^2; #Squared sine of angle between V and P1: sinp1q := Multiply(Transpose(P2),Vn)^2; #Compute the principal curvatures k1,k2 k1 := solve((k*E-L)*(k*G-N) - (k*F-M)^2,k)[1]; k2 := solve((k*E-L)*(k*G-N) - (k*F-M)^2,k)[2]; #Now we have everything set up to compute e_a using (2) #Our error function e_a_2 := simplify(cosp1q*sinp1q*(k2-k1)^2); #We can simplify e_a (2) e_a_2s := factor(expand(numer(e_a_2))/expand(denom(e_a_2))); ############################################# #Now we compute e_a using (4) B := CrossProduct(V, NNormal); e_a_4 := (Multiply(Multiply(Transpose(NNormal),J),B))^2 / Norm(V,2)^4; ############################################# #We compare e_a_4 and e_a_3s, it should be zero to complete the proof simplify(e_a_2s - e_a_4);