現場の点群をUnityに入れると、なぜ位置がずれるのか

建設会社でSEをやっています。
社内システムの面倒を見つつ、現場のDX化を進めています。
（70mのLANケーブル持って現場入ったり、重機の配線いじったりもしています）

過去に平面と点群の数学の話を2本書きました

今回はその続きというより、その手前の話です。
数学をやる以前に、そもそもデータが壊れていたら意味がないという話です。

そして今回、その対策を入れたLASインポータをOSSとして公開しました。人生初のOSS公開です。

MITで商用利用可能です。測量座標をどれだけ引き寄せたか（アンカー値）と座標系をアセットに保存するので、Unityで扱った結果を元の測量座標に戻せます。
ここが既存インポータとの一番の違いです。

以下はVirtualShizuoka様よりいただいた点群データをUnityに反映したものです。
こんな感じになります。

この記事で扱うこと

なぜずれるのか
float32の有効桁と、測量座標の大きさの相性の悪さ
どう直したか
点群の中心をアンカーにして原点へ引き寄せる
作ってみてわかったこと
LOD設計、日本の公開点群の座標系事情、LAS仕様の罠

使い方

Unityのメニューから Window > Package Manager を開く
左上の + から Install package from git URL... を選び、以下を貼り付ける

https://github.com/KazukiHarada-Seikyo/UnityLasImporter.git

.las ファイルを Assets フォルダにドラッグ＆ドロップする

インポート時に自動で4段階のLODメッシュに変換され、シーンに置けば距離に応じて切り替わります。

なぜ現場の点群をUnityに入れると、位置がずれるのか

現場でスキャンした点群をUnityに入れる。
ありそうな作業ですよね。私もやろうとしたことあります。

ところが、こういう経験はないでしょうか。

カメラを動かすと点群がプルプル震える
平面フィッティングの結果が妙にガタつく
ミリ単位で測ったはずなのに、Unity上ではセンチ単位でしか合わない

「点群だからこんなものか」で片付けてしまいがちですが、
割と簡単な話で対策ができてしまいます。

float32は「7桁しか覚えられない電卓」

Unityの座標は float32 です。float32 は数字を約7桁しか保持できません。

大事なのは、小数点の位置は関係ないということです。全体で7桁。

平面直角座標系の座標値で、ひたすらに遠い点の座標は10万m台になることがあります。

例えば

165375.05

これを7桁の電卓に入れると、こうなります。

1 6 5 3 7 5 . 0   ← ここまでで7桁
                5 ← 入らない

整数部分で6桁を使い切ってしまい、小数側に1桁しか残りません。

10進数で7桁だけ入る箱を考えるとわかりやすいです。
「刻み」とは、一番右の桁が小数点第何位の精度で表せるかということです。

4500.000の整数部を増やしてみましょう。

数	整数部の桁数	刻み
4500.000	4	0.001
45000.00	5	0.01
450000.0	6	0.1

同じ7桁でも、整数部が増えるほど小数部（精度）が削られるのがわかりますね。
Unityでは7桁までしか保存できないということから、

1mm精度で測ったデータが、Unityに入れた瞬間にcm精度になってしまう。

ということになります。

2進数なので、刻みが半端な数になる

2進数でもう少しコンピュータサイエンスらしく見てみましょう。

コンピュータは10進数ではなく2進数で桁を数えます。なので刻みは 0.01 や 0.001 のような綺麗な数ではなく、1を2で割り続けた数になります。

1  0.5  0.25  0.125  0.0625  0.03125  0.015625 ...

2進の桁を、整数と小数で取り合う

さっき「7桁」と書きましたが、正確には 2進で24桁 です。
2²⁴ = 16,777,216 で、10進に直すと約7桁になる、という話です。

そしてこの24桁を、整数部分と小数部分で取り合います。

165,000 m の場合

165,000 は 2¹⁷ = 131,072 と 2¹⁸ = 262,144 の間にあります。
つまり整数部分に 18桁 必要です。

24 − 18 = 6桁

小数に残るのは6桁だけ。2進の小数の桁は上から

1/2、1/4、1/8、1/16、1/32、1/64 ...

一番細かい刻みが 1/64 = 0.015625 m = 15.625 mm になります。

4,500 m の場合

2¹² = 4,096 と 2¹³ = 8,192 の間なので、整数部分は 13桁。

24 − 13 = 11桁

小数に11桁残るので、一番細かい刻みは 1/2048 = 0.00048828125 m = 0.488 mm。

まとめると

座標	整数部に必要な桁	小数に残る桁	刻み
4,500 m	13	11	0.488 mm
165,000 m	18	6	15.625 mm

整数部の差は5桁。2進は1桁で2倍なので、刻みの差は 2⁵ = 32倍です。

刻みが 15.625 mm という半端な数になるのは、2で割り続けた数しか使えないからです。
0.01 や 0.001 のような綺麗な数は、2進数では作れません。

「1mmを保てるのは16,777mまで」と私が勘違いしていた話

どういうことかというと、

float32 は2進で24桁。2²⁴ = 16,777,216 です。
だから「1mm刻みを1677万個ならべられる」→「16,777 m まで1mm精度を保てる」。

そう思っていました。間違いでした。

小数側に何桁必要か？から考える

1mm を保つには、小数側の一番細かい刻みが 1mm 以下でなければいけません。
2進数の小数の桁を、上から順に並べてみます。

小数の桁	大きさ
1桁目	1/2 = 500 mm
2桁目	1/4 = 250 mm
...	...
9桁目	1/512 = 1.95 mm
10桁目	1/1024 = 0.98 mm ← ここで初めて1mmを切る

小数側に最低 10桁 必要だとわかります。

24桁のうち10桁を小数に取られるので、整数に使えるのは 14桁 だけ。

2進の14桁で表せる最大の数は 2¹⁴ = 16,384 です。

16,384 で何が起きるか

16,384 になった瞬間、整数部分に 15桁 必要となります。

すると小数に残るのは 24 − 15 = 9桁。つまり 1/512 = 1.95 mm精度 となってしまうわけです。

なぜ 16,777 ではないのか

グラフにすると一目でわかります。実際の刻みは階段です。
2のべき乗のところで、2倍ずつ跳ね上がる。

一方、私が思い描いていた 16,777 という数字は、点線のようになだらかに落ちることを
前提にしていました。「距離が2倍になれば刻みも2倍」を、連続的に当てはめた数字です。

でも実際は、桁は整数個ずつしか動きません。
「整数に 14.4 桁」という配分は存在しません。24マスの箱を整数と小数で分け合っていて、
棚は1マス単位でしか動かせないからです。

だからコンピュータ上では境目は必ず2のべき乗になります。16,777 は階段のどの角にも当たらないため、刻みの基準になる数ではないのです。

実際に見てみる

理論だけでは信じられないので、テスト用のLASファイルを自分で書き出しました。

座標を 165,375 m 付近に置く
50mm間隔の格子を作る
LAS側のスケールは 0.001（mm単位で保持）

これを読み込んで、隣り合う点の間隔を測ります。

何が起きているか

座標 165,000m 付近で float が乗れる目盛りは 15.625mm 刻みでしたね。

ここに 50mm 間隔で点を並べようとすると、どうなるか。

50mm は目盛り 3.2 個ぶんです。3個でも4個でもない。
なので点は、一番近い目盛りに吸着します。

行	本来の位置	必要な目盛り数	実際に乗る目盛り	実際の位置	前行との差
0	0 mm	0	0	0 mm	—
1	50 mm	3.2	3	46.875 mm	46.875
2	100 mm	6.4	6	93.750 mm	46.875
3	150 mm	9.6	10	156.250 mm	62.500
4	200 mm	12.8	13	203.125 mm	46.875
5	250 mm	16.0	16	250.000 mm	46.875

3行目を見てください。9.6 は 10 に切り上がります。
それまで足りない側にズレていたのが、ここで初めて多い側に転びます。
だからこの1行だけ、目盛り4個ぶん = 62.5mm 進みます。

そして 5行目で 16.0 とぴったり整数に戻る。
だから5行で1周期になります。

隣の点との間隔
正しい値	50.000 mm
そのまま float 化	46.875 / 46.875 / 62.500 / 46.875 / 46.875 の繰り返し

X軸では違う壊れ方をする

同じファイルのX軸は座標 4,500m 付近で、目盛りは 0.488mm 刻みです。

50mm は 102.4 個ぶん。以下の表のようになります。

隣の点との間隔
正しい値	50.000 mm
そのまま float 化	49.805 / 50.293 / 49.805 / 50.293 / 49.805 の繰り返し

こちらも 102.4 × 5 = 512 で整数に戻るので、やはり5行で1周期。

周期は 50mm × 5 = 250mm なので、30m四方なら 30m ÷ 250mm = 120本の縞が立つ計算になります。
実際に描画したら、細かい縞が見えました。

対策：点群の中心をアンカーにして、原点へ引き寄せる

やることは単純です。

点群の範囲の中心をアンカー（基準点）に決めて、全点からその値を引く。

こうすると原点からの最大距離が、点群の半径まで縮みます。

	座標	誤差
アンカーなし・Y軸	165,000 m	±7.8 mm
アンカーなし・X軸	4,500 m	±0.24 mm
アンカーあり	15 m	±0.0005 mm

整数部分が席を空けたから、小数がその席に座れた、という話です。

さきほどのテストファイルでも、中心を引いてから float 化すると、隣の点との間隔は全行 50.000 mm ぴったりに戻りました。

重要なのは、引いた量を捨てないことです。多くのインポータはこれを黙ってやるか、まったくやりません。どちらの場合も、測量座標やロボット座標に戻す方法が残りません。今回作ったものは、アンカーの値と座標系をアセットに保存しています。

外した予想の話：重いのは描画だと思っていた

ここからは、作りながらわかったことの話です。

点が多いから描画が重い、だからLODが必要と考えましたが、
実測すると違いました。Core Ultra 5（内蔵GPU）で計測した結果です。

なお表の「マス目」は、点を間引くときの格子の間隔です。マス目2mなら、2m四方に代表点1つまで残します。

マス目	点数	FPS
2 m	約22万	180
1 m	約86万	145
0.5 m	345万	99

点数が15倍違うのに、フレーム時間は1.8倍しか違いません。
しかも 2m のときは1タイルでも3タイルでも 180 FPS のまま張り付いていました。点数によらずFPSが一定ということは、描画以外の固定コストが支配的だということです。

描画は全然ボトルネックではありませんでした。

本当の制約はこちらです。

	実測	100タイルなら
読み込み時間	3枚（3456万点）で6.5秒	約3.5分
ディスク	1枚300MB	30GB
描画	49万点で180FPS	余裕

遅いのは描画ではなく読み込みで、しかもPlayするたびに待つことになる。

だから設計をこう変えました。実行時にLODを計算するのではなく、Unityのインポート時に複数の粗さを作ってアセットに保存する。

Unityのインポート機構は変換結果をキャッシュするので、6.5秒かかるのは最初の1回だけです。そして読み込みを速くする仕組みと、LODの仕組みが、同じ一つの作業になりました。

日本の公開点群データは、座標系を名乗らない

LASには座標系情報を入れる領域（VLR）があります。
EPSGコードを読めば、軸の並びまでわかるはずでした。

日本の平面直角座標系（EPSG 6669-6687 ほか）は、定義上 Northing, Easting の順です。一方 LAS の仕様と世界の大半の投影座標系は Easting, Northing の順。ここが食い違うと、縦横が入れ替わって鏡写しのような表示になります。

そこで実データを開いてみたら——

手元の2ファイル、どちらもVLRが0個でした。

愛知の航空レーザ測量データ
VIRTUAL SHIZUOKA の車載計測データ

どちらも配布ページには「平面直角座標系第○系」と明記されているのに、ファイル自体は何も名乗っていません。 人間向けの説明にはあって、機械向けのデータには無い。

なので、こういう設計にしました。

EPSGが読めたら、軸順の推奨値を出す
ただし手動で上書きできる
読めなかったら、推測せず「不明」と表示して既定のまま

自動判定はあくまで参考程度です。
「たぶんこうです、違ったら直してください」が現場では正しい態度だと思っています。

ちなみに静岡のデータについては、座標値から緯度経度を逆算して伊豆半島東岸に落ちることと、図郭（点群を配布するときの地図の区画単位）の縦横比が3:4であることの両方から、X=東・Y=北（LAS標準）だと推測できました。
座標系が書いてなくても、座標値がヒントを教えてくれることがあります。

LODの切り替えをどこで決めるか

LODは「画面上で点の間隔が何ピクセルに見えるか」で切り替えています。

距離 $d$、マス目 $g$、垂直画角 $\theta$、画面高 $h$ px のとき、画面上の点間隔は

s = \frac{g \cdot h}{2d \tan(\theta / 2)} \ \text{px}

1ピクセルを下回るほど細かくしても画面には出ません。必要な細かさは距離に反比例するだけなので、遠景をいくら広げても点数はあまり増えません。

実装で1つ踏んだところがあります。最初は毎フレーム Camera.main を見て判定していたのですが、編集モードではシーンビューのカメラが反映されず、距離が固定されたままでした。

正しくは、それぞれのカメラが描き始める直前に、そのカメラ用の段階を選ぶ。

RenderPipelineManager.beginCameraRendering += OnBeginCameraRendering; // URP/HDRP
Camera.onPreCull += Evaluate;                                        // Built-in

LODは「毎フレーム1回決めるもの」ではなく「カメラごとに決めるもの」でした。

タイルをまたぐときの工夫

LASは重いので、図郭単位で分割配布されています。複数枚を並べると、2つ問題が起きます。

アンカーが揃わない

タイルごとに自分の中心を引くと、隣り合うタイルが同じ場所に重なります。
↓
頂点座標はタイル中心からの相対位置（小さい数）、タイル同士の位置関係は GameObject の Transform（大きい数）に持たせました。

マス目が揃わない

インポータは1ファイルずつしか見られないので、隣に何があるか知りません。
↓
マス目の原点を floor(座標 / マス目) × マス目に丸めました。こうすると隣を知らなくても格子線が全部同じ位置に来ます。 継ぎ目に二重や抜けが出ません。

LASの仕様で引っかかったところ

せっかくなので、実装で踏んだ細かい罠も置いておきます。

LAS 1.4 で点数の場所が変わる

8バイトの別フィールドに入り、従来の位置には 0 が書かれます。古い読み方だと「点が0個のファイル」に見えて、何も表示されずに終わります。

色の位置がフォーマットで動く

XYZ は全形式で先頭12バイトですが、色は形式2なら20バイト目、形式3なら28バイト目、形式7なら30バイト目。間にGPS時刻が挟まるかどうかで8バイトずれます。

テスト用に、位置から色を作ったLASを書き出して確認しました。
赤が横方向、緑が奥方向。形式2と形式3で同じ絵が出れば、
位置の切り替えが効いている証拠になります。

1点のサイズが仕様より大きいことがある

メーカー独自の情報が後ろに付く場合です。1点ぶん読み進める幅（歩幅）は必ずヘッダの値を使い、自分で計算してはいけません。

色が8bitか16bitか、ファイルによって違う

仕様では1チャンネル16bit（0〜65535）ですが、実際には 0〜255 を16bitの箱に入れているファイルが多くあります。判定を誤ると真っ黒か真っ白になります。

形式0〜5では分類バイトの上位3ビットが別の意味のフラグ

潰さずに読むと分類番号が化けます。

できなかったこと

LAZ（圧縮形式）に対応していません。
技術ではなくライセンスの問題です。C#からLAZを読む現実的な実装は LGPL ですが、UnityのIL2CPPビルドは静的リンクになるため、LGPLの条件を満たすのが困難です。MITで配布する以上、本体には含められません。LAZを展開してから使ってください。

標高で色を塗ると、タイル間で色がつながりません。
インポートが1ファイル単位なので、そのタイル自身の高さ範囲で正規化しています。RGBを持つデータでは起きません。

実行時（エディタ外）の読み込みには対応していません。
エディタのインポート経由のみです。

1000万点超のストリーミングはやっていません。
今回の設計だと全点をメモリに載せて間引くので、桁が上がると破綻します。
Potree系のオクツリー方式に頼るしかありません。

こんな使い方もできそう

作りながら思ったのですが、座標を揃えた点群がUnityに載るというのは、それだけで応用の幅があります。

3Dハザードマップ
他で計算した結果を、座標を合わせて点群の上に重ねて見せるのがUnityの得意な領域だと思っています。点群データの上に他で計算したハザードマップの情報を描画することで3Dハザードマップを作製することができます。
Unityなので、例えばカメラを人の目の高さに置いて、自分の家の玄関から見て、水位2mでどこまで沈むかのような見せ方ができます。
ドローンの飛行経路
障害物を避けて通る経路を引くこともできるかと思います。
設計モデルと違って、点群には電線やクレーン、仮設足場も入っています。
図面に無いものにぶつからない経路が引けるはずです。
地形の差分
同じ場所を2つの時期でスキャンして重ねることもできます。
造成の前後、地滑りの前後、災害の前後。座標が揃っているので引き算するだけで、どこがどれだけ動いたかを立体で視覚的に把握することができます。

おわりに

LASの仕様書を1バイトずつ読むところから始めて、公開できるパッケージになるまで、丸二週間かかりました。

MITライセンスです。Issue も Pull Request も歓迎します。
基礎的な質問から技術的な相談まで鋭意対応いたします。

感想・いいねお待ちしております。
ではでは皆さま、今日も一日、ご安全に！