読者です 読者をやめる 読者になる 読者になる

TadaoYamaokaの日記

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

SSDによる物体検出を試してみた

SSD TensorFlow Windows

先日の日記でYOLOv2による物体検出を試してみたが、YOLOと同じくディープラーニングで物体の領域検出を行うアルゴリズムとしてSSD(Single Shot MultiBox Detector)がある。

YOLOv2の方が精度が高いとYOLOv2の論文に書かれているが、SSDの精度も高いようなので試してみた。

オリジナルのSSDの実装は、Caffeが使用されているが、WindowsでビルドできるCaffeとバージョンが異なるものが使用されており、トライしてみたがビルドがうまくいかなかった。

Python系のフレームワークの実装がないか調べたところ、TensorFlowの実装が2つとKerasの実装が1つとChainerの実装が1つ見つかった。

TensorFlowの実装
ssd_tensorflow
SSD-Tensorflow
Kerasの実装
ssd_keras
Chainerの実装
chainer-SSD

個人的にはChainerを使い慣れているが、他に比べて更新が古くてメンテされていないようだったので、熱心に更新されていたTensorFlowの実装を試すことにした。

github.com

生のTensorFlowは記述が冗長なので、あまりさわりたくないが、この実装ではTF-Slimが使用されており少しは読みやすくなっている。


学習済みモデルとデモも配布されていたので、README.mdの「SSD minimal example」の通り実行すると、問題なく動作した。

f:id:TadaoYamaoka:20170317233452p:plain

実行環境には、WindowsにインストールしたTensorFlowを使用した。


これを使って、先日のYOLOv2で試した動画で、物体検出を試してみた。

www.nicovideo.jp

そもそもアニメが学習データに入っていないので比較することには意味がないかもしれないが、YOLOv2と比べて検知漏れや物体でないところでの誤検知、キャラクターがaeroplaneになることが多い。

実行速度は、FPS:45.8であった。YOLOv2よりは遅いがリアルタイムに近い速度が出ている。


なお、TensorFlowの実装には、動画に対して検出するコードは含まれていなかったので、自分で作成した。
以下に作成した動画対応のコードを示しておく。

SSD minimal example」の最後の「In [11]」の部分を以下のコードに置き換える。

VOC_LABELS = {
    0: 'none',
    1: 'aeroplane',
    2: 'bicycle',
    3: 'bird',
    4: 'boat',
    5: 'bottle',
    6: 'bus',
    7: 'car',
    8: 'cat',
    9: 'chair',
    10: 'cow',
    11: 'diningtable',
    12: 'dog',
    13: 'horse',
    14: 'motorbike',
    15: 'person',
    16: 'pottedplant',
    17: 'sheep',
    18: 'sofa',
    19: 'train',
    20: 'tvmonitor',
}

colors = [(random.randint(0,255), random.randint(0,255), random.randint(0,255)) for i in range(len(VOC_LABELS))]

def write_bboxes(img, classes, scores, bboxes):
    """Visualize bounding boxes. Largely inspired by SSD-MXNET!
    """
    height = img.shape[0]
    width = img.shape[1]
    for i in range(classes.shape[0]):
        cls_id = int(classes[i])
        if cls_id >= 0:
            score = scores[i]
            ymin = int(bboxes[i, 0] * height)
            xmin = int(bboxes[i, 1] * width)
            ymax = int(bboxes[i, 2] * height)
            xmax = int(bboxes[i, 3] * width)
            cv2.rectangle(img, (xmin, ymin), (xmax, ymax),
                                 colors[cls_id],
                                 2)
            class_name = VOC_LABELS[cls_id]
            cv2.rectangle(img, (xmin, ymin-6), (xmin+180, ymin+6),
                                 colors[cls_id],
                                 -1)
            cv2.putText(img, '{:s} | {:.3f}'.format(class_name, score),
                           (xmin, ymin + 6),
                           cv2.FONT_HERSHEY_PLAIN, 1,
                           (255, 255, 255))


import time

vid = cv2.VideoCapture('path/to/movie')
if not vid.isOpened():
    raise IOError(("Couldn't open video file or webcam. If you're "
    "trying to open a webcam, make sure you video_path is an integer!"))

vidw = vid.get(cv2.CAP_PROP_FRAME_WIDTH)
vidh = vid.get(cv2.CAP_PROP_FRAME_HEIGHT)
fps = vid.get(cv2.CAP_PROP_FPS)

fourcc = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
out = cv2.VideoWriter('path/to/output.avi', int(fourcc), fps, (int(vidw), int(vidh)))

prev_time = time.time()
frame_cnt = 0
while True:
    retval, img = vid.read()
    if not retval:
        print("Done!")
        break
    
    rclasses, rscores, rbboxes =  process_image(img)
    write_bboxes(img, rclasses, rscores, rbboxes)
    out.write(img)
    frame_cnt += 1

curr_time = time.time()
exec_time = curr_time - prev_time
print('FPS:{0}'.format(frame_cnt/exec_time))

vid.release()
out.release()

WindowsにCaffeをインストールする

Caffe Windows

以前にWindowsでCaffeをビルドしてGPUで実行する方法について書きましたが、内容が古くなったので書き直します。

以前に日記を書いた時点では、ビルド済みバイナリが配布されていなかったため、自分でビルドを行いましたが、ビルド済みバイナリが配布されるようになりました。

公式からリンクされているGitHubWindowsブランチからダウンロードできます。

github.com

Visual Studio 2013用と2015用、GPU対応、Pythonのバージョンの組み合わせで、5種類のビルド済みバイナリが用意されています。

ここでは、
Visual Studio 2015, CUDA 8.0, Python 3.5: Caffe Release
をダウンロードした前提で記述します。
他のバージョンでもおそらく同様の手順で大丈夫です。

環境

インストール方法

ビルド済みバイナリのダウンロード

WindowsブランチのREADME.mdのPrebuilt binariesに記載されているリンクから環境に合わせたバイナリをダウンロードします。

アーカイブ展開

ダウンロードしたcaffe.zipを任意のディレクトリに展開します。

インストールは以上です。

MNISTサンプルの実行

公式の手順に従って、MNISTサンプルを実行します。

リポジトリ取得

caffe.zipにはサンプルが含まれていないので、Githubリポジトリを取得します。

任意のディレクトリで、

git clone https://github.com/BVLC/caffe.git

を実行して、リポジトリを取得します。

ディレクトリのコピー

git cloneしたローカルリポジトリから、以下のディレクトリをビルド済みのcaffeを展開したディレクトリに、コピーします。

  • data
  • examples
MNISTデータダウンロード

以下の手順で実行するシェルは、bashのシェルになっているので、Windowsで実行するにはMinGWなどでbashをインストールします。

bashを起動し、Caffeのディレクトリに移動し、

./data/mnist/get_mnist.sh

を実行し、MNISTのデータをダウンロードします。
wgetとgunzipが使われています。不足しているコマンドがあれば、MinGWのパッケージマネージャでインストールしてください。

データ準備

examples\mnist\create_mnist.shを編集します。

「BUILD=build/examples/mnist」→「BUILD=bin」に修正
「convert_mnist_data.bin」→「convert_mnist_data.exe」に修正(2か所)
編集後、

./examples/mnist/create_mnist.sh

を実行します。

訓練実行

examples\mnist\train_lenet.shを編集します。

「./build/tools/caffe」→「./bin/caffe」に修正
編集後、

./examples/mnist/train_lenet.sh

を実行します。

成功すれば以下のように出力されます。

$ ./examples/mnist/train_lenet.sh
I0315 23:17:26.000556 14060 caffe.cpp:219] Using GPUs 0
I0315 23:17:26.162506 14060 caffe.cpp:224] GPU 0: GeForce GTX 1080
I0315 23:17:26.425354 14060 common.cpp:36] System entropy source not available, using fallback algorithm to generate seed instead.
I0315 23:17:28.068279 14060 solver.cpp:44] Initializing solver from parameters:
test_iter: 100
test_interval: 500
base_lr: 0.01
display: 100
max_iter: 10000
lr_policy: "inv"
gamma: 0.0001
power: 0.75
momentum: 0.9
weight_decay: 0.0005
snapshot: 5000
snapshot_prefix: "examples/mnist/lenet"
solver_mode: GPU
device_id: 0
net: "examples/mnist/lenet_train_test.prototxt"
train_state {
  level: 0
  stage: ""
}
I0315 23:17:28.069281 14060 solver.cpp:87] Creating training net from net file: examples/mnist/lenet_train_test.prototxt
I0315 23:17:28.071815 14060 net.cpp:296] The NetState phase (0) differed from the phase (1) specified by a rule in layer mnist
I0315 23:17:28.071815 14060 net.cpp:296] The NetState phase (0) differed from the phase (1) specified by a rule in layer accuracy
I0315 23:17:28.072289 14060 net.cpp:53] Initializing net from parameters:
(省略)
I0315 23:19:28.510016 14060 solver.cpp:219] Iteration 9700 (90.204 iter/s, 1.1086s/100 iters), loss = 0.0021487
I0315 23:19:28.510545 14060 solver.cpp:238]     Train net output #0: loss = 0.00214882 (* 1 = 0.00214882 loss)
I0315 23:19:28.511018 14060 sgd_solver.cpp:105] Iteration 9700, lr = 0.00601382
I0315 23:19:29.621973 14060 solver.cpp:219] Iteration 9800 (90.028 iter/s, 1.11077s/100 iters), loss = 0.00816679
I0315 23:19:29.622474 14060 solver.cpp:238]     Train net output #0: loss = 0.0081669 (* 1 = 0.0081669 loss)
I0315 23:19:29.622474 14060 sgd_solver.cpp:105] Iteration 9800, lr = 0.00599102
I0315 23:19:30.711405 14060 solver.cpp:219] Iteration 9900 (91.8658 iter/s, 1.08854s/100 iters), loss = 0.00674243
I0315 23:19:30.711905 14060 solver.cpp:238]     Train net output #0: loss = 0.00674254 (* 1 = 0.00674254 loss)
I0315 23:19:30.712406 14060 sgd_solver.cpp:105] Iteration 9900, lr = 0.00596843
I0315 23:19:31.783195 14060 solver.cpp:448] Snapshotting to binary proto file examples/mnist/lenet_iter_10000.caffemodel
I0315 23:19:31.795729 14060 sgd_solver.cpp:273] Snapshotting solver state to binary proto file examples/mnist/lenet_iter_10000.solverstate
I0315 23:19:31.805755 14060 solver.cpp:311] Iteration 10000, loss = 0.00296754
I0315 23:19:31.805755 14060 solver.cpp:331] Iteration 10000, Testing net (#0)
I0315 23:19:32.516114  7412 data_layer.cpp:73] Restarting data prefetching from start.
I0315 23:19:32.544718 14060 solver.cpp:398]     Test net output #0: accuracy = 0.9909
I0315 23:19:32.544718 14060 solver.cpp:398]     Test net output #1: loss = 0.0277917 (* 1 = 0.0277917 loss)
I0315 23:19:32.545192 14060 solver.cpp:316] Optimization Done.
I0315 23:19:32.545692 14060 caffe.cpp:260] Optimization Done.

YOLOv2による物体検出を試してみた

Darknet Windows YOLOv2

前回の日記WindowsにインストールしたDarknetを使ってYOLOv2による物体検出を試してみました。

Darknetの学習済みモデルを使用して、ニコニコ動画の上位にあった動画に対して行ってみました。

こちらの動画です。

www.nicovideo.jp


キャラクターがだいたいpersonと表示されています。
たまに、dogになったり、birdになったりcowになったりするのはある意味正しい?

検出速度はFPS=55.3でしたので、かなり高速です。

今回使用した学習済みモデルは、VOCという写真のデータセットで学習しているので、アニメ向け特化したモデルを学習すればさらに精度が上がると思われます。

DarknetをWindowsにインストールする

Darknet Windows YOLOv2

以前に行ったマンガの顔パーツ検出では、顔の領域の検出にdlibのHOG+SVMを使った方法で矩形領域を抽出してから、DCNNで顔パーツ検出を行っていた。
HOG+SVMである程度マンガから顔の領域を抽出できるが、少なからず検出漏れがあった。
また、スライディングウィンドウを使用しているため、処理時間がかかる問題もある。

領域の検出自体もディープラーニングで行うことができれば精度が上がるかもしれない。
ディープラーニングを使用した領域の検出にどのような方法があるか調べたら、YOLOv2というアルゴリズムが現時点で最も精度が高く高速で動作するようだ。

YOLOv2の実装は、YOLOv2の作者自身によるC言語で記述されたDarknetというフレームワークが使われている。

pjreddie.com

ビルド環境はLinux向けになっており、Windowsで試すにはプロジェクトの修正が必要になる。

すでにWindows向けにポーティングされていないか調べたら、フォークされたリポジトリがあった。

github.com

これをありがたく使用させていただき、Windowsでビルドを行った。

ビルドは、README.mdの通りに行えば問題なく成功した。

以下に実施した手順を示す。

環境

ビルド手順

Gitリポジトリのclone
git clone https://github.com/AlexeyAB/darknet.git
Visual Studioでプロジェクトを開く

git cloneしたディレクトリ配下の
darknet\build\darknet\darknet.sln
Visual Studioで開く。

ビルド構成設定

ビルド構成を「Release/x64」に変更する。

OpenCVのパス設定

README.mdに書かれているとおり、OpenCV 2.4.9をダウンロードして任意のディレクトリに展開する。
https://sourceforge.net/projects/opencvlibrary/files/opencv-win/2.4.9/opencv-2.4.9.exe/download

プロジェクトのプロパティ->VC++ディレクトリから、OpenCVを展開したディレクトリ配下のopencv\build\includeをインクルードディレクトリに追加する。

ライブラリディレクトリに、OpenCVを展開したディレクトリ配下のopencv\build\x64\vc12\libを追加する。

なお、ディレクトリ名がvc12になっているが、Visual Studio 2015でも問題なく使用できた。

ビルド

ソリューションのビルドを実行する。

成功すると、
darknet\build\darknet\x64
にdarknet.exeができる。

学習済みモデルダウンロード

README.mdに書かれているURLから学習済みモデルをダウンロードして、darknet.exeと同じディレクトリに格納する。
モデルファイルは、4数種類ある。
モデルごとに使用メモリサイズ、精度、速度に違いがある。

とりあえずサンプル実行に必要な
yolo-voc.weights
をダウンロードした。
http://pjreddie.com/media/files/yolo-voc.weights

サンプル実行

画像ファイルからオブジェクトを検出するサンプルを試す。

コマンドプロンプトを起動する。

環境変数PATHにOpenCVのパスを設定する。

set path=%path%;H:\src\opencv\build\x64\vc12\bin

OpenCVをH:\srcに展開した場合の例
※システムのプロパティで環境変数PATHに追加してもよい。その場合、設定後はコマンドプロンプトの再起動が必要。

darknet.exeのあるディレクトリに移動する。

cd darknet\build\darknet\x64

サンプルを実行する。

darknet_voc.cmd

以下のような出力がされた後、画像ファイルのパスの入力プロンプトが表示される。

layer     filters    size              input                output
    0 conv     32  3 x 3 / 1   416 x 416 x   3   ->   416 x 416 x  32
    1 max          2 x 2 / 2   416 x 416 x  32   ->   208 x 208 x  32
    2 conv     64  3 x 3 / 1   208 x 208 x  32   ->   208 x 208 x  64
    3 max          2 x 2 / 2   208 x 208 x  64   ->   104 x 104 x  64
    4 conv    128  3 x 3 / 1   104 x 104 x  64   ->   104 x 104 x 128
    5 conv     64  1 x 1 / 1   104 x 104 x 128   ->   104 x 104 x  64
    6 conv    128  3 x 3 / 1   104 x 104 x  64   ->   104 x 104 x 128
    7 max          2 x 2 / 2   104 x 104 x 128   ->    52 x  52 x 128
    8 conv    256  3 x 3 / 1    52 x  52 x 128   ->    52 x  52 x 256
    9 conv    128  1 x 1 / 1    52 x  52 x 256   ->    52 x  52 x 128
   10 conv    256  3 x 3 / 1    52 x  52 x 128   ->    52 x  52 x 256
   11 max          2 x 2 / 2    52 x  52 x 256   ->    26 x  26 x 256
   12 conv    512  3 x 3 / 1    26 x  26 x 256   ->    26 x  26 x 512
   13 conv    256  1 x 1 / 1    26 x  26 x 512   ->    26 x  26 x 256
   14 conv    512  3 x 3 / 1    26 x  26 x 256   ->    26 x  26 x 512
   15 conv    256  1 x 1 / 1    26 x  26 x 512   ->    26 x  26 x 256
   16 conv    512  3 x 3 / 1    26 x  26 x 256   ->    26 x  26 x 512
   17 max          2 x 2 / 2    26 x  26 x 512   ->    13 x  13 x 512
   18 conv   1024  3 x 3 / 1    13 x  13 x 512   ->    13 x  13 x1024
   19 conv    512  1 x 1 / 1    13 x  13 x1024   ->    13 x  13 x 512
   20 conv   1024  3 x 3 / 1    13 x  13 x 512   ->    13 x  13 x1024
   21 conv    512  1 x 1 / 1    13 x  13 x1024   ->    13 x  13 x 512
   22 conv   1024  3 x 3 / 1    13 x  13 x 512   ->    13 x  13 x1024
   23 conv   1024  3 x 3 / 1    13 x  13 x1024   ->    13 x  13 x1024
   24 conv   1024  3 x 3 / 1    13 x  13 x1024   ->    13 x  13 x1024
   25 route  16
   26 reorg              / 2    26 x  26 x 512   ->    13 x  13 x2048
   27 route  26 24
   28 conv   1024  3 x 3 / 1    13 x  13 x3072   ->    13 x  13 x1024
   29 conv    125  1 x 1 / 1    13 x  13 x1024   ->    13 x  13 x 125
   30 detection
Loading weights from yolo-voc.weights...Done!
Enter Image Path: 

画像ファイルのパスを入力してEnterを押すと、コンソールに以下のような結果が出力され、検出された領域がウィンドウに表示される。

****.jpg: Predicted in 0.187000 seconds.
boat: 11%
person: 12%
person: 11%
person: 69%
person: 28%

 SRC output_video = 0000000000000000

 cvCreateVideoWriter, DST output_video = 000001D281CC12F0

 cvWriteFrame

また、結果が画像として、darknet.exeのあるディレクトリにpredictions.jpgというファイル名で保存される。

ウィンドウを閉じると次の画像のパスを入力するプロンプトが表示される。
Ctrl+Cで終了する。

Free-Form Deformationでテクスチャを変形する

Python PyOpenGL ベジエ曲面 FFD

前回までの日記で、OpenGLでテクスチャ(2D画像)を平面に描画する方法と、ベジエ曲面を描画する方法について記述した。
この2つを組み合わせて、テクスチャの変形を行う。

テクスチャの変形には、自由形状変形(FFD: Free-Form Deformation)という手法を使う。
FFDは、ベジエ曲面を用いて座標系を変形する。

FFDによる座標系(u, v)の変形は、制御点を\mathbf{P_{ij}}とすると、以下の式で表される。
\displaystyle
\begin{eqnarray*}
\mathbf{X}(u, v) &=& \sum_{i=0}^{m} \sum_{j=0}^{n} B_i^m(u) B_j^n(v) \mathbf{P_{ij}},\;\;\; 0 \leq u \leq 1, 0 \leq v \leq 1 \\
B_i^n(t) &=& {}_nC_i t^i(1 - t)^{n-i}
\end{eqnarray*} \tag{1}
これは、前回の日記で説明したベジエ曲面の式と同じである。

FFDを使用してテクスチャを変形するには、テクスチャを格子状の頂点で構成される面に貼り付けた後、頂点をFFDで変形する。
なお、テクスチャを貼り付ける面は格子状でなくてもよい。この記事では簡単なため格子状とした。

Live2Dの曲面デフォーマーは、このFFDを使用していると思われる。
これによりLive2Dと同じような変形が一部実現できる。

Pythonで記述すると以下のようになる。

import sys
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
from PIL import Image
import numpy as np
import scipy.misc as scm

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

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))

def load_texture():
    img = Image.open("sample2.png")
    w, h = img.size
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, w, h, 0, GL_RGBA, GL_UNSIGNED_BYTE, img.tobytes())

def display():
    #P = np.array([
    #    [[ 0. , 0. ], [ 0. , 1/3], [ 0. , 2/3], [ 0. , 1. ]],
    #    [[ 1/3, 0. ], [ 1/3, 1/3], [ 1/3, 2/3], [ 1/3, 1. ]],
    #    [[ 2/3, 0. ], [ 2/3, 1/3], [ 2/3, 2/3], [ 2/3, 1. ]],
    #    [[ 1. , 0. ], [ 1. , 1/3], [ 1. , 2/3], [ 1. , 1. ]],
    #    ])
    P = np.array([
        [[ 0. , 0. ], [ 0. , 0.3], [ 0. , 0.6], [ 0. , 1. ]],
        [[ 0.3, 0.1], [ 0.3, 0.4], [ 0.3, 0.8], [ 0.3, 1.1]],
        [[ 0.6, 0.1], [ 0.7, 0.4], [ 0.7, 0.7], [ 0.6, 1. ]],
        [[ 1. , 0. ], [ 1.1, 0.3], [ 1.1, 0.6], [ 1. , 1. ]],
        ])

    X = np.zeros((11 * 11, 2))
    T = np.zeros((11 * 11, 2)) # texture
    for i in range(11):
        for j in range(11):
            u = i / 10.0
            v = j / 10.0
            X[i * 11 + j] = bezier_patches(3, 3, u, v, P) # FFD

            T[i * 11 + j] = [i / 10.0, 1.0 - j / 10.0]

    V = np.zeros((11 * 11, 3))
    for i in range(11):
        for j in range(11):
            V[i * 11 + j] = [X[i * 11 + j, 0] - 0.5, X[i * 11 + j, 1] - 0.5, 0]

    IDX = np.zeros((10, 11 * 2))
    for i in range(10):
        for j in range(11):
            IDX[i, j * 2] = i * 11 + j
            IDX[i, j * 2 + 1] = (i + 1) * 11 + j

    glClear(GL_COLOR_BUFFER_BIT)

    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
    glEnable(GL_BLEND)

    glEnableClientState(GL_VERTEX_ARRAY)
    glEnableClientState(GL_TEXTURE_COORD_ARRAY)

    glEnable(GL_TEXTURE_2D)

    glVertexPointerf(V)
    glTexCoordPointerf(T)

    for idx in IDX:
        glDrawElementsui(GL_TRIANGLE_STRIP, idx)

    glDisable(GL_TEXTURE_2D)

    glDisable(GL_BLEND)

    glColor3f(1.0, 1.0, 1.0)
    glPointSize(2.0)
    glBegin(GL_POINTS)
    for p in P.reshape((-1, 2)):
        glVertex2f(p[0] - 0.5, p[1] - 0.5)
    glEnd()

    glFlush()

def main():
    glutInit(sys.argv)
    glutInitDisplayMode(GLUT_RGBA)
    glutInitWindowSize(300, 300)
    glutCreateWindow(b"FFDSample1")
    glutDisplayFunc(display)
    glClearColor(0.0, 0.0, 1.0, 1.0)

    load_texture()
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)

    glutMainLoop()

if __name__ == '__main__':
    main()

1.0×1.0の大きさの面を10分割してテクスチャを貼り付けている。
4×4の制御点を等間隔の格子から少しずらしている。
その制御点を使用して、FFDで、テクスチャを貼り付けた頂点を変形している。
テクスチャを貼り付けた面の中心が(0, 0)になるように(-0.5, -0.5)だけ移動して表示している。

実行結果
  • 変形前

f:id:TadaoYamaoka:20170309225604p:plain

  • 変形後

f:id:TadaoYamaoka:20170309225659p:plain

白い点は制御点を示している。

参考文献

Pythonでベジエ曲面を描く

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

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

Pythonでベジエ曲線を描く

Python ベジエ曲線

Free-Form Deformation(FFD)を使用して2D画像の変形を行おうと試みているが、FFDはベジエ(Bezier)曲面で表される座標系を使用して変形する。
そこで、ベジエ曲線を復習を兼ねてmatplotlibで描画してみる。

ベジエ曲線を2次元に拡張するとベジエ曲面となる。
まずは、単純なベジエ曲線を描画する。

ベジエ曲線の定義

ベジエ曲線は次式で表される曲線である。
\displaystyle
\begin{eqnarray*}
P(t) &=& \sum_{i=0}^{n} B_i^n(t)\mathbf{q_i},\;\;\; 0 \leq t \leq 1 \\
B_i^n(t) &=& {}_nC_i t^i(1 - t)^{n-i}
\end{eqnarray*} \tag{1}

{}_nC_iは二項係数で、\binom{n}{i}と表記する場合もあり、以下の式で展開される。
\displaystyle
{}_nC_i = \binom{n}{i} = \frac{n!}{i!(n-i)!} \tag{2}

B_i^n(t)は、バーンスタイン(Bernstein)多項式と呼ばれる。

\mathbf{q_i}は制御点で、nはベジエ曲線の次数である。
ベジエ曲線は、始点\mathbf{q_0}と終点\mathbf{q_n}を通る。

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(n, t, q):
    p = np.zeros(2)
    for i in range(n + 1):
        p += bernstein(n, i, t) * q[i]
    return p

q = np.array([[0, 0], [1, 1], [4, -1], [5, 0]], dtype=np.float)

list = []
for t in np.linspace(0, 1, 100):
    list.append(bezier(3, t, q))
P = np.array(list)

plt.plot(P.T[0], P.T[1])
plt.plot(q.T[0], q.T[1], '--o')
plt.show()

二項係数は、scipy.misc.combを使用して計算している。

上記のコード例では、制御点の次数は3で、制御点(0, 0)、(1, 1)、(4, -1)、(5, 0)によって表されるベジエ曲線を描画している。

なお、list.appendの部分をリスト内包を使ってもう少しスマートに書くと、

P = np.array([bezier(3, t, q) for t in np.linspace(0, 1, 100)])

のように記述できる。

bezier関数でΣの計算にforを使用しているが、以下のようにリスト内包と行列の積を使用しても記述できる。

np.dot([bernstein(n, i, t) for i in range(n + 1)], q)

numpyの行列は、横ベクトルなので、行列の積を計算する際は間違わないように注意が必要である。

実行結果

f:id:TadaoYamaoka:20170307220018p:plain