概要

360度画像(パノラマ画像, Equirectangular)を平面に展開するpythonツールです。雑な解説を載せますが、詳しいことはソースコードを読んでください。

解説

平面展開(__call__)

360度画像の指定方向の像を画像ファイルとして出力する。

最終的にcv2.remapを使うため、展開結果の画像の座標(i, j)が元の360度画像のどの座標(x, y)に対応しているかということを表す行列を作ることから始める。

いきなり360度画像の座標に割り当てるのは難しいので、最初は直径1の球面上の座標(x, y, z)に割り当てることから始める。「像の座標」と「視点中心からの角度」の対応関係は複数の射影方式(正射影・立体射影など)があり、それぞれmapping functionに従っている(英語版wikipediaなど参照)。いずれかのmapping functionを使って座標ごとの「視点中心からの角度」を算出し、z軸正の方向の球面上の点を割り当てる。

次に視野の角度を考慮して割り当てた点全てを一斉に回転させる。そして最後に回転後の座標から緯度(lon)経度(lat)を求め、対応する360度画像上の座標を求めて終了。

逆展開(inverse)

平面展開の逆で、平面展開された画像と同じサイズのレイヤー(透過pngを想定)を読み込み、元の360度画像に重なる形で引き戻す。

考え方は簡単で、平面展開の逆を行うだけ。最初に球面上の座標(x, y, z)がわかり、最後にレイヤー上の座標(x, y)が求められる。途中、視界の裏側の点を一様に無限遠点に取るという操作で不要な範囲を取り除いている(ソースコードを読んでください)。

コード

必要なpythonパッケージと360度画像は各自手に入れてください。360度画像はEquirectangularと検索すると出てくると思います。

from abc import ABCMeta, abstractmethod
import tkinter
import tkinter.filedialog

import numpy as np
import cv2


class PanoramaProjector(metaclass=ABCMeta):
    @staticmethod
    def parse_mapping_style(mapping_style):
        mode = 1
        if isinstance(mapping_style, str):
            if mapping_style == 'perspective':
                mode = 1
            elif mapping_style == 'stereographic':
                mode = 1 / 2
            elif mapping_style == 'equidistant':
                mode = 0
            elif mapping_style == 'equisolid_angle':
                mode = -1 / 2
            elif mapping_style == 'orthographic':
                mode = -1
        elif isinstance(mapping_style, tuple):
            f, n = mapping_style
            if f == 'tan':
                mode = 1 / n
            elif f == 'id':
                mode = 0
            elif f == 'sin':
                mode = -1 / n
        else:
            mode = mapping_style

        return mode

    def __call__(self, panorama, lon=0, lat=0, fov=105, width=800, height=600, mapping_style=1, fov_mode=True):
        lon = np.radians(lon)
        lat = np.radians(lat)
        if fov_mode:
            fov = np.radians(fov)

        [panorama_height, panorama_width, _] = panorama.shape
        spmap = self._spmap(lon, lat, fov, width, height, mapping_style=mapping_style, fov_mode=fov_mode)
        pxmap = self._pxmap(spmap, panorama_width, panorama_height)
        return cv2.remap(panorama,
                         pxmap[0].astype(np.float32),
                         pxmap[1].astype(np.float32),
                         cv2.INTER_CUBIC,
                         borderMode=cv2.BORDER_WRAP,
                         )

    @staticmethod
    def r(angle, mapping_style=1):
        # projection fomula of distance
        if mapping_style == 0:
            return angle
        elif mapping_style > 0:
            return np.tan(mapping_style * angle) / mapping_style
        else:
            return np.sin(mapping_style * angle) / mapping_style

    @staticmethod
    def r_inv(d, mapping_style=1):
        if mapping_style == 0:
            inv = d
        elif mapping_style > 0:
            inv = np.arctan(mapping_style * d) / mapping_style
        else:
            inv = np.arcsin(mapping_style * d) / mapping_style

        if isinstance(inv, np.ndarray):
            inv[inv > np.pi / 2] = np.float('inf')

        return inv

    def _spmap(self, lon, lat, fov, width, height, mapping_style=1, fov_mode=True):
        xyz = self._3dmap(fov, width, height, mapping_style=mapping_style, fov_mode=fov_mode)

        # x軸周りに -(pi/2-lat) 回転
        rx = np.array([[1, 0, 0],
                       [0, np.sin(lat), np.cos(lat)],
                       [0, -np.cos(lat), np.sin(lat)]])
        # z軸周りに -(pi/2-lon) 回転
        rz = np.array([[np.sin(lon), np.cos(lon), 0],
                       [-np.cos(lon), np.sin(lon), 0],
                       [0, 0, 1]])

        xyz = xyz.reshape([height * width, 3]).T
        xyz = np.dot(rx, xyz)
        xyz = np.dot(rz, xyz)

        # 球面上の座標(x, y, z)を緯度lonと経度latに変換する
        lon_map = np.arctan2(xyz[1], xyz[0])  # 東経
        lat_map = np.arcsin(xyz[2])  # 北緯

        lon_map = lon_map.reshape([height, width])
        lat_map = lat_map.reshape([height, width])

        return lon_map, lat_map

    def _3dmap(self, fov, width, height, mapping_style=1, fov_mode=True):
        # 画像の座標(i, j)から回転前の球面上の座標(x, y, z)への対応を作る

        # 像の範囲を求める
        if fov_mode:
            # fovが視野角を意味するモード
            # 視野角fovから(像の対角線の長さ)/2=r(fov/2)に変換
            fov = self.r(fov / 2, mapping_style=mapping_style)
        a = fov / np.sqrt(width ** 2 + height ** 2)
        vw = width * a  # z=1平面上の幅の半分
        vh = height * a

        # 画像の座標(i, j)から像の中の座標への対応
        x = np.tile(np.linspace(-vw, vw, width), [height, 1])
        y = np.tile(np.linspace(-vh, vh, height), [width, 1]).T

        # 球面上の座標に変換するための倍率を計算
        # self._r_inv(a2) が像の中心からの角度になっている
        a2 = np.sqrt(x ** 2 + y ** 2)
        a2 = np.sin(self.r_inv(a2, mapping_style=mapping_style)) / a2

        # 球面上の座標(x, y, z)を求める
        x *= a2
        y *= a2
        z = np.sqrt(1 - x ** 2 - y ** 2)

        return np.stack((x, y, z), axis=2)  # shape = height, width, 3

    @abstractmethod
    def _pxmap(self, sperical_map, panorama_width, panorama_height):
        pass

    def inverse(self, perspective, lon=0, lat=0, fov=105, panorama_width=8192, panorama_height=4096,
                mapping_style=1, fov_mode=True):
        lon = np.radians(lon)
        lat = np.radians(lat)
        if fov_mode:
            fov = np.radians(fov)

        [height, width, _] = perspective.shape
        pxmap_inv = self._pxmap_inv(panorama_width, panorama_height)
        x, y = self._spmap_inv(pxmap_inv, lon, lat, fov, width, height,
                               mapping_style=mapping_style, fov_mode=fov_mode)

        panorama = cv2.remap(perspective,
                             x.astype(np.float32),
                             y.astype(np.float32),
                             cv2.INTER_CUBIC,
                             borderMode=cv2.BORDER_TRANSPARENT)

        return panorama

    def _spmap_inv(self, pxmap_inv, lon, lat, fov, width, height, mapping_style=1, fov_mode=True):
        xyz, panorama_width, panorama_height = self._3dmap_inv(pxmap_inv, lon, lat)

        if fov_mode:
            fov = self.r(fov / 2, mapping_style=mapping_style)
        a = fov / np.sqrt(width ** 2 + height ** 2)
        vw = width * a
        vh = height * a

        # 裏側を取り除く
        f = xyz[2] <= 0  # filter
        xyz[:, f] = [[float('inf')], [float('inf')], [-1]]

        a2 = self.r(np.arccos(xyz[2]), mapping_style=mapping_style) / np.sqrt(xyz[0] ** 2 + xyz[1] ** 2)
        xyz *= a2

        xyz[0] += vw
        xyz[1] += vh

        xyz /= a*2

        x = xyz[0].reshape([panorama_height, panorama_width])
        y = xyz[1].reshape([panorama_height, panorama_width])

        return x, y

    def _3dmap_inv(self, pxmap_inv, lon, lat):
        lon_map, lat_map = pxmap_inv
        panorama_height = lon_map.shape[0]
        panorama_width = lat_map.shape[1]

        # x = sin(pi/2 - lat) * cos(lon) = cos(lat) * cos(lon)
        # y = sin(pi/2 - lat) * sin(lon) = cos(lat) * sin(lon)
        # z = cos(pi/2 - lat)            = sin(lat)
        x = np.cos(lat_map) * np.cos(lon_map)
        y = np.cos(lat_map) * np.sin(lon_map)
        z = np.sin(lat_map)
        xyz = np.stack((x, y, z), axis=2)

        # x軸周りに (pi/2-lat) 回転
        rxi = np.array([[1,           0,            0],
                        [0, np.sin(lat), -np.cos(lat)],
                        [0, np.cos(lat),  np.sin(lat)]])
        # z軸周りに (pi/2-lon) 回転
        rzi = np.array([[np.sin(lon), -np.cos(lon), 0],
                        [np.cos(lon),  np.sin(lon), 0],
                        [          0,            0, 1]])

        xyz = xyz.reshape([panorama_height * panorama_width, 3]).T
        xyz = np.dot(rzi, xyz)
        xyz = np.dot(rxi, xyz)

        return xyz, panorama_width, panorama_height

    @abstractmethod
    def _pxmap_inv(self, panorama_width, panorama_height):
        pass


# equirectangular
class EquirecProjector(PanoramaProjector):
    def _pxmap(self, spmap, panorama_width, panorama_height):
        i = panorama_width * (1 / 2 - spmap[0] / 2 / np.pi)
        j = panorama_height * (1 / 2 - spmap[1] / np.pi)

        return i, j

    def _pxmap_inv(self, width, height):
        lon_map = np.tile(np.linspace(np.pi, -np.pi, width), [height, 1])
        lat_map = np.tile(np.linspace(np.pi/2, -np.pi/2, height), [width, 1]).T

        return lon_map, lat_map


class PanoramaViewer(metaclass=ABCMeta):
    def __init__(self, image_path, projector, lon=0, lat=0, fov=105, width=800, height=600, mapping_style='perspective'):
        self._image = cv2.imread(image_path, cv2.IMREAD_COLOR)

        self._projector = projector

        self._lon = lon
        self._lat = lat
        self._width = width
        self._height = height
        self._mapping_style = projector.parse_mapping_style(mapping_style)
        self._d = self._projector.r(np.radians(fov) / 2, mapping_style=self._mapping_style)
        self._stride = np.degrees(self._projector.r_inv(self._d / 10))

    def __call__(self):
        if self._mapping_style > 0:
            p = 1 - np.cos(self._mapping_style * np.pi / 2)
        elif self._mapping_style < 0:
            p = -1 + np.cos(self._mapping_style * np.pi / 2)
        else:
            p = 0

        # f:p -> style
        # arccosに従って滑らかに遷移させる
        def f(p):
            if p > 0:
                return 2 / np.pi * np.arccos(1 - p)
            elif p < 0:
                return -2 / np.pi * np.arccos(1 + p)
            else:
                return 0

        while True:
            image = self._projector(self._image, self._lon, self._lat, self._d, self._width, self._height,
                                    mapping_style=self._mapping_style, fov_mode=False)
            cv2.imshow("image", image)
            key = cv2.waitKeyEx(0)

            if key == 2424832:  # left
                self._lon += self._stride
            elif key == 2490368:  # up
                self._lat += self._stride
                if self._lat > 90:
                    self._lat = 90
            elif key == 2555904:  # right
                self._lon -= self._stride
            elif key == 2621440:  # down
                self._lat -= self._stride
                if self._lat < -90:
                    self._lat = -90
            elif key == ord('z'):
                self._d /= 1.1
                self._stride = np.degrees(self._projector.r_inv(self._d / 10))
            elif key == ord('x'):
                self._d *= 1.1
                self._stride = np.degrees(self._projector.r_inv(self._d / 10))
            elif key == ord('m'):
                p -= 0.1
                if p < -1:
                    p = -1
                self._mapping_style = f(p)
            elif key == ord('n'):
                p += 0.1
                if p > 1:
                    p = 1
                self._mapping_style = f(p)
            elif key == ord('s'):
                self.output(image=image)
            elif key == ord('i'):
                root = tkinter.Tk()
                root.withdraw()
                layer_path = tkinter.filedialog.askopenfilename(filetypes=[('PNG', '*.png')])
                if layer_path is not '':
                    self.inverse(layer_path)
            elif cv2.getWindowProperty('img', cv2.WND_PROP_AUTOSIZE) == -1:
                cv2.destroyAllWindows()
                return

    def output(self, image=None, filename=None):
        if image is None:
            mapping_style = self._projector.parse_mapping_style(self._mapping_style)
            image = self._projector(self._image, self._lon, self._lat, self._d, self._width, self._height,
                                    mapping_style=mapping_style, fov_mode=False)
        if filename is None:
            filename = 'save/{}_{}_{}_{}_{}.jpg'.format(self._lon, self._lat, self._d, self._width, self._height)
        cv2.imwrite(filename, image)

    def inverse(self, layer_path, image=None, filename=None):
        if image is None:
            perspective = cv2.imread(layer_path, cv2.IMREAD_UNCHANGED)
            [panorama_height, panorama_width, _] = self._image.shape

            mapping_style = self._projector.parse_mapping_style(self._mapping_style)
            image = self._projector.inverse(
                perspective, self._lon, self._lat, self._d, panorama_width, panorama_height,
                mapping_style=mapping_style, fov_mode=False)
        if filename is None:
            filename = 'save/{}_{}_{}_{}_{}.png'.format(self._lon, self._lat, self._d, self._width, self._height)
        cv2.imwrite(filename, image)


if __name__ == '__main__':
    projector = EquirecProjector()
    viewer = PanoramaViewer('resource/panorama.png', projector)
    viewer()