点群のファイル形式，特にCOPCについて理解する

# きっかけなど

点群データを使うにあたって，最初に目にした／手にしたファイル形式はLAS（拡張子`.las`）。でも，どうやらQGISではCOPC（拡張子`.copc.laz`）を使っている。ドラッグアンドドロップで読み込むと，勝手に変換される。
ファイルを見比べてみると，サイズが圧縮されている。でも，それだけではなくて，扱いやすいように工夫されているらしい。
使うにしても，中身が分からないのは気持ち悪いので，調べてみることにする。

なお，自分で実装する必要や意思はないので，概要がわかれば良いことにする。

![qgis1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3617623/50819821-a4d6-5f9e-4ec1-77cc552f9c8b.png)
画像：QGISで点群を読み込んだところ。赤枠は直ぐに表示されるが，表示するための処理(`.las`から`.copc.laz`への変換，画面下のプログレスバー。)が終わらないと，その中身は表示されない。

# 要約

- COPC形式を知るにはLAZ形式を知る必要があり，LAZ形式を知るにはLAS形式を知る必要がある。
- LAS形式は点群データの基本。ヘッダーおよびデータ本体からなる。
    - ヘッダーには規定のメタデータを指すPUBLIC HEADER BLOCKと認証ユーザーによって定義された情報VARIABLE LENGTH RECORDSがある。
    - データ本体はいくつかのフォーマットが用意されている。いまは5以前のフォーマットは使わない。
- LAZ形式はLAS形式の上位互換
    - 可逆（lossless），非プログレッシブ（non-progressive），ストリーミング（streaming），順序保持（order preserving），ランダムアクセス（random access）の圧縮
    - ファイルデータはLAS形式の10%くらいになる
    - LASと比べると，構造はそのまま，メタデータは圧縮されず，データ本体のみの圧縮。
- COPC形式はLAZ形式に，データを読み込むための有用なヘッダー情報を付与している
    - ユーザー定義のメタ情報VARIABLE LENGTH RECORDSに情報を付与することによって，LAZ形式のまま，八分木による効率的なデータ呼び出しが実現されている

# LAS形式（`.las`）

公式の仕様はココ。LIDAR（または他の）点群データレコードを含むオープンフォーマットであることが謳われている。
https://www.ogc.org/standard/las/

記事作成時の最新仕様はLAS 1.4

## 構成

Table 1に，構成が示されている。

- PUBLIC HEADER BLOCK
    - 点の数や境界などのデータ
- VARIABLE LENGTH RECORDS (VLR)
    - オプショナル。
    - 投影情報，メタデータ，波形パケット情報，ユーザー・アプリケーション・データなど。65,535バイトが上限。
- POINT DATA RECORDS
- EXTENDED VARIABLE LENGTH RECORDS (EVLR)
    - オプショナル。
    - VLRに書ききれない情報。

## ヘッダー（PUBLIC HEADER BLOCK）

次の情報が含まれる（Table 3参照）。
メタデータが含まれる一般的なヘッダー。

Item|Format|Size|Required
---|---|---|---
File Signature (“LASF”)|char[4]|4 bytes|*
File Source ID|unsigned short|2 bytes|*
Global Encoding|unsigned short|2 bytes|*
Project ID - GUID data 1|unsigned long|4 bytes
Project ID - GUID data 2|unsigned short|2 byte
Project ID - GUID data 3|unsigned short|2 byte
Project ID - GUID data 4|unsigned char[8]|8 bytes
Version Major|unsigned char|1 byte|*
Version Minor|unsigned char|1 byte|*
System Identifier|char[32]|32 bytes|*
Generating Software|char[32]|32 bytes|*
File Creation Day of Year|unsigned short|2 bytes|*
File Creation Year |unsigned short|2 bytes|*
Header Size|unsigned short|2 bytes|*
Offset to point data|unsigned long|4 bytes|*
Number of Variable Length Records|unsigned long|4 bytes|*
Point Data Record Format |unsigned char|1 byte|*
Point Data Record Length |unsigned short|2 bytes|*
Legacy Number of point records|unsigned long|4 bytes|*
Legacy Number of points by return|unsigned long|[5]|20 bytes|*
X scale factor|double|8 bytes|*
Y scale factor|double|8 bytes|*
Z scale factor|double|8 bytes|*
X offset|double|8 bytes|*
Y offset|double|8 bytes|*
Z offset|double|8 bytes|*
Max X|double|8 bytes|*
Min X|double|8 bytes|*
Max Y|double|8 bytes|*
Min Y|double|8 bytes|*
Max Z|double|8 bytes|*
Min Z|double|8 bytes|*
Start of Waveform Data Packet Record|unsigned long long|8 bytes|*|
Start of first Extended Variable Length Record|unsigned long long|8 bytes|*
Number of Extended Variable Length Records| unsigned long |4 bytes|*|
Number of point records|unsigned long long|8 bytes|*|
Number of points by return|unsigned long long [15]|120 bytes|*

ここで，`Point Data Record Format`は，LAS 1.4仕様では0～10の値をとり，点群データのフォーマットを表す。この値によって，データ形式の格納仕様が変わる（後述）。

また，座標のオフセットや最大値がヘッダーで指定されていること（つまり点群の位置情報が直接的に緯度経度（浮動小数点）で表されないこと）について，LASZip論文に重要な指摘がある。

>ある点のx,y,z座標を格納するのに，浮動小数点表現は整数表現よりも精度が高いという誤解がよくある。しかしそうではない。
>LiDARの点の座標はx-y平面に一様に分布しており，ほぼ同じ精度で点が取得されている。高い精度でポイントを保存する必要があるようなエリアは存在しない。
>しかし，浮動小数点フォーマットは，一様な分布の数値を保存するようには設計されていない。浮動小数点表現で数値を格納するということは，その数値の値によって精度が異なることを意味し，ゼロに近ければ近いほど精度は高くなる。点座標の保存に浮動小数点フォーマットを使うということは，ある点（原点）を中心にデータサンプルの間隔がどんどん正確になっていくことを意味する。
>座標を緯度経度で直接保存すると（つまり浮動小数点表現で保存すると），わずか0.5メートルの精度しか保てない。ただし，原点がバウンディングボックスの中央に来るように座標から一定のオフセットを引くと，原点付近のサンプルは信じられないほどの精度で保存される。
>LiDARポイントの座標を保存するのに適切なフォーマットは，適切にスケーリングされオフセットされた整数である。

## VLR

LASの管理団体によって認証されたユーザーが，任意のメタデータを付与できるスペースとして，VLRが準備されている。VLRは，次のヘッダーから始まる（Table 6参照）。

Item|Format|Size|Required
---|---|---|---
Reserved|unsigned short|2 bytes|
User ID|char[16]|16 bytes|*
Record ID|unsigned short|2 bytes|*
Record Length After Header|unsigned short|2 bytes|*
Description|char[32]|32 bytes

各フィールドについて，概要は次の通り
- `User ID`: LASの管理団体によって認証されたユーザーのID（たとえば`copc`）。任意のユーザーIDが登録できるわけではない。
- `Record ID`: 各ユーザーによって管理されている，0から65536までの値。
- `Record Length After Header`: このヘッダーに続くデータの長さ。

## 点群データ (POINT DATA RECORDS)

`Point Data Record Format`の値によって，格納データの仕様が異なる。
ざっくりいえば，Format 0～5が比較的古いフォーマット，6～10が比較的新しいフォーマット。

### Format 0～5

Format 0（基本となる形式）に含まれるフィールドは次の通り。（Table 7参照）

Item|Format|Size|Required
---|---|---|---
X|long|4 bytes|*
Y|long|4 bytes|*
Z|long|4 bytes|*
Intensity|unsigned short|2 bytes
Return Number |3 bits (bits 0 – 2)|3 bits|*
Number of Returns (given pulse) |3 bits (bits 3 – 5) |3 bits|*
Scan Direction Flag|1 bit (bit 6)|1 bit|*
Edge of Flight Line|1 bit (bit 7)|1 bit|*
Classification|unsigned char|1 byte|*
Scan Angle Rank (-90 to +90) – Left side|char|1 byte|*
User Data|unsigned char|1 byte|
Point Source ID|unsigned short|2 bytes|*

- Format 1では，Format 0に，GPS Time (Double, 8 bytes)が追加される（Table 10参照）。
    - Format 4では，Format 1に，Wave packetに関するフィールド（詳細略）が追加される（Table 13参照）。
- Format 2では，Format 0に，Red, Green, Blue (unsigned short, 各2 bytes)が追加される（Table 11参照）。
- Format 3では，Format 0に，GPS Time (Double, 8 bytes), Red, Green, Blue (unsigned short, 各2 bytes)が追加される（Table 12参照）。
    - Format 5では，Format 3に，Wave packetに関するフィールド（詳細略）が追加される（Table 13参照）。


### Format 6～10

Format 6（基本となる形式）に含まれるフィールドは次の通り（Table 15参照）。
Format 0と比較すると，Return NumberおよびNumber of Returnsのビット数が1大きくなっている等の改善が行われている。

Item|Format|Size|Required
---|---|---|---
X|long|4 bytes|*
Y|long|4 bytes|*
Z|long|4 bytes|*
Intensity|unsigned short|2 bytes
Return Number |4 bits (bits 0 – 3)|4 bits|*
Number of Returns (given pulse) |4 bits (bits 4 – 7) |4 bits|*
Classification Flags|4 bits (bit 0 - 3)|4 bits|
Scan Channel|2 bits (bit 4 - 5)|2 bits|*
Scan Direction Flag|1 bit (bit 6)|1 bit|*
Edge of Flight Line|1 bit (bit 7)|1 bit|*
Classification|unsigned char|1 byte|*
Scan Angle |short|1 bytes|*
Point Source ID|unsigned short|2 bytes|*
GPS Time|double|8 bytes|*

- Format 7では，Format 6に，Red, Green, Blue (unsigned short, 各2 bytes)が追加される（Table 18参照）。
    - Format 8では，Format 7に，Near InfraRed channel (unsigned short, 2 bytes)が追加される（Table 19参照）。
    - Format 10では，Format 7に，Wave packetに関するフィールド（詳細略）が追加される（Table 21参照）。
- Format 9では，Format 6に，Wave packetに関するフィールド（詳細略）が追加される（Table 20参照）。

## EVLRs

次の情報が含まれる。Table 22参照。

Item|Format|Size|Required
---|---|---|---
Reserved|unsigned short|2 bytes|
User ID|char[16]|16 bytes|*
Record ID|unsigned short|2 bytes|*
Record Length After Header|unsigned short|8 bytes|*
Description|char[32]|32 bytes

# LAS Zip形式（`.laz`）

`.laz`はlasのzip圧縮を指しているらしい。公式はココ。

https://laszip.org/

元となった論文は上記リンク先にも紹介されているが，リンク切れ。おそらくコレのこと。

https://www.cs.unc.edu/~isenburg/lastools/download/laszip.pdf

（圧縮技術に詳しくないのできわめてラフなまとめだが，）LASzipコンプレッサーの特徴は次の通り。

- 可逆圧縮（lossless）
- 非プログレッシブ（non-progressive）
    - プログレッシブ形式では，解凍が完了する前であっても，低解像度バージョンのデータが表示される
- ストリーミング（streaming）
    - 点群データの圧縮時に，データの一部を読み込んで圧縮ファイルを出力できる。解凍時には，データの一部を読み込んで解凍ファイルを出力できる。
- 順序保持（order preserving）
- ランダムアクセス（random access）

LASzipはヘッダ（PUBLIC HEADER BLOCK）やVLRを圧縮せず，それらは単にLASからLAZにそのままコピーされる。点群データ（POINT DATA RECORDS）の圧縮については，フィールド毎に適切な手法が選択される。
この様にして，ファイルサイズを元ファイルの7～10%にしつつ，圧縮ファイルの全てを解凍する必要なしに，点群をアプリケーションに読み込むことができる。

# COPC形式（`.copc`）

ようやくこれにたどり着いた。公式はコチラ。

https://copc.io/

![copc-vlr-chunk-table-illustration.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3617623/cf74c257-8ed0-4211-a640-186b40ab2375.png)


ただ，introductionを見ると，EPT形式から派生した形式の様なので，まずこれを確認する。


## EPT形式

公式はココ。

https://entwine.io/en/latest/entwine-point-tile.html

点群データのためのシンプルで柔軟な八分木(octree)ベースのストレージフォーマットらしい。
八分木になっているのは，3次元空間データだから。タテ2分割×ヨコ2分割×タカサ2分割で8分割になる。

次のようなファイル構成になっている。
```
├── ept.json
├── ept-data
│   └── 0-0-0-0.laz
├── ept-hierarchy
│   └── 0-0-0-0.json
└── ept-sources
    ├── list.json
    └── 0.json
```

### `ept.json`

ヘッダーに相当し，次のメタデータが含まれる。なお，COPC形式では利用されない。

- bounds
- boundsConforming
- dataType
- hierarchyType
- points
- schema
- span
- srs
- version

### `ept-data`

タイルマップに似た形式で点群ファイルが保持される。
タイルマップでは，`Zoom-X-Y`という形式だったが，EPT形式では`Depth-X-Y-Z`という形式になる。また，タイルマップのデータが排他的なのに対し，EPT形式では相加的なデータである。つまり，ルートデータ（`0-0-0-0`）から表示したい細部データまで，データが加算されていく。

データ`Depth-X-Y-Z`において，ルートボリュームは常に`0-0-0-0`で表され，空間解像度は`span`に基づき定まる。
深度Depthの各ノードには，それぞれ8つの子ノードが存在し，深度が1増えるごとに空間が2分割されていく。

### `ept-hierarchy`

`ept-hierarchy`には，`ept-data`にどの様なデータファイルが存在し，どの様な子構造を持ち，それぞれに何個の点が含まれているかが`json`ファイルで保存される。

### `ept-sources`

LAZデータについてさらにメタ情報がある際に，`ept-sources`に保存される。


## COPC形式の概略

まず，EPTと比較した際に，次の差異が挙げられている
- `ept.json`がない。LASファイルのヘッダーとVLRが利用される。
- `ept-sources.json`がサポートされない。
- LAZ点形式のみをサポートし、バイナリ点配置をサポートしない。
- 点データだけをLAZとして保存する。（EPTは，LASヘッダーとVLRを含む完全なLAZファイルが使用される。）

肝心なのは，LAZファイルを含む複数のファイルで構成されていたEPT形式と比較して，COPC形式では必要となる全てのメタ情報をLAZ本体に入れてしまったということ。
メタ情報の格納には，VLR領域が用いられ，`infoVLR`および`hierarchyVLR`が利用される。

COPC形式が利用可能なアプリケーション（QGISなど）でファイルを開くと，VLR領域の情報が利用され，八分木による効率的なデータ管理ができるようになる。

# おわりに

平面位置情報ならR-treeを使ったflatgeobuf，3次元情報ならoctreeを使ったcopcが，これからの標準の様ですね。学ぶことが多くて大変。このメモが未来の自分やその他の人々の助けになれば。

なお，上記には不正確な内容が多く含まれると思われるので，信用せず，利用目的に応じて原典に当たってください。