前回までの記事で、世界地図から地球上の座標を計算し、それをステレオ投影によって平面上に投影し、新しい画像を得ることにした。これまでの手法では画像の中心は南極か北極にしかできなかったが、今回は画像の中心を地球上の任意の点にすることに挑戦したい。
もう一度、以下の画像を見てほしい。ステレオ投影では、球の上から光を当てて平面に投影するので、球の下が中心になることがわかる。よって、地球を回転させて中心にしたい点が下に来るようにしてから投影すれば、その点が中心となる画像が得られる。
前回、世界地図の画像から始めて最終的な画像に投影するのではなく、最終的な画像の座標から逆算して元の座標を得ることで、隙間のない画像が得られることを学んだ。今回も、この逆算の手法を使う。
まず、最終的な画像の座標が対応する、ステレオ投影前の球上の座標を求める。この球の下端は、最終的な画像の中心にしたい点であり、南極はどこかほかのところにある。そしてこの球を、南極が下に、北極が上に、そして下端が元の地球上の位置に来るように回転させる。そうすると地球が正しい方向を向くので、元の世界地図上の座標が得られる。
言葉で説明するのが難しいので、絵をかいたらわかるかな？
このことを数式で表す。ただし、今回は少し違う座標系を用いる。今までは緯度と経度を用いて球上の座標を表現していたが、今回は球面座標系というのを用いる。球面座標系は数学においてより一般的に用いられている座標系であり、3次元直交座標へシンプルな式で変換できる公式が存在する。そして今回、回転には直交座標と回転行列を用いるので、球面座標系を使うのが便利である。
球面座標系を、地球の中心を原点、地球の半径を1とし、x軸の正の方向が本初子午線の方向を向くように定める。中心にしたい点の緯度と経度をそれぞれ\(lat\)、\(long\)とすると、その点の球面座標は、
\[
(\phi, \theta) = \left( \frac{\pi}{180} (90 - lat), \frac{\pi}{180} long \right)
\]
である。球の下端をその中心にしたい点に来るように回転させるには、球をx軸を中心に反時計回りに\(\pi-\phi\)、z軸を中心に時計回りに\(\frac{\pi}{2} - \theta\)（反時計回りに\(\theta - \frac{\pi}{2}\)）回転させる。
回転には直交座標系を使う。球面座標\( (\phi, \theta)\)は、以下の公式によって直交座標に変換できる。
\[
(x, y, z) = (\sin\phi \cos\theta, \sin\phi \sin\theta, \cos\phi)
\]
x軸中心に反時計回りに\(\theta\)回転させる回転行列は、
\[
M_x = \begin{pmatrix}
1 & 0 & 0 \\
0 & \cos\theta & -\sin\theta \\
0 & \sin\theta & \cos\theta\\
\end{pmatrix}
\]
z軸中心は
\[
M_z = \begin{pmatrix}
\cos\theta & -\sin\theta & 0 \\
\sin\theta & \cos\theta & 0 \\
0 & 0 & 1\\
\end{pmatrix}
\]
\(\theta\)に適当な値を代入して、これらの回転行列を利用することで、地球を回転させることができる。
回転させた後は、直交座標をふたたび球面座標に戻す。
\[
(\phi, \theta) = (\arccos z, \arctan2(y, x))
\]
最後にこれを緯度と経度に直し、適当にスケールすることで、元の世界地図上の座標が得られる。
これが必要な数学のすべてだ。これらをプログラムで実装する。
まず必要な変数を定義する。これらは以前の記事ですでに説明してある。
d = 2000 k = .80 img_original = np.array(Image.open("world_map.jpg"))
次に、いくつかのヘルパー関数を定義する。to_cartesianは球面座標を直交座標に、to_sphericalは直交座標を球面座標に変換する関数で、rotateは座標のリストvを、x軸中心にphi、z軸中心にtheta回転させる関数である。
def to_cartesian(phi, theta): x = np.sin(phi) * np.cos(theta) y = np.sin(phi) * np.sin(theta) z = np.cos(phi) return x, y, z def to_spherical(x, y, z): phi = np.arccos(z) theta = np.arctan2(y, x) return phi, theta def rotate(v, phi, theta): rotation_x = [ [1, 0, 0], [0, np.cos(phi), -np.sin(phi)], [0, np.sin(phi), np.cos(phi)] ] rotation_z = [ [np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1] ] rotation = np.matmul(rotation_z, rotation_x) return np.einsum('ij,klj->kli', rotation, v)
最後に、上で述べたプロセスをすべてcreate_stereographic_projection関数で実行する。
def create_stereographic_projection(img, d, k, lat, long): h = img.shape[0] w = img.shape[1] x, y = np.meshgrid(np.arange(0, d), np.arange(0, d), indexing='xy') factor = -2 * np.tan(k * np.pi / 2) / d x0 = factor * (x - d/2) y0 = factor * (y - d/2) r = np.sqrt(np.square(x0) + np.square(y0)) phi = 2 * np.arctan(1/r) theta = np.arctan2(y0, x0) # rotate phi0 = np.pi/180 * (90 - lat) theta0 = np.pi/180 * long x, y, z = to_cartesian(phi, theta) v = np.dstack((x, y, z)) v = rotate(v, np.pi - phi0, theta0 - np.pi/2) x = v[:,:,0] y = v[:,:,1] z = v[:,:,2] phi, theta = to_spherical(x, y, z) x1 = w/(2*np.pi) * (theta + np.pi) y1 = h/np.pi * phi x1 = np.int_(np.round(x1)) y1 = np.int_(np.round(y1)) x1 = np.maximum(x1, 0) y1 = np.maximum(y1, 0) x1 = np.minimum(x1, w - 1) y1 = np.minimum(y1, h - 1) return img[y1, x1]
これで必要な道具はそろった。ではさっそく動作を確認してみよう。
今回は、東京（北緯36度、東経140度）を中心に、地球の半分をステレオ投影した画像を作ってみる。
img = create_stereographic_projection(img_original, 2000, .5, 36, 140) # Japan plt.imshow(img)
確かに、日本が中心となっている！
プロジェクトのまとめ
以上が、Flat Earthプロジェクトの全貌である。まず、世界地図上のピクセル一つ一つに対して地球上の座標を計算し、それに北極から光を当ててステレオ投影することによって、新たな画像を得た。しかし、中心から離れるにつれてピクセルが散乱してしまうという問題が発生した。これは、最終的な画像のピクセルの座標から、元の世界地図上の座標を逆算することによって解決できた。さらに、地球をうまく回転させて座標を計算することによって、画像の中心を南極と北極だけでなく、任意の点にできるようにした。
このプロジェクトからはいくつかのことを学んだ。まず技術面では、NumPyの行列計算についての理解を深め、よりきれいなコードで行列を操ることができるようになった。この経験は、これからPythonで機械学習のプログラミングをするうえで役に立つだろう。また、数学を駆使して事象を数式で表し、それをプログラムで実装することによって、思い通りの結果が得られることにも強い感動を覚えた。数学を使って何かを作るという経験は初めてなので、非常にいい刺激になった。
久しぶりに熱中できるプロジェクトができたし、結果も完璧といっていいほど思い通りのものになったのでとてもうれしい。また何か面白いアイデアが思いついたら挑戦したい。
