Python
画像処理
ComputerVision
Jupyter

「実践コンピュータビジョン」の全演習問題をやってみた。詳細編 第4章 カメラモデルと拡張現実感

まえおき

先日投稿した「実践コンピュータビジョン」の全演習問題をやってみたの詳細です。なお、あくまで「やってみた」であって、模範解答ではありません。「やってみた」だけなので、うまくいかなかった、という結果もあります。

この記事は「実践コンピュータビジョン」(オライリージャパン、Jan Erik Solem著、相川愛三訳)の演習問題に関する記事です。この記事を書いた人間と書籍との間には何の関係もありません。著者のJan氏と訳者の相川氏にはコンピュータビジョンについて学ぶ機会を与えて頂いたことに感謝します。

実践コンピュータビジョンの原書("Programming Computer Vision with Python")の校正前の原稿は著者のJan Erik Solem氏がCreative Commons Licenseの元で公開されています。。

回答中で使用する画像の多く、また、sfm.py, sift.py, stereo.pyなどの教科書中で作られるファイルはこれらのページからダウンロードできます。
オライリージャパンのサポートページ
Jan Erik Solem氏のホームページ
相川愛三氏のホームページ

ここで扱うプログラムはすべてPython2.7上で動作を確認しています。一部の例外を除きほとんどがJupyter上で開発されました。また必要なライブラリーを適宜インストールしておく必要があります。Jupyter以外の点は、上記の書籍で説明されています。

結構な分量なので章毎に記事をアップしています。今回は第4章です。

4章 カメラモデルと拡張現実感

ピンホールカメラモデルを中心としたカメラモデルと、ヴァーチャルリアリティへの応用です。

4.1.3のfactorの定義が教科書のままだとうまくいきませんでした。

4.5.1 Modify the example code for the motion in Figure 4.2 to transform the points instead of the camera. You should get the same plot. Experiment with different transformations and plot the results.

Figure4.2のコードをカメラの代わりに点を動かすように変えろ

私の回答

https://github.com/moizumi99/CVBookExercise/blob/master/Chapter-4/CV%20Book%20Ch4%20Ex1.ipynb

import camera
from PIL import Image
from numpy import *
from pylab import *
import numpy as np

points = loadtxt('house.p3d').T
points = vstack((points, ones(points.shape[1])))

r = 0.05*np.random.rand(3)
rot = camera.rotation_matrix(r)
P = hstack((eye(3), array([[0],[0],[-10]])))
cam = camera.Camera(P)
x = cam.project(points)
figure()
plot(x[0], x[1], 'k.')
show()
figure()
for t in range(20):
    cam.P = dot(cam.P, rot)
    x = cam.project(points)
    plot(x[0], x[1], 'k.')
show()

image.png

image.png

P = hstack((eye(3), array([[0],[0],[-10]])))
cam = camera.Camera(P)
p = points.copy()
figure()
for t in range(20):
    p = dot(rot, p)
    x = cam.project(p)
    plot(x[0], x[1], 'k.')
show()

image.png

上の図がカメラを動かしたもの(教科書の例)、下の図が逆に点の方を動かしたものです。同じになっているようです。あってるようですね。

4.5.2 Some of the Oxford multi-view datasets have camera matrices given. Compute the camera positions for one of the sets an plot the camera path. Does it match with what you are seeing in the images?

私の回答

オックスフォードのデータベースはこちらです。
http://www.robots.ox.ac.uk/~vgg/data/data-mview.html
今回はここからMerton College Iをダウンロードして使いました。

https://github.com/moizumi99/CVBookExercise/blob/master/Chapter-4/CV%20Book%20Ch4%20Ex%202.ipynb

なお、教科書内で作成したライブラリcamera.py1を使います。

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import camera

Pmatrix = [loadtxt('house.'+str(i).zfill(3)+'.P').T for i in range(10)]

tarray = []
for P in Pmatrix:
    c = camera.Camera(P.T)
    K, R, t = c.factor()
    tarray.append(list(t))
points = array(tarray).T
figure()
plot(points[0], points[1], 'k.')
show()

image.png

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(points[0], points[1], points[2])
show()

image.png

これが軌道らしいのですが、正しいのかなんなのかわかりません。そこで特徴点がカメラ視点に対して正常に移動するか見てみます。

points = loadtxt('house.p3d').T
points = vstack((points, ones(points.shape[1])))
for i in range(len(Pmatrix)):
    P = Pmatrix[i].T
    cam = camera.Camera(P)
    x = cam.project(points)
    figure()
    gray()
    im1 = Image.open('house.'+str(i).zfill(3)+'.pgm')
    imshow(im1)
    plot(x[0],x[1], 'k.')
    show()

image.png

image.png

image.png

特徴点がカメラ行列によって他の画像に正しくマッピングされているようです。

4.5.3 Take some images of a scene with a planar marker or object. Match features to a full frontal image to compute the pose of each image’s camera location. Plot the camera trajectory and the plane of the marker. Add the feature points if you like.

マーカーを使って、正面からとった画像と、他の画像の特徴点を対応付けて、カメラ位置を計算しましょう。カメラの軌道を計算しろ

私の回答

ニンテンドー3DSについてきたマーカーを使いましたが、あまりうまくいきませんでした。

https://github.com/moizumi99/CVBookExercise/blob/master/Chapter-4/CV%20Book%20Ch%204%20Ex%203.ipynb

from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import camera
import homography
import sift

def my_calibration(sz):
    row, col = sz
    fx = 3600.0*col/3120
    fy = 3545.0*row/4160
    K = diag([fx, fy, 1])
    K[0,2] = 0.5*col
    K[1,2] = 0.5*row
    return K

def cube_points(c, wid):
    """ Creates a list of points for plotting
    a cube with plot. (the first 5 points are
    the bottom square, some sides repeated). """

    p=[]
    #bottom
    p.append([c[0]-wid,c[1]-wid,c[2]-wid])
    p.append([c[0]-wid,c[1]+wid,c[2]-wid])
    p.append([c[0]+wid,c[1]+wid,c[2]-wid])
    p.append([c[0]+wid,c[1]-wid,c[2]-wid])
    p.append([c[0]-wid,c[1]-wid,c[2]-wid]) # same as first to close plot

    #top
    p.append([c[0]-wid,c[1]-wid,c[2]+wid])
    p.append([c[0]-wid,c[1]+wid,c[2]+wid])
    p.append([c[0]+wid,c[1]+wid,c[2]+wid])
    p.append([c[0]+wid,c[1]-wid,c[2]+wid])
    p.append([c[0]-wid,c[1]-wid,c[2]+wid]) # same as first to close plot

    #vertical sides
    p.append([c[0]-wid,c[1]-wid,c[2]+wid])
    p.append([c[0]-wid,c[1]+wid,c[2]+wid])
    p.append([c[0]-wid,c[1]+wid,c[2]-wid])
    p.append([c[0]+wid,c[1]+wid,c[2]-wid])
    p.append([c[0]+wid,c[1]+wid,c[2]+wid])
    p.append([c[0]+wid,c[1]-wid,c[2]+wid])
    p.append([c[0]+wid,c[1]-wid,c[2]-wid])

    return array(p).T

imname0 = 'marker_small0.jpg'
sift.process_image(imname0, 'im0.sift')
l0, d0 = sift.read_features_from_file('im0.sift')

im0 = array(Image.open(imname0))
l1 = []
d1 = []
for i in range(l0.shape[0]):
    x, y, s, w = l0[i]
    if (325<x<1285 and 160<y<1760):
        l1.append(l0[i]) 
        d1.append(d0[i])
l1 = array(l1)
d1 = array(d1)
print im0.shape
figure()
imshow(im0)
sift.plot_features(im0, l1, circle=True)
show()

正面からの画像
image.png

sift.write_features_to_file('marker_small0.sift', l1, d1)
imname1 = 'marker_small0.jpg'
imname2 = 'marker_small2.jpg'
l1, d1 = sift.read_features_from_file('marker_small0.sift')
sift.process_image(imname2, 'im2.sift')
l2, d2 = sift.read_features_from_file('im2.sift')
im1 = array(Image.open(imname1))
im2 = array(Image.open(imname2))
print im1.shape
print im2.shape
figure()
subplot(1, 2, 1)
imshow(im1)
sift.plot_features(im1, l1, circle=True)
subplot(1, 2, 2)
imshow(im2)
sift.plot_features(im2, l2, circle=True)
show()

image.png

matches = sift.match_twosided(d1, d2)
ndx = matches.nonzero()[0]
fp = homography.make_homog(l1[ndx,:2].T)
ndx2 = [int(matches[i]) for i in ndx]
tp = homography.make_homog(l2[ndx2, :2].T)

figure(figsize=(16, 16))
gray()
sift.plot_matches(im1, im2, l1, l2, matches)
show()

image.png

model = homography.RansacModel()
H = homography.H_from_points(fp, tp)
# camera calibration
K = my_calibration((2000, 1600)) # this is the size of the image
print "K"
print K


# 3D points at plane z=0 with sides of length 0.4
box = cube_points([0, 0, 0.2], 0.2)
print 'box'
print box.T

#project bottom square in first image
cam1 = camera.Camera( hstack((K, dot(K, array([[0], [0],[-1]])))))
print "cam1.P"
print cam1.P

#fist points are the bottom square
box_cam1 = cam1.project(homography.make_homog(box[:, :5]))

# use H to transfer points to the second image
box_trans = homography.normalize(dot(H, box_cam1))

# compute second camera matrix from cam1 and H
cam2 = camera.Camera(dot(H, cam1.P))
A = dot(linalg.inv(K), cam2.P[:, :3])
A = array([A[:,0],A[:,1],cross(A[:,0],A[:,1])]).T
cam2.P[:, :3] = dot(K, A)
print "cam2.P"
print cam2.P

#project with the second camera
box_cam2 = cam2.project(homography.make_homog(box))
im1 = array(Image.open(imname1))
im2 = array(Image.open(imname2))

figure()
imshow(im1)
plot(box_cam1[0, :], box_cam1[1, :], linewidth=3)
figure()
imshow(im2)
plot(box_trans[0,:], box_trans[1,:], linewidth=3)
show()

image.png

image.png

まともに使える画像が1個しかとれなかったので、カメラ軌道は描画できなかったのが残念です

4.5.5 Take a look at the online documentation for .obj model files and see how to use textured models. Find a model (or create your own) and add it to the scene.

.objファイルでのテクスチャの使い方を調べて、画像に追加しなさい。

私の回答

https://github.com/moizumi99/CVBookExercise/blob/master/Chapter-4/CV%20Book%20Ch%204%20Ex%205.ipynb

pythonのコードは書籍の4.5のものと同じなのでブログ上では省きます
新たにwood.mtlファイルを追加します。

wood.mtl
newmtl wood
  illum 0
  Ka 1.0 1.0 1.0
  Kd 1.0 1.0 1.0
  Ks 0.0 0.0 0.0
  d 1.0
  map_Kd wood.png

wood.pngは次のページからダウンロードして色合いを調整したものです。
http://www.oyonale.com/modeles.php

さらに、toyplane.obj内の以下の部分を

mtllib toyplane.mtl
g PlaneRigCenter_mesh
#usemtl material1
usemtl lightblue

以下のように変更します。

mtllib wood.mtl
g PlaneRigCenter_mesh
#usemtl material1
usemtl wood

他にもusemtl lightblueとなっている部分はすべてusemtl woodにします。

結果は
image.png

無事テクスチャがつきました

4章終わり

これで4章終わりです。そろそろ中盤