================ BEZIER 1 ================ Mapleを用いた実習を行います． Mapleを起動してください． 適当にプログラムの中の数字などを変更して試してみてください． 以下の説明の中で＊＊＊＊＊ではさまれた部分は Mapleに入力する文字列です．自分で入力あるいはコピーしてから Enter(Return) キーで実行します． 上から順に，実行するようにしてください． 上で定義した関数を後で利用すること があります． 関数のグラフを描くにはパッケイジ　plots と plottools を用います． ＊＊＊＊＊＊ with(plots): with(plottools): ＊＊＊＊＊＊ ベジエ曲線を描くには、plotを用います． これは、y=f(x)の形になるもの以外も描く必要があるためです． plotを用いて放物線を描いてみましょう． ＊＊＊＊＊＊ plot([-1 + 2*t, (-1 + 2*t)^2, t=0..1]); ＊＊＊＊＊＊ 縦と横の比を１対１にするには，1:1のボタンを押すか， プロット（Plot）メニューの絶対軸（Scaling Constrained）を選んで下さい． plotを用いて円（の一部）を描いてみましょう． ＊＊＊＊＊＊ plot([(1-t^2)/(1+t^2), (2*t)/(1+t^2), t=0..1]); ＊＊＊＊＊＊ 平面上にベジエ曲線を描くための準備として、まず次を入力してください． ＊＊＊＊＊ lerp := proc(p,q,r,s,t) local R; R:= (s - t)/(s - r)*p + (t - r)/(s - r)*q; end: ＊＊＊＊＊ ＊＊＊＊＊ badecas:=proc(L,r,s) local i, b; 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); od; b; end: ＊＊＊＊＊ ＊＊＊＊＊ 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: ＊＊＊＊＊ ＊＊＊＊＊ beziercurve:=proc(L, u) local n, b, 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)]: end: ＊＊＊＊＊ ＊＊＊＊＊ bezierplot:=proc(L, r, s) local bb; bb:=badecas(L, 0, 1): plot([bb[1],bb[2], t=r..s],color=blue, thickness=2); end: ＊＊＊＊＊ では、放物線を描く例から始めましょう. ＊＊＊＊＊ cpoly:=[[-1, 1], [0, -1], [1, 1]]; bezierplot(cpoly, 0,1); ＊＊＊＊＊ 少し曲線をのばしてみましょう.（tの範囲を変える） ＊＊＊＊＊ bezierplot(cpoly, -1, 2); ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[-1, 0], [0, 2], [1, 0]]; bezierplot(cpoly, 0,1); ＊＊＊＊＊ 問題１ 次の答えを別紙に記入せよ． ＊＊＊＊＊ beziercurve(cpoly, t); ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[1, 0], [1, 1], [0, 1]]; bezierplot(cpoly, 0,1); ＊＊＊＊＊ ＊＊＊＊＊ beziercurve(cpoly, t); ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[0, -1], [2, 0], [0, 1]]; AA:=plot(cpoly, color=black, thickness=2): display(AA, bezierplot(cpoly, 0,1)); ＊＊＊＊＊ ＊＊＊＊＊ beziercurve(cpoly, t); ＊＊＊＊＊ ド・カステリョのアルゴリズムと関係をみましょう． ＊＊＊＊＊ cpoly:=[[-1, 1], [0, -1], [1, 1]]; bezier:=badecas(cpoly,0,1): ppoint:=k->pointplot(eval(bezier, t=k/64) , color=black, thickness=2): BB:=decas(cpoly, 0,1): frame:=k-> seq(plot(eval(BB[i],t=k/64), thickness=2), i=1..nops(cpoly)-1): AA:=plot(cpoly, color=black, thickness=2): frame(0):=AA: p:=NULL: for k from 0 to 64 do p:=p, display(AA,frame(k), ppoint(k), plot([bezier[1],bezier[2], t=0..1],color=blue, thickness=2) ): end do: display(p,insequence=true, axes=none); ＊＊＊＊＊ cpolyの値を変えてみましょう．cpoly:=[[1, 0], [1, 1], [0, 1]];　に変えてみましょう． ＊＊＊＊＊ cpoly:=[[1, 0], [1, 1], [0, 1]]; bezier:=badecas(cpoly,0,1): ppoint:=k->pointplot(eval(bezier, t=k/64), color=black, thickness=2): BB:=decas(cpoly, 0,1): frame:=k-> seq(plot(eval(BB[i],t=k/64), thickness=2), i=1..nops(cpoly)-1): AA:=plot(cpoly, color=black, thickness=2): frame(0):=AA: p:=NULL: for k from 0 to 64 do p:=p, display(AA,frame(k), ppoint(k), plot([bezier[1],bezier[2], t=0..1],color=blue, thickness=2) ): end do: display(p,insequence=true, axes=none); ＊＊＊＊＊ 3次曲線を考えましょう．cpolyの値を変えてみましょう． ＊＊＊＊＊ cpoly:=[[3,2], [-1, -3], [-1, 3], [3,-2]]; bezier:=badecas(cpoly,0,1): ppoint:=k->pointplot(eval(bezier, t=k/64) , color=black, thickness=2): BB:=decas(cpoly, 0,1): frame:=k-> seq(plot(eval(BB[i],t=k/64), thickness=2), i=1..nops(cpoly)-1): AA:=plot(cpoly, color=black, thickness=2): frame(0):=AA: p:=NULL: for k from 0 to 64 do p:=p, display(AA,frame(k), ppoint(k), plot([bezier[1],bezier[2], t=0..1],color=blue, thickness=2) ): end do: display(p,insequence=true, axes=none); bezier; ＊＊＊＊＊ つぎ場合のcpolyと比較してみましょう． ＊＊＊＊＊ cpoly:=[[3,2], [-1, -2], [-1,2], [3,-2]]; ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[3,2], [-1, -1], [-1, 1], [3,-2]]; ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[0, 0], [1, 2], [4, 4], [6, 0]]; ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[0, 0], [2, 2], [4, -2], [6, 0]]; ＊＊＊＊＊ 高次のベジエ曲線の場合を考えましょう． ＊＊＊＊＊ cpoly:=[[-1,0], [-1,2], [2,3], [4,1], [3, -1]]; ＊＊＊＊＊ 問題２ 次の答えを別紙に記入せよ． ＊＊＊＊＊ cpoly:=[[-1,-1],[-3/5,1],[-1/5,-1],[1/5,1], [3/5,-1], [1,1]]; beziercurve(cpoly, t);　 ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[-1, 0], [1, 2], [-1, 2], [0, 3], [1, 2], [1, 0], [-1, 2], [-1, 0]]; ＊＊＊＊＊ ＊＊＊＊＊ cpoly:=[[2652/15625, 200512416/244140625], [-106132/234375, -9917022752/3662109375], [-1342484/1640625, -19394337248/25634765625], [-528372/546875, 1869010208/1220703125], [-219388/234375, 8206418272/3662109375], [-251836/328125, 11238002464/8056640625], [-273468/546875, -317292768/3759765625], [-40612/234375, -9268160096/8056640625], [40612/234375, -9268160096/8056640625], [273468/546875, -317292768/3759765625], [251836/328125, 11238002464/8056640625], [219388/234375, 8206418272/3662109375], [528372/546875, 1869010208/1220703125], [1342484/1640625, -19394337248/25634765625], [106132/234375, -9917022752/3662109375], [-2652/15625, 200512416/244140625]]; plot([beziercurve(cpoly,t)[1],beziercurve(cpoly,t)[2], t=0..1], color=blue, thickness=2); ＊＊＊＊＊ ＊＊＊＊＊ AA:=plot(cpoly, color=black, thickness=2): p:=NULL: for j from 0 to 64 do p:=p, display(AA, plot([beziercurve(cpoly,t)[1],beziercurve(cpoly,t)[2], t=-0.001..j/64], color=blue, thickness=2)); end do: display(p,insequence=true, axes=none); ＊＊＊＊＊ ベジエ曲線の次数上げを考えましょう．　 ＊＊＊＊＊　 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, 1], [0, -1], [1, 1]]; ＊＊＊＊＊ ＊＊＊＊＊ plot(cpoly, color=blue, thickness=2); ＊＊＊＊＊　 問題３ 次の答えを別紙に記入せよ． ＊＊＊＊＊ elevt(cpoly, 4); ＊＊＊＊＊ ＊＊＊＊＊ plot(elevt(cpoly, 4), color=red, thickness=2); ＊＊＊＊＊ 3次曲線の場合の次数上げ ＊＊＊＊＊ cpoly:=[[3,2], [-1, -2], [-1, 2], [3,-2]]; ＊＊＊＊＊ ＊＊＊＊＊ a1:=plot(elevt(cpoly, 1), color=red, thickness=2): a2:=plot(elevt(cpoly, 2), color=blue, thickness=2): display(a1,a2); ＊＊＊＊＊ つぎ場合と比較してみましょう． ＊＊＊＊＊　 cpoly:=[[3,2], [-1, -3], [-1, 3], [3,-2]];　 ＊＊＊＊＊　 ＊＊＊＊＊　 cpoly:=[[3,2], [-1, -1], [-1, 1], [3,-2]]; ＊＊＊＊＊　 次に，ベジエ曲線の細分割を考えましょう． 数回細分割するとベジエ曲線にほとんど一致することを確かめましょう． ＊＊＊＊＊ 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]]; subdivstep2(cpoly, 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]]; ＊＊＊＊＊