CGとCADの数理

Author

gentam

Status
archived
Date
2018-04-23
Category

lecturenote

Tag

cg, lecture

Updated

2018-07-22

CGとCADの数理 / Geometric Modeling and Computer Graphics (中野亜希人) の授業ノート

Contents

シラバス

科目概要

1963年にアイバン・サザランドが開発したSketchpadはコンピュータと人間との対話方式 を革新的に変えた、CG/CADシステム(以下、CG/CAD)の先駆けです。コンピュータの技術 的発展と共にCG/CADは、現実世界の光や剛体、流体の物理シミュレーションを視覚化する 媒体に留まらず、仮想世界のゲームやアート、映画において独自の表現を発明してきまし た。本講義「CGとCADの数理」では、このような高度な応用への深化を考慮しながら、2/3 次元CG/CADの数理的な理解と汎用的なアルゴリズムによる実践を主題とします。

主題と目標/授業の手法など

1963年にアイバン・サザランドが開発したSketchpadはコンピュータと人間との対話方式 を革新的に変えた、CG/CADシステム(以下、CG/CAD)の先駆けです。コンピュータの技術 的発展と共にCG/CADは、現実世界の光や剛体、流体の物理シミュレーションを視覚化する 媒体に留まらず、仮想世界のゲームやアート、映画において独自の表現を発明してきまし た。本講義「CGとCADの数理」では、このような高度な応用への深化を考慮しながら、2/3 次元CG/CADの数理的な理解と汎用的なアルゴリズムによる実践を主題とします。

プログラミング言語の文法よりもCG/CADの数理的な理解と汎用的なアルゴリズムによる実 践に注力するため、講義では直感的なプログラミング文法を持つオープンソースの開発環 境 Processing を使用します。プログラミングの初心者の方も歓迎します。

講義は、毎回の冒頭でプログラムの雛形を配布し、スライドによるCG/CADの数理の解説と 担当教員によるライブコーディングを往来して進めます。

教材・参考文献: 鳥谷浩志, 千代倉弘明. 3次元CADの基礎と応用, 1995, 228p.

第1回 CG/CADの歴史

第1回では、研究や作品例を紹介しながらCG/CADの歴史を概観し、その誕生と隆盛、今 後の展開に目を向けます。後半に、本講義を履修する上での注意事項、本講義の概要、 成績の付け方、Processingの導入方法を説明します。

インロトダクション

  • この授業の目的は映画や最新のゲームで見られるようなリッチCGを最新の技術を使って 作れるようになることではない.

  • CGの基礎理論を数式から始めて実装していき理解することが趣旨.よって出力される絵 はシンプルなものしかない.

歴史

第2回 プリミティブ形状

CG/CADにおけるプリミティブ形状は、頂点と稜線と面で構成されたポリゴンで構成され ます。後半では、ポリゴンを組み合わせて正多面体を描画します。

Processingでは line(x0, y0, x1, y2); で線が引けてしまうが,それを座標x,yのピ クセルの色を設定する(set(x, y, color))だけで実現するにはどうすればいいか?

ブレゼンハムのアルゴリズム を使って描写する.

Bresenham's line algorithm

与えられた始点と終点の間に連続した点を置き,近似的な直線を引くためのアルゴリズ ム.コンピュータのディスプレイに直線を描画するのによく使われ,整数の加減算とビ ットシフトのみで実装できる.コンピュータグラフィックスの分野の最初期のアルゴリ ズムの1つである.

直感的な説明: 直線の傾きが1であれば,x0,y0からx,yを1ずつインクリメントしながら点 を打っていけばよく,傾きが1/2の場合は,xを2回インクリメントする度にyを1インクリ メントすればいい.これを,直線が次のピクセルと交わる点が中心より上か下かによって 判定して進めていく.

このアルゴリズムの基本的な形では \(x \ge 0, y \ge 0\) かつ \(\frac{\varDelta y}{\varDelta x} < 1\) という制限がある.

bresenham.pde

class CVertex {
  int x, y;
  CVertex () { x = 0; y = 0; }
}

CVertex v0, v1;

size(400, 400);
background(40);

v0 = new CVertex();
v0.x = 20;
v0.y = 20;
v1 = new CVertex();
v1.x = 200;
v1.y = 100;

int dx = v1.x - v0.x;
int dy = v1.y - v0.y;
int E = -dx; // 2
int x, y;

x = v0.x;
y = v0.y;

for( ; x < v1.x; x++ ) {
  set(x, y, color(255));
  E += 2 * dy; // 3
  if(E >= 0) { // 1
    y++;
    E -= 2 * dx; // 4
  }
}

演習課題: 基本アルゴリズムを全方位へ線を引くことができるように拡張する.

最もシンプルに考えれば,線を引く目標点(x1, y1)の符号と線の傾きによって,グラフを 8等分の領域に考えて場合分けをしてやればよい.

それを,

ことによって解決するのが以下のコード

bresenham-opt1.pde

void draw () {
  background(40);
  v1.x = mouseX;
  v1.y = mouseY;

  int dx = abs(v1.x - v0.x);
  int dy = abs(v1.y - v0.y);
  if ((steep = dx < dy)) { // swap x and y for v1
    int tmp = v1.x;
    v1.x = v1.y;
    v1.y = tmp;
    dx = abs(v1.x - v0.x);
    dy = abs(v1.y - v0.y);
  }

  int x, y;
  int sx = (v0.x < v1.x) ? 1 : -1;
  int sy = (v0.y < v1.y) ? 1 : -1;
  int E = dy - dx/2;

  for (x = v0.x, y = v0.y; x != v1.x; x += sx) {
    if (steep) set(y, x, color(255));
    else       set(x, y, color(255));
    E += dy;
    if (E > 0) {
      y += sy;
      E -= dx;
    }
  }
}

さらに \(\varDelta y\) を負の値, \(\varDelta x\) を正の値として共通の変 数に足していけば,傾きに関係なくx,y両軸それぞれを誤差の蓄積に合わせて進めていく ことができる.こういうことをすっと思いつけるようになりたい... (Wikipedia参考)

bresenham-opt2.pde

void draw () {
  background(40);
  v1.x = mouseX;
  v1.y = mouseY;
  int x = v0.x,
      y = v0.y;

  int dx =  abs(v1.x - v0.x),
      dy = -abs(v1.y - v0.y);
  int sx = (v0.x < v1.x) ? 1 : -1,
      sy = (v0.y < v1.y) ? 1 : -1;
  int err = dx + dy,
      E;

  for (;;) {
    set(x, y, color(255));
    if (x == v1.x && y == v1.y) break;
    E = err;
    if (E > dy) { err += dy; x += sx; }
    if (E < dx) { err += dx; y += sy; }
  }
}

Note

Processing IDEの編集画面で書くのは疲れるので,vim-processing <https://github.com/sophacles/vim-processing> を入れて nnoremap <Leader>m :make<CR> などと設定しておくと,vimで編集 → <Leader>m で結果を確認.とい う快適なワークフローができる.(Processing.app > Tools > Install "processing-java" が必要)

第3回 ベクトル

3次元空間の位置は一般に(x,y,z)の3次元ベクトルで表されます。第3回では、線形代数 の基礎としてベクトルについて概説します。その理解を深めるため、CG/CADへの利用を 念頭とした3次元ベクトルの四則演算を実装します。

vector.pde

四則演算

class CVector {
  float x, y, z;
  CVector () {
    x = 0.0; y = 0.0; z = 0.0;
  }

  void add (CVector Q, CVector result) {
    result.x = x + Q.x;
    result.y = y + Q.y;
    result.z = z + Q.z;
  }

  void sub (CVector Q, CVector result) {
    result.x = x - Q.x;
    result.y = y - Q.y;
    result.z = z - Q.z;
  }

  void mult (float n, CVector result) {
    result.x = x * n;
    result.y = y * n;
    result.z = z * n;
  }

  void div (float n, CVector result) {
    // if (n == 0) {...}
    result.x = x / n;
    result.y = y / n;
    result.z = z / n;
  }
  ...
}

内積 (Inner Product)

float mag () {
  return (float) sqrt(x*x + y*y + z*z);
}

void normalize (CVector result) {
  float m = mag();
  if (m != 0)
    div(m, result);
}

float dot (CVector v) {
  return x*v.x + y*v.y + z*v.z;
}

射影ベクトル (Projection)

Pを別ベクトルQに対して2つの並行な成分 \(proj_{q}P\) (x軸に並行な成分), \(perp_{q}P\) (y軸に並行な成分) に分解したい.

  • \(proj_{q}P = \frac{P \cdot Q}{|Q|^2}Q\)

  • \(|P|\,\cos\alpha\)

  • \(proj_{q}P + perp_{q}P = P\) より \(perp_{q}P = P - proj_{q}P\)

void proj (CVector Q, CVector result) {
  float d = dot(Q);
  float Q2 = Q.x*Q.x + Q.y*Q.y + Q.z*Q.z;
  Q.mult(d/Q2, result);
}

void perp (CVector Q, CVector result) {
  proj(Q, result);
  sub(result, result);
}

第4回 ベジエ曲線

パラメトリックな曲線表現としてまずベジエ曲線を解説します。Bernstein既定関数や 、その二項定理からの導出、閉包性といったベジエ曲線の数理を解説します。

bezier-cubic.pde

制御点と制御ポリゴン

  • 制御点 control point

  • 制御ポリゴン control polygon

  • 2次ベジェ曲線 Quadratic Bézier curve

  • 3次ベジエ曲線 Cubic Bézier Curve

3次ベジエ曲線は,4つの制御点からなる. ベジエ曲線では必ず,「制御点の数-1」字 数のベジエ曲線が生成される. また,生成されたベジエ曲線は必ず制御ポリゴンの内側にあり,これを凸包性という.

processing では bezier(x0, y0, x1, y1, x2, y2, x3, y3); で簡単に書くことがで きる.この中では何が行われているのかを理解する.

ベルンシュタイン基底関数

Bernstein Basis Function

n次のベジエ曲線は以下の式で計算される:

\(\mathbf{R}(t) = \sum_{i=0}^{n} \mathbf{B}_i^n (t) P_i\) ただし \(0 \le t \le 1\)

ここで, \(\mathbf{B}_i^n(t) = {n \choose i} (1 - t)^{n-i} t^i\) がベルンスタイン基底関数である.

  • n次 → (n+1)個の制御点

  • \(\binom n i\)\(_nC_i\) と同じ

順番にtを当てはめて行くと,このベルンスタイン基底関数がやっていることは単純な二 項定理であることがわかる.

これをコードで実現すると:

class BezierCurve {

  PVector P0, P1, P2, P3;
  PVector[] R;
  int tn;

  BezierCurve () {
    P0 = new PVector(...);
    P1 = ...
    tn = 100;
    R = new PVector[tn];
  }

  void draw () {
    int   i, tt;
    float t = 0.0;
    float ts = 1.0 / tn;
    float B30t, B31t, B32t, B33t;

    ...

    noFill();
    stroke(255);
    for(tt = 0; tt < tn; tt++) {
        B30t =     (1-t) * (1-t) * (1-t);
        B31t = 3 * (1-t) * (1-t)         * t;
        B32t = 3 * (1-t)                 * t * t;
        B33t =                             t * t * t;
        R[tt] = new PVector();
        R[tt].x = B30t * P0.x +
                  B31t * P1.x +
                  B32t * P2.x +
                  B33t * P3.x;
        R[tt].y = B30t * P0.y +
                  B31t * P1.y +
                  B32t * P2.y +
                  B33t * P3.y;
      if (tt != 0)
        line(R[tt-1].x, R[tt-1].y, R[tt].x, R[tt].y);
      t += ts;
    }
  }
}

上記から分かるように,実際の線は line() 関数のみで描写されている.(そして, line関数にはブレゼンハムのアルゴリズムを使うこともできる!)

3次のベジエ曲線ならprocessingの組み込み関数を使えばいいが,上に書いたプログ ラムを拡張すれば4次以上のベジエ曲線を書くことができる.

課題:

  1. 最後の点が繋がっていないバグをなおす

  2. 4次, 5次のベジエ曲線を書けるようにする

  3. N次のベジエ曲線を書けるようにする

Hint: 係数はパスカルの三角形になっている

  • 2次のベジエ曲線の係数は, 1, 2, 1

  • 3次のベジエ曲線の係数は, 1, 3, 3, 1

  • 4次のベジエ曲線の係数は,1, 4, 6, 4, 1

バグ修正とN次のベジエ曲線の実装をしたもの:

bezier-curve.pde

操作説明:

  • 左クリックで新しい制御点を追加

  • 右クリックで制御点を削除

  • ドラッグ&ドロップで制御点を移動

実際に曲線を書く部分の関数.係数はbinominal関数を別に作って求めている:

void bezierN (int n) {
    noFill();
    stroke(255, 255, 255);
    int tt;
    float t = 0.0;
    float ts = 1.0 / tn;

    for (tt = 0; tt < tn ; tt++) {
      float sx = 0,
            sy = 0, b;
      for (int k = 0; k <= n; k++) {
        b = binomial(n, k) * pow((1-t), n-k) * pow(t, k);
        sx += W.get(k).x * b;
        sy += W.get(k).y * b;
      }
      R[tt].set(sx, sy);

      if (tt != 0)
        line(R[tt-1].x, R[tt-1].y, R[tt].x, R[tt].y);
      t += ts;
    }
    line(R[tn-1].x, R[tn-1].y, W.get(n).x, W.get(n).y);
  }

実装する上で次のページを参考にした.直感的には De Casteljau のアルゴリズムを使っ た再帰的なバージョンの方が分かりやすかったので,時間があったらそっちも実装したい.

参考:

第5回 曲率

曲線の曲がり具合の指標となる曲率の計算を、ベジエ曲線を題材に学習します。具体的 には、ベジエ曲線の一階微分、二階微分の計算結果から弧長や動標構を求め、曲率円を 描画します。

curvature.pde

ベジェ曲線上の各点における接線と垂線が表示されている.これは曲線の微分をとること で求められる.

曲率 Curvature

正則な曲線がパラメータtで表されている時,曲率 \(K(t)\) は曲線の一階,二階 微分を使って以下の式で求められる:

\begin{equation*} K(t) = \frac{R'(t).x R''(t).y - R''(t).x R'(t).y}{(R'(t).x R'(t).x + R'(t).y R'(t).y)^{3/2}} \end{equation*}

\(a^{3/2} = \sqrt{a^3} = a\sqrt a\) となるので,これをコードに落とすと:

// Curvature K
K[tt] =  (D1[tt].x * D2[tt].y - D2[tt].x * D1[tt].y) /
         (sqrt(D1[tt].x * D1[tt].x + D1[tt].y * D1[tt].y) *
              (D1[tt].x * D1[tt].x + D1[tt].y * D1[tt].y)
         );

課題: 接線や法線を使ったアルゴリズムの拡張:

曲率円 (Circle of Curvature)

曲率が大きくなる(強く曲がる)ほど半径が小さくなる,曲率が小さくなる(直線に近くな る)ほど,円の半径は大きくなる.

曲率半径は曲率の逆数になる \(R = \frac{1}{K(t)}\)

// Circle of Curvature
stroke(255, 255, 0);
ellipseMode(RADIUS);
temp.x = e2[tt].x; temp.y = e2[tt].y;
temp.normalize();
float KR = 1.0 / K[tt];
temp.mult(KR);
ellipse(R[tt].x + temp.x, R[tt].y + temp.y, KR, KR);

curvature-circle.pde

(これで合ってるかは不明)

第6回 次数上げと接続

より局所的な曲線形状の操作には、既にある曲線の形状を変えずに、曲線の自由度を上 げる必要があります。つまり、曲線の形状を変えずに制御点を増やす必要があります。 後半では、複数への分割やその接続のアルゴリズムも解説します。

次数上げ

degree-elevation.pde

操作説明: "u" を押してる間,次数を上げる

2次のベジェ曲線の形を変えずに次数を上げて,3次のベジェ曲線にすることを考える.

\begin{equation*} \binom 3 0 t^0 (1-t)^3 Q_0 + \binom 3 1 t^1 (1-t)^2 Q_1 + \binom 3 2 t^2 (1-t)^1 Q_2 + \binom 3 3 t^3 (1-t)^0 Q_3 \end{equation*}

それぞれを展開,整理すると以下のようになる.同じ線なので端点の \(Q0, Q3\) は当然変わらない.

  • \(Q_0 = P_0\)

  • \(Q_1 = {P_0 + 2P_1 \over 3}\)

  • \(Q_2 = {2P_1 + P_2 \over 3}\)

  • \(Q_3 = P_2\)

これをコードに落としている部分:

// inside BezierCurve2 class
void degreeElevation (BezierCurve3 b) {
  b.P0.x = P0.x;
  b.P0.y = P0.y;
  b.P1.x = (  P0.x + 2*P1.x)/3;
  b.P1.y = (  P0.y + 2*P1.y)/3;
  b.P2.x = (2*P1.x +   P2.x)/3;
  b.P2.y = (2*P1.y +   P2.y)/3;
  b.P3.x = P2.x;
  b.P3.y = P2.y;
}

分割

split.pde

操作説明: "d" を押してる間,曲線を分割する

de Casteljau's Algorithm

直感的には「点をつないだ線を内分し,その内分した点を繋いだ線を内分する」という作 業を最初は制御点から始めて,再帰的におこなっていく. 最後,1つに収束した点がベジェ曲線上にある. つまり,内分する点(の割合)tを0から1まで動かしていきながら点を打つとベジェ曲線が 描ける.

\(P^r_i (t) = (1-t)P^{r-1}_i(t) + tP^{r-1}_{i+1}(t)\)

ここで

  • \(r = 1, \,\ldots, \,n\)

  • \(i = 0, \,\ldots, \,n-r\)

  • \(r = 1\) のとき \(i = 3-1\):

    • \(P^1_0 (t) = (1-t) P^0_0 (t) + tP^0_1 (t)\)

    • \(P^1_1 (t) = (1-t) P^0_1 (t) + tP^0_2 (t)\)

    • \(P^1_2 (t) = (1-t) P^0_2 (t) + tP^0_3 (t)\)

  • \(r = 2\,\ldots\)

パラメータtで分割(split)するということは,前半と後半のベジェ曲線に分けるというこ と.

De Casteljau's Algorithm | A Primer on Bézier Curves <https://pomax.github.io/bezierinfo/#decasteljau>

課題: インタラクティブ(キーボードやマウス操作)な次数上げや分割を実装する

bezier-decasteljau.pde

操作説明:

  • 制御点を左クリックで追加/右クリックで削除

  • ドラッグ&ドロップで制御点を移動

  • 'd' を押している間,曲線をを分割 (アルゴリズムの過程を表示)

第7回 B-Spline曲線

ベジエ曲線は、制御点の移動が曲線全体に影響を与えますが、B-Spline曲線はこの影 響は局所的なためより細かな形状操作が可能です。第7回では、B-Spline基底関数やノ ットベクトルについて学習し、B-Spline曲線の数理を理解します。

b-spline.pde

次回やるNURBSの基本となる: Non-Uniform Rational B-Spline

B-Spline曲線 (B-Spline Curve)

B-スプライン曲線(B-spline curve)とは、与えられた複数の制御点から定義される滑 らかな曲線である。区分多項式により表現されているため、一部を変更しても曲線全体 に影響は及ばない等の性質がある。ベジェ曲線とともに、コンピュータグラフィックス の世界で広く利用されている。なお、B-splineはBasis spline(Basis=基底)の省略 形である。基本的に曲線は制御点を通らない。(from Wikipedia <https://ja.wikipedia.org/wiki/B-スプライン曲線>)

\(k = 次数+1\)\(n+1\) 個の制御点 \(P_i\) で定義されるB-Spline曲 線は以下の式で表される:

\begin{align*} \mathbf{R}(t) = \sum_{i=0}^{n} N_i^k (t) P_i \\ (0 \le t \le t_{max}) \end{align*}

パラメータtの範囲 \(t_{max}\)\(n-k+2\). ここで \(N\) がB-スプライン基底関数(B-Spline basis function):

\begin{equation*} N_i^k (t) = (t-u_i) \frac{N_i^{k-1}(t)}{u_{i+k-1} - u_i} + (u_{i+k}-t)\frac{N_{i+1}^{k-1}(t)}{u_{i+k} - u_{i+1}} \end{equation*}

\(k = 1\) のとき:

\begin{equation*} N_i^1 = \begin{cases} 1 & \quad\text{if$\quad u_i \le t < u_{i+1}$} \\ 0 & \quad\text{otherwise} \end{cases} \end{equation*}

Nの定義にNが含まれているため,k=1になるまで再帰的に計算をする関数として実装する ことができる.(が授業のプログラムではそれぞれを別に計算している)

ここで \(u\)ノット(knot) と呼ばれている.このノットの作り方によって, open (nonperiodic),close(periodic) なB-Spline曲線を作ることができる.

ノットは \(\mathbf{u}_j = [0, 0, 0, 1, 2, 2, 2]\) (k=3のとき) のように成分を 列挙する形で書かれ,ノットベクトル (Knot Vector) と呼ばれることもある.

\begin{equation*} u_j = \begin{cases} 0, & j < k \\ j-k+1, & k \le j \le n \\ n-k+2, & n < j \end{cases} \end{equation*}

授業のコードでそのまま対応する部分:

int knot(int j, int n, int k) {
  int uj = 0;
  if      (j < k)            uj = 0;
  else if (k <= j && j <= n) uj = j-k+1;
  else if (j > n)            uj = n-k+2;
  return uj;
}

課題: B-Spline曲線の改良.より高次やN次のB-Spline曲線,または分割など.

N次のB-Spline曲線をN個の制御点から描けるようにしたもの:

b-spline-n.pde

操作説明:

  • 制御点を左クリックで追加/右クリックで削除

  • ドラッグ&ドロップで制御点を移動

  • より高い次数の曲線を'u'で追加,'d'で削除

第8回 NURBS曲線

第7回のB-Spline曲線を発展させ、制御点が重み(ウェイト)を持つNon-Uniform Rational B-Splines(NURBS)を解説します。

nurbs.pde

操作説明: キーボードの'w'/'W'で,P1の重みを増減

NURBS

Non-Uniform Rational B-Spline (非一様有理Bスプライン) は、曲線や曲面を生成する ためにコンピュータグラフィックスで一般的に採用される数学的モデルである。その柔 軟性と正確性からモデリング用の形状にも、解析的な用途にも向いている。 (Wikipedia <https://ja.wikipedia.org/wiki/NURBS>)

復習

  • ベジェ曲線: 制御点は曲線全体に影響を与える

  • B-Spline曲線: 曲線の部分制御ができる

(1変数の)有理関数とは: \(\mathbf{R}(t) = \frac{A(t)}{B(t)}\) が多項式(Polynomial) の関数.

ベジェ曲線の有理化

3次ベジェ曲線:

\begin{equation*} R(t) = \sum_{i=0}^3 B_i^3 (t) P_i \end{equation*}

直感的には,特定のtの値に対して,各制御点が持つ重みの割合がベルンシュタイン基底 関数 \(B_i^k(t)\) によって決められる. e.g. \(P_0: 61\%, P_1: 33\%, P_2: 6\%, P_3:0\%\) を足した位置.(これを混ぜ合 わせ関数と呼ぶ)

  • ここでもし重みの総和が1ではないとき,総和が1であるとして割合を計算してしまうと 当然正しくない

  • 正しい割合を得るためには: 重みの総和 = 1 に戻したい

    • よって 重みの総和で,それぞれの重みを割る(重みの正規化)↓

有理3次ベジェ曲線 (Rational Cubic Bezier Curve):

\begin{equation*} R(t) = \frac{ \sum_{i=0}^3 B_i^3 (t) w_i P_i } { \sum_{i=0}^3 B_i^3 (t) w_i } \end{equation*}
  • \(w_i\) は制御点にかかっている重み(weight)のこと.

  • \(\sum_{i=0}^3 B_i^3 (t) w_i\): 重みの総和.これで割ることで上で説明した重 みの正規化を行っている.

非一様有理B-Spline曲線 (Non-Uniform Rational B-Spline Curve)

上で行った3次ベジェ曲線の有理化と同じことを B-Spline曲線 で行う → NURBS曲線:

\begin{equation*} R(t) = \frac{ \sum_{i=1}^{k} N_i^n w_i P_i } { \sum_{i=1}^{k} N_i^n w_i } \end{equation*}
  • \(w_i\): 制御点にかかっている重み(weight)

  • \(\sum_{i=1}^k N_i^n w_i\) は重みの総和.正規化を行う.(正規化係数)

重みが全て1のときはただのB-Spline曲線になる.(重みつけされていないのと同じなので)

課題:

有理3次ベジェ曲線を書く.制御点に重みをつける.

bezier-weighted.pde

操作説明:

  • 左クリックで制御点追加,右クリックで削除

  • 制御点をドラッグ&ドロップで移動

  • 制御点をクリックして選択 → 'w'/'W' キーでウェイトを増減

第9回 パラメトリックな曲面

第4回, 第5回で既習したベジエ曲線を曲面に発展させます。

bezier-surface-bicubic.pde

操作説明: キーボードの'w'/'W'で,重みを増減

パラメトリックな曲線(2次元)を曲面(3次元)へ発展.

Bezier,B-Spline,NURBS曲線それぞれに曲面がある.今日の授業ではベジェ曲面をやる .

双3次ベジェ曲面

(Bicubic Bezier Surface)

式の形としてはベジェ曲線の式の基底関数とシグマ計算が二重になったもの:

\begin{equation*} S(u,v) = \sum_{i=0}^3 \sum_{j=0}^3 B_i^3 (u) B_j^3 (v) P_{ij} \end{equation*}

この曲面も前回と同じように有理化することができる:

双3次有利ベジェ曲面

(Bicubic Rational Bezier Surface)

\begin{equation*} S(u,v) = \frac{ \sum_{i=0}^3 \sum_{j=0}^3 B_i^3 (u) B_j^3 (v) P_ij } { \sum_{i=0}^3 \sum_{j=0}^3 B_i^3 (u) B_j^3 (v) } \end{equation*}

この式も混ぜ合わせ関数.直感的な説明としては,曲面上の各点は各制御点から,特定の 割合(全体が100%になる)で引っ張り合った均衡点になっている.

課題: 有利ベジェ曲面の改良

  • より高次(3x3, 3x4, 4x3, 4x4, ...)のベジェ曲面.

  • NxMのベジェ曲面や接続ができるとすごい.

第10回 行列

ベクトルをアフィン変換(移動、回転、拡大・縮小)する為に、CG/CADでは4x4の同次 座標表現が多用されます。第7回では、行列の未経験者を念頭にした解説も加えて、行 列によるベクトル操作を解説します。

matrix.pde

操作説明: クリックで開始/停止

アフィン変換 (Affine Transformation)

行列の掛け算によって,移動,回転,拡大・縮小すること

単位行列 (Unit Matrix)

対角成分は1,それ以外は全て0の正方行列.掛けても何も変換しない

\begin{equation*} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*}

移動 (Transformation)

\begin{align*} v_o[0] &= v_i[0] + x \\ v_o[1] &= v_i[1] + y \\ v_o[2] &= v_i[2] + z \\ v_o[3] &= 1 \end{align*}
\begin{equation*} (v_o[0], v_o[1], v_o[2], 1) = (v_i[0], v_i[1], v_i[2], 1) \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ x & y & z & 1 \end{bmatrix} \end{equation*}

ここで掛けている方を移動行列という. これをコードに落とすと:

void matrix_translate(float[][] M, float x, float y, float z) {
  init_matrix(M);
  M[3][0] = x;
  M[3][1] = y;
  M[3][2] = z;
}

回転 (Rotation)

  • Z軸周りに回転させる場合:

    \begin{align*} v_o[0] &= \cos\theta \ v_i[0] - \sin\theta \ v_i[1] \\ v_o[1] &= \sin\theta \ v_i[0] + \cos\theta \ v_i[1] \\ v_o[2] &= v_i[2] \\ v_o[3] &= 1 \end{align*}
    \begin{equation*} (v_o[0], v_o[1], v_o[2], 1) = (v_i[0], v_i[1], v_i[2], 1) \begin{bmatrix} \cos\theta & \sin\theta & 0 & 0 \\ -\sin\theta & \cos\theta & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*}
  • Y軸周り:

    \begin{align*} v_o[0] &= \sin\theta \ v_i[2] - \cos\theta \ v_i[0] \\ v_o[1] &= v_i[1] \\ v_o[2] &= \cos\theta \ v_i[2] - \sin\theta \ v_i[0] \\ v_o[3] &= 1 \end{align*}
    \begin{equation*} (v_o[0], v_o[1], v_o[2], 1) = (v_i[0], v_i[1], v_i[2], 1) \begin{bmatrix} \cos\theta & 0 & -\sin\theta & 0 \\ 0 & 1 & 0 & 0 \\ \sin\theta & 0 & \cos\theta & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*}
  • X軸周り:

    \begin{align*} v_o[0] &= v_i[0] \\ v_o[1] &= \cos\theta \ v_i[1] - \sin\theta \ v_i[2] \\ v_o[2] &= \sin\theta \ v_i[1] + \cos\theta \ v_i[2] \\ v_o[3] &= 1 \end{align*}
    \begin{equation*} (v_o[0], v_o[1], v_o[2], 1) = (v_i[0], v_i[1], v_i[2], 1) \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos\theta & \sin\theta & 0 \\ 0 & -\sin\theta & \cos\theta & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*}

これらの操作をコードにしている部分:

// rotation
void matrix_rotate(float[][] M, char axis, float theta) {
  init_matrix(M);
  if (axis == 'x' || axis == 'X') {
    M[1][1] =  cos(theta);
    M[1][2] =  sin(theta);
    M[2][1] = -sin(theta);
    M[2][2] =  cos(theta);
  } else if (axis == 'y' || axis == 'Y') {
    M[2][2] =  cos(theta);
    M[2][0] =  sin(theta);
    M[0][2] = -sin(theta);
    M[0][0] =  cos(theta);
  } else if (axis == 'z' || axis == 'Z') {
    M[0][0] =  cos(theta);
    M[0][1] =  sin(theta);
    M[1][0] = -sin(theta);
    M[1][1] =  cos(theta);
  }
}

拡大・縮小 (Scale)

\begin{align*} v_o[0] &= a \cdot v_i[0] \\ v_o[1] &= b \cdot v_i[1] \\ v_o[2] &= c \cdot v_i[2] \\ v_o[3] &= 1 \end{align*}
\begin{equation*} (v_o[0], v_o[1], v_o[2], 1) = (v_i[0], v_i[1], v_i[2], 1) \begin{bmatrix} a & 0 & 0 & 0 \\ 0 & b & 0 & 0 \\ 0 & 0 & c & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*}
// scale
void matrix_scale(float[][] M, float x, float y, float z) {
  init_matrix(M);
  M[0][0] = x;
  M[1][1] = y;
  M[2][2] = z;
}

行列の操作にすると何が良いのか:

例えばある図形に対して,拡大,回転,移動という操作を行いたい場合は,先にそれらの 拡大,回転,移動行列を掛け合わせた下の様な行列を求めておく.そうすれば座標に対し ては一回の行列掛け算で求めることができる.

\(\begin{bmatrix} 0 & -2 & 0 & 0 \\ 3 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 5 & 5 & 0 & 1 \end{bmatrix}\)

しかも行列計算はGPUなどが最適化しているため高速化ができる.

参考: 3x3行列のアフィン変換の図解

https://upload.wikimedia.org/wikipedia/commons/2/2c/2D_affine_transformation_matrix.svg

(By Cmglee, CC BY-SA 3.0)

最終課題の告知

得点:

  1. Processingの関数を用いずに独自実装するほど配点は高い

  2. 独創性の高いアルゴリズムや作品ほど配点は高い

  3. 提出者は最終回で発表する

第11回 3次元物体の表示

3次元CG/CADへ展開するために、3次元座標系(右手座標系)を導入し、ワールド座標系 、カメラ座標系、スクリーン座標系への変換、陰面処理を解説します。

(早慶戦の休講の影響で飛ばされた)

第12回 レイトレーシング

レイトレーシングいう、Ray(光線)をTracing(追跡する)事で3次元物体を描画する古 典的なレンダリング手法を解説します。

Ray Tracing

CGアニメーションの制作行程:

  • モデリング → テクスチャ → (リギング) → アニメーション → レンダリング → ポストプロダクション

この中でレイトレーシングは,レンダリングよりの部分で使われる.

  • スクリーン上のPixelから視点へ入射してくる光線(Ray)を逆方向に追跡(Trace)する手法

  • どうすれば set() 関数だけで(OpenGLなどを使わず)3Dの球体を描くことができるか?

今日授業でやるのは 平行投影 = Z軸と平行な光線を想定する. 視点からの光線にもできる.

衝突判定

  • \(P(Px, Py)\) を通り,Z軸に平行な直線と球との交点(基本2つある)を求め,手前 の交点の色が見える

  • 球は \(x^2+y^2+(z-z_o)^2 = r^2\,(z_o: \text{center of sphere})\) で表される.

  • \((P_x, P_y, P_z)\) を通り,単位ベクトル \(E(E_x, E_y, E_z)\) に平 行な直線上の点は:

    \begin{equation*} \left\{\begin{array}{} x = E_x t + P_x \\ y = E_y t + P_y \\ z = E_z t + P_z \end{array}\right. \end{equation*}

    と表される.(\(t\) はパラメータ)

  • 今回はZ軸と平行な直線なので \(E(E_x, E_y, E_z)\)\(E = (0, 0, 1)\)

    • なお \(E(E_x, E_y, E_z)\) のように方向を表す単位ベクトルを 方向余弦 ということがある.

  • したがって,

    \begin{equation*} \left\{\begin{array}{} x = 0\ t + P_x &= P_x \\ y = 0\ t + P_y &= P_y \\ z = 1\ t + 0 &= t \end{array}\right. \end{equation*}
    • つまり \(P_x, P_y\) によってZ軸に並行な直線が決まり,後はパラメータ \(t\) によって直線上の点が移動することになる

  • \(P_x^2 + P_y^2 + (t - Z_o)^2 = r^2 \leftrightarrow\\ t^2 - (2 Z_o)t + (P_x^2 + P_y^2 + Z_o^2 - r^2) = 0\)

  • ここから \(t\) の値は解の公式より,

    • \(t = Z_o \pm \sqrt{r^2 - P_x^2 - P_y^2}\)

    • \(r^2 - P_x^2 - P_y^2 < 0\) の時,交点は存在しない.

画面の全てのピクセルごとに RayTrace() が呼ばれて,色を返す.

raytracing-0.pde

color RayTrace(float px, float py) {
  float t = RayCast(px, py);
  return (t == -1) ? color(255*Bg[0], 255*Bg[1], 255*Bg[2])
                   : color(255*K[0], 255*K[1], 255*K[2]);
}

float RayCast(float px, float py) {
  float d = R*R - px*px - py*py;
  float t = (d > 0) ? (Zo - sqrt(d))
                    : -1;
  return t;
}
r/lec-cgcad/raytracing-0.png

この段階では,光が当たるか当たらないかの判断しかしていないので,のっぺりしていて 球体ではなく面にしか見えない.球体らしさは陰影によって生まれる.

→ Shadeをつけることをシェーディングという

  • shade: 物体の色が暗くなる陰

  • shadow: 物体が作る輪郭のある影

Lambert's cosine law

ランバートの余弦則
  • 拡散反射光の明るさは見る方向によって変わらない

  • 明るさは,面の法線(L)と入射光(L)のなす角のθの余弦に比例

\begin{equation*} L \cdot N = |L| |N| \cos\theta = L_x \cdot N_x + L_y \cdot N_y + L_z \cdot N_z \end{equation*}

内積が0以上の場合と0未満の場合に分かれる

  • 内積が0より小さい場合は,シェーディングの結果は計算しない

  • 0以上の場合は,光の拡散反射係数に内積の結果をかけて,物体の放射輝度をかける

    • 拡散反射係数: 要は物体が何色かということ. float K[] = {1.0,0.0,0.0}; は RGBのRを100%反射してGBは反射しない = color(255, 0, 0) という意味.

raytracing-1.pde

color Shading(float[] N) {
  float LN = innerProduct(L, N);
  return (LN < 0) ? color(0,0,0)
                  : color(255*K[0]*(dif[0]*LN),
                          255*K[1]*(dif[1]*LN),
                          255*K[2]*(dif[2]*LN));
}
r/lec-cgcad/raytracing-1.png

Ambient Light

上の方法だと,光が当たっていない部分が真っ黒になってしまうが,現実の世界ではそう はならない.(もっとも,宇宙空間では物体に当たらなかった光は反射してこないので, その場合は上の方が正しいが)

→ 光が直接は当たってないが見える部分を, 環境光 によって再現する.

raytracing-2.pde

color Shading(float[] N) {
  float LN = innerProduct(L, N);
  return (LN < 0) ? color(255 * K[0] * amb[0], // 光源が反対にある場合にも
                          255 * K[1] * amb[1], // 環境光によって色をつける
                          255 * K[2] * amb[2])
                  : color(255 * K[0] * (dif[0]*LN + amb[0]),
                          255 * K[1] * (dif[1]*LN + amb[1]),
                          255 * K[2] * (dif[2]*LN + amb[2]));
}
r/lec-cgcad/raytracing-2.png

参考になりそう:

課題: シェーディングの改良

光源や球の増加,三角形のシェーディングなど.

  • ランバートモデル

  • フォンモデル (プラスチックっぽい)

  • ブリンモデル (金属っぽい)


Phongの反射モデル: raytracing-phong.pde

r/lec-cgcad/raytracing-phong.png

点pを視点ベクトルVから見たときの反射光は以下の式で求められる:

\begin{equation*} R_{specular}(p, V) = k_{specular}(p)\ (R \cdot V)^S\ L(p) \end{equation*}

(これは物理的には間違っているらしい)

  • \(k_{specular}\): 鏡面反射係数

  • \(R\): 反射ベクトル (\(L\) が表面で完全反射されたもの)

    \begin{align*} R &= \text{reflect}(L, N) \\ &= 2(L \cdot N) N - L \end{align*}
    • \(N\): 表面上の点における法線ベクトル

  • \(S\): 鏡面反射指数

  • \(L\): 光源ベクトル

直感的には,完全に反射された \(R\)\(L\) に反射係数の \(k_{specular}\) をかけたものになり, \(R\) との差が大きくなるほど反射光 は弱くなっていく(内積 \(R \cdot V\) が小さくなっていく)と理解できる.

参考:


最終課題の採点基準:

  • 技術的/アルゴリズム的な難しさ

  • 審美的

  • プレゼンテーション

  • 努力

  • 独創性

第13回 シェーディング

CG/CADでは、影と陰は区別されます。第9回では後者の、光源が物体に落とす陰影モデ ルについて、具体的にランバートモデル、グーローモデル、フォンモデルの数理を解説 します。

r/lec-cgcad/shadowing.png

Shadowing (影付け)

前回で,物体の陰(shade)は表現できたが,物体が作る影(shadow)も実現したい.

Shadowingの原理は非常にシンプル: ある点から光源が...

  • 見える → そこは影ではない

  • 見えない → そこは影

この判定を実現するには,以下2つの交点を求めればいい.

  • 交点 \(P (P_x, P_y, P_z)\) を通り,方向余弦 \(L (L_x, L_y, L_z)\) を持 つ直線.(つまり点から光源へ向かう直線)

    \begin{equation*} \left\{\begin{array}{} x &= L_x t + P_x \\ y &= L_y t + P_y \\ z &= L_z t + P_z \end{array}\right. \end{equation*}
  • 半径 \(R\) で中心が \((X_o, Y_o, Z_o)\) の球

    \begin{equation*} (x-X_o)^2 + (y-Y_o)^2 + (z-Z_o)^2 = r^2 \end{equation*}

これを \(t\) について整理して \(at^2 + bt + c = 0\) の形にすると:

\begin{align*} &(L_x + P_x + X_o)^2 + (L_y + P_y + Y_o)^2 + (L_z + P_z + Z_o)^2 = r^2 \\ \to\ &L_xt^2 + 2tL_x(P_x - X_o) + (P_x - X_o)^2 + \\ &L_yt^2 + 2tL_y(P_y - Y_o) + (P_y - Y_o)^2 + \\ &L_zt^2 + 2tL_z(P_z - Z_o) + (P_z - Z_o)^2 = r^2 \\ \to\ &(L_x^2 + L_y^2 + L_z^2) t^2 + \\ &2 \left\{ L_x(P_x - X_o) + L_y(P_y - Y_o) + L_z(P_z - Z_o) \right\} t + \\ &(P_x - X_o)^2 + (P_y - Y_o)^2 + (P_z - Z_o)^2 - r^2 = 0 \end{align*}

ここから,

\begin{equation*} \left\{\begin{array}{} a &= (L_x^2 + L_y^2 + L_z^2) \\ 2b &= L_x(P_x - X_o) + L_y(P_y - Y_o) + L_z(P_z - Z_o) \\ c &= (P_x - X_o)^2 + (P_y - Y_o)^2 + (P_z - Z_o)^2 - r^2 \end{array}\right. \end{equation*}

後は前回と同じように解の公式を使って,交点の存在を確かめることができる. 解が存在するということは,その点は影の中にあるため block = true になっている.

\begin{equation*} t = \frac{-b \pm \sqrt{b^2 - ac}}{a} \end{equation*}

shadowing.pde

float Shadowing(float[] P) {
  boolean block = false;
  float t0;

  for (int i = 0; i < OBJ_NUM; i++) {
    float a, b, c, D;
    if (primitives[i] instanceof Sphere) {
      Sphere s = new Sphere();
      s = (Sphere)primitives[i];
      a = L[0]*L[0] + L[1]*L[1] + L[2]*L[2];
      b = (L[0] * (P[0] - s.Xo) +
           L[1] * (P[1] - s.Yo) +
           L[2] * (P[2] - s.Zo));
      c = (P[0] - s.Xo) * (P[0] - s.Xo) +
          (P[1] - s.Yo) * (P[1] - s.Yo) +
          (P[2] - s.Zo) * (P[2] - s.Zo) - s.R*s.R;
      D = b*b - a*c;
      if (D >= 0) {
        t0 = -b - sqrt(D);
        if (t0 > 0) {
          block = true; break;
        }
      }
    }
  }
  return (block) ? 0 : 1;
}

ただこの実装では,影かどうかの2値しか判断していないのでエッジがギザギザになって いる.そこで,もっとリアリスティックにするために Soft shadow などの技術が開発さ れている.

第14回 最終発表会

講義内容を振り返るともに、今後のCG/CADの動向について触れます。

最終課題として,NxN次のベジェ曲面を描くことができ,各制御点の位置と重みをマウス とキーボードで自由に操作ができるというプログラムを書いた. が,発表会の場には間に合わなかったのでデバッグ中のものを見せるという恥ずかしいこ とになってしまった (課題の提出期限には間に合った)

マウスを3D空間上の平面へマップする関数 unProjectPoint は <http://d.hatena.ne.jp/kougaku-navi/20160102/p1> を参考にした. が下では,Processing.jsで動かすためにマウス関連の処理を消している.

操作方法:

bezier-surface.pde