Krassimira Vlachkova: Testing for coincidence two irreducible tensor product 
Bezier (TPB) surfaces of degree (n,m), notebook with examples   
Created by Wolfram Mathematica 10.0                    
Acknowledgments. This work was supported by Sofia University Science Fund
Grant No.48/2016

Input: Irreducible TPB surfaces S1 and S2 of same degree (n,m) given by their control polygons 
Output: (i) S1 and S2 are different;
                 (ii) S1 and S2 are disjoint;
               (iii) S1 and S2 have coincident part.Report its control points.

Clear[b,p];
(*Input the degrees n and m, and the control points of the two Bezier surfaces b and p.*)
(*EXAMPLES*)

(*EXAMPLE 1*)
(*n=2;
m=4;
b={{{0,0,0},{0,7/10,-1/2},{0,2,-4/5}},{{1,0,1},{1,1,0},{1,1,1}},
{{2,0,3/4},{2,1,6/5},{2,2,3/4}},
{{3,-1,3/20},{3,1,2/5},{3,2,1/5}},{{4,0,1/5},{4,2,1/4},{4,3,3/4}}
};

p={{{(2/3),(10649/32400),(179489/648000)},{(1/12),(21053/64800),-(43963/1296000)},{-(1/2),(6961/25920),-(1495087/2592000)},{-(13/12),(49961/259200),-(1416983/1036800)},{-(5/3),(1233973/2592000),-(22960591/10368000)}},{{(2/3),(15121/25920),(4897/32400)},{(1/12),(29017/51840),-(51089/259200)},{-(1/2),(53281/103680),-(9503/12960)},{-(13/12),(98281/207360),-(287029/207360)},{-(5/3),(311761/414720),-(915767/518400)}},{{(2/3),(1073/1296),(9689/103680)},{(1/12),(4243/5184),-(61339/207360)},{-(1/2),(4439/5184),-(373243/414720)},{-(13/12),(21235/20736),-(276875/165888)},{-(5/3),(176833/103680),-(3733231/1658880)}}}
*)


(*EXAMPLE 2*)
(*n=2;
m=4;
b={{{0,0,0},{1/2,7/10,-1/2},{1,2,-4/5}},
{{1,0,1},{1,1,0},{1,1,1}},
{{2,0,3/4},{2,1,6/5},{2,2,3/4}},
{{1,-1,3/20},{3,1,2/5},{2,2,1/5}},{{4,0,1/5},{4,2,1/4},{4,3,3/4}}
}

p={{{(1283/784),(61/560),(8933/15680)},{(1849/1568),(331/1568),(19993/31360)},{(2327/3136),(1067/3136),(4259/8960)},{-(35/128),(775/6272),-(27887/125440)},{-(6709/12544),(2567/12544),-(332491/250880)}},{{(1439/672),(1747/1120),(3557/6720)},{(1867/1344),(591/448),(4447/13440)},{(1991/2688),(1087/896),-(7123/26880)},{(2755/5376),(1975/1792),-(67673/53760)},{(8303/10752),(3687/3584),-(88939/107520)}},{{(19/12),(34/15),(889/1440)},{(4/3),(91/48),(2219/2880)},{(61/48),(27/16),(4879/5760)},{(5/6),(425/192),-(4471/11520)},{(463/192),(653/96),-(135203/23040)}}}*)


(*EXAMPLE 3*)
(*n=2;
m=3;
b={{{0,0,0},{0,0.7,-0.5},{0,2,-.8}},
{{1,0,1},{1,1,0},{1,1,1}},
{{2,0,.75},{2,1,1.2},{2,2,0.75}},
{{3,-1,0.15},{3,1,0.4},{3,2,0.2}}};

p={{{(1/2),0.32759259259259255`,0.18815740740740738`},{-(1/12),0.30879629629629624`,-0.18903240740740734`},{-(2/3),0.24689814814814826`,-0.8024884259259258`},{-(5/4),0.2736990740740721`,-1.6089386574074083`}},{{(1/2),0.5766203703703705`,0.05164351851851845`},{-(1/12),0.5466435185185182`,-0.3502893518518518`},{-(2/3),0.5024884259259259`,-0.9192418981481483`},{-(5/4),0.5533275462962958`,-1.4934403935185163`}},{{(1/2),0.8252314814814805`,-0.017766203703703687`},{-(1/12),0.8292824074074079`,-0.4684317129629626`},{-(2/3),0.9042245370370356`,-1.119719328703705`},{-(5/4),1.2187789351851863`,-1.8351634837962938`}}};*)





(*EXAMPLE 4*)

(*n=4;
m=2;
b={{{0,0,0},{1,0,1},{0,0,2},{1,0,3},{2,2,2}},
{{1,1,0},{4,1,1},{1,1,2},{2,1,4},{1,1,3}},
{{0,2,0},{1,2,1},{0,2,2},{2,2,0},{0,2,1}}};

p={{{(2959/2401),(4803/4802),(5457/9604)},{(11803/9604),(1203/4802),(11133/19208)},{-(1688/2401),-(9579/19208),(21597/38416)}},{{(2823/1372),(1031/1029),(7153/4116)},{(5583/2744),(1049/4116),(14957/8232)},{-(81/112),-(502/1029),(27913/16464)}},{{-(187/882),(449/441),(1228/441)},{-(43/126),(521/1764),(1441/441)},{(311/3528),-(341/882),(317/126)}},{{(51/14),(221/189),(587/189)},{(407/126),(509/756),(989/189)},{(1123/504),(211/378),(338/189)}},{{-(367/162),(209/81),(16/27)},{(209/162),(1361/324),-(16/9)},{(3851/648),(1519/162),-(16/27)}}};*)


(*EXAMPLE 5*)
(*n=3; m=2;
b={{{0,0,0},{3/2,0,−3},{3/2,0,3},{3,0,0}},
{{0,3/2,−3},{4/3,4/3,0},{7/6,5/3,−10 /3},{3,1,1}},{{0,3/2,3},{5/3,7/6,10/3},{11/6,11/6,0},{3,2,1}}(*,{{0,3,0},{1,3,−1},{2,3,−1},{3,3,0}}*)
};

p={{{(18/25),(331/750),-(4657/4500)},{(187/250),-(97/375),-(9097/9000)},{(437/500),-(293/200),(16883/18000)}},{{(203/200),(527/1200),-(883/900)},{(1267/1200),-(619/2400),-(1927/1800)},{(119/96),-(1397/960),(331/720)}},{{(197/160),(71/160),-(41/60)},{(1231/960),-(253/960),-(61/120)},{(2903/1920),-(559/384),(13/10)}},{{(93/64),(85/192),-(53/144)},{(193/128),-(103/384),(13/288)},{(449/256),-(365/256),(1195/576)}}};*)





(*EXAMPLE 6*)
n=2;
m=3;
b={{{0,0,0},{0,7/10,-1/2},{0,2,-4/5}},
{{1,0,1},{1,1,0},{1,1,1}},
{{2,0,3/4},{2,1,6/5},{2,2,3/4}},{{3,-1,3/20},{3,1,2/5},{3,2,1/5}}
};
(*b={{{0,0,0},{0,0.7,-0.5},{0,2,-0.8}},
{{1,0,1},{1,1,0},{1,1,1}},
{{2,0,0.75},{2,1,1.2},{2,2,0.75}},{{3,-1,0.15},{3,1,0.4},{3,2,0.2}}
}*)

p={{{(3/2),(347/1960),(4821/7840)},{(3/4),(233/784),(8301/15680)},{0,(307/1568),(321/31360)},{-(3/4),(557/3136),-(60063/62720)}},{{(3/2),(2351/1680),(667/1680)},{(3/4),(839/672),-(223/3360)},{0,(1531/1344),-(6233/6720)},{-(3/4),(2831/2688),-(13051/13440)}},{{(3/2),(727/360),(259/360)},{(3/4),(253/144),(1183/1440)},{0,(587/288),(17/720)},{-(3/4),(3037/576),-(23279/5760)}}};



b=Transpose[b]
p

(*Step 1.Compute vectors ρ_i,ρ_i^{1,0},and ρ_i^{0,1},i=1,2.*)
ro1=Sum[Sum[(-1)^(n+m-i-j)* Binomial[n,i]* Binomial[m,j]* b[[i+1,j+1]],{j,0,m}],{i,0,n}];
ro110=Sum[Sum[(-1)^(n+m-i-j)*Binomial[n,i]* Binomial[m,j]*(i-n)*b[[i+1,j+1]],{j,0,m}],{i,0,n}]/n;
ro101=Sum[Sum[(-1)^(n+m-i-j) Binomial[n,i]* Binomial[m,j]* (j-m)*b[[i+1,j+1]],{j,0,m}],{i,0,n}]/m;
ro2=Sum[Sum[(-1)^(n+m-i-j) Binomial[n,i] *Binomial[m,j]* p[[i+1,j+1]],{j,0,m}],{i,0,n}];
ro210=Sum[Sum[(-1)^(n+m-i-j) Binomial[n,i]* Binomial[m,j] *(i-n)*p[[i+1,j+1]],{j,0,m}],{i,0,n}]/n;
ro201=Sum[Sum[(-1)^(n+m-i-j) Binomial[n,i]* Binomial[m,j] *(j-m)*p[[i+1,j+1]],{j,0,m}],{i,0,n}]/m;

(*Step 2. Check the conditions of Lemma1.
      2.1 If ro1 and ro2 are non-collinear then return (i): S1 and S2 are different. *)

(If[Norm[Cross[ro1,ro2]]<0.0000001,Goto[l1],Goto[l2]];
Label[l2];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l1];)
kappa=(ro2.ro1)/(ro1.ro1);

(* 2.2 If either ro1,ro2,ro110,ro210, or ro1,ro2,ro101,ro201 are non-coplanar then return (i): S1 and S2 are different; else the main nonlinear system has either unique, or no solution. *)

mat1={{ro1[[1]],ro1[[2]],ro1[[3]]},{ro110[[1]],ro110[[2]],ro110[[3]]},{ro210[[1]],ro210[[2]],ro210[[3]]}};
(If[Det[mat1]==0.,Goto[l3],Goto[l4]];
Label[l4];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l3];)
mat2={{ro1[[1]],ro1[[2]],ro1[[3]]},{ro101[[1]],ro101[[2]],ro101[[3]]},{ro201[[1]],ro201[[2]],ro201[[3]]}};
(If[Det[mat2]==0.,Goto[l5],Goto[l6]];
Label[l6];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l5];)


(*Step 3. Solve the two linear systems. If any of them is inconsistent
then return (i);else denote their unique solutions by ((SuperStar[a]),(SuperStar[b])) and ((SuperStar[c]),(SuperStar[d]))*)

Clear[a1,b1,c1,d1]; 
sol1=Solve[Rationalize[
FullSimplify[b1*({1,0,0}.ro210)-a1*({1,0,0}.(ro2+ro210))==kappa*({1,0,0}.ro110 )]&&
FullSimplify[b1*({0,1,0}.ro210)-a1*({0,1,0}.(ro2+ ro210))==kappa*({0,1,0}.ro110 )]&&
FullSimplify[b1*({0,0,1}.ro210)-a1*({0,0,1}.(ro2+ro210))==kappa*({0,0,1}.ro110 )
]],{a1,b1}]


sol2=Solve[Rationalize[
FullSimplify[d1*({1,0,0}.ro201)-c1*({1,0,0}.(ro2+ro201))==kappa*({1,0,0}.ro101 )]&&
FullSimplify[d1*({0,1,0}.ro201)-c1*({0,1,0}.(ro2+ro201))==kappa*({0,1,0}.ro101 )]&&
FullSimplify[d1*({0,0,1}.ro201)-c1*({0,0,1}.(ro2+ro201))==kappa*({0,0,1}.ro101 )
]],{c1,d1}]

a1=sol1[[All,1,2]];
b1=sol1[[All,2,2]];
c1=sol2[[All,1,2]];
d1=sol2[[All,2,2]];
a1=a1[[1]];
b1=b1[[1]];
c1=c1[[1]];
d1=d1[[1]];

(*Step 4.*)
(If[Abs[(((b1-a1))^n)*(((d1-c1))^m)-kappa]<0.0000001,Goto[l7],Goto[l8]];
Label[l8];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l7];)


(*Step 5. Compute the control polygon of the transformed tensor product Bezier surface cc using blossoming principle.*)


Clear[u1,u2,u3,u4,u5,v1,v2,v3,v4,v5];
Clear[pwr,bern,bzSurf,g1,g2,u,v,s,t];
(*Handle Indeterminate condition*)
pwr[x_,y_]:=If[x==0&&y==0,1,x^y];
bern[n_,i_,u_]:=Binomial[n,i]pwr[u,i]pwr[1-u,n-i];
bzSurf[u_,v_]:=Sum[bern[n,i,u] b[[i+1,j+1]] bern[m,j,v],{i,0,n},{j,0,m}]






bzSurf[u,v]=Expand[Simplify[bzSurf[u,v]]];
bzSurf1=bzSurf[u,v][[1]];
bzSurf2=bzSurf[u,v][[2]];
bzSurf3=bzSurf[u,v][[3]];

Clear[pp1,pp2,pp3,cc];
pp=Array[{0,0,0}&,{n+1,m+1}];
cc=Array[{0,0,0}&,{n+1,m+1}];
pp1=Array[0&,{n+1,m+1}];
pp2=Array[0&,{n+1,m+1}];
pp3=Array[0&,{n+1,m+1}];

For[i=1,i≤n,i++,For[j=1,j≤m,j++,pp1[[i+1,j+1]]=Coefficient[bzSurf1,(u^i)* (v^j)]]];
For[i=1,i≤n,i++,{v=0;pp1[[i+1,1]]=Coefficient[bzSurf1,(u^i)],Clear[v]}];
For[j=1,j≤m,j++,{u=0;pp1[[1,j+1]]=Coefficient[bzSurf1,(v^j)];Clear[u]}];
pp1[[1,1]]=CoefficientList[bzSurf[u,v][[1]],{u,v}][[1,1]];

For[i=1,i≤n,i++,For[j=1,j≤m,j++,pp2[[i+1,j+1]]=Coefficient[bzSurf2,(u^i)* (v^j)]]];
For[i=1,i≤n,i++,{v=0;pp2[[i+1,1]]=Coefficient[bzSurf2,(u^i)],Clear[v]}];
For[j=1,j≤m,j++,{u=0;pp2[[1,j+1]]=Coefficient[bzSurf2,(v^j)];Clear[u]}];
pp2[[1,1]]=CoefficientList[bzSurf[u,v][[2]],{u,v}][[1,1]];

For[i=1,i≤n,i++,For[j=1,j≤m,j++,pp3[[i+1,j+1]]=Coefficient[bzSurf3,(u^i)* (v^j)]]];
For[i=1,i≤n,i++,{v=0;pp3[[i+1,1]]=Coefficient[bzSurf3,(u^i)],Clear[v]}];
For[j=1,j≤m,j++,{u=0;pp3[[1,j+1]]=Coefficient[bzSurf3,(v^j)];Clear[u]}];
pp3[[1,1]]=CoefficientList[bzSurf[u,v][[3]],{u,v}][[1,1]];

For[i=0,i≤n,i++,For[j=0,j≤m,j++,cc[[i+1,j+1]]={pp1[[i+1,j+1]],pp2[[i+1,j+1]],pp3[[i+1,j+1]]}]];
cc;

For[nu=0,nu≤n,nu++,For[mu=0,mu≤m,mu++,
pp[[nu+1,mu+1]]=Sum[Sum[cc[[i+1,j+1]]/(Binomial[n,i]*Binomial[m,j])*First[(Sum[(Sum[Binomial[nu,k]*Binomial[n-nu,i-k]*Binomial[mu,r]*Binomial[m-mu,j-r]*a1^k*b1^{i-k}*c1^r*d1^{j-r},{r,Max[0,j+mu-m],Min[j,mu]}]),{k,Max[0,i+nu-n],Min[i,nu]}])],{j,0,m}],{i,0,n}]]]

(*Compare the control points of surface p and the transformed surface pp. If they coincide (up to eight different configurations ot the control points) then surfaces p and pp are coincident.*)
p
pp

(*Plot the input surfaces b and p, and the transformed surface pp obtained by subdivision from b*)
f1=BezierFunction[p];
f2=BezierFunction[b];
f3=BezierFunction[pp];
g1=ParametricPlot3D[f1[u,v],{u,0,1},{v,0,1},Mesh->None,PlotStyle->Directive[Specularity[White,20],Opacity[0.8]],ColorFunction->"Rainbow",(*ColorFunction->Function[{x,y,z},Hue[z]],*)Boxed->False,Axes->None];
g1=Show[g1,Graphics3D[{EdgeForm[],Black,Thickness[0.005],BezierCurve[First[p],SplineDegree->m],BezierCurve[Last[p],SplineDegree->m],Black,Thickness[0.005],BezierCurve[p[[All,1]],SplineDegree->n],BezierCurve[p[[All,-1]],SplineDegree->n],Antialiasing->True}],PlotRange->All,Boxed->False, Axes->False];
g2=ParametricPlot3D[f2[u,v],{u,0,1},{v,0,1},Mesh->None,PlotStyle->Directive[Specularity[White,20],Opacity[0.8]],ColorFunction->"Rainbow",(*ColorFunction->Function[{x,y,z},Hue[z]],*)Boxed->False,Axes->None];
g2=Show[g2,Graphics3D[{EdgeForm[],Black,Thickness[0.005],BezierCurve[First[b],SplineDegree->m],BezierCurve[Last[b],SplineDegree->m],Black,Thickness[0.005],BezierCurve[b[[All,1]],SplineDegree->n],BezierCurve[b[[All,-1]],SplineDegree->n],Antialiasing->True}],PlotRange->All,Boxed->False, Axes->False];
g3=ParametricPlot3D[f3[u,v],{u,0,1},{v,0,1},Mesh->None,PlotStyle->Directive[Specularity[White,20],Opacity[0.8]],ColorFunction->"Rainbow",(*ColorFunction->Function[{x,y,z},Hue[z]],*)Boxed->False,Axes->None];
g3=Show[g3,Graphics3D[{EdgeForm[],Black,Thickness[0.005],BezierCurve[First[pp],SplineDegree->m],BezierCurve[Last[pp],SplineDegree->m],Black,Thickness[0.005],BezierCurve[pp[[All,1]],SplineDegree->n],BezierCurve[pp[[All,-1]],SplineDegree->n],Antialiasing->True}],PlotRange->All,Boxed->False, Axes->False];
Show[g1,g2]
Show[g1,g3]