LAZ形式の点群データがLASやCSVと比較してなぜ軽量なのか、をまとめてみる

みなさん点群データ使ってますか！！！
最近は神奈川県も点群データを出してくれていますね！目指せ全国制覇！

点群データに関してはいろんな視点で記事を書いているので見てみてください！

https://qiita.com/nokonoko_1203/items/a8d58c12c31faf55a7f5

https://qiita.com/nokonoko_1203/items/ddd14090d4141e0476f1

https://qiita.com/nokonoko_1203/items/c6e2f47da74f92ec3666

点群データは（機材はそれなりに高価ですが）ドローンを飛ばすだけで高域かつ人があまり立ち入れない場所の地形データが入手できます。
そのため基本的には測量に利用されていますが、**歴史的建造物の3Dアーカイブ**に利用されたり、**iPhoneでも取得できたり**と利活用の幅が広まっています！

自動運転でも積極的に利用されていたと記憶していますが、ただ、その**ファイルサイズの大きさから意外と取り回しが難しく**、他の手法に転換せざるを得なかった〜なんていう話も聞きました。

それもそのはずで、例えば測量成果として利用するためには1m2あたり、4点以上の密度を持つ必要があり、国土基本図図郭のレベル500の領域（300m × 400m）で**1000万ポイント**を超えることもあります。

![Pasted image 20241028213121.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/b9f59ce5-1176-a73a-d723-3b519caa53f7.png)
（G空間情報センターのダウンロード画面（矩形がレベル500の範囲）: https://www.geospatial.jp/ckan/dataset/kanagawa-2022-pointcloud ）

点群データはXYZの座標値の他にRGBのカラー情報やそれ以外の情報も大量に保持することができます。
このためこの程度の領域でも100MBを超えてきます。
（というか、**平気で500MB**あったりします）

なので、配布する場合でも利用する場合でも基本的にはLASファイルを圧縮した**「LAZ」ファイル**の方が良いです！
**ファイルサイズも1/10**になったりします。

ただ、なぜかあまり利用されていないように見えるので、今回はなんでLAZファイルの方が良いのか！っていう点を熱く語っていこうかと思います！
ついでに、変換手法も載せときます！

## CSVの点群はとても重たい

それなりの頻度で、点群データがCSVなどのテキスト形式で保持されていることがあります。
テキストファイルは人間にとって可読性が高いので、配布されているとすぐにデータが確認できてわかりやすい！

…ということはあまりなく、**ファイルサイズが大きすぎてエディターで開くのにすら苦戦**したり、**ヘッダー行が存在せず仕様の記載もないため、どの列がどのデータなのかよくわからない**といったケースがかなり多いです。

ファイルサイズの比較ですが基本的には以下のようになります。

```
CSV > LAS >>>> LAZ
```

CSVとLASについては、いろんな事情（後述）でCSVとほぼ変わらなかったりしますが、CSVはテキストファイル、LASはバイナリファイルのため**読み込みはLASファイルの方が高速**です。

CSVはなぜファイルサイズが大きくなるのか、について具体例を挙げていきます！

### なぜCSVは大きくなるのか

CSVは**テキスト形式**でデータを保存します。
このため、**各データが文字列として記録**されるため、**1文字1バイト**消費しますし、数値を文字列に変換する際、**余分な文字（例えば小数点やカンマ、改行コードなど）が含まれ**、ファイルサイズが大きくなります。

広く公開されている点群データは測量業務の一環として計測されるため、公共測量で利用される「平面直角座標系」で座標値を格納することが多いです。

平面直角座標系では、単位がメートルで表されます。このため、32ビット浮動小数点数（f32）で十分な精度を確保できます。RGBの色情報は8ビットの符号なし整数（u8）で表現できますが、経緯度で座標値を保持している場合、f32では精度が不足するため注意が必要です。

### ファイルサイズの推定

例として、XYZRGBのデータ構成で、XYZがf32、RGBがu8の場合、1000万ポイントの点群データを考えてみましょう。
各ポイントは以下のようになります：

- **XYZ座標**:
    - X座標: 「123456.7890」→ 約11文字
    - Y座標: 「1234567.890」→ 約11文字
    - Z座標: 「12345.67890」→ 約11文字
- **RGB値**:
    - R値: 「255」→ 3文字
    - G値: 「255」→ 3文字
    - B値: 「255」→ 3文字
- **区切り文字**: カンマ5個 → 5文字
- **改行コード**: 1文字（環境によっては2文字）

各値は文字数にするとおよそ：

- XYZ座標: 11 + 11 + 11 = **33文字**
- RGB値: 3 + 3 + 3 = **9文字**
- 区切り文字と改行コード: 5 + 1 = **6文字**
- **合計**: 34 + 9 + 6 = **49文字**

カンマや改行コードも含めると、1ポイントあたり約40〜50バイトになります。
したがって、1000万ポイントでは約400〜500MB、場合によってはそれ以上のサイズになります。

- 1文字あたり1バイトとして、1ポイント = 約49バイト
- 1000万ポイントの場合: 49バイト/ポイント × 1,000万ポイント = **約490MB**

CSVは以下のような理由でファイルサイズが増大していることがわかりました。

- 数値の文字列化: CSVでは、数値データも文字列として保存します。これは、数値をバイナリ形式で保存する場合に比べて、各数値が占めるバイト数が増加します。例えば、32ビットの浮動小数点数（f32）はバイナリ形式では4バイトですが、文字列では桁数に応じてバイト数が増加します
- 区切り文字の存在: CSVでは、フィールド間をカンマやタブで区切ります。これらの区切り文字もデータサイズに含まれ、全体のファイルサイズを増加させます
- 改行コードの使用: 各レコードの終わりには改行コード（`\n`や`\r\n`）が含まれます。これもファイルサイズの増加要因です

## LASファイルのメリット

前述の通り、CSVでは1000万ポイントを保持すると500MB程度になってしまうだろうということが予想つきます。
さらに1000万ポイントあっても、300m × 400m程度の範囲なので、より広域で利用しようと思うと取り扱いがとても難しいです。

ではLASファイルはどうでしょうか？

### ScaleとOffsetの利用

LASでは、ヘッダー部分にスケール（scale）とオフセット（offset）の情報を格納します。
これにより、XYZ座標を整数値として保持できます。具体的には、座標値をスケールで割り、オフセットを引いた後、整数化します。
これにより、データの精度を保ちながら、数値を小さくし、保存するビット数を削減できます。

### バイナリ形式の利点

LASはバイナリファイルであり、データをそのままの形式（例えばf32）で保存できます。
テキスト形式のように数値を文字列に変換する必要がないため、ファイルサイズを大幅に抑えることができます。

**ファイルサイズの削減例**:

- **1ポイントあたりのバイト数**:
    - XYZ座標: int32 × 3 → 4バイト × 3 = **12バイト**
    - RGB値: uint8 × 3 → 1バイト × 3 = **3バイト**
    - **合計**: 12バイト + 3バイト = **15バイト**
- **1000万ポイントの場合**: 15バイト/ポイント × 1,000万ポイント = **約150MB**

CSV形式の約490MBと比較すると、LAS形式では約150MBと、大幅にファイルサイズが削減されます。

## 注意点

ただし、**実際にはCSV形式とLAS形式で似たようなファイルサイズになっている**ことがあります。
（実際に、神奈川県の点群データでは、TXT（CSV）をダウンロードしても、LASでダウンロードしてもほぼ変わりませんでした。）
これにはちゃんと理由があります。
LASにはポイントフォーマットと呼ばれる「**必ず保有しているべきカラム一覧**」が定められており、測量機器などの都合により該当する属性を取得できなかったとしても、必ずデータとして含んでいなければいけない情報があります。
（つまり、0やnullなどの列が多いということ）
例えばLAS v1.4のPoint Data Record Format 0ではXYZの他にClassificationやPoint Source IDなど、必須項目があり、最低でも1ポイント24バイト消費します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/203944/030ea355-5d3f-e278-806b-422cd0fcbdbe.png)
（LAS Specification 1.4 - R14: https://www.asprs.org/wp-content/uploads/2019/03/LAS_1_4_r14.pdf ）

これに対し、CSVはスキーマが定義されていないため、とにかく**ファイルサイズを減らすことを目的とし、XYZ以外の情報を持たない！**という選択をすることができます。
そうなると、データ次第ではCSVのファイルサイズがLASに匹敵するレベルまで削減されることもあります。
（当然、情報量が足りないので純粋な比較ではないですが）

ファイルサイズのみを比較する場合、この点は考慮する必要があるでしょう。

## LAZファイルのメリット

では、LAZファイルはどのようになっているのでしょうか？
LAZファイルの仕様は以下に記載があるため、見ていきましょう！
https://downloads.rapidlasso.de/doc/LAZ_Specification_1.4_R0.pdf

と思って見ていたんですが、英語がめちゃくちゃ苦手かつ技術力不足により、概要しかわかりませんでした。
**※間違っている説明も多いと思うので、以下の情報はあくまで参考情報として受け取ってもらえると助かります。仕様書をお読みの際には、ぜひコメントいただけると嬉しいです！**

### LAZの圧縮プロセス

LAZ形式では、データの特性によって動的にエンコーディング手法が選定されるようでした。
整数値は整数エンコーダーが利用され、分類などの場合はシンボルエンコーダーが利用されるとかそう言った手法らしいのですが、詳しいことは読み取れなかったので、複数あるエンコーダーの中身を紹介していきます。

#### シンボルエンコーダー

シンボルエンコーダーは、データをシンボル（記号）として扱い、頻繁に出現するシンボルに対して短いコード、あまり出現しないシンボルに対しては長いコードを割り当てることで圧縮します。これにより、データ全体を短縮し、効率的な圧縮を実現します。
ハフマンエンコーディングによく似ています。

例えば、データが「100, 200, 300」しか出てこないのであれば「A、B、C」の3つのシンボルを割り当てることができます。出現頻度が「A > B > C」の順である場合、シンボル「A」に対しては短いコード、「C」に対しては長いコードを使います。Aにはより短いビット数のシンボル「0」が割り当てられ、Bは「1」、Cは「10」といったイメージです。
圧縮されたデータストリームを復号する際には、割り当てられた範囲に基づいて各シンボルを正確に復元します。
（データストリームについても、仕様書に記載がありました。）

#### ビットシンボルエンコーダー

ビットシンボルエンコーダーは、2つのシンボル（0と1）を圧縮するための特化型エンコーダーです。このエンコーダーは特にビットデータを圧縮するために使われ、例えば整数の圧縮において最初のビットを符号化する際に使われます。
ランレングスエンコーディングに似ていますが、こちらは2つのシンボルにしか利用されません。

例えば、ビット0がビット1よりも多く出現する場合、ビット0に対して短いコード、ビット1に対して長いコードを割り当てます。

#### 整数圧縮エンコーダー

整数圧縮エンコーダーは、データ間の差分を圧縮します。具体的には、連続するデータポイントの間の値の差分を符号化することで、整数の変動が少ない場合に効率的な圧縮を行います。
デルタエンコーディングに似ています。

このエンコーダーは、ある整数値とその予測値（通常は前のデータポイント）との差分（=残差）を符号化します。この差分が小さい場合に、圧縮効率が向上します。符号化された整数は、最も上位のビットがシンボルエンコーダーによって圧縮され、残りのビットはRawエンコーダーで保存されます。
（Rawエンコーダーは何もしないエンコーダーらしいです。）

例えば、連続する整数値のデータがあり、前の値が100で現在の値が102であるとします。整数圧縮エンコーダーは、差分（=2）だけを符号化することで、データを効率的に圧縮します。

### 事前処理による更なる圧縮

上記のエンコーダーの組み合わせにより、データが1/5~1/10程度まで圧縮されますが、さらに[空間充填曲線（space-filling curve）](https://ja.wikipedia.org/wiki/%E7%A9%BA%E9%96%93%E5%85%85%E5%A1%AB%E6%9B%B2%E7%B7%9A)により事前に点群をソートしておくことで更なる圧縮を行うことができます。

pdalを利用すると、LASからLAZに変換する際に**モートンオーダー**を選択することができます。
- 参考
    - [八分木（モートンオーダー）を使ってエリアを分割して処理負荷を軽減する](https://edom18.hateblo.jp/entry/2017/07/28/083153)
    - [filters.mortonorder](https://pdal.io/en/latest/stages/filters.mortonorder.html)


点群データは、**現実世界をスキャンして得るデータ**という特性上、データに計測時間を格納することができますが、**計測時間が近いデータは物理的にも近い座標**にあり、尚且つ**データの差分が少ない**ことが多いです。
このため、空間充填曲線よりも、デルタエンコーディングが効きやすく、圧縮効率が良くなるそうです。
なので、計測時間を持っているデータであれば、モートンオーダーではなく、**計測時間で事前にソート**してあげると良いでしょう！
（PDALには存在しないので、何かしら別のツールを利用する必要があります。）
- 参考
    - https://groups.google.com/g/lastools/c/OipcNpx1HoM
    - https://gis.stackexchange.com/questions/330019/further-compression-of-lidar-data-besides-laszip

### エンコードの流れ

結果として、以下のような流れでエンコーダーが協調して動作しているように見えます。

- そもそもLAS自体の仕様として、scaleとoffsetをヘッダーに格納することでXYZ座標は整数値で保持できるようになっている
- 空間情報を考慮してソートする
- カラムごとに予測モデルを構築し、予測値を算出する
    - ソートしているので、予測モデルの精度が上がる
- ソートしているので、予測結果と実際の値の残差が小さくなる
- 残差が小さく整数値なので、特定の出現頻度が高くなり、シンボルエンコーダーが効きやすい
- （同一シンボルも並びやすいので、ランレングスエンコーディングをするとさらに圧縮されるのではーとか思ったけど、やってなさそう？）

## 点群圧縮の具体例

すごいざっくりですが、こんなLASデータがあるとします。

|ID|X整数値|Y整数値|Z整数値|
|---|---|---|---|
|1|0|0|0|
|2|10|5|5|
|3|20|10|10|
|4|30|15|15|
|5|40|20|20|
|6|50|25|25|
|7|60|30|30|
|8|70|35|35|
|9|80|40|40|
|10|90|45|45|

なんらかの予測モデルを構築し、デルタエンコーディングで残差を求めます。

|ID|X整数値|前の値との差分|
|---|---|---|
|1|0|N/A|
|2|10|10 - 0 = **10**|
|3|20|20 - 10 = **10**|
|4|30|30 - 20 = **10**|
|5|40|40 - 30 = **10**|
|6|50|50 - 40 = **10**|
|7|60|60 - 50 = **10**|
|8|70|70 - 60 = **10**|
|9|80|80 - 70 = **10**|
|10|90|90 - 80 = **10**|

ここまで綺麗に並んでいると残差は全て0になります。

|ID|予測値（X）|実際のX整数値|残差（X）|
|---|---|---|---|
|1|N/A|0|N/A|
|2|0 + 10 = 10|10|10 - 10 = **0**|
|3|10 + 10 = 20|20|20 - 20 = **0**|

この値を、X座標として格納しておきます。

シンボルエンコーダーを通すと、残差0の出現回数が9回で最も多いので（というか先頭以外は全て0）、0という数値を0という最小のビットに割り当てることができます。

そうするとX座標は「**残差0が9回出現する**」という情報のみを格納しておくことで、元の数値をデコードすることが可能です。

めちゃ圧縮されましたね！！

## PDALを使って実際に圧縮してみる

今回は冒頭で記載した神奈川県のデータを使っていこうと思います！
場所はどこでも良いので割愛します！
https://www.geospatial.jp/ckan/dataset/kanagawa-2022-pointcloud

今回は「数値地形図_オリジナルデータ（TXT形式）」を取得してきました。
このデータは以下のようになっており、カラム名がついておらず、どのデータがなんの情報なのかわかりませんね。

```text
1,-19000.00,-59937.56,2.29,2
2,-19000.00,-59859.98,13.80,1
3,-19000.00,-59729.79,12.51,1
4,-19000.00,-59476.85,12.62,1
5,-19000.00,-59357.90,2.22,2
...
```

今回は、「数値地形図_オリジナルデータ（LAS形式）」も一緒にダウンロードして比較することで、XYZの情報であるということがわかりました。
（この辺りが、CSVのめんどくさいところですね…！）

大容量なのでテキストエディターでもまともに編集できないため、ターミナルからvimを利用して以下のように編集しました。

```text
id,X,Y,Z,Classification
1,-19000.00,-59937.56,2.29,2
2,-19000.00,-59859.98,13.80,1
3,-19000.00,-59729.79,12.51,1
4,-19000.00,-59476.85,12.62,1
5,-19000.00,-59357.90,2.22,2
...
```

このファイルのサイズは「295.3MB」でした。

ではPDALを利用してLAS形式に変換しましょう！

```sh
pdal translate 09ld9534_org.txt 09ld9534_org.las --writer las --reader text
```

出力されたファイルは「291.4MB」でした。
CSVと比較してほとんど圧縮されていないですね。
前項で語った通り、LASのポイントフォーマットに従って、CSVには存在していなかった情報が付与されてしまっていることがわかります。

```sh
❯ pdal info -p 1 09ld9534_org.las
{
  "file_size": 291414183,
  "filename": "09ld9534_org.las",
  "now": "2024-10-29T13:27:11+0900",
  "pdal_version": "2.8.0 (git-version: Release)",
  "points":
  {
    "point":
    {
      "Blue": 0,
      "Classification": 1,
      "EdgeOfFlightLine": 0,
      "GpsTime": 0,
      "Green": 0,
      "Intensity": 0,
      "KeyPoint": 0,
      "NumberOfReturns": 1,
      "Overlap": 0,
      "PointId": 1,
      "PointSourceId": 0,
      "Red": 0,
      "ReturnNumber": 1,
      "ScanAngleRank": 0,
      "ScanChannel": 0,
      "ScanDirectionFlag": 0,
      "Synthetic": 0,
      "UserData": 0,
      "Withheld": 0,
      "X": -19000,
      "Y": -59859.98,
      "Z": 13.8
    }
  },
  "reader": "readers.las"
}
```

では、LAZ形式に変換していきましょう。

```sh
pdal translate 09ld9534_org.las 09ld9534_org.laz
```

ファイルサイズは「**29.2MB**」まで圧縮されました！
**1/10**になっていますね！！！

では最後にモートンオーダーでソートした上で圧縮して見ましょう！

```sh
pdal translate 09ld9534_org.las 09ld9534_org_mortonorder.laz --filter mortonorder
```

ファイルサイズは「24.1MB」になりました！
良いですね！

## 終わりに

ということで、今回はLAZファイルにフォーカスを当てて、その仕様を探っていきました。
特にLAZファイルの圧縮方式については間違っているところも多いと思いますので、「優しく」指摘していただけるとありがたいです。
LAZについては日本語での情報がとても少ないので、情報共有できると嬉しいです！！！