【B02】ベジェ曲面
この記事はアドカレに参加しています。
ベジェ曲線
n次のベジェ曲線は制御点$P$が与えられたとき、シグマを用いて以下のような式で表わせるのでした。
p=\sum_{i=0}^{n} {}_nC_i B_{i}^n(t)p_{i}
このとき、
\begin{align}
B_{i}^n(t)&=t^i(1-t)^{n-i} \\
{}_nC_i&=\frac{n!}{i!(n-i)!}
\end{align}
となります。
ベジェ曲面
さて、ベジェ曲線のように制御点を用いて面も変形できると便利ですよね。
制御点を用いて面の関数を変形できるものに、ベジェ曲面というものがあります。
ベジェ曲面は外側の四点(図ではp00,p03,p30,p33)を通る面の関数であり、外側の四点の間にある制御点で面の形状を決めます。
また、図で言えば横に四本、縦に四本のベジェ曲線を更に面として補間しているのが分かると思います。ベジェ曲線は黄色の線です。ベジェ曲線は三次なので、4つの点の組み合わせは、
p00,p01,p02,p03、
p10,p11,p12,p13、
p20,p21,p22,p23、
p30,p31,p32,p33、
p00,p10,p20,p30、
p01,p11,p21,p31、
p02,p12,p22,p32、
p03,p13,p23,p33
の八組です。
式でみてみましょう。
横方向の位置を表す引数uと縦方向の位置を表す引数vを用いて、ベジェ曲面は以下のような関数で表現できます。
p=f(u,v) (0 \leqq u \leqq 1,0 \leqq v \leqq 1)
図では、
(u,v)=(0,0)でp00
(u,v)=(1,0)でp03
(u,v)=(0,1)でp30
(u,v)=(1,1)でp33
のようになります。
f(x)を求めるには、まず横方向にベジェ曲線として補間をします。横方向のベジェ曲線をn次としたとき、
p_{j}=\sum_{i=0}^{n} {}_nC_i B_{i}^n(u)p_{j,i}
のようになります。
更に、縦方向のベジェ曲線をm次としたとき、横方向のベジェ曲線はm+1本あることになります。横方向にあるm+1本のベジェ曲線を、縦方向に補間します。
\begin{align}
p&=\sum_{j=0}^{m} {}_mC_j B_{j}^m(v) p_{j} \\
p&=\sum_{j=0}^{m} {}_mC_j B_{j}^m(v) (\sum_{i=0}^{n} {}_nC_i B_{i}^n(u) p_{j,i})
\end{align}
シグマは、係数を外に出したり中に入れたりできるので、式を整理します。
\begin{align}
p&=\sum_{j=0}^{m} {}_mC_j B_{j}^m(v) \sum_{i=0}^{n} {}_nC_i B_{i}^n(u) p_{j,i} \\
p&=\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) p_{j,i} \\
\end{align}
はい。n×m次元のベジェ曲面を求めることができました。
3×3のベジェ曲面
n×m次元のベジェ曲面は以下のように表すことができたのでした。
p=\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m {}_nC_i B_{i}^n p_{j,i} \\
ここで、3×3次元のベジェ曲面について考えてみましょう。
三次のベジェ曲線は、
\begin{align}
T&=\begin{bmatrix} t^3 & t^2 & t & 1 \end{bmatrix} \\
A&=\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix} \\
P&=\begin{bmatrix} p0 \\ p1 \\ p2 \\ p3 \end{bmatrix}
\end{align}
として、
p=TAP
のように表すことができるのでした。
この行列を使用して3×3のベジェ曲面についてみていきましょう。
まずは横方向に補間します。
\begin{align}
P_{j}&=\begin{bmatrix} p_{j,0} \\ p_{j,1} \\ p_{j,2} \\ p_{j,3} \end{bmatrix} \\
U&=\begin{bmatrix} u^3 & u^2 & u & 1 \end{bmatrix}
\end{align}
として、
p_{j}=UAP_{j}
とすることができます。縦方向にベジェ曲線はm+1本あるので、全て書くと以下のようになります。
\begin{align}
p_{0}&=UAP_{0} \\
p_{1}&=UAP_{1} \\
p_{2}&=UAP_{2} \\
p_{3}&=UAP_{3}
\end{align}
次に縦方向に補間します。
\begin{align}
P_{y}&=\begin{bmatrix} p_{0} \\ p_{1} \\ p_{2} \\ p_{3} \end{bmatrix} \\
V&=\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\end{align}
として、
\begin{align}
p&=VAP_{y} \\
p&=
\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} p_{0} \\ p_{1} \\ p_{2} \\ p_{3} \end{bmatrix} \\
p&=
\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} UAP_{0} \\ UAP_{1} \\ UAP_{2} \\ UAP_{3} \end{bmatrix}
\end{align}
ここで、
\begin{align}
UAP&=P^T A U^T \\
\begin{bmatrix} u^3 & u^2 & u & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} p_{j,0} \\ p_{j,1} \\ p_{j,2} \\ p_{j,3} \end{bmatrix}
&=
\begin{bmatrix} p_{j,0} & p_{j,1} & p_{j,2} & p_{j,3} \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} u^3 \\ u^2 \\ u \\ 1 \end{bmatrix}
\end{align}
のように置き換えることができるので、
\begin{align}
p&=
\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} UAP_{0} \\ UAP_{1} \\ UAP_{2} \\ UAP_{3} \end{bmatrix} \\
p&=
\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} P_{0}^T A U^T \\ P_{1}^T A U^T \\ P_{2}^T A U^T \\ P_{3}^T A U^T \end{bmatrix} \\
p&=
\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} P_{0}^T \\ P_{1}^T \\ P_{2}^T \\ P_{3}^T \end{bmatrix} A U^T \\
p&=
\begin{bmatrix} v^3 & v^2 & v & 1 \end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
P_{0,0} & P_{0,1} & P_{0,2} & P_{0,3} \\
P_{1,0} & P_{1,1} & P_{1,2} & P_{1,3} \\
P_{2,0} & P_{2,1} & P_{2,2} & P_{2,3} \\
P_{3,0} & P_{3,1} & P_{3,2} & P_{3,3}
\end{bmatrix}
\begin{bmatrix}
-1 & 3 & -3 & 1 \\
3 & -6 & 3 & 0 \\
-3 & 3 & 0 & 0 \\
1 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} u^3 \\ u^2 \\ u \\ 1 \end{bmatrix} \\
\end{align}
制御点の行列を
P=
\begin{bmatrix}
P_{0,0} & P_{0,1} & P_{0,2} & P_{0,3} \\
P_{1,0} & P_{1,1} & P_{1,2} & P_{1,3} \\
P_{2,0} & P_{2,1} & P_{2,2} & P_{2,3} \\
P_{3,0} & P_{3,1} & P_{3,2} & P_{3,3}
\end{bmatrix}
として、
p=VAPAU^T
これが3×3のベジェ曲面になります。
ベジェ曲面を座標系に落とし込む
n次のベジェ曲面の式は以下のようになるのでした。
p=\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) p_{j,i}
ですが、これは一次元の座標系での式となります。
二次元平面でのベジェ曲線は以下の式で表すことができます。
\begin{bmatrix} x \\ y \end{bmatrix}=
\begin{bmatrix}
\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) x_{j,i}
\\
\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) y_{j,i}
\end{bmatrix}
同様にして、三次元空間では以下のようになります。
\begin{bmatrix} x \\ y \\ z \end{bmatrix}=
\begin{bmatrix}
\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) x_{j,i}
\\
\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) y_{j,i}
\\
\sum_{j=0}^{m} \sum_{i=0}^{n} {}_mC_j B_{j}^m(v) {}_nC_i B_{i}^n(u) z_{j,i}
\end{bmatrix}
プログラム例
三次のベジェ曲面のプログラム例です。
BezierLine_M.hppではベジェ曲線を求める関数を定義しておきます。
BezierLinePは制御点から求めます。
BezierLineCは三次関数の係数から求めます。
BezierCoefficientは制御点から三次関数の係数を求めます。
/*
BezierLine_M.hpp
*/
#pragma once
//return p[0]*s^3 + t(3*p[1]*s^2 + t(3*p[2]*s + t*p[3]))
double BezierLineP(double* p, double t) {
double s = 1 - t;
return *(p + 0) * s * s * s + t * (3 * *(p + 1) * s * s + t * (3 * *(p + 2) * s + t * *(p + 3)));
//(3*s*s*x1+(3*s*x2+t)*t)*t
}
//return ((at+b)t+c)t+d
double BezierLineC(double* c, double t) {
return ((*(c + 0) * t + *(c + 1)) * t + *(c + 2)) * t + *(c + 3);
}
//p->c
void BezierCoefficient(double* p, double* c) {
*(c + 0) = -1 * *(p + 0) + 3 * *(p + 1) - 3 * *(p + 2) + 1 * *(p + 3);
*(c + 1) = 3 * *(p + 0) - 6 * *(p + 1) + 3 * *(p + 2);
*(c + 2) = -3 * *(p + 0) + 3 * *(p + 1);
*(c + 3) = 1 * *(p + 0);
}
BezierSurface_M.hppではベジェ曲面を求めるクラスBezierSurface_Mを定義しています。
コンストラクタでは16個の制御点から横方向のベジェ曲線の三次関数の係数を求めています。これは、三次関数の係数から値を求めたほうが計算回数が少なくて済むからです。(BezierLineCは掛け算3回、BezierLinePは掛け算11回)
関数runは、ベジェ曲面の関数です。
横方向の補間の結果を配列aに入れておき、縦方向に補間して値を返しています。
式で表すと、
p=VAPAU^T
という式で$PA=C$として、
p=VACU^T
というように計算しているということです。
実際の変数に当てはめて考えてみます。
横方向の三次関数の係数を
\begin{align}
C_{0}&=A \begin{bmatrix} p_{0,0} \\ p_{0,1} \\ p_{0,2} \\ p_{0,3} \end{bmatrix}
\\
C_{1}&=A \begin{bmatrix} p_{1,0} \\ p_{1,1} \\ p_{1,2} \\ p_{1,3} \end{bmatrix}
\\
C_{2}&=A \begin{bmatrix} p_{2,0} \\ p_{2,1} \\ p_{2,2} \\ p_{2,3} \end{bmatrix}
\\
C_{3}&=A \begin{bmatrix} p_{3,0} \\ p_{3,1} \\ p_{3,2} \\ p_{3,3} \end{bmatrix}
\end{align}
として、
\begin{align}
a_{0}&=UC_{0} \\
a_{1}&=UC_{1} \\
a_{2}&=UC_{2} \\
a_{3}&=UC_{3} \\
return&=VA \begin{bmatrix} a_{0} \\ a_{1} \\ a_{2} \\ a_{3} \end{bmatrix}
\end{align}
/*
BezierSurface_M.hpp
*/
#pragma once
#include "BezierLine_M.hpp"
class BezierSurface_M {
double c[4][4];
public:
BezierSurface_M(double* p) {
BezierCoefficient(p + 0, c[0]);
BezierCoefficient(p + 4, c[1]);
BezierCoefficient(p + 8, c[2]);
BezierCoefficient(p + 12, c[3]);
}
double run(double u, double v) {
double a[4];
a[0] = BezierLineC(c[0], u);
a[1] = BezierLineC(c[1], u);
a[2] = BezierLineC(c[2], u);
a[3] = BezierLineC(c[3], u);
return BezierLineP(&a[0], v);
}
};
むすび
ベジェ曲面についてみてきました。行列しんどいですね。