あかすくぱるふぇ

同人サークル「あかすくぱるふぇ」のブログです。

三次元幾何

Pythonで、SpaceCarvingのためのボクセルデータ投影をやってみました。
Pythonだと本当にすっきりと書けていいですね。
#coding:utf-8
import numpy as np
import math
from PIL import Image

# ユーザー定義変数
voxelWidth = 50 # ボクセルの各次元の個数
spaceSize = 1000.0 # ボクセル群が埋める空間の大きさ[mm]
radius = 250.0 # 球の半径[mm]
imageShape = np.array((640, 640)) # 画像解像度
viewAngleDeg = 90.0 # 画角[deg]
cameraPos = np.array([0.0, 0.0, -500.0]) # カメラ位置[mm]

# 自動算出変数
voxelSize = spaceSize / voxelWidth # ボクセル1つの大きさ[mm]
principlePoint = imageShape / 2.0 # 主点
viewAngleRad = viewAngleDeg / 180.0 * math.pi # 画角[rad]
focalLengthPix = imageShape / 2 / math.tan(viewAngleRad / 2.0) # 焦点距離[pix]

# 中心からの距離を要素として持つボクセルデータを生成
print 'calc distance'
voxelNo = np.arange(voxelWidth**3) # ボクセル番号。座標計算に用いる。
coords3d = np.zeros((3, voxelWidth**3)) # 三次元座標。各列がボクセル、各行がX, Y, Z。
coords3d[0,:] = (voxelNo % voxelWidth - voxelWidth / 2 + 0.5) * voxelSize # X座標
coords3d[1,:] = (voxelNo / voxelWidth % voxelWidth - voxelWidth / 2 + 0.5) * voxelSize # Y座標
coords3d[2,:] = (voxelNo / voxelWidth / voxelWidth - voxelWidth / 2 + 0.5) * voxelSize # Z座標
distance = np.sqrt(np.sum((coords3d * coords3d), axis=0)) # 各ボクセルの中心からの距離

# 球のボクセルデータを生成。内部は1。外部は0
print 'make voxel data'
voxels = np.zeros(voxelWidth**3)
voxels[distance < radius] = 1
del distance

# 射影行列の生成
print 'make projection matrix'
intrinsic = np.array([[focalLengthPix[0], 0.0, principlePoint[0]],\
[0.0, focalLengthPix[1], principlePoint[1]],\
[0.0, 0.0, 1.0]])
extrinsic = np.hstack((np.eye(3), -cameraPos.reshape(1, len(cameraPos)).transpose()))
proj = np.dot(intrinsic, extrinsic)

# 投影
print 'projection'
coords3d = np.vstack((coords3d, np.ones(voxelWidth**3)))
coords2d = np.dot(proj, coords3d)
coords2d[:2,:] /= coords2d[2,:]
del coords3d

# 球内部の投影座標を抽出
print 'extract coord'
intersectCoord = np.array(coords2d[:2, (voxels == 1) *\
(coords2d[0,:] >= 0) *\
(coords2d[0,:] < imageShape[0]) *\
(coords2d[1,:] >= 0) *\
(coords2d[1,:] < imageShape[1])], dtype=np.int64)
del coords2d, voxels

# 画像生成
print 'make image'
image = np.zeros(imageShape)
for i in range(intersectCoord.shape[1]):
image[intersectCoord[0,i], intersectCoord[1,i]] = 255

# 画像表示
Image.fromarray(image).show()
以下、投影結果の画像です。
voxelWidth=400で、メモリ6GBくらい使います。
全ボクセルの三次元座標を行列として確保してるのが痛い。

・voxelWidth=100
100

・voxelWidth=200
200


・voxelWidth=300
300

・voxelWidth=400
400

「カメラ位置姿勢」と言った場合、普通は「世界座標系におけるカメラの位置と姿勢」のことを指す。

・カメラ姿勢とビュー変換
カメラ姿勢は、姿勢行列で表す場合と回転軸+回転角で表す場合とがある。
姿勢行列で表す場合、ビュー変換は、姿勢行列の逆行列をかけることで成される。
回転軸+回転角で表す場合、ビュー変換は、回転軸の周りを"マイナス回転角"することで成される。

・透視投影の式Xc=RXw+Tについて
回転行列Rは姿勢行列ではなく、姿勢行列の逆行列である。
Tは世界座標系におけるカメラ位置ではなく、ビュー座標系におけるカメラ位置であり、世界座標系におけるカメラ位置に対してRをかけた結果である。

・基底変換との関係
回転行列Rは、各行にビュー座標系の基底行列を並べたものと解釈できる。
上記同様、これは姿勢行列ではなく、姿勢行列の逆行列である。

・OpenGLによる変換
glRotate();glTranslate();の順番で変換する。
OpenGLの処理はスタックで実現されるので、glTranslate(世界座標系におけるカメラ位置)→glRotate()の順で実行される。
この順番で実行しないと、世界座標系におけるカメラ位置に対して回転行列Rがかけられない。

・glGetFloatv()による確認
glRotatef(-90, 0, 1, 0); glTranslatef(-3, 0, 0);だと、Tは(0, 0, -3, 1)になる。
glRotatef(180, 0, 1, 0); glTranslatef(-3, 0, 0);だと、Tは(3, 0, 0, 1)になる。
同じglTranslatef(-3, 0, 0)でも、Tは異なる。

ジンバルロックは以下のサイトで体験できます。
http://wonderfl.net/c/uTHS

ジンバルの軸には親子関係があります。
上記サイトの例では、Z軸周りに回転すると、Y軸とX軸も回転します。
Y軸周りに回転すると、X軸も回転します。
X軸周りに回転しても、他の軸は回転しません。

この親子関係の元で、Y軸周りに回転すると、
X軸はY軸周りに回転するが、Z軸は回転しないため、X軸とZ軸が重なります。
これがジンバルロックが起こる理由です。

この親子関係は、オイラー角における回転の順番(ヨー・ピッチ・ロール)と対応しています。
例えば、ヨーとロールについて考えると、ヨーの方が先に回転するので、
ヨーの回転はロールの回転軸に影響します。
しかし、ロールはヨーの後に回転するので、ロールの回転はヨーの回転軸には影響しません。

では、どんな時にジンバルロックが問題になるのでしょうか?
回転の順番が決まっているのであれば、Y軸周りに回転(ピッチ)した後に、
Z軸周りの回転(ヨー)とX軸周りの回転(ロール)が一緒になってしまっても問題はないはずです。
Y軸周りに回転(ピッチ)した時には、もうZ軸周りの回転(ヨー)は済んでるわけですから。

問題は、2つの方向AとBの間を補間して、少しずつ動かそうとした時に起こります。
少しずつ動かした結果、ある時刻において、Y軸周りの回転が90°になったとします。
その次の時刻に動かせる方向が、X軸周りとZ軸周りとで同じになってしまうのです。

以上が、ジンバルロックが起こる理由と問題になる状況の説明です。

↑このページのトップヘ