Krassimira Vlachkova: Testing for coincidence two irreducible tensor product 
Bezier (TPB) surfaces of degrees (n,m) and (n+m,n+m), notebook with examples   
Created by Wolfram Mathematica  10.0                    
Acknowledgments. This work was supported by Sofia University Science Fund
Grant No.10-30-130/2017

Input: Integers n and m; Irreducible TPB surfaces S1 and S2 of degrees (n,m) and (n+m,n+m), respectively,   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},{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={{{(3/2),-(1/4),(93/160)},{(337/192),-(113/960),(327/640)},{(329/160),(73/1200),(2947/6400)},{(391/160),(51/160),(331/800)},{(35/12),(83/120),(811/2400)},{(41/12),(7/6),(7/24)},{4,(7/4),(29/80)}},{{(91/64),(47/960),(407/640)},{(5/3),(167/960),(2737/4608)},{(1841/960),(4673/14400),(8041/14400)},{(2129/960),(317/600),(20009/38400)},{(619/240),(593/720),(177/400)},{(71/24),(173/144),(71/192)},{(41/12),(5/3),(23/60)}},{{(611/480),(341/1200),(11689/19200)},{(1481/960),(6001/14400),(4463/7200)},{(717/400),(3367/6000),(8957/14400)},{(1979/960),(17737/24000),(14761/24000)},{(2843/1200),(8909/9000),(8859/16000)},{(643/240),(941/720),(19/40)},{(61/20),(203/120),(1091/2400)}},{{(33/32),(359/800),(1489/3200)},{(1273/960),(237/400),(20957/38400)},{(7619/4800),(1183/1600),(9671/16000)},{(591/320),(2897/3200),(82871/128000)},{(3403/1600),(9049/8000),(29893/48000)},{(2297/960),(847/600),(21401/38400)},{(431/160),(279/160),(209/400)}},{{(47/60),(337/600),(263/1200)},{(87/80),(2573/3600),(5323/14400)},{(1627/1200),(7781/9000),(71003/144000)},{(7733/4800),(1637/1600),(7181/12000)},{(15/8),(3701/3000),(11239/18000)},{(2027/960),(21503/14400),(16871/28800)},{(377/160),(2153/1200),(10537/19200)}},{{(7/12),(13/20),-(1/15)},{(7/8),(193/240),(203/1440)},{(91/80),(3413/3600),(4747/14400)},{(443/320),(437/400),(6479/12800)},{(1571/960),(18511/14400),(17153/28800)},{(89/48),(489/320),(13757/23040)},{(395/192),(579/320),(91/160)}},{{(1/2),(17/20),-(9/20)},{(3/4),(59/60),-(1/5)},{(59/60),(219/200),(71/1200)},{(193/160),(959/800),(1041/3200)},{(691/480),(1621/1200),(9673/19200)},{(105/64),(499/320),(55/96)},{(29/16),(29/16),(37/64)}}};*)

(*Example 2*)
n=2; m=3;
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}}
};


p={{{(5/4),-(1/8),(27/40)},{(117/80),(7/400),(31/50)},{(527/320),(257/1600),(3821/6400)},{(59/32),(5/16),(381/640)},{(83/40),(21/40),(391/800)},{(9/4),(3/4),(23/80)}},{{(93/80),(67/400),(127/200)},{(571/400),(663/2000),(1017/1600)},{(331/200),(1989/4000),(5221/8000)},{(601/320),(5331/8000),(22183/32000)},{(423/200),(89/100),(2557/4000)},{(91/40),(9/8),(399/800)}},{{(299/320),(581/1600),(2957/6400)},{(979/800),(537/1000),(8471/16000)},{(471/320),(5701/8000),(19203/32000)},{(273/160),(111/125),(11107/16000)},{(3109/1600),(8931/8000),(22287/32000)},{(67/32),(109/80),(389/640)}},{{(57/80),(81/160),(303/1600)},{(1609/1600),(5479/8000),(2027/6400)},{(503/400),(689/800),(1767/4000)},{(2387/1600),(8261/8000),(18963/32000)},{(69/40),(5029/4000),(1043/1600)},{(599/320),(2417/1600),(3893/6400)}},{{(11/20),(31/50),-(17/200)},{(83/100),(799/1000),(7/80)},{(1713/1600),(7719/8000),(8471/32000)},{(207/160),(1117/1000),(7643/16000)},{(607/400),(2643/2000),(963/1600)},{(133/80),(627/400),(487/800)}},{{(1/2),(17/20),-(9/20)},{(3/4),(51/50),-(49/200)},{(77/80),(37/32),-(17/1600)},{(371/320),(2021/1600),(361/1280)},{(109/80),(567/400),(397/800)},{(3/2),(13/8),(93/160)}}};


(*EXAMPLE 3*)
(*n=2; m=3;
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}}
};


p={{{-(1/2),-(11/20),(11/20)},{(21/40),-(193/200),(703/400)},{(1241/960),-(1657/1920),(23009/19200)},{(93367/69120),-(114839/115200),(683189/1382400)},{(5429/8640),-(89119/57600),(62351/691200)},{-(22187/27648),-(68791/27648),-(58051/552960)}},{{(19/60),-(7/15),(139/200)},{(39/40),-(397/2000),(2609/4000)},{(82517/57600),-(67817/288000),(730007/1152000)},{(66233/43200),-(11819/24000),(123833/259200)},{(372203/345600),-(5335481/5184000),(4779601/20736000)},{-(307/3456),-(330347/172800),-(1399/14400)}},{{(197/240),(1187/2400),-(2749/14400)},{(11549/9600),(70039/144000),(14533/64000)},{(93227/57600),(25537/86400),(2411099/3456000)},{(3956999/2073600),(1288969/31104000),(92916557/124416000)},{(6793/3840),-(2034671/5184000),(2592319/6912000)},{(240709/276480),-(65171/55296),-(158267/614400)}},{{(351/320),(5011/4800),(1223/19200)},{(3503/2400),(5219/6000),(20401/28800)},{(230143/115200),(159571/192000),(470963/460800)},{(285163/115200),(1252991/1728000),(3088163/3456000)},{(598643/230400),(53589/128000),(123251/512000)},{(9833/5120),-(7573/25600),-(226697/307200)}},{{(361/240),(613/400),(337/600)},{(1501/800),(17527/12000),(128513/144000)},{(727/300),(81181/48000),(422873/576000)},{(22207/7680),(944341/576000),(24861/51200)},{(2007/640),(12583/9600),-(17019/64000)},{(3369/1280),(3533/6400),-(7339/5120)}},{{(17/12),(727/360),(259/360)},{(749/480),(5059/2400),(1379/2400)},{(2317/1280),(46621/19200),(14159/76800)},{(7751/3840),(40921/19200),(2597/9216)},{(889/384),(14951/9600),-(7771/38400)},{(267/128),(429/640),-(717/512)}}};*)







b=Transpose[b]
p

(*Step 1.Compute vectors ρ,ρ^{1,0},ρ^{0,1}, and the corressponding eight finite diferences for the boundary curves of p.*)
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;


delta1=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*p[[i+1,1]],{i,0,n+m}];
delta2=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*p[[1,i+1]],{i,0,n+m}];
delta3=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*p[[i+1,n+m+1]],{i,0,n+m}];
delta4=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*p[[n+m+1,i+1]],{i,0,n+m}];


delta110=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*(i-n-m)*p[[i+1,1]],{i,0,n+m}];
delta201=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*(i-n-m)*p[[1,i+1]],{i,0,n+m}];
delta310=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*(i-n-m)*p[[i+1,n+m+1]],{i,0,n+m}];
delta401=Sum[(-1)^(n+m-i)*Binomial[n+m,i]*(i-n-m)*p[[n+m+1,i+1]],{i,0,n+m}];


(*Step 2. Check the conditions of Lemma2.
      2.1 If ro1 and any of delta1, delta2, delta3, and delta4 are non-collinear then return (i): S1 and S2 are different. *)

(If[Norm[Cross[ro1,delta1]]<0.0000001,Goto[l1],Goto[l2]];
Label[l2];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l1];
If[Norm[Cross[ro1,delta2]]<0.0000001,Goto[l21],Goto[l22]];
Label[l22];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l21];
If[Norm[Cross[ro1,delta3]]<0.0000001,Goto[l31],Goto[l32]];
Label[l32];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l31];
If[Norm[Cross[ro1,delta4]]<0.0000001,Goto[l41],Goto[l42]];
Label[l42];Print["Surfaces S1 and S2 are different."];
Exit[];
Label[l41];
)
kappa1=(delta1.ro1)/(ro1.ro1);
kappa2=(delta2.ro1)/(ro1.ro1);
kappa3=(delta3.ro1)/(ro1.ro1);
kappa4=(delta4.ro1)/(ro1.ro1);



(*Step 3. Solve the nonlinear system and find the solution a1,a2,b1,b2,c1,c2,d1,d2.*)
Clear[a1,a2,b1,b2,c1,c2,d1,d2];
solution=Solve[Rationalize[
(((b1-a1))^(n-1))*(((b2-a2))^(m-1))*({n*ro110[[1]],m*ro101[[1]],n*ro1[[1]],m*ro1[[1]]}.{b2-a2,b1-a1,a1*(b2-a2),a2*(b1-a1)})==delta110[[1]]&&
(((b1-a1))^(n-1))*(((b2-a2))^(m-1))*({n*ro110[[2]],m*ro101[[2]],n*ro1[[2]],m*ro1[[2]]}.{b2-a2,b1-a1,a1*(b2-a2),a2*(b1-a1)})==delta110[[2]]&&
(((b1-a1))^(n-1))*(((b2-a2))^(m-1))*({n*ro110[[3]],m*ro101[[3]],n*ro1[[3]],m*ro1[[3]]}.{b2-a2,b1-a1,a1*(b2-a2),a2*(b1-a1)})==delta110[[3]]&&
(((d1-a1))^(n-1))*(((d2-a2))^(m-1))*({n*ro110[[1]],m*ro101[[1]],n*ro1[[1]],m*ro1[[1]]}.{d2-a2,d1-a1,a1*(d2-a2),a2*(d1-a1)})==delta201[[1]]&&
(((d1-a1))^(n-1))*(((d2-a2))^(m-1))*({n*ro110[[2]],m*ro101[[2]],n*ro1[[2]],m*ro1[[2]]}.{d2-a2,d1-a1,a1*(d2-a2),a2*(d1-a1)})==delta201[[2]]&&
(((d1-a1))^(n-1))*(((d2-a2))^(m-1))*({n*ro110[[3]],m*ro101[[3]],n*ro1[[3]],m*ro1[[3]]}.{d2-a2,d1-a1,a1*(d2-a2),a2*(d1-a1)})==delta201[[3]]&&
(((c1-d1))^(n-1))*(((c2-d2))^(m-1))*({n*ro110[[1]],m*ro101[[1]],n*ro1[[1]],m*ro1[[1]]}.{c2-d2,c1-d1,d1*(c2-d2),d2*(c1-d1)})==delta310[[1]]&&
(((c1-d1))^(n-1))*(((c2-d2))^(m-1))*({n*ro110[[2]],m*ro101[[2]],n*ro1[[2]],m*ro1[[2]]}.{c2-d2,c1-d1,d1*(c2-d2),d2*(c1-d1)})==delta310[[2]]&&
(((c1-d1))^(n-1))*(((c2-d2))^(m-1))*({n*ro110[[3]],m*ro101[[3]],n*ro1[[3]],m*ro1[[3]]}.{c2-d2,c1-d1,d1*(c2-d2),d2*(c1-d1)})==delta310[[3]]&&
(((c1-b1))^(n-1))*(((c2-b2))^(m-1))*({n*ro110[[1]],m*ro101[[1]],n*ro1[[1]],m*ro1[[1]]}.{c2-b2,c1-b1,b1*(c2-b2),b2*(c1-b1)})==delta401[[1]]&&
(((c1-b1))^(n-1))*(((c2-b2))^(m-1))*({n*ro110[[2]],m*ro101[[2]],n*ro1[[2]],m*ro1[[2]]}.{c2-b2,c1-b1,b1*(c2-b2),b2*(c1-b1)})==delta401[[2]]&&
(((c1-b1))^(n-1))*(((c2-b2))^(m-1))*({n*ro110[[3]],m*ro101[[3]],n*ro1[[3]],m*ro1[[3]]}.{c2-b2,c1-b1,b1*(c2-b2),b2*(c1-b1)})==delta401[[3]]&&
kappa1==(((b1-a1))^n)*(((b2-a2))^m)&&
kappa2==(((d1-a1))^n)*(((d2-a2))^m)&&
kappa3==(((c1-d1))^n)*(((c2-d2))^m)&&
kappa4==(((c1-b1))^n)*(((c2-b2))^m)],{a1,a2,b1,b2,c1,c2,d1,d2},Reals,Method->{"UseSlicingHyperplanes"->False},VerifySolutions->True,WorkingPrecision->16
]

a1=solution[[All,1,2]];
a2=solution[[All,2,2]];
b1=solution[[All,3,2]];
b2=solution[[All,4,2]];
c1=solution[[All,5,2]];
c2=solution[[All,6,2]];
d1=solution[[All,7,2]];
d2=solution[[All,8,2]];

a1=a1[[1]];
b1=b1[[1]];
c1=c1[[1]];
d1=d1[[1]];
a2=a2[[1]];
b2=b2[[1]];
c2=c2[[1]];
d2=d2[[1]];

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



Clear[pwr,bern,bzSurf,g1,g2,u,v,w,s,t];
(*Handle Indeterminate condition*)
pwr[x_,y_]:=If[x==0&&y==0,1,x^y];
bern[n_,i_,w_]:=Binomial[n,i]pwr[w,i]pwr[1-w,n-i];
bzSurf[s_,t_]:=Sum[bern[n,i,s] b[[i+1,j+1]] bern[m,j,t],{i,0,n},{j,0,m}]


s=(1-u)*(1-v)*a1+(1-u)*v*d1+u*v*c1+u*(1-v)*b1;
t=(1-u)*(1-v)*a2+(1-u)*v*d2+u*v*c2+u*(1-v)*b2;


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


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

For[i=1,i≤n+m,i++,For[j=1,j≤n+m,j++,pp1[[i+1,j+1]]=Coefficient[bzSurf1,(u^i)* (v^j)]]];
For[i=1,i≤n+m,i++,{v=0;pp1[[i+1,1]]=Coefficient[bzSurf1,(u^i)],Clear[v]}];
For[j=1,j≤n+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+m,i++,For[j=1,j≤n+m,j++,pp2[[i+1,j+1]]=Coefficient[bzSurf2,(u^i)* (v^j)]]];
For[i=1,i≤n+m,i++,{v=0;pp2[[i+1,1]]=Coefficient[bzSurf2,(u^i)],Clear[v]}];
For[j=1,j≤n+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+m,i++,For[j=1,j≤n+m,j++,pp3[[i+1,j+1]]=Coefficient[bzSurf3,(u^i)* (v^j)]]];
For[i=1,i≤n+m,i++,{v=0;pp3[[i+1,1]]=Coefficient[bzSurf3,(u^i)],Clear[v]}];
For[j=1,j≤n+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+m,i++,For[j=0,j≤n+m,j++,cc[[i+1,j+1]]={pp1[[i+1,j+1]],pp2[[i+1,j+1]],pp3[[i+1,j+1]]}]];
cc;

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

(*Step 5. 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

(*Visualize the graphics of:  1). the two input surfaces S1 and S2; 2). S1 and the transformed surface obtained from S1. Since they coincide we must see only one surface.  *)
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->n+m],BezierCurve[Last[p],SplineDegree->n+m],Black,Thickness[0.005],BezierCurve[p[[All,1]],SplineDegree->n+m],BezierCurve[p[[All,-1]],SplineDegree->n+m],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->n+m],BezierCurve[Last[pp],SplineDegree->n+m],Black,Thickness[0.005],BezierCurve[pp[[All,1]],SplineDegree->n+m],BezierCurve[pp[[All,-1]],SplineDegree->n+m],Antialiasing->True}],PlotRange->All,Boxed->False, Axes->False];
Show[g1,g2]
Show[g1,g3]