TadaoYamaokaの日記

山岡忠夫 Home で公開しているプログラムの開発ネタを中心に書いていきます。

Pythonでベジエ曲面を描く

前回の日記Pythonベジエ曲線を描いたが、今回はベジエ曲面(Bezier patches)を描く。

ベジエ曲面は、以下の式で表される。
\displaystyle
P(u, v) = \sum_{i=0}^{m} \sum_{j=0}^{n} B_i^m(u) B_j^n(v) \mathbf{q_{ij}},\;\;\; 0 \leq u \leq 1, 0 \leq v \leq 1
\tag{1}

B_i^n(t)は、前回の日記で説明したバーンスタイン(Bernstein)多項式である。

\mathbf{q_{ij}}は、平面上に格子状に配置した制御点で、m, nは、それぞれu軸、v軸方向の制御点の次数である。
曲面の四隅は、制御点の四隅を通る。

Pythonのコード例

ベジエ曲面をPythonで記述すると以下のようになる。

import scipy.misc as scm
import numpy as np
import matplotlib.pyplot as plt

def bernstein(n, i, t):
    return scm.comb(n, i) * t**i * (1 - t)**(n-i)

def bezier_patches(m, n, u, v, q):
    p = np.zeros(2)
    for i in range(m + 1):
        for j in range(n + 1):
            p += bernstein(m, i, u) * bernstein(n, j, v) * q[i, j]
    return p

q = np.array([
    [[0, 0], [0, 1], [0, 3], [0, 5]],
    [[1, 1], [1, 2], [1, 4], [1, 6]],
    [[2, 1], [2, 3], [2, 5], [2, 6]],
    [[4, 1], [4, 3], [4, 4], [4, 5]],
    [[5, 0], [6, 1], [6, 3], [5, 5]],
    ])

list1 = []
for u in np.linspace(0, 1, 10):
    list2 = []
    for v in np.linspace(0, 1, 10):
        list2.append(bezier_patches(4, 3, u, v, q))
    list1.append(list2)
P = np.array(list1)

for i in range(10):
    plt.plot(P[i].T[0], P[i].T[1], color='g')
    P_T = np.array([P[j, i] for j in range(10)])
    plt.plot(P_T.T[0], P_T.T[1], color='g')

for i in range(5):
    plt.plot(q[i].T[0], q[i].T[1], '--o', color='b')
for j in range(4):
    q_T = np.array([q[i, j] for i in range(5)])
    plt.plot(q_T.T[0], q_T.T[1], '--o', color='b')

plt.show()

次数m=4, n=3の制御点を配置して、ベジエ曲面を横軸、縦軸方向にそれぞれ10本の線にして描画している。

関数bezier_patchesのfor文は、リスト内包と行列の積を使って記述すると以下のようになる。

def bezier_patches(m, n, u, v, q):
    return np.dot([bernstein(n, j, v) for j in range(n + 1)], np.tensordot([bernstein(m, i, u) for i in range(m + 1)], q, axes=1))
実行結果

f:id:TadaoYamaoka:20170308214747p:plain

緑の線がベジエ曲面、青の点線が制御点の格子を表している。