(* 以下をコピーして mathematica のファイルに張付けて(ペースト)ください *) (* ベジエ曲線・曲面を描く準備です *) bernstein[n_, r_][t_] := If[n == r, If[r == 0, n!/(r!(n - r)!), n! t^r/(r!(n - r)!)], If[r == 0, n!(1 - t)^(n - r)/(r!(n - r)!), n!(1 - t)^(n - r)t^r/(r!(n - r)!)]] bezcurve[b_][t_] := Sum[b[[i]] bernstein[Length[b] - 1, i - 1][t], {i, 1, Length[b]}] ratbezcurve[w_, b_][t_] := (Sum[w[[i]] b[[i]] bernstein[Length[b] - 1, i - 1][t], {i, 1, Length[b]}])/ (Sum[w[[i]] bernstein[Length[b] - 1, i - 1][t], {i, 1, Length[b]}]) poly[b_] := Graphics[Table[Line[{b[[i]], b[[i + 1]]}], {i, 1, Length[b] - 1}]] pts[b_] := Graphics[{PointSize[0.02], Table[Point[b[[i]]], {i, 1, Length[b]}]}] bezsurface[b_][u_, v_] := Sum[Sum[b[[i + 1]][[j + 1]]*bernstein[Length[b] - 1, i][u]*bernstein[Length[b[[1]]] - 1, j][v], {i, 0, Length[b] - 1}], {j, 0, Length[b[[1]]] - 1}] ratbezsur[w_, b_][u_, v_] := (Sum[Sum[w[[i + 1]][[j + 1]] b[[i + 1]][[j + 1]]*bernstein[Length[b] - 1, i][u]* bernstein[Length[b[[1]]] - 1, j][v], {i, 0, Length[b] - 1}], {j, 0, Length[b[[1]]] - 1}])/(Sum[Sum[w[[i + 1]][[j + 1]]* bernstein[Length[b] - 1, i][u] bernstein[Length[b[[1]]] - 1, j][v], {i, 0, Length[b] - 1}], {j, 0, Length[b[[1]]] - 1}]) beznet[b_] := Module[{tmp1, tmp2}, tmp1 = Table[Graphics3D[Line[Table[b[[i + 1]][[j + 1]], {j, 0, Length[b[[1]]] - 1}]]], {i, 0, Length[b] - 1}]; tmp2 = Table[Graphics3D[Line[Table[b[[i + 1]][[j + 1]], {i, 0, Length[b] - 1}]]], {j, 0, Length[b[[1]]] - 1}]; Join[tmp1, tmp2]] bezpts[b_] := Graphics3D[{PointSize[0.015], Flatten[Table[Point[b[[i + 1]][[j + 1]]], {i, 0, Length[b] - 1}, {j, 0, Length[b[[1]]] - 1}], 1]}] Needs["Graphics`Colors`"] beziercurve[b_, color_]:=ParametricPlot[bezcurve[ b][t]//Evaluate,{t,0,1},AspectRatio\[Rule]Automatic,Axes\[Rule] True,PlotRange\[Rule]All,PlotStyle\[Rule]{Thickness[0.01],color}] bezpoly[b_]:=Show[poly[b],AspectRatio->Automatic, Axes->False] rationalbeziercurve[b_,w_, color_]:=ParametricPlot[ratbezcurve[w,b][t]// Evaluate,{t,0,1},AspectRatio\[Rule] Automatic,Axes\[Rule]True,PlotRange\[Rule] All,PlotStyle\[Rule]{Thickness[0.01],color}] ********************* (* それでは放物線の例から始めましょう *) bb1={{-1,1},{0,-1},{1,1}} AA1=beziercurve[bb1, Red] BB1=bezpoly[bb1] Show[AA1,BB1] bb2={{-1,0},{0,2},{1,0}} AA2=beziercurve[bb2, Blue] BB2=bezpoly[bb2] Show[AA2,BB2] bb3={{0,-1},{2,0},{0,1}} AA3=beziercurve[bb3, Blue] BB3=bezpoly[bb3] Show[AA3,BB3] (* 4点のベジエ曲線の例 *) cc1={{2,-2},{0,2},{0,-2},{2,2}} AA4=beziercurve[cc1, Blue] BB4=bezpoly[cc1] Show[AA4,BB4] cc2={{2,-2},{0,3},{0,-3},{2,2}} AA5=beziercurve[cc2, Red] BB5=bezpoly[cc2] Show[AA5,BB5] cc3={{2,-2},{0,1},{0,-1},{2,2}} AA6=beziercurve[cc3, Red] BB6=bezpoly[cc3] Show[AA6,BB6] cc4={{0,0},{4/3,-4/3},{4/3,4/3},{0,0}} AA7=beziercurve[cc4, Blue] BB7=bezpoly[cc4] Show[AA7,BB7] ******************* (* 有理ベジエ曲線の例1 *) bb1={{-1,1},{0,-1},{1,1}} ww1={1,2,1} AA1=rationalbeziercurve[bb1,ww1, Red] BB1=bezpoly[bb1] Show[AA1,BB1] bb1={{-1,1},{0,-1},{1,1}} ww2={1,1/2,1} AA2=rationalbeziercurve[bb1,ww2, Red] BB2=bezpoly[bb1] Show[AA2,BB2] bb3={{1,0},{1,1},{0,1}} ww3={1,1,2} AA3=rationalbeziercurve[bb3,ww3, Blue] BB3=bezpoly[bb3] Show[AA3,BB3] (* 有理ベジエ曲線の例2 *) bb4={{0,0},{-2,0},{-(2/3),-(8/3)},{4,-( 8/3)},{4,8/3},{-(2/3),8/3},{-2,0},{0,0}} ww4={1,1/7,1/7,3/35,3/35,1/7,1/7,1} AA4=rationalbeziercurve[bb4,ww4, Blue] BB4=bezpoly[bb4] Show[AA4,BB4] bb5={{0,0},{0,2},{-3, 1/2},{-3,-(15/2)},{20/3,-5},{20/3,5},{-3,15/2},{-3,-(1/2)},{0,-2},{0,0}} ww5={1,1/9,1/9,1/21,1/21,1/21,1/21,1/9,1/9,1} AA5=rationalbeziercurve[bb5,ww5, Blue] BB5=bezpoly[bb5] Show[AA5,BB5] ************************ (* ベジエ曲面を描く *) wireframe[surface_, {u_,u0_,u1_}, {v_,v0_,v1_}, opts__]:= Module[{plottmp, grtmp}, plottmp = ParametricPlot3D[surface//Evaluate, {u, u0, u1}, {v, v0, v1}, DisplayFunction->Identity, opts]; grtmp = plottmp/. (Polygon[pts_] :> Line[Append[pts, First[pts]]]); Show[grtmp, DisplayFunction->$DisplayFunction]] ********************* cc1 = {{{0, 0, 0}, {0, 1, 0}}, {{1, 0, 0}, {1, 1, 1}}} Show[beznet[cc1], bezpts[cc1], Axes -> True, PlotRange -> All] ParametricPlot3D[bezsurface[cc1][u,v]//Evaluate,{u,0,1},{v,0,1}, ViewPoint->{2, -1.5, 1}] wireframe[bezsurface[cc1][u,v]//Evaluate,{u,0,1},{v,0,1}, ViewPoint->{2,- 1.5, 1/2}] ParametricPlot3D[bezsurface[cc1][u, v] // Evaluate, {u, -1, 1}, {v, -1, 1}, ViewPoint -> {2, 1.5, 1/2}] wireframe[bezsurface[cc1][u, v] // Evaluate, {u, -1, 1}, {v, -1, 1}, ViewPoint -> {3, 2, 1}] ********************* ********************* cc2 = {{{-1, -1, 0}, {-1, 0, 2}, {-1, 1, 0}}, {{0, -1, -2}, {0, 0, 0}, {0, 1, -2}}, {{1, -1, 0}, {1, 0, 2}, {1, 1, 0}}} Show[beznet[cc2], bezpts[cc2], Axes -> True, PlotRange -> All, ViewPoint -> {1, 1.5, 3}] ParametricPlot3D[bezsurface[cc2][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, ViewPoint -> {1, 1.5, 1}] wireframe[bezsurface[cc2][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, ViewPoint -> {1, 1.5, 1}] ********************* ********************* cc3={{{0,0,0},{1,0,1},{2,0,0}},{{0,1,1},{1,1,2},{2,1,1}},{{0,2,0},{1,2,1},{2, 2,0}}} Show[beznet[cc3],bezpts[cc3],ViewPoint->{2,1, 2},PlotRange\[Rule]All] wireframe[bezsurface[cc3][u,v]//Evaluate,{u,0,1},{v,0,1}, ViewPoint->{1,1.5, 1}] Show[ParametricPlot3D[bezsurface[cc3][u,v]//Evaluate,{u,0,1},{v,0,1}, DisplayFunction->Identity],beznet[cc3],bezpts[cc3], ViewPoint->{2,1, 1}, PlotRange->All, DisplayFunction->$DisplayFunction] ********************* ********************* cc4 = {{{6, 0, 0}, {6, 6, 0}, {0, 6, 0}}, {{6, 0, 2}, {6, 6, 2}, {0, 6, 2}}, {{4, 0, 2}, {4, 4, 2}, {0, 4, 2}}} Show[beznet[cc4], bezpts[cc4], ViewPoint -> {1, 3, 2}, PlotRange -> All] Show[ParametricPlot3D[bezsurface[cc4][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, DisplayFunction -> Identity], beznet[cc4], bezpts[cc4], ViewPoint -> {1, 3, 2}, PlotRange -> All, DisplayFunction -> $DisplayFunction] ********************* ********************* cc5 = {{{1, 0, 0}, {1, 0, 1}, {0, 0, 1}}, {{1, 1, 0}, {1, 1, 1}, {0, 0, 1}}, {{0, 1, 0}, {0, 1, 1}, {0, 0, 1}}} Show[beznet[cc5], bezpts[cc5], ViewPoint -> {1, 3, 2}, PlotRange -> All] Show[wireframe[bezsurface[cc5][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, ViewPoint -> {1, 3, 2}], beznet[cc5], bezpts[cc5]] ********************* ********************* (* 有理ベジエ曲面を描く *) cc4 = {{{6, 0, 0}, {6, 6, 0}, {0, 6, 0}}, {{6, 0, 2}, {6, 6, 2}, {0, 6, 2}}, {{4, 0, 2}, {4, 4, 2}, {0, 4, 2}}} ww4 = {{1,1,2},{1,1,2},{2,2,4}} Show[ParametricPlot3D[ratbezsur[ww4, cc4][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, DisplayFunction -> Identity], beznet[cc4], bezpts[cc4], ViewPoint -> {-2, 1, 1}, PlotRange -> All, DisplayFunction -> $DisplayFunction] Show[ParametricPlot3D[ ratbezsur[ww4, cc4][u, v] // Evaluate, {u, -1, 1}, {v, -1, 1}, DisplayFunction -> Identity], beznet[cc4], bezpts[cc4], ViewPoint -> {1, -3, 2}, PlotRange -> All, DisplayFunction -> $DisplayFunction] Show[ParametricPlot3D[ ratbezsur[ww4, cc5][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, DisplayFunction -> Identity], beznet[cc5], bezpts[cc5], ViewPoint -> {2, -4, 1}, PlotRange -> All, DisplayFunction -> $DisplayFunction] Show[ParametricPlot3D[ ratbezsur[ww4, cc5][u, v] // Evaluate, {u, -1, 1}, {v, -1, 1}, DisplayFunction -> Identity], beznet[cc5], bezpts[cc5], ViewPoint -> {1, -2, 0}, PlotRange -> All, DisplayFunction -> $DisplayFunction] ********************* ********************* (** スプーンを描いてみましょう。行列として、制御点(ベジエ点)を入力します。 まずスプーンの先を描きましょう。**) cc = {{{4, 0, 0}, {4, 1/2, 0}, {3, 3/2, 0}, {3/2, 3/2, 0}, {1/2, 2, 0}, {0, 1/4, 0}}, {{4, 0, 0}, {4, 1/2*1/8, -1/3}, {3, 3/2*1/4, -1/2}, {3/2, 3/2*1/2, -1/2}, {1/2, 2*1/2, -1}, {0, 1/4, 0}}, {{4, 0, 0}, {4, 0, -1/3}, {3, 0, -1/2}, {3/2, 0, -1/2}, {1/2, 0, -1}, {0, 0, 0}}, {{4, 0, 0}, {4, -1/2*1/8, -1/3}, {3, -3/2*1/4, -1/2}, {3/2, -3/2*1/2, -1/2}, {1/2, -2*1/2, -1}, {0, -1/4, 0}}, {{4, 0, 0}, {4, -1/2, 0}, {3, -3/2, 0}, {3/2, -3/2, 0}, {1/2, -2, 0}, {0, -1/4, 0}}} Show[ beznet[cc], bezpts[cc], Axes -> True, PlotRange -> All, ViewPoint -> {2,3, 3}] (* <{{-0.1,4.1}, {-2,2},{-1,0}}, PlotPoints->{30,50}, ViewPoint->{1, -1.5, 1.3}], beznet[cc], bezpts[cc]] (**** つぎにスプーンの柄の部分を描きましょう。***) bb1 ={{{0, 1/8, 0}, {-1, 5/16, 1}, {-6, 6/16, 2}}, {{0, -1/8, 0}, {-1, -5/16, 1}, {-6, -6/16, 2}}} WWW = ParametricPlot3D[bezsurface[bb1][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {8, 16}, Boxed -> False] WWW1 = Show[wireframe[bezsurface[bb1][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {6, 40}, Boxed -> False], beznet[bb1], bezpts[bb1], PlotRange->All] (***** スプーン全体です。***) ZZZ = Show[WWW, VVV, PlotRange->All,ViewPoint->{1, -2, 1}] Show[beznet[cc], bezpts[cc], beznet[bb1], bezpts[bb1], PlotRange->All, ViewPoint->{10, 20, 10}] ZZZ1 = Show[WWW1, VVV1, PlotRange->All,ViewPoint->{10, 20, 10}] ZZZ1 = Show[WWW, VVV, PlotRange->All,ViewPoint->{10, 20, 10}] (*** 次に、ワイングラスを描きましょう。まず円形にするために重みを決めます。***) w = { {1, 1, 2}, {1, 1, 2}, {1, 1, 2}} dd1 = {{ {-3, 0, 12}, {-3, 3, 12}, {0, 3, 12}}, { {-6, 0, 6}, {-6, 6, 6}, {0, 6, 6}}, { {-3/2, 0, 6}, {-3/2, 3/2, 6}, {0, 3/2, 6}}} (***** ベジエネットを描きます。****) P1 = Show[ beznet[dd1], bezpts[dd1], Axes -> True, PlotRange -> All] QQ1 = ParametricPlot3D[ ratbezsur[w, dd1][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {20, 20}] Show[P1, QQ1] Q1 = ParametricPlot3D[ ratbezsur[w, dd1][u, v] // Evaluate, {u, 0, 1}, {v, -1, 1}, PlotPoints -> {20, 20}] dd2 = {{ {3, 0, 12}, {3, 3, 12}, {0, 3, 12}}, { {6, 0, 6}, {6, 6, 6}, {0, 6, 6}}, { {3/2, 0, 6}, {3/2, 3/2, 6}, {0, 3/2, 6}}} P2 = Show[ beznet[dd2], bezpts[dd2], Axes -> True, PlotRange -> All] QQ2 = ParametricPlot3D[ ratbezsur[w, dd2][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {20, 20}] Show[P2, QQ2, ViewPoint -> {2., -1., -2}] Q2 = ParametricPlot3D[ ratbezsur[w, dd2][u, v] // Evaluate, {u, 0, 1}, {v, -1, 1}, PlotPoints -> {20, 20}] dd3 = {{ {-3/2, 0, 6}, {-3/2, 3/2, 6}, {0, 3/2, 6}}, { {-1/2, 0, 6}, {-1/2, 1/2, 6}, {0, 1/2, 6}}, {{-1/2, 0, 7/2}, {-1/2, 1/2, 7/2}, {0, 1/2, 7/2}}} P3 = Show[ beznet[dd3], bezpts[dd3], Axes -> True, PlotRange -> All] Q3 = ParametricPlot3D[ ratbezsur[w, dd3][u, v] // Evaluate, {u, 0, 1}, {v, -1, 1}, PlotPoints -> {10, 20}] dd4 = {{ {3/2, 0, 6}, {3/2, 3/2, 6}, {0, 3/2, 6}}, { {1/2, 0, 6}, {1/2, 1/2, 6}, {0, 1/2, 6}}, {{1/2, 0, 7/2}, {1/2, 1/2, 7/2}, {0, 1/2, 7/2}}} P4=Show[ beznet[dd4], bezpts[dd4], Axes -> True, PlotRange -> All] Q4 = ParametricPlot3D[ ratbezsur[w, dd4][u, v] // Evaluate, {u, 0, 1}, {v, -1, 1}, PlotPoints -> {10, 20}] dd5 = {{{-1/2, 0, 7/2}, {-1/2, 1/2, 7/2}, {0, 1/2, 7/2}}, { {-1/2, 0, 1}, {-1/2, 1/2, 1}, {0, 1/2, 1}}, { {-7/2, 0, 1}, {-7/2, 7/2, 1}, {0, 7/2, 1}}} P5 = Show[ beznet[dd5], bezpts[dd5], Axes -> True, PlotRange -> All] Q5 = ParametricPlot3D[ ratbezsur[w, dd5][u, v] // Evaluate, {u, 0, 1}, {v, -1, 1}, PlotPoints -> {10, 20}] dd6 = {{{1/2, 0, 7/2}, {1/2, 1/2, 7/2}, {0, 1/2, 7/2}}, { {1/2, 0, 1}, {1/2, 1/2, 1}, {0, 1/2, 1}}, { {7/2, 0, 1}, {7/2, 7/2, 1}, {0, 7/2, 1}}} P6 = Show[ beznet[dd6], bezpts[dd6], Axes -> True, PlotRange -> All] Q6 = ParametricPlot3D[ ratbezsur[w, dd6][u, v] // Evaluate, {u, 0, 1}, {v, -1, 1}, PlotPoints -> {10, 20}] Show[Q1, Q2, Q3, Q4, Q5, Q6] Show[ beznet[dd1], bezpts[dd1], beznet[dd2], bezpts[dd2], beznet[dd3], bezpts[dd3], beznet[dd4], bezpts[dd4], beznet[dd5], bezpts[dd5], beznet[dd6], bezpts[dd6], Axes -> True, PlotRange -> All, Boxed->False]; (***** マウスの形(?) ******) mm1 ={ { {2, 10, 0}, {2, 6, 0}, {1, 4, 0}, {4, 0, 0}, {2, 0, 0}, {0, 0, 0},{-2, 0, 0}, {-4, 0, 0}, {-1, 4, 0}, {-2, 6, 0}, {-2, 10, 0}}, {{2, 10, 0}, {2, 10.2, 0},{2, 10.5, 0}, {1, 10.8, 0}, {0, 11, 0}, {0, 11.3, 0}, {0, 11, 0}, {-1, 10.8, 0}, {-2, 10.5, 0}, {-2, 10.2, 0},{-2, 10, 0} }} ****** ****** P1 = Show[ beznet[mm1], bezpts[mm1], Boxed -> False, PlotRange -> All] ****** ****** B1=ParametricPlot3D[bezsurface[mm1][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {20, 20}, Boxed -> False] ****** ****** mm2 ={ { {2, 10, 0}, {2, 6, 0}, {1, 4, 0}, {4, 0, 0}, {2, 0, 0},{0, 0, 0}, {-2, 0, 0}, {-4, 0, 0}, {-1, 4, 0}, {-2, 6, 0}, {-2, 10, 0}}, { {2, 10, 1}, {2, 6, 1}, {1, 4, 1}, {4, 2, 5}, {2, 1,5},{0, 1, 6}, {-2,1, 5}, {-4, 2, 5}, {-1, 4, 1}, {-2, 6, 1}, {-2, 10, 1}}, {{2, 10, 1}, {2, 10.2, 1},{2, 10.5, 1}, {1, 10.8, 1}, {0, 11, 1}, {0, 11.3, 1}, {0, 11, 1}, {-1, 10.8, 1}, {-2, 10.5, 1}, {-2, 10.2, 1},{-2, 10, 1} }} ****** ****** P2 = Show[ beznet[mm1], beznet[mm2], bezpts[mm2], Boxed -> False, PlotRange -> All] ****** ****** B2=ParametricPlot3D[bezsurface[mm2][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {20, 20}, Boxed -> False, ViewPoint->{2.656, 1.919, 0.845}] ***** mm3={{{2, 10, 0}, {2, 10.2, 0},{2, 10.5, 0}, {1, 10.8, 0}, {0, 11, 0}, {0, 11.3, 0}, {0, 11, 0}, {-1, 10.8, 0}, {-2, 10.5, 0}, {-2, 10.2, 0},{-2, 10, 0} }, {{2, 10, 1}, {2, 10.2, 1},{2, 10.5, 1}, {1, 10.8, 1}, {0, 11, 1}, {0, 11.3, 1}, {0, 11, 1}, {-1, 10.8, 1}, {-2, 10.5, 1}, {-2, 10.2, 1},{-2, 10, 1} }} ****** P2 = Show[ beznet[mm1], beznet[mm2], beznet[mm3], bezpts[mm2], Boxed -> False, PlotRange -> All] ****** B3=ParametricPlot3D[bezsurface[mm3][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {5, 16}, Boxed -> False, ViewPoint->{2.656, 1.919, 0.845}] Show[ B1, B2, B3, ViewPoint->{3, 3.252, 1.659}] BB1=wireframe[bezsurface[mm1][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {20, 20}, Boxed -> False] BB2=wireframe[bezsurface[mm2][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {20, 20}, Boxed -> False] BB3=wireframe[bezsurface[mm3][u, v] // Evaluate, {u, 0, 1}, {v, 0, 1}, PlotPoints -> {6, 20}, Boxed -> False] Show[B1, BB2, BB3, ViewPoint->{2, 5, 2}]