#まえおき
先日投稿した「実践コンピュータビジョン」の全演習問題をやってみたの詳細です。なお、あくまで「やってみた」であって、模範解答ではありません。「やってみた」だけなので、うまくいかなかった、という結果もあります。
この記事は「実践コンピュータビジョン」（オライリージャパン、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.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()
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()
上の図がカメラを動かしたもの（教科書の例）、下の図が逆に点の方を動かしたものです。同じになっているようです。あってるようですね。
###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をダウンロードして使いました。
なお、教科書内で作成したライブラリ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()
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()
これが軌道らしいのですが、正しいのかなんなのかわかりません。そこで特徴点がカメラ視点に対して正常に移動するか見てみます。
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()
特徴点がカメラ行列によって他の画像に正しくマッピングされているようです。
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についてきたマーカーを使いましたが、あまりうまくいきませんでした。
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()
正面からの画像
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()
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()
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()
まともに使える画像が１個しかとれなかったので、カメラ軌道は描画できなかったのが残念です
###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ファイルを追加します。
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にします。
結果は
無事テクスチャがつきました
##４章終わり
これで４章終わりです。そろそろ中盤