前回の記事では、世界地図から地球上の座標を計算し、それを平面にステレオ投影して、新たな画像を作り出すことに成功した。しかし、新たな画像上では、元の画像のピクセルが散乱してしまうことも確認した。ピクセルの散乱を防ぎ、すべてのピクセルが色を持つようにするには、どうしたらいいだろうか。
今までの手法では、元の世界地図から出発し、地図上の各点を投影することによって新たな画像を得た。このとき、有限の点を無限の平面に投影するという過程で、ピクセルが散乱してしまった。だから、散乱を防ぐためには、有限を有限に投影することが必要である。
これは、最終的な画像から出発し、画像上の各点を元の世界地図上に投影することで解決できる。最終的な画像の各ピクセルが、世界地図上のどの点からくるのかを逆にたどり、その世界地図上の点の色を、そのピクセルの色とするのだ。これによって、有限の大きさの画像を、有限の大きさの世界地図に投影することになるので、ピクセルの散乱がなくなり、すべてのピクセルが色を持つようになる。
この方法には、もう一つ利点がある。前回のやり方だと、地球のより大きい面積をカバーするには、画像の大きさを大きくする必要があったが、この方法では、任意の大きさの画像で、地球上の任意の範囲を画像に含めることができる。理由は、あとで説明する。
画像に含めたい表面積の割合をkとする。例えば、k=0.25は、地球の4分の1、つまり南極から南緯45°までの範囲、k=0.5は地球の半分、つまり南極から赤道までの範囲、k=0.75は地球の4分の3、つまり南極から北緯45°までの範囲を表す。
最終的な画像は、緯度がラジアンで \(-\frac{\pi}{2}\) から \(k\pi-\frac{\pi}{2}\)までの範囲を含む。緯度が \(k\pi-\frac{\pi}{2}\)となる点は、単位球から平面に投影されたときに原点から\(\tan(\frac{k\pi}{2})\)離れるので、最終的な画像のピクセルを、原点を中心とし、一片の長さが\(2\tan(\frac{k\pi}{2})\)の正方形の中に納まるようにスケールする。
最終的な画像上での座標が\((x, y)\)となる点を、その正方形内の点\((x',y')\)に移るようにスケールすると、
$$
\left( x', y' \right) = \left( \frac{2\tan \left( \frac{k\pi}{2} \right)}{d} \left( x - \frac{d}{2} \right), \frac{2\tan \left( \frac{k\pi}{2} \right)}{d} \left( y - \frac{d}{2} \right) \right)
$$
y座標にはマイナス符号がつくべきだが、前回も示した通り、結局投影すると鏡像になり、それを直すためにyに-1をかけるので、yはそのままでいい。
次に、平面上の点を単位球上の点に投影する。前回導いた公式を逆に用いると、経度\(\lambda\)、緯度\(\phi\)は、
$$
\lambda = \arctan2 \left(y', x' \right)
$$
$$
\phi = 2 \arctan \left(\sqrt{x'^2 + y'^2} \right) - \frac{\pi}{2}
$$
\(\lambda\)の範囲は\(-\pi\)から\(\pi\)なので、\(\arctan\)ではなく\(\arctan2\)を用いていることに注意する。
最後に、元の世界地図上の座標\((x'',y'')\)は、再び公式を逆に用いて、
$$
x'' = \frac{w}{2\pi}(\lambda + \pi)
$$
$$
y'' = \frac{h}{\pi} \left( -\phi + \frac{\pi}{2} \right) = h - \frac{2h}{\pi} \arctan \left( \sqrt{x'^2 + y'^2} \right)
$$
\(x''\)と\(y''\)は整数になるとは限らないので、その座標を四捨五入した点のRGB値を、最終的な画像の点\((x,y)\)に代入する。
このことをプログラムであらわしてみよう。
まず、必要な変数を定義する。dは最終画像の大きさ、kは画像に含める地球の面積の割合、imgは最終画像のピクセルを表す行列である。
d = 2000 k = .80 img_original = np.array(Image.open("world_map.jpg")) h = img_original.shape[0] w = img_original.shape[1] img = np.empty((d, d, 3), dtype=int) 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) x1 = w/(2*np.pi) * (np.arctan2(y0, x0) + np.pi) y1 = h - 2*h/np.pi * np.arctan(np.sqrt(np.square(x0) + np.square(y0)))
座標を四捨五入し、四捨五入によって世界地図からはみ出た点を世界地図上に戻すと、
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)
最後に、世界地図上に投影された点のRGB値を最終的な画像上の点に代入し、画像を表示すると、
img[y, x] = img_original[y1, x1] plt.imshow(img)
すべてのピクセルが埋め尽くされた、きれいな画像が得られた！kの値を変えることによって、画像に含める範囲を変えることができる。また、元の世界地図を逆さまにすることによって、北極を中心とした画像も作れる。
前回同様、緯度、経度ではなく、球面座標を利用する方法も考えた。このやり方だと、単純な公式が利用できるので、導かれる公式がよりきれいな形をしている。しかし、緯度、経度を利用したほうが直感的にわかりやすいので、今回はそちらを利用した。球面座標を利用する方法は、GitHub上に公開してあるので、参照してほしい。
最後に、今までのプロセスを関数にして、再利用可能にする。
def create_stereographic_projection(img, d, k, center='S'): h = img.shape[0] w = img.shape[1] if center == 'N': img = np.flip(np.flip(img, 0), 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) x1 = w/(2*np.pi) * (np.arctan2(y0, x0) + np.pi) y1 = h - 2*h/np.pi * np.arctan(np.sqrt(np.square(x0) + np.square(y0))) 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]
img = create_stereographic_projection(img_original, 2000, .5, 'S') plt.imshow(img)
これでプロジェクトの目的は達成できたが、次回はもう一歩踏み込んで、画像の中心を南極と北極以外にもできるようにする。北極の場合は画像をひっくり返すだけでよかったが、それ以外を中心とするとなると少し複雑になる。
