laspyのインスタンスをOpen3Dのgeometryに変換し、jupyter上で可視化

## 初めに

点群データに限らずjupyterでデータ処理しながら、可視化を行うという作業はよくあると思いますが、**点群データ（特にlasデータ）をjupyter上で可視化**しようとすると、実はちょっと手間です。

- **Open3Dでは点群データを読み取るメソッドが存在しているが、lasは未対応**
- **pyntcloudではインスタンスをOpen3Dの点群インスタンスにコンバートする機能があるが、pyntcloud自体に、las読み取りのバグが存在する**
- **leafmapではサクッと読み取れるが、ビューワーのカスタマイズができないため、高度（Z座標）を使って勝手に色付けされる**

など、痒いところに手が届かない感じです。

この記事では、jupyterを使ってlaspyで読み取った点群データをうまいことnumpyに処理させながら、Open3Dのインスタンスを作成し、**ビューワーをカスタマイズ（色味・ポイントの大きさ・背景色など）する方法**を伝えていきます。

## データのダウンロードと加工

今回は以下の記事で作成したLAZ形式（LASの圧縮形式）のデータを利用していきます。

[古いlasデータをLAStoolsで最新版にアップグレードして、さらにLAZに圧縮にしてファイルサイズを1/4以下にする](https://qiita.com/nokonoko_1203/items/3d7078b6aff4084784cc)

が、お手持ちのlasデータでも問題ありません。

手頃なlasデータをお持ちでない方は、以下のサイトから沢山ダウンロードできますので活用してみてください。

https://www.geospatial.jp/ckan/dataset/shizuoka-2019-pointcloud/resource/d5e98a7b-f15c-45b0-bf40-0287f5b1de68

## laspyで読み込み

データは上記サイトからダウンロードし、LAZに圧縮し、さらにバージョンを上げた`08OF1348_1_4.laz`を利用していきますが、お手持ちのファイル名に適宜読み替えてください。

また、jupyterをインストールしていない！というかたはこちらを参考にしてみてください。

- [MacでGISデータ分析を始めるためにサクッとAnacondaとjupyter labをインストールしてみる](https://qiita.com/nokonoko_1203/items/3b5ceb911af6b22a449f)
- [jupyter labで最低限これだけは入れとけっていう拡張機能の紹介！](https://qiita.com/nokonoko_1203/items/edd010616a74468ea531)
- [Pythonで地理空間情報（GIS）やるために必要なパッケージ全部入りの「geospatial」が便利すぎた](https://qiita.com/nokonoko_1203/items/55faf4455f7986080e59)

まずは必要なモジュールをインポートします。

```python
import laspy as lp
import numpy as np
import open3d as o3d
import pandas as pd
```

laspyにファイルを読み込ませます。

[古いlasデータをLAStoolsで最新版にアップグレードして、さらにLAZに圧縮にしてファイルサイズを1/4以下にする](https://qiita.com/nokonoko_1203/items/3d7078b6aff4084784cc)　を参考にされた方はインストール済みかと思いますが、laspyはそのままではLAZ形式を取り扱うことはできず、**拡張版をpypiからインストールして利用する**必要がありますのでご注意ください。

```python
file_path = "08OF1348_1_4.laz"
laz = lp.read(file_path)
```

lazの情報を見てみましょう。

```python
laz
# <LasData(1.4, point fmt: <PointFormat(7, 0 bytes of extra dims)>, 10647410 points, 0 vlrs)>
```

ポイントが持っているデータを確認してみましょう。
point format 7という形式なので、xyzとrgbを持ってることが確認できます。

```python
list(laz.point_format.dimension_names)
# ['X',
#  'Y',
#  'Z',
#  'intensity',
#  'return_number',
#  'number_of_returns',
#  'synthetic',
#  'key_point',
#  'withheld',
#  'overlap',
#  'scanner_channel',
#  'scan_direction_flag',
#  'edge_of_flight_line',
#  'classification',
#  'user_data',
#  'scan_angle',
#  'point_source_id',
#  'gps_time',
#  'red',
#  'green',
#  'blue']
```

## numpyの配列を見てみる

読み込んだlazのXYZやRGBなどプロパティは、numpyのndarrayになっています。

lasデータは以下の記事に記載の通り、**XYZには整数しか格納**できません。
後ほど諸々いじっていきますが、正しい座標を取得するには**scaleやoffsetを考慮する必要がある**ことを頭に入れておきましょう。

[iPhoneでも使えるようになったLiDARの標準ファイル形式「LAS」ってどんなデータなの？を、PDAL/Laspyを使って調べてみる。](https://qiita.com/nokonoko_1203/items/8b8568545148acfe11d3)

そのうえで、X座標を見てみると以下のように整数値が格納されているのが確認できます。

```python
laz.X
# array([55200070, 55200090, 55200090, ..., 55599975, 55599985, 55599995],
#       dtype=int32)
```

Xだけ取り出すと、当然1次元の配列になっていることがわかります。

```python
laz.X.ndim
# 1
```

形状を確認するとn行1列の配列に見えますが、**numpyでは1次元配列は行と列の区別がないらしく、行ベクトルとして扱われる上、転置しても形状は変わりません**。

```python
laz.X.shape
# (10647410,)

laz.X.T.shape
# (10647410,)
```

余談ですが、n行1列の列ベクトルを明示的に示したいときはreshapeを利用して2次元配列にする必要があります。

```python
laz.X.reshape((-1, 1))
# array([[55200070],
#        [55200090],
#        [55200090],
#        ...,
#        [55599975],
#        [55599985],
#        [55599995]], dtype=int32)

laz.X.reshape((-1, 1)).shape
# (10647410, 1)
```

## scaleとoffsetを考慮して再計算し、XYZを持ったOpen3Dのインスタンスを作成する

Open3Dのインスタンスを作るには、**Open3Dの`Vector3dVector`という形式に合わせる必要がある**ので、形状を修正していきましょう。
とはいえ難しいことではなく、n行3列の形式に結合・変更してあげれば良いだけです。

numpyのstackメソッドには同じ形状を持った複数の配列を、配列に入れて渡すことで、入力よりも1次元多い配列が返されます。

例えばstackを行方向（axis=0）で使うと垂直に繋がるので3行n列になります。
前述の通り、1次元配列は行と列の区別がなく、行ベクトルであるため、垂直方向に要素が結合され、3行n列の2次元配列ができました。

```python
points = np.stack([laz.X, laz.Y, laz.Z], axis=0)
points.shape
# (3, 10647410)
```

ちなみにvstackというメソッドは縦（行）方向に結合するので、stack(axis=0)と同様の結合ができます。

```python
points = np.vstack([laz.X, laz.Y, laz.Z])
points.shape
# (3, 10647410)
```

列方向（axis=1）でstackを使うと当然n行3列の2次元配列ができます。

```python
points = np.stack([laz.X, laz.Y, laz.Z], axis=1)
points.shape
# (10647410, 3)
```

stackをaxis=0で利用した場合や、vstackを使った2次元配列に対して、transposeを利用すると、axis=1を利用した時と同じ結果を得られます。

```python
points = np.stack([laz.X, laz.Y, laz.Z], axis=0).transpose()
points.shape
# (10647410, 3)
```

似たようなメソッドにhstackというものがありますが、これは横（列）方向に結合するので、1次元配列をaxis=1でstackするのとはちょっとイメージが異なり、今回の場合は1次元配列のまま、配列長が伸びます。

```python
points = np.hstack([laz.X, laz.Y, laz.Z])
points.shape
# (31942230,)
```

ごちゃごちゃと話しましたが結果的に、**XYZの1次元配列3つに対して、stackをaxis=1で利用してn行3列の配列を作って**いきます。

```python
points = np.stack([laz.X, laz.Y, laz.Z], axis=1)
```

ただし、これではscaleとoffsetが考慮されていないので、座標が整数のままです。
これを修正していきましょう。

以下のようにして、lazのヘッダーが持つscalesとoffsetsを考慮して座標を計算していきましょう。

```python
data = pd.DataFrame(las.points.array)
data.columns = (x.lower() for x in data.columns)

data.loc[:, ["x", "y", "z"]] *= las.header.scales
data.loc[:, ["x", "y", "z"]] += las.header.offsets
```

その後、再計算されたポイントデータを利用して再度変数を作ります。

```python
points = np.stack([data.x, data.y, data.z], axis=1)
```

Open3Dの点群インスタンスを作成し、そこにデータを格納していきます。

```python
geom = o3d.geometry.PointCloud()
geom.points = o3d.utility.Vector3dVector(points)
```

これで点群を可視化できる基盤ができました！
実際にビューワーで確認してみましょう。

```python
o3d.visualization.draw_geometries([geom])
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/3afc1a5a-7182-4134-9cca-d131e4681909.png)

うーん。RGBを取り込んでいない上、デフォルトでは高度ごとに勝手に色分けされて表示されるため、色々微妙ですね。
また、このビューワー、**閉じるのがめっちゃめんどくさくて、起動するのに利用したインタープリターを強制的に停止させてあげないとなかなか閉じてくれません**。

ひとまずビューワーを改善していきましょう。

## Open3Dのビューワーをカスタマイズ

こんな関数を定義してあげましょう。

```python
import open3d.visualization.gui as gui
import open3d.visualization.rendering as rendering


def visualize(data):
    app = gui.Application.instance
    app.initialize()

    vis = o3d.visualization.O3DVisualizer("Open3D visualizer", 1024, 768)
    vis.show_settings = True
    vis.add_geometry("data", data)

    vis.reset_camera_to_default()

    app.add_window(vis)
    app.run()
```

詳しい設定方法は以下ですが、これでこのようなことが出来るようになっています。
[http://www.open3d.org/docs/release/python_api/open3d.visualization.O3DVisualizer.html](http://www.open3d.org/docs/release/python_api/open3d.visualization.O3DVisualizer.html)

- **1024 * 768の解像度で表示**
- **`Open3D visualizer`という名称でビューワーを表示**
- **点群のサイズなど変更できるメニューを追加**

字面だけでは分かりづらいと思いますので、実際に見てみましょう！

```python
visualize(geom)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/9ea9d066-a279-a5fc-b08d-5c390297e933.png)

だいぶリッチになりましたね！
ただ、肝心の色情報がないので、そちらを追加していきましょう。

## RGBをインスタンスに付与

lasの仕様として、RGBはuint16じゃないといけなかったはずですが、たまに守られていないデータもあるっぽいので、正規化していきます。
uint8だった場合は値に256を掛けていきます。

```python
color_dtypes = [data[column_name].dtype for column_name in ["red", "green", "blue"]]
color_dtype = color_dtypes[0]

if color_dtype == "uint8":
    data.loc[:, ["red", "green", "blue"]] *= 256
    data = data["points"].astype({"red": "uint16", "green": "uint16", "blue": "uint16"})
```

その後pointと同様に2次元配列にまとめて、Vector3dVectorに変換して先ほど作成したgeometryに格納していきましょう。
このとき、colorの範囲に制限（0~1？）があるので、RGB値を65535で割っていきましょう。
先ほど正規化したおかげで、ちゃんとしたカラーで表示されるようになるはずです。

```python
colors = np.stack([data.red, data.green, data.blue], axis=1)

geom.colors = o3d.utility.Vector3dVector(colors / 65535)
```

では可視化していきましょう。

```python
visualize(geom)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/4315c309-a703-1c6f-eede-7b1a1bea7563.png)

色もしっかり表示されましたね！

## 終わりに

ということで、jupyterとlaspyとOpen3Dを利用して、ただ表示するだけではなく、Open3Dのgeometryの作り方を学びながら、かつビューワーをカスタマイズして表示させていきました。

これでデータ処理しながら点群を可視化できるようになったと思いますので、みなさんどんどん点群を活用していきましょう！