> #######################################################
> #
> # Semi-sharp Creases on Subdivision Curves and Surfaces
> #
> # SGP 2014
> #
> # Supplementary Maple file
> #
> # Input: my_deg - degree
> #
> # Output: mtrS - change of basis matrix
> #         NsubS - subdivision matrix
> #
> # Jiri Kosinka, 13.6.2014
> #
> ######################################################## 
> 
> restart:
> with(LinearAlgebra): with(plots): with(plottools): with(CurveFitting):
> interface(rtablesize=20):
> 
> # input degree d>=2
> my_deg:=3:

> # the Oslo algorithm
> 
> findInt:=proc(kn,refTau,j,tau)
>  local i:
>  for i from 1 to kn do
>   if ((tau[i] <= refTau[j]) and (refTau[j] < tau[i + 1])) then
>    return(i):
>   end if:
>  end do:
>  print("Knot vector problem!"):
> end proc:
> 
> findPoi:=proc(refTau,rp1,i,j,tau,plg,k)
>  local r,p1,p2,p12,poi1,poi2,tmp1,tmp2:
>  global p:
>  r:=rp1-1:
> 
>  if (r > 0) then
>   p1 := refTau[j + k - r] - tau[i]:
>   p2 := tau[i + k - r] - refTau[j + k - r]:
>   if p1 <> 0 then
>    poi1 := findPoi(refTau, r, i, j,tau,plg,k):
>   end if:
>   if p2 <> 0 then
>    poi2 := findPoi(refTau, r, i - 1, j,tau,plg,k):
>   end if:
>   p12 := p1 + p2:
>   tmp1 := (p1 / p12) * poi1:
>   tmp2 := (p2 / p12) * poi2:
>   return(tmp1 + tmp2):
>  else
>   return (plg[i]):
>  end if:
> end proc:
> 
> oslo:=proc(deg, tau, refTau, plg)
>  local q, mu, j, poi, refPlg, numv, k:
>  numv:=nops(plg):
>  k := deg + 1:
>  q := nops(refTau):
>  for j from 1 to q - k do
>   mu := findInt(k + numv, refTau, j, tau):
>   poi := findPoi(refTau, k, mu, j, tau,plg,k):
>   refPlg[j]:=poi:
>  end do:
>  return([seq(refPlg[j],j=1..q-k)]):
> end proc:

> # construct the subdivision matrix
> 
> Msize:=2*my_deg-1:         
> size:=ceil((3*my_deg-4)/2):
> p:=2*(size+my_deg+1):         
> plg:=[seq([seq(0,i=1..p)],j=1..p)]:
> for i from 1 to p do
>  plg[i][i]:=2^my_deg:
> end do:
> my_tau:=[seq(-2*(size+1),i=1..my_deg+1),seq(-2-2*size+2*i,i=1..size),seq(0,i=1..my_deg),seq(2*i,i=1..size),seq(2*(size+1),i=1..my_deg+1)]:
> my_refTau:=[seq(-2*(size+1),i=1..my_deg+1),seq(-2-2*size+i,i=1..size*2+1),seq(0,i=1..my_deg),seq(i,i=1..size*2+1),seq(2*(size+1),i=1..my_deg+1)]:
> Os:=oslo(my_deg,my_tau,my_refTau,plg):
> matice:=convert(Os,Matrix):
> f:=ceil(ColumnDimension(matice)/2):
> ff:=ceil(RowDimension(matice)/2):
> Bsub:=matice[ff-my_deg..ff+my_deg,f-my_deg..f+my_deg]:

> # construct the selection matrix
> 
> my_tau:=[seq(-2*(size+1),i=1..my_deg+1),seq(-2-2*size+2*i,i=1..size),seq(0,i=1..1),seq(2*i,i=1..size),seq(2*(size+1),i=1..my_deg+1)]:
> my_refTau:=[seq(-2*(size+1),i=1..my_deg+1),seq(-2-2*size+2*i,i=1..size),seq(0,i=1..my_deg),seq(2*i,i=1..size),seq(2*(size+1),i=1..my_deg+1)]:
> Os:=oslo(my_deg,my_tau,my_refTau,plg):
> matice:=convert(Os,Matrix)/2^my_deg:
> f:=ceil(ColumnDimension(matice)/2):
> ff:=ceil(RowDimension(matice)/2):
> mtrr:=matice[ff-my_deg..ff-my_deg+2*my_deg,ff-my_deg..ff-my_deg+my_deg+1]:
> m[ceil(RowDimension(mtrr)/2)]:=1:
> for i from 1 to ceil(my_deg/2)+1 do
>  m[i]:=0:
> end do:
> 
> for i from 1 to floor(RowDimension(mtrr)) do
>  m[RowDimension(mtrr)+1-i]:=m[i]:
> end do:
> V:=Vector([seq(m[i],i=1..RowDimension(mtrr))]):
> mtr:=<mtrr[1..RowDimension(mtrr), 1..ceil(ColumnDimension(mtrr)/2)] | V | mtrr[1..RowDimension(mtrr), ceil(ColumnDimension(mtrr)/2)+1..ColumnDimension(mtrr)]>:

> # construct the subdivision matrix with control vectors
> 
> my_tau:=[seq(-2*(size+1),i=1..my_deg+1),seq(-2-2*size+2*i,i=1..size),seq(0,i=1..1),seq(2*i,i=1..size),seq(2*(size+1),i=1..my_deg+1)]:
> my_refTau:=[seq(-2*(size+1),i=1..my_deg+1),seq(-2-2*size+i,i=1..size*2+1),seq(0,i=1..1),seq(i,i=1..size*2+1),seq(2*(size+1),i=1..my_deg+1)]:
> Os:=oslo(my_deg,my_tau,my_refTau,plg):
> matice:=convert(Os,Matrix):
> f:=2*my_deg:
> ff:=my_deg:
> Nsubb:=matice[f..f+my_deg+1,ff..ff+my_deg+1]:
> V:=Vector([seq(0,i=1..RowDimension(Nsubb))]):
> Nsubbb:=Transpose(<Transpose(Nsubb)[1..RowDimension(Nsubb), 1..ceil(ColumnDimension(Nsubb)/2)] | V | Transpose(Nsubb)[1..RowDimension(Nsubb), ceil(ColumnDimension(Nsubb)/2)+1..ColumnDimension(Nsubb)]>):
> for i from 1 to floor(RowDimension(Nsubbb)/2)-1 do
>  t[RowDimension(Nsubb)+2-i]:=t[i]:
> end do:
> V:=Vector([seq(t[i],i=1..RowDimension(Nsubb)+1)]):
> Nsub:=<Nsubbb[1..RowDimension(Nsubbb), 1..ceil(ColumnDimension(Nsubbb)/2)] | V | Nsubbb[1..RowDimension(Nsubbb), ceil(ColumnDimension(Nsubbb)/2)+1..ColumnDimension(Nsubbb)]>:

> # solve the resulting system
> 
> eq:=Multiply(Bsub,mtr)-Multiply(mtr,Nsub):
> all:=seq(seq(eq[i,j],i=1..RowDimension(eq)),j=1..ColumnDimension(eq)):
> sol:=solve({all}):
> nops([sol]);
> if nops([sol])=1 then
>  sl[1]:=sol:
> else
>  sl:=sol:
> end if:
> # alternatively, use Groebner basis computation
> 
> lst:=simplify({all}):
> lst2:=[op(lst[2..nops(lst)])]:
> nops(lst2):
> ind:=op(indets(lst2)):
> with(Groebner):
> b:=Basis(lst2[1..nops(lst2)],plex(ind)):

> # find the unique crease solution
> 
> for i from 1 to nops([sol]) do
>  if subs(sl[i],m[my_deg])<>1 then
>   mtrS:=subs(sl[i],mtr);
>   print(mtrS);
>   NsubS:=subs(sl[i],Nsub);
>   print(NsubS);
>  end if:
> end do:
>