# Partial MAPLE code for the paper # "Minimal bi-6 $G^2$ completion of bicubic spline surfaces" # by K. Karciauskas and J. Peters # For n=6, as an example, # the following Maple code resolves all the G1 and G2 constraints # by assigning some of the bi-6 BB-coefficients of the cap # in terms of # 1. input quad net: a[] (B-spline coeff's), E (extraordinary node) # a2- a5- a6-- * # a3- a6- E a6+ # a4 a5 a6 a5+ # a1 a2 a3 a4+ # 2. tau # # This allows providing a[] and E (the quad net of a challenge geometry) # and choosing tau # (repeatedly to optimize the observed highlight lines) # and # applying functionals (eg F5) to obtain the remaining # degrees of freedom = set P in step 4 of page 7 # # In the solution presented in the paper, we optimized tau # and chose the functional by generating surfaces # for every choice of input quad mesh from the obstacle course and # for every functional (there are only a few functionals), # We rejected functionals yielding poorly distributed # highlight lines and performed a form of binary search on tau # starting with tau=1. ################ Notation ################ # n, c, d[2], tau -- as in the paper # mu1 -- dot mu, mu2 -- ddot mu, mu3 -- dddot mu, # mut -- tilde mu, mub -- bar mu # p -- BB-coefficients of bi-6 patches, but indexing is different: # p[i,j] in the paper, but here is p[6-i,6-j] ################## TBdata: tensor border from bi-3 ####### # here p[i,j] the BB-coefficients of input deg3 tensor-border, # output -- reparameterizied tensor-border = a part of bi-6 patches TBdata:=unapply( [[p[0,0], 1/2*p[0,0]+1/2*p[1,0], 1/5*p[0,0]+1/5*p[2,0]+3/5*p[1,0], 1/20*p[0,0]+9/20*p[1,0]+9/20*p[2,0]+1/20*p[3,0], 1/5*p[3,0]+1/5*p[1,0]+3/5*p[2,0], 1/2*p[2,0]+1/2*p[3,0], p[3,0]], [1/2*p[0,0]+1/2*p[0,1], 1/4*p[0,0]+1/4*p[1,0]+1/4*p[0,1]+1/4*p[1,1], 1/10*p[0,0]+1/10*p[0,1]+3/10*p[1,0]+1/10*p[2,0]+3/10*p[1,1]+1/10*p[2,1], 1/40*p[0,0]+9/40*p[1,0]+9/40*p[2,0]+1/40*p[3,0]+1/40*p[0,1]+9/40*p[1,1 ]+9/40*p[2,1]+1/40*p[3,1], 1/10*p[1,0]+3/10*p[2,0]+1/10*p[3,0]+1/10*p[1,1]+3/10*p[2,1]+1/10*p[3,1 ], 1/4*p[2,0]+1/4*p[3,0]+1/4*p[2,1]+1/4*p[3,1], 1/2*p[3,0]+1/2*p[3,1]], [3/5*p[0,1]+1/5*p[0,0]+1/5*p[0,2], 1/10*p[1,2]+1/10*p[0,0]+1/10*p[1,0]+3/10*p[0,1]+3/10*p[1,1]+1/10*p[0,2 ], 3/25*p[1,2]+1/25*p[2,2]+1/25*p[0,0]+3/25*p[1,0]+1/25*p[2,0]+3/25*p[0,1 ]+9/25*p[1,1]+3/25*p[2,1]+1/25*p[0,2]+1/50*p[0,1]*d[2]-1/25*d[2]*p[0,0 ]+1/50*d[2]*p[1,0], 1/600*(18*d[2]*p[2,0]*tau+3*p[0,1]*tau*d[2]+2*p[0,1]*tau*c-12*p[1,0]*c *tau+10*p[0,0]*c*tau-6*p[0,0]*d[2]*tau+27*p[1,1]*d[2]*tau+12*p[1,0]*c+ 54*p[1,2]*tau+162*p[2,1]*tau+18*p[0,1]*tau+162*p[1,1]*tau+54*p[2,0]* tau+6*p[3,0]*tau+6*p[0,0]*tau+54*p[1,0]*tau+54*p[2,2]*tau+6*p[3,2]*tau + 18*p[3,1]*tau-42*d[2]*p[1,0]*tau+6*p[0,2]*tau-12*p[0,0]*c)/tau, 1/150*(6*p[1,2]*tau+18*p[2,2]*tau+6*p[3,2]*tau+6*p[1,0]*tau+18*p[2,0]* tau+6*p[3,0]*tau+18*p[1,1]*tau+54*p[2,1]*tau+18*p[3,1]*tau+9*p[2,1]*d[ 2]*tau+2*p[1,1]*tau*c+3*p[1,1]*d[2]*tau-5*d[2]*p[1,0]*tau+6*p[1,0]*c* tau-8*p[1,0]*c-10*d[2]*p[2,0]*tau-8*p[2,0]*c*tau+8*p[2,0]*c+3*d[2]*p[3, 0]*tau)/tau, 1/60*(6*p[2,2]*tau+6*p[3,2]*tau+6*p[2,0]*tau+6*p[3,0]*tau+18*p[2,1]* tau+18*p[3,1]*tau+3*p[3,1]*d[2]*tau+2*p[2,1]*tau*c+3*p[2,1]*d[2]*tau-4 *d [2]*p[2,0]*tau+2*p[2,0]*c*tau-4*p[2,0]*c-2*d[2]*p[3,0]*tau-4*p[3,0]*c* tau+4*p[3,0]*c)/tau, 1/5*p[3,2]+1/5*p[3,0]+3/5*p[3,1]+1/15*p[3,1]*c+1/10*p[3,1]*d[2]-1/15*p [3,0]*c-1/10*d[2]*p[3,0]]], seq(seq(p[i,j],i=0..3),j=0..2),c,d[2],tau): ########## BtoBB ################ # conversion from B-spline P to bi-3 BB-form b BtoBB:=proc(P) local i,j,b,p; for i from 0 to 3 do for j from 0 to 3 do p[i,j]:=P[j+1][i+1]: od:od: b[0,0]:=(16*p[1,1]+4*(p[1,0]+p[2,1]+p[1,2]+p[0,1])+p[0,0]+p[2,0]+p[2,2 ]+p[0,2])/36: b[3,0]:=(16*p[2,1]+4*(p[2,0]+p[3,1]+p[2,2]+p[1,1])+p[1,0]+p[3,0]+p[3,2 ]+p[1,2])/36: b[0,3]:=(16*p[1,2]+4*(p[1,1]+p[2,2]+p[1,3]+p[0,2])+p[0,1]+p[2,1]+p[2,3 ]+p[0,3])/36: b[3,3]:=(16*p[2,2]+4*(p[2,1]+p[3,2]+p[2,3]+p[1,2])+p[1,1]+p[3,1]+p[3,3 ]+p[1,3])/36: b[1,0]:=4/9*p[1,1]+2/9*p[2,1]+1/9*(p[1,0]+p[1,2])+1/18*(p[2,0]+p[2,2]) : b[2,0]:=2/9*p[1,1]+4/9*p[2,1]+1/18*(p[1,0]+p[1,2])+1/9*(p[2,0]+p[2,2]) : b[1,3]:=4/9*p[1,2]+2/9*p[2,2]+1/9*(p[1,1]+p[1,3])+1/18*(p[2,1]+p[2,3]) : b[2,3]:=2/9*p[1,2]+4/9*p[2,2]+1/18*(p[1,1]+p[1,3])+1/9*(p[2,1]+p[2,3]) : b[0,1]:=4/9*p[1,1]+2/9*p[1,2]+1/9*(p[0,1]+p[2,1])+1/18*(p[0,2]+p[2,2]) : b[0,2]:=2/9*p[1,1]+4/9*p[1,2]+1/18*(p[0,1]+p[2,1])+1/9*(p[0,2]+p[2,2]) : b[3,1]:=4/9*p[2,1]+2/9*p[2,2]+1/9*(p[1,1]+p[3,1])+1/18*(p[1,2]+p[3,2]) : b[3,2]:=2/9*p[2,1]+4/9*p[2,2]+1/18*(p[1,1]+p[3,1])+1/9*(p[1,2]+p[3,2]) : b[1,1]:=4/9*p[1,1]+2/9*p[2,1]+2/9*p[1,2]+1/9*p[2,2]: b[2,1]:=4/9*p[2,1]+2/9*p[1,1]+2/9*p[2,2]+1/9*p[1,2]: b[1,2]:=4/9*p[1,2]+2/9*p[1,1]+2/9*p[2,2]+1/9*p[2,1]: b[2,2]:=4/9*p[2,2]+2/9*p[2,1]+2/9*p[1,2]+1/9*p[1,1]: [seq([seq(b[i,j],i=0..3)],j=0..3)]: # output end: sf:=k->k mod n: ################### MAIN ###################### # # see the paper "Minimal bi-6 $G^2$ completion of bicubic spline surfaces" # for steps (S1),(S2),(S3),(S4) n:=6: c:=cos(2*Pi/n): # simplifies to 1/2 # (S1) ---- define central quadratic for s from 0 to n-1 do # central point p[s][6,6]:=P00: od: p[0][5,6]:=TAN[0]: p[0][6,5]:=TAN[1]: p[0][4,6]:=TWO[0]: p[0][5,5]:=TWO[1]: p[0][6,4]:=TWO[2]: # --- glue together the quadratics of the sectors C2 Eq. (6) for s from 0 to n-2 do p[sf(s+1)][5,6]:=p[s][6,5]: p[sf(s+1)][4,6]:=p[s][6,4]: p[sf(s+1)][6,5]:=factor(-p[s][5,6]+2*c*p[s][6,5]+2*(1-c)*p[s][6,6]): p[sf(s+1)][5,5]:=factor(5/3*p[s][6,4]*c+1/3*c*p[s][6,6]-2*p[s][6,5]*c- p[s][5,5]+2*p[s][6,5]): p[sf(s+1)][6,4]:= factor(24/5*tau*c*p[s][5,6]-24/5*tau*c^2*p[s][6,5]-24/5*tau*c*p[s][6,6 ]+24/5*tau*c^2*p[s][6,6]-4*p[s][5,6]+44/5*p[s][6,5]*c+ 4*p[s][6,6]-4*c*p[s][6,6]+p[s][4,6]+4*p[s][6,4]*c^2+4/5*c^2*p[s][6,6]- 24/5*p[s][6,5]*c^2-24/5*p[s][5,5]*c): od: # (S2) assignment of coefficients with subscript (in paper) 21: # circulant system mu1[0] := 2-mu1[1]-mu1[2]-mu1[3]-mu1[4]-mu1[5]: mu1[1] := 2/5*c: mu1[2] := -c+2: mu1[3] := c: mu1[4] := -2/5*c: mu1[5] := 1/15*c: for r from 0 to n-1 do lyg[2*r+1]:=p[sf(r+1)][4,5]-(-p[r][5,4]+ mu1[0]*p[r][6,6]+mu1[1]*p[r][6,5]+mu1[2]*p[r][6,4]+ mu1[3]*p[r][6,2]+mu1[4]*p[r][6,1]+mu1[5]*p[r][6,0]): od: mu2[0] := -mu2[1]-mu2[2]-mu2[3]-mu2[4]-mu2[5]-mu2[6]-mu2[7]-mu2[8]: mu2[1] := -4/5*c*tau^2: mu2[2] := 24/5*c^2*tau+4/5*c^2-4*c+4/5*c^2*tau^2+4-24/5*c*tau: mu2[3] := 24/5*c*tau-4: mu2[4] := -2/3*c*(-11+6*c*tau+3*c): mu2[5] := -4*c: mu2[6] := 2*c^2: mu2[7] := -4/5*c^2: mu2[8] := 2/15*c^2: for r from 0 to n-1 do lyg[2*r]:=p[sf(r+1)][5,4]-(p[r][4,5]+ mu2[0]*p[r][6,6]+mu2[1]*p[r][5,6]+mu2[2]*p[r][6,5]+mu2[3]*p[r][5,5]+ mu2[4]*p[r][6,4]+mu2[5]*p[r][5,4]+ mu2[6]*p[r][6,2]+mu2[7]*p[r][6,1]+mu2[8]*p[r][6,0]): od: # solve equations (7) and assign dependent variables spr1:=solve({seq(lyg[r]=0,r=0..2*n-1)}, {seq(p[r][5,4],r=0..n-1),seq(p[r][4,5],r=1..n-1)}): assign(spr1); S1err := [seq(factor(lyg[r]),r=0..2*n-1)]; # (S2) assignment of coefficients with subscript (in paper) 22: # circulant system mu3[0] := -mu3[1]-mu3[2]-mu3[3]-mu3[4]-mu3[5]-mu3[6]-mu3[7]-mu3[8]: mu3[1] := 4/25*tau^3*c: mu3[2] := -4/25*c*(6*c*tau^2-6*tau^2-5+3*c*tau+tau^3*c): mu3[3] := -24/25*c*tau^2: mu3[4] := -2*c+4/5*c^2*tau^2+6/5*c^2*tau+4-12/5*c*tau: mu3[5] := 12/5*c*tau-4: mu3[6] := -2/5*c*(3*c*tau-5): mu3[7] := 4/25*c*(3*c*tau-5): mu3[8] := -2/75*c*(3*c*tau-5): mu3[9] := -3/5*c/tau: mu3[10] := -6/125*c/tau^2*(-5+c*tau): mu3[11] := 1/125*c*(-5+6*c*tau)/tau^3: for r from 0 to n-1 do lyc[r]:=factor(p[sf(r+1)][4,4]-(p[r][4,4]+ mu3[0]*p[r][6,6]+mu3[1]*p[r][5,6]+mu3[2]*p[r][6,5]+mu3[3]*p[r][5,5]+ mu3[4]*p[r][6,4]+mu3[5]*p[r][5,4]+ mu3[6]*p[r][6,2]+mu3[7]*p[r][6,1]+mu3[8]*p[r][6,0] +mu3[9]*(p[r][4,2]-p[sf(r+1)][2,4])+mu3[10]*(p[r][4,1]-p[sf(r+1)][1,4] )+mu3[11]*(p[r][4,0]-p[sf(r+1)][0,4]))): od: # solve equations (8) and assign dependent variables spr2:=solve({seq(lyc[r]=0,r=1..n-1)},{seq(p[r][4,4],r=1..n-1)}): assign(spr2); for r from 1 to n-1 do p[r][4,4]:=factor(p[r][4,4]): od: S2err := [seq(factor(lyc[r]),r=1..n-1)]; # show 1:n-1 is zero sal:=factor(lyc[0]): # not yet zero ########## d2 (tau) see (12,13) ################################## d[2]:=factor(1/2*(3*tau^2+tau^2*c-3*tau-2*tau*c+1)/tau^2): # ---Border data ----------- for s from 0 to n-1 do tam[s]:=BtoBB([ #[a1[s],a2[s],a3[s],a4[sf(s+1)]], #[a4[s],a5[s],a6[s],a5[sf(s+1)]], #[a3[sf(s-1)],a6[sf(s-1)],E,a6[sf(s+1)]], #[a2[sf(s-1)],a5[sf(s-1)],a6[sf(s-2)],0]]): [0,0,0,0], [0,0,0,0], [0,0,0,0], [0,0,0,0]]): # 16th node = e.o.position is not used (many-valued for n>4), # hence is set to 0. od: for s from 0 to n-1 do # shift for TBdata for i from 0 to 3 do for j from 0 to 3 do pp[s][i,j]:=tam[s][j+1][i+1]: od:od:od: # get boundary dependent data via TBdata for s from 0 to n-1 do tad[s]:=TBdata(seq(seq(pp[s][i,j],i=0..3),j=0..2),c,d[2],tau): tak[s]:=TBdata(seq(seq(pp[s][j,i],i=0..3),j=0..2),c,d[2],tau): od: for s from 0 to n-1 do for i from 0 to 6 do for j from 0 to 2 do p[s][i,j]:=factor(tad[s][j+1][i+1]): p[s][j,i]:=factor(tak[s][j+1][i+1]): od:od:od: factor(sal); # now remaining 22 equation is satisfied # reduced degree sector curve Eq (5) for r from 0 to n-1 do p[r][6,3]:=3/4*p[r][6,4]+3/4*p[r][6,2]+1/20*p[r][6,6]-3/10*p[r][6,5]-3 /10*p[r][6,1]+1/20*p[r][6,0]: od: for r from 0 to n-1 do # glue sectors for i from 0 to 6 do p[sf(r+1)][i,6]:=p[r][6,i]: od:od: # (S3) assignment of coefficients with subscript (in paper) 31, 32 (9) mut[0] := 1-mut[1]-mut[2]-mut[3]-mut[4]-mut[5]-mut[6]: mut[1] := -11/20: mut[2] := 1/8*(-10+11*c)/c: mut[3] := 5/4/c: mut[4] := -1/4*c+1/8: mut[5] := 3/10*c-1/20: mut[6] := -1/20*c+1/120: mut[7] := 3/8/tau: mut[8] := 3/100/tau^2*(c*tau-5): mut[9] := -1/200*(-5+6*c*tau)/tau^3: mut[10] := -5/16/c: for r from 0 to n-1 do p[r][5,3]:=factor( mut[0]*p[r][6,6]+mut[1]*p[r][6,5]+mut[2]*p[r][6,4]+mut[3]*p[r][5,4]+ mut[4]*p[r][6,2]+mut[5]*p[r][6,1]+mut[6]*p[r][6,0] +mut[7]*(p[r][4,2]-p[sf(r+1)][2,4])+mut[8]*(p[r][4,1]-p[sf(r+1)][1,4]) +mut[9]*(p[r][4,0]-p[sf(r+1)][0,4]) +mut[10]*(p[r][4,4]-p[sf(r+1)][4,4])): p[sf(r+1)][3,5]:=factor( mut[0]*p[r][6,6]+mut[1]*p[r][6,5]+mut[2]*p[r][6,4]+mut[3]*p[sf(r+1)][4 ,5]+ mut[4]*p[r][6,2]+mut[5]*p[r][6,1]+mut[6]*p[r][6,0] -mut[7]*(p[r][4,2]-p[sf(r+1)][2,4])-mut[8]*(p[r][4,1]-p[sf(r+1)][1,4]) -mut[9]*(p[r][4,0]-p[sf(r+1)][0,4]) -mut[10]*(p[r][4,4]-p[sf(r+1)][4,4])): od: mub[0] := -mub[1]-mub[2]-mub[3]-mub[4]-mub[5]-mub[6]: mub[1] := 1: mub[2] := -5/2*(-2+c)/c: mub[3] := -5/c: mub[4] := 5/2: mub[5] := -1: mub[6] := 1/6: mub[7] := 3/10/tau*(c*tau-5): mub[8] := 3/125/tau^2*(-20*c*tau+25+c^2*tau^2): mub[9] := -1/250*(-45*c*tau+25+6*c^2*tau^2)/tau^3: mub[10] := 5/4/c: for r from 0 to n-1 do p[sf(r+1)][3,4]:=factor(p[r][4,3]+ mub[0]*p[r][6,6]+mub[1]*p[r][6,5]+mub[2]*p[r][6,4]+mub[3]*p[r][5,4]+ mub[4]*p[r][6,2]+mub[5]*p[r][6,1]+mub[6]*p[r][6,0] +mub[7]*(p[r][4,2]-p[sf(r+1)][2,4])+mub[8]*(p[r][4,1]-p[sf(r+1)][1,4]) +mub[9]*(p[r][4,0]-p[sf(r+1)][0,4]) +mub[10]*(p[r][4,4]-p[sf(r+1)][4,4])): od: ######### smoothness verification across sector boundaries ################### # # verify : is G1 -------- for s from 0 to n-1 do F[s]:=0: G[s]:=0: for i from 0 to 6 do for j from 0 to 2 do F[s]:=F[s]+p[s][6-j,6-i]*binomial(6,i)*binomial(6,j)*(1-u)^(6-i)*u^i*( 1-v)^(6-j)*v^j: G[s]:=G[s]+p[sf(s+1)][6-i,6-j]*binomial(6,i)*binomial(6,j)*(1-u)^(6-i) *u^i*(1-v)^(6-j)*v^j: od:od: Fu[s]:=factor(unapply(diff(F[s],u),v)(0)): Fv[s]:=factor(unapply(diff(F[s],v),v)(0)): Gv[s]:=factor(unapply(diff(G[s],v),v)(0)): Fuu[s]:=factor(unapply(diff(F[s],u,u),v)(0)): Fuv[s]:=factor(unapply(diff(F[s],u,v),v)(0)): Fvv[s]:=factor(unapply(diff(F[s],v,v),v)(0)): Gvv[s]:=factor(unapply(diff(G[s],v,v),v)(0)): od: for s from 0 to n-1 do G1[s]:=factor(Fv[s]+Gv[s]-2*c*(1-u)^2*Fu[s]): od: G1err := [seq(G1[s],s=0..n-1)]; # verify : is G2 -------- A:=24*c*(tau-1)/(1-u+tau*u): for s from 0 to n-1 do G2[s]:=factor(-Gvv[s]+Fvv[s]-4*c*(1-u)^2*Fuv[s]+4*c^2*(1-u)^4*Fuu[s] + c*(1-u)^3*(A*u-8*c-A)*Fu[s]+A*(1-u)^2*Fv[s]): od: G2err := [seq(G2[s],s=0..n-1)]; # collect free values (other than tau) "p[s][i,j], s=0:5, i,j=0:6 contain the patches in terms of Pfree"; # renaming to be consistent with the paper p[0][4,4] := P220: p[0][4,5] := P210: for r from 0 to n-1 do p[r][4,3] := P23[r]; od: for r from 0 to n-1 do p[r][3,3] := P33[r]; od: for r from 0 to n-1 do for i from 0 to 6 do for j from 0 to 6 do p[r][i,j] := simplify(p[r][i,j]): od: od: od: Pfree := [ P00, TAN[0], TAN[1], TWO[0], TWO[1], TWO[2], P220, P210, seq(P23[r],r=0..n-1), seq(P33[r],r=0..n-1)];