================
BEZIER5（ベジエ曲面）
================
Mathematicaを用いた練習を行います。

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

以下の説明の中で
＊＊＊＊＊
ではさまれた部分を
Mathematica に入力して下さい。自分で入力するか、コピー＆ペーストして下さい。
ペーストしたものをShift+Enter(Return) キーで実行します。
上から順に，すべて実行するようにしてください。
上で定義した関数を後で利用すること があります。

ベジエ曲面を描くための準備として、まず次を入力してください。
＊＊＊＊＊＊
<<Graphics`Animation` 
(*  <<RealTime3D` *)
＊＊＊＊＊＊ 
少し長いですが、ベジエ曲面、有理ベジエ曲面を描いたり、
ベジエ点・重みを求めるための次のプログラムを
入力（コピーしてペースト）して下さい。

＊＊＊＊＊＊
combination[n_,r_]:=  n!/(r! (n-r)!)

Bernstein[n_,r_,t_]:= combination[n,r] t^r (1-t)^(n-r)

Bernsteinlist[n_,t_]:=Table[Bernstein[n,j,t],{j,0,n}]

Beziertransformation[n_]:= 
Table[Flatten[{Table[0,{j}],Table[combination[i,j]/combination[n,j],{i, j,n}]}],{j,0,n}]

Bezierpoints[polynomiallist_,DegreesElevation___]:=
  Transpose[
    Table[
      Table[Coefficient[polynomiallist[[j]],t,i],{i,0,
    DegreesElevation+Max[Exponent[polynomiallist,t]]}].
        Beziertransformation[DegreesElevation+Max[Exponent[polynomiallist,
            t]]],
      {j,1,Length[polynomiallist]}
      ]
    ]


Projectiveparameterlist[rationalfunctionlist_]:=
  Expand[
    Cancel[
      Flatten[
        {
          Apply[PolynomialLCM,
    Denominator[rationalfunctionlist]]*rationalfunctionlist,Apply[
        PolynomialLCM,Denominator[rationalfunctionlist]]
          }
        ]
      ]
    ]

GenerateRationalControlpoints[lists_]:=  
Table[Delete[lists[[i]]/Last[lists[[i]]],Length[lists[[i]]]],{i,1,Length[lists]}]

GenerateRationalWeights[lists_]:=  Table[Last[lists[[i]]],{i,1,Length[lists]}]

Reflection[list_]:={{1,0,0},{0,1,0},{0,0,-1}}.list

Rotation[angle_,list_]:=
{{Cos[angle],-Sin[angle],0},{Sin[angle],Cos[angle],0}, {0,0,1}}.list

Beziersurface[lists_,u_,v_]:= 
 Bernsteinlist[Dimensions[lists][[2]]-1,v].(Bernsteinlist[Dimensions[lists][[1]]-1,u].lists)

Beziersurfacegraphics[lists_,options___]:=
  Show[
    ParametricPlot3D[Evaluate[Beziersurface[lists,u,v]],{u,0,1},{v,0,1},
      AxesLabel -> {"x","y","z"},AspectRatio -> Automatic,options],
    Graphics3D[
      {
        Table[Line[Transpose[lists][[i]]],{i,1,Dimensions[lists][[2]]}],
        Table[Line[lists[[j]]],{j,1,Dimensions[lists][[1]]}]
        }
      ]
    ]

ColorBeziersurfacegraphics[lists_,color_,options___]:=
  Show[
    ParametricPlot3D[Evaluate[Beziersurface[lists,u,v]],{u,0,1},{v,0,1},
      AxesLabel -> {"x","y","z"},AspectRatio -> Automatic,options],
    Graphics3D[
      {
        {color,
          Table[Line[Transpose[lists][[i]]],{i,1,Dimensions[lists][[2]]}],
          Table[Line[lists[[j]]],{j,1,Dimensions[lists][[1]]}]
          }
        }
      ]
    ]

WeightBezierlistsforsurface[lists_,weightlist_]:=
Table[ lists[[i,j]]*weightlist[[i,j]],{i,1,Dimensions[weightlist][[1]]},{j,1,Dimensions[weightlist][[2]]}]

RationalBeziersurface[lists_,weightlists_,u_,v_]:= 
Beziersurface[WeightBezierlistsforsurface[lists,weightlists],u,v]/Beziersurface[weightlists,u,v]

RationalBeziersurfacegraphics[lists_,weightlists_,options___]:=
  Show[
    ParametricPlot3D[
      Evaluate[RationalBeziersurface[lists,weightlists,u,v]]
      ,{u,0,1},{v,0,1},
      AxesLabel -> {"x","y","z"},
      AspectRatio -> Automatic,
      options
      ],
    Graphics3D[
      {
        Table[Line[Transpose[lists][[i]]],{i,1,Dimensions[lists][[2]]}],
        Table[Line[lists[[j]]],{j,1,Dimensions[lists][[1]]}]
        }
      ]
    ]

ColorRationalBeziersurfacegraphics[lists_,weightlists_,color_,options___]:=
  Show[
    ParametricPlot3D[
      Evaluate[RationalBeziersurface[lists,weightlists,u,v]]
      ,{u,0,1},{v,0,1},
      AxesLabel -> {"x","y","z"},
      AspectRatio -> Automatic,
      options
      ],
    Graphics3D[
      {
        {color,
          Table[Line[Transpose[lists][[i]]],{i,1,Dimensions[lists][[2]]}],
          Table[Line[lists[[j]]],{j,1,Dimensions[lists][[1]]}]
          }
        }
      ]
    ]
    
    Beziercoefficients[polynomiallist_, i_] :=
  Transpose[Beziertransformation[Max[Exponent[polynomiallist, u]]]].
    Table[
      If[p + q == 0, polynomiallist[[i]] /. u -> 0 /. v -> 0,
        Coefficient[polynomiallist[[i]], (u^p)*(v^q)] /. u -> 0 /. v -> 0],
      {p, 0, Max[Exponent[polynomiallist, u]]},
      {q, 0, Max[Exponent[polynomiallist, v]]} ].
    Beziertransformation[Max[Exponent[polynomiallist, v]]]

Beziersurfacepoints[polynomiallist_] := Simplify[
    Partition[
      Transpose[
        Table[
          Flatten[Beziercoefficients[polynomiallist, i], 1],
          {i, 1, Length[polynomiallist]}]
        ],
      Max[Exponent[polynomiallist, v]] + 1
      ]
    ]

RationalBeziersurfaceControlpoints[rationalfunctionlist_] :=
  Simplify[
    Partition[
      GenerateRationalControlpoints[
        Transpose[
          Table[
            Flatten[
      Beziercoefficients[Projectiveparameterlist[rationalfunctionlist], i], 
        1],
            {i, 1, Length[Projectiveparameterlist[rationalfunctionlist]]}]
          ]
        ],
      Max[Exponent[Projectiveparameterlist[rationalfunctionlist], v]] + 1
      ]
    ]
    
 RationalBeziersurfaceWeights[rationalfunctionlist_] :=
  Simplify[
    Partition[
      GenerateRationalWeights[
        Transpose[
          Table[
            Flatten[
      Beziercoefficients[Projectiveparameterlist[rationalfunctionlist], i], 
        1],
            {i, 1, Length[Projectiveparameterlist[rationalfunctionlist]]}]
          ]
        ],
      Max[Exponent[Projectiveparameterlist[rationalfunctionlist], v]] + 1
      ]
    ]

ColorBeziernet[lists_, color_] :=
  Show[
    Graphics3D[
      {
        {color,
          Table[Line[Transpose[lists][[i]]], {i, 1, Dimensions[lists][[2]]}],
          Table[Line[lists[[j]]], {j, 1, Dimensions[lists][[1]]}]
          }
        }, Axes -> True, AxesLabel -> {"x", "y", "z"}
      ]
    ]
＊＊＊＊＊＊

双曲放物面のベジエ点を求めてみましょう。
＊＊＊＊＊＊
SS1[u_,v_]:={u,v,u*v} 
Beziersurfacepoints[SS1[2u-1,2v-1]]
＊＊＊＊＊＊

＊＊＊＊＊＊
ColorBeziersurfacegraphics[ Beziersurfacepoints[SS1[2u-1,2v-1]],  Blue] 
＊＊＊＊＊＊ 
＊＊＊＊＊＊
ColorBeziernet[Beziersurfacepoints[SS1[2u-1,2v-1]],  Red]
＊＊＊＊＊＊

放物線から生成された曲面のベジエ点を求めてみましょう。
＊＊＊＊＊＊
SS2[u_,v_]:={u,v,u^2} 
Beziersurfacepoints[SS2[2u-1,2v-1]]
＊＊＊＊＊＊

＊＊＊＊＊＊
ColorBeziersurfacegraphics[ Beziersurfacepoints[SS2[2u-1,2v-1]],  Blue] 
＊＊＊＊＊＊ 
＊＊＊＊＊＊
ColorBeziernet[Beziersurfacepoints[SS2[2u-1,2v-1]],  Red]
＊＊＊＊＊＊

放物面のベジエ点を求めてみましょう。
＊＊＊＊＊＊
SS3[u_,v_]:={u,v,u^2 +v^2} 
Beziersurfacepoints[SS3[2u-1,2v-1]]
＊＊＊＊＊＊

＊＊＊＊＊＊
ColorBeziersurfacegraphics[ Beziersurfacepoints[SS3[2u-1,2v-1]],  Blue] 
＊＊＊＊＊＊ 
＊＊＊＊＊＊
ColorBeziernet[Beziersurfacepoints[SS3[2u-1,2v-1]],  Red]
＊＊＊＊＊＊


異なった表し方の双曲放物面のベジエ点を求めてみましょう。
＊＊＊＊＊＊
SS4[u_,v_]:={u, v, u^2 - v^2} 
Beziersurfacepoints[SS4[2u-1,2v-1]]
＊＊＊＊＊＊

＊＊＊＊＊＊
ColorBeziersurfacegraphics[ Beziersurfacepoints[SS4[2u-1,2v-1]],  Blue] 
＊＊＊＊＊＊ 
＊＊＊＊＊＊
ColorBeziernet[Beziersurfacepoints[SS4[2u-1,2v-1]],  Red]
＊＊＊＊＊＊

ベジエ点を与えて、ベジエ曲面を描いてみましょう。
＊＊＊＊＊＊
cpoly={{{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}}}  
ColorBeziersurfacegraphics[ cpoly,  Blue, PlotRange->All] 
ColorBeziernet[cpoly,  Red] 
＊＊＊＊＊＊
ベジエ曲面のベルンシュタイン多項式による表現を求めておく。
＊＊＊＊＊＊
Beziersurface[cpoly,u,v] 
%
＊＊＊＊＊＊

同じように、次のベジエ点のベジエ曲面を描いてみましょう。
また、ベジエ曲面のベルンシュタイン多項式による表現を求ましょう。
＊＊＊＊＊＊
cpoly ={{{-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}}} 
＊＊＊＊＊＊

＊＊＊＊＊＊
cpoly ={{{2,0, 2},{4, 0, 2}, {4, 0, 0}}, 
{{2, 2, 2}, {4, 4, 2}, {4, 4, 0}},  {{0, 2, 2}, {0, 4, 2}, {0, 4, 0}}}  
＊＊＊＊＊＊

＊＊＊＊＊＊
cpoly = {{{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}}} 
＊＊＊＊＊＊

＊＊＊＊＊＊
cpoly ={{{-22,-22, 0},{10, -18, 16}, {10, 18, 16}, {-22, 22, 0}},
 {{-18, 10, -16},{-22/3,-22/3, 0}, {-22/3,22/3, 0}, {-18, -10, -16}},  
{{18, 10, -16},{22/3,-22/3, 0}, {22/3,22/3, 0}, {18, -10, -16}},
{{22,-22, 0},{-10, -18, 16}, {-10, 18, 16}, {22, 22, 0}}} 
＊＊＊＊＊＊

ベジエ点と重みを与えて、有理ベジエ曲面を描いてみましょう。

８分の１球面を描いてみましょう
＊＊＊＊＊＊
cpoly = {{{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}}} 
weight ={{1,1,2}, {1,1,2}, {2,2,4}}
＊＊＊＊＊＊
＊＊＊＊＊＊
ColorRationalBeziersurfacegraphics[cpoly, weight , Blue]
ColorBeziernet[cpoly,  Red] 
＊＊＊＊＊＊
＊＊＊＊＊＊
RationalBeziersurface[cpoly,weight,u,v]//Simplify
＊＊＊＊＊＊

トーラスの一部を描いてみましょう。
＊＊＊＊＊＊
cpoly1 = {{{8, 0, 0},{8, 0, 2}, {6, 0, 2}}, 
{{8, 8, 0},{8, 8, 2}, {6, 6, 2}},  {{0, 8, 0},{0, 8, 2}, {0, 6, 2}}} 
weight1 ={{1,1,2}, {1,1,2}, {2,2,4}} 
＊＊＊＊＊＊
＊＊＊＊＊＊
A1=ColorRationalBeziersurfacegraphics[cpoly1, weight1, Blue]
ColorBeziernet[cpoly1,  Red] 
＊＊＊＊＊＊
＊＊＊＊＊＊
RationalBeziersurface[cpoly1, weight1, u, v]//Simplify
＊＊＊＊＊＊
＊＊＊＊＊＊
cpoly2 = {{{4, 0, 0},{4, 0, 2}, {6, 0, 2}}, 
{{4, 4, 0},{4, 4, 2}, {6, 6, 2}},  {{0, 4, 0},{0, 4, 2}, {0, 6, 2}}} 
weight2 ={{1,1,2}, {1,1,2}, {2,2,4}} 
＊＊＊＊＊＊
＊＊＊＊＊＊
A2=ColorRationalBeziersurfacegraphics[cpoly2, weight2 , Blue]
ColorBeziernet[cpoly2,  Red] 
＊＊＊＊＊＊
＊＊＊＊＊＊
RationalBeziersurface[cpoly2, weight2, u, v]//Simplify
＊＊＊＊＊＊
＊＊＊＊＊＊
Show[A1, A2] 
＊＊＊＊＊＊

１葉の双曲面の一部を描いてみましょう．
＊＊＊＊＊
cpoly = {{{-3,0, -2*2^(1/2)},{-3, 3, -2*2^(1/2)}, {0, 3, -2*2^(1/2)}}, 
{{-1/3, 0, 0}, {-1/3, 1/3, 0}, {0, 1/3, 0}},  
{{-3,0, 2*2^(1/2)},{-3, 3, 2*2^(1/2)}, {0, 3, 2*2^(1/2)}}}
weight ={{1,1,2}, {3,3,6}, {1,1,2}};
＊＊＊＊＊＊
ColorRationalBeziersurfacegraphics[cpoly, weight, Blue]
ColorBeziernet[cpoly,  Red] 
＊＊＊＊＊＊
＊＊＊＊＊＊
RationalBeziersurface[cpoly, weight, u, v]//Simplify
＊＊＊＊＊＊


いろいろなベジエ点を与えて，ベジエ曲面を作って見ましょう 。
何かおもしろそうなものができたら、教えて下さい。