================
BEZIER2　(次数上げ，細分割）
================
Mapleを用いた実習を行います．
Mapleを起動してください．


適当にプログラムの中の数字などを変更して試してみてください．

以下の説明の中で＊＊＊＊＊ではさまれた部分は
Mapleに入力する文字列です．自分で入力あるいはコピーし
Enter(Return) キーで実行します．
上から順に，すべて実行するようにしてください．
上で定義した関数を後で利用すること 
があります．

関数のグラフを描くには Maple のパッケージ　plots と plottools
を用います．

＊＊＊＊＊＊
with(plots): with(plottools):
＊＊＊＊＊＊

平面上にベジエ曲線を描くための準備として、まず次を入力してください．

＊＊＊＊＊
 bezier1:=proc(L, clr, thk, r,s) local n, b, bb, k;
 n:=nops(L)-1; b:= [sum(n!/(k!*(n-k)!)*(u^k)*((1-u)^(n-k))*L[k+1][1], k=0..n), 
sum(n!/(k!*(n-k)!)*(u^k)*((1-u)^(n-k))*L[k+1][2],k=0..n)]:
 bb:=plot([b[1], b[2], u=r..s], color=clr, thickness=thk);
 end:
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[[-1,0], [-1,2], [2,3], [4,1], [3, -1],[2,-1]];
＊＊＊＊＊

＊＊＊＊＊
bezier1(cpoly, red,3,0, 1);
＊＊＊＊＊


＊＊＊＊＊
AA:=plot(cpoly, color=black, thickness=2):
p:=NULL:
for j from 1 to 64 do
p:=p, display(AA, bezier1(cpoly,blue, 2, 0, j/64) ):
end do: 
display(p,insequence=true, axes=none);
＊＊＊＊＊

ベジエ曲線の次数上げを考えましょう．
＊＊＊＊＊
lerp := proc(p,q,r,s,t)  local R;
 R:= (s - t)/(s - r)*p + (t - r)/(s - r)*q;  end:
＊＊＊＊＊

＊＊＊＊＊
elevt:=proc(L, rr) local j, i, b, bb ; b:=L; for i from 1 to rr do 
bb:=seq( lerp(b[j], b[j+1], 0,1,(nops(L)-j+i-1)/(nops(L)+i-1) ), j=1..nops(L)-2+i);
b:=[ L[1], bb, L[nops(L)] ]; od; b;
end;
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[[-1, 0], [1, 2], [-1, 2], [0, 3], [1, 2], [1, 0], [-1, 2], [-1, 0]];
＊＊＊＊＊

＊＊＊＊＊
plot(cpoly,  color=blue, thickness=2);
＊＊＊＊＊

＊＊＊＊＊
elevt(cpoly, 4);
＊＊＊＊＊

＊＊＊＊＊
plot(elevt(cpoly, 4), color=red, thickness=2);
＊＊＊＊＊


＊＊＊＊＊
p:=NULL:
for j from 1 to 64 do
p:=p, display(plot(elevt(cpoly, 1), color=black, thickness=2), bezier1(elevt(cpoly, 1),blue, 2, 0, j/64) ):
end do:
display(p,insequence=true, axes=none);
＊＊＊＊＊


＊＊＊＊＊
BB1:=bezier1(cpoly,blue, 1, 0, 1):
BB2:=bezier1(elevt(cpoly,1),red,4, 0, 1):
display(BB2,  BB1);
＊＊＊＊＊


3次曲線の場合の次数上げ
＊＊＊＊＊
cpoly:=[[3,2], [-1, -3], [-1, 3], [3,-2]];
p:=NULL:
for j from 1 to 64 do
p:=p, display(plot(elevt(cpoly, 1), color=black, thickness=2), bezier1(elevt(cpoly, 1),blue, 2, 0, j/64) ):
end do:
display(p,insequence=true, axes=none);
＊＊＊＊＊

つぎ場合と比較してみましょう．

＊＊＊＊＊
cpoly:=[[3,2], [-1, -2], [-1,2], [3,-2]];
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[[3,2], [-1, -1], [-1, 1], [3,-2]];
＊＊＊＊＊

高次のベジエ曲線の場合の次数上げを考えましょう．
次数を上げて行くとどうなりますか．
＊＊＊＊＊
cpoly:=[[-1,0], [-1,2], [2,3], [4,1], [3, -1]];
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[[0,0],[-3,0],[-3,1],[-1,1], [-1,3], [-3,3]];
＊＊＊＊＊


次に，ベジエ曲線の細分割を考えましょう．
数回細分割するとベジエ曲線にほとんど一致することを確かめましょう．

次の項目を入力したとき，選択
function definition
remmber table assignment
が出たときは，maple11では
remmber table assignmentを選ぶ．

＊＊＊＊＊
decas:=proc(L,r,s) local i,j, b, cc,c;
 b:=L;
 for i from 1 to nops(L)-1 do 
 b:=seq( lerp(b[j], b[j+1], r,s,t), j=1..nops(L)-i);
 cc(i):=b;  
 od; 
 c:=seq([cc(i)],i=1..nops(L)-1);
 end:
＊＊＊＊＊

＊＊＊＊＊
lsubdivstep:=proc(L,r,s) local  lpoly ;
 lpoly:=[L[1], seq(eval(decas(L, r,s)[i][1],t=(r+s)/2), i=1..nops(L)-1)];
 end;
＊＊＊＊＊

＊＊＊＊＊
rsubdivstep:=proc(L,r,s) local  rpoly; 
 rpoly:=[seq(eval(decas(L, r,s)[nops(L)-i][i],t=(r+s)/2), i=1..nops(L)-1), L[nops(L)]];
 end;
＊＊＊＊＊


＊＊＊＊＊
subdivstep2:=proc(L, n) local i,j, b, c,d, dpoly; 
 b[1]:=L; c[1]:=lsubdivstep(b[1], 0, 1); d[1]:=rsubdivstep(b[1], 0, 1);  
 for j from 1 to n do
 b[2*i-1]:=c[i]; b[2*i]:=d[i];
 for i from 1 to 2^j do
 c[i]:=lsubdivstep(b[i], 0, 1);   d[i]:=rsubdivstep(b[i], 0, 1); od;
 for i from 1 to 2^j do
 b[2*i-1]:=c[i]; b[2*i]:=d[i];
 od: od:
 dpoly:=seq(b[j], j=1..2^(n));
 end:
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[[-1, 1], [0, -1], [1, 1]];
subdivstep2(cpoly, 2);
＊＊＊＊＊

＊＊＊＊＊
AA:=plot(cpoly, color=black, thickness=2);
subdivstep2(cpoly, 2);
display(AA, plot([subdivstep2(cpoly, 2)],color=red, axes=none, thickness=2));
＊＊＊＊＊

＊＊＊＊＊
display(AA, plot([subdivstep2(cpoly, 4)], color=red, axes=none, thickness=2));
＊＊＊＊＊

次の場合の細分割を考え，数回細分割するとベジエ曲線にほとんど一致することを確かめましょう．
＊＊＊＊＊
cpoly:=[[0,0],[-3,0],[-3,1],[-1,1], [-1,3], [-3,3]];
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[[-1, 0], [1, 2], [-1, 2], [0, 3], [1, 2], [1, 0], [-1, 2], [-1, 0]];
＊＊＊＊＊


次に，空間内のベジエ曲線の細分割を考えましょう．
数回細分割するとベジエ曲線にほとんど一致することを確かめましょう．

＊＊＊＊＊
cpoly:=[ [-1, 1, -1], [-1/3,-1/3,1], [1/3,-1/3,-1], [1,1,1]];
dpoly:=subdivstep2(cpoly, 1):
＊＊＊＊＊

＊＊＊＊＊
frame:=k-> spacecurve(dpoly[k], color=blue, thickness=2):
AA:=spacecurve(cpoly, color=red, thickness=2):
frame(0):=AA:
p:=NULL:
for k from 0 to 2^1 do
p:=p, display(AA,frame(k)):
 end do:               
display(p,axes=boxed);
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[ [-1, 1, -1], [-1/3,-1/3,1], [1/3,-1/3,-1], [1,1,1]];
dpoly:=subdivstep2(cpoly, 4):
frame:=k-> spacecurve(dpoly[k], color=blue, thickness=2):
AA:=spacecurve(cpoly, color=red, thickness=2):
frame(0):=AA:
p:=NULL:
for k from 0 to 2^4 do
p:=p, display(AA,frame(k)):
 end do:               
display(p,axes=boxed);
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[ [2, 0, 0], [2,2,0], [2,2,2], [0,0,2], [2,0,0]];
＊＊＊＊＊

＊＊＊＊＊
cpoly:=[ [2, 0, 0], [2,2,0], [2,2,2], [0,0,2], [-1,2,0]];
＊＊＊＊＊