equirectangular画像をcubemap画像に変換する方法

## equirectangular画像とcubemap画像
equirectangular画像は360度カメラで撮影したデータとして、よく用いられるフォーマットです。下画像のように360度の景色を1枚の写真に収めることができます。
![original_equirectangular.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/663246/efb30e01-7df2-e9dc-3505-4c60075ba627.png)

ただ、equirectangular画像は真上と真下のピクセルが必要以上に引き伸ばされて**歪みが大きくなる**という問題があります。この歪みを解消するには、360度画像で撮影した景色を球とすると、それを囲むような立方体を考えて、その6面に投影するという方法があります。そのような6面画像は**cubemap画像**(下画像)と呼ばれています。
![cube_map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/663246/1bfb6177-5d69-9803-4774-cf8df2d9b359.png)

今回はequirectangular画像からcubemap画像への変換をPythonで実装してみました。データセットとしては、Youtuberの[Young360](https://www.youtube.com/watch?v=gsQcw1NZDR4)さんが公開しているデータを使用しています。

## Algorithm
6枚の画像を作成し、最後に1枚の画像にまとめることでcubemapを作成します。今回は、真下画像の作成を例として説明します。同じ考え方で他の面の画像も作成することができます。

### 1. 撮像面を作成する
真下画像の場合は、z軸方向のベクトルを固定して球の真下に撮像面を配置します。撮像面のサイズは、出力画像サイズとなります。

### 2. 3次元ベクトルの緯度・経度を計算する
球の中心から1で作成した撮像面の1ピクセルごとへの3次元ベクトルの緯度・経度を計算します。下画像のφが緯度でθが経度です。
![plane.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/663246/8e949fc7-2010-d18b-4741-9912d28744ea.jpeg)

緯度を計算するには、z成分と3次元ベクトルのスカラーのarccosを計算します。

```math
\phi=\arccos\frac{z}{\sqrt{x^2+y^2+z^2}}
```

![phi.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/663246/1d9146a1-1be6-c6aa-1d7f-7c2089c95051.jpeg)

経度を計算するには、x, y成分のarctanを計算します。

```math
\theta=\arctan\frac{y}{x}
```

![theta.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/663246/10e458cc-5de6-8d66-2f1f-547f5403d7c5.jpeg)

### 3. 緯度・経度を用いて、出力画像のピクセルごとに、equirectangular画像(元画像)上の位置を算出する

2で緯度と経度がわかったので、その値を用いて、元画像(equirectangular画像)ではどのピクセルになるかを計算します。
equirectangular画像では横方向には360°、縦方向には180°の情報があるので、前のステップで得た緯度と経度を、横方向と縦方向の全体の角度(2πとπ)で正規化した後、2方向の全体のピクセルを掛けることで、equirectangular画像上でのピクセル値が得られます。

```math
output_x[pixel]=\frac{\theta}{2\pi}\times input_{width}
```

```math
output_y[pixel]=\frac{\phi}{\pi}\times input_{height}
```

あとは、このピクセル値を使ってopencvのremapなどで変換すれば、真下画像を作成できます。

```python3
def get_theta(x, y):
    if y < 0:
        theta = (-1) * np.arctan2(y, x)
    else:
        theta = 2 * math.pi - np.arctan2(y, x)

    return theta


def create_equirectangler_to_bottom_and_top_map(input_w, input_h, output_sqr, z):
    map_x = np.zeros((output_sqr, output_sqr))
    map_y = np.zeros((output_sqr, output_sqr))
    for row in tqdm(range(output_sqr)):
        for col in range(output_sqr):
            x = row - output_sqr / 2.0
            y = col - output_sqr / 2.0

            rho = np.sqrt(x * x + y * y + z * z)
            norm_theta = get_theta(x, y) / (2 * math.pi)
            norm_phi = (math.pi - np.arccos(z / rho)) / math.pi
            ix = norm_theta * input_w
            iy = norm_phi * input_h

            if input_w <= ix:
                ix = ix - input_w
            if input_h <= iy:
                iy = iy - input_h

            map_x[row, col] = ix
            map_y[row, col] = iy

    return map_x, map_y
```

## All code
6面全てを変換して、cubemapを作成したコードは[こちら](https://github.com/Hiroaki-K4/3d_computer_vision/blob/main/equirectangular_to_cubemap/equirectangular_to_cubemap.py)です。

## Dataset
- [Insta360 X3 vs X2 Camera Review | DOWNLOAD REAL Footage](https://www.youtube.com/watch?v=gsQcw1NZDR4)
- [Dataset in Google Drive](https://drive.google.com/drive/folders/1HnwAKZAkJ8ZKL8S-IBMWmDN71RWK3ll8)

## Referrence
- [Converting to/from cubemaps](http://paulbourke.net/panorama/cubemaps/)