[Unity] スマホ1台で「実寸の点群データ」をつくる
点群データを皆様ご存知でしょうか?
先月、GoogleCloudNext Tokyoに行きまして、GoogleMapに関する講演を聴講していました。そこで点群データの紹介があったのですが、かなり詳しく説明されていたのでそろそろ時代がくるんじゃないかと感じています。
点群とは、現実の物や場所を3次元座標を持った点の集まりとして写し取ったデータのことを言います。スマホで撮った写真からこれを起こす技術(フォトグラメトリ)は昔からあるのですが、写真だけでは実際の大きさが決まらないという根本的な壁があります。
今回は、その壁をスマホの加速度計で越え、実寸と地図座標の付いた点群を作ります。点群を作成する理論から書いていますので、是非参考にしていってください。
これまで3本、見る側の話を書きました。
今回は作る側です。
先に結論を置いておきます。
| 実測対象 | 実測精度 |
|---|---|
| 形と大きさ(相対) | 1%以内。173cmの巻尺を三角測量して172.6cm |
| 置き場所と向き(絶対) | 水平2.4m / 方位4.7度 |
| 出来上がった点群 | 6,221,638点(密ステレオまで通した場合) |
形を出すだけならよいのですが、やはり画像から向きや座標を決め切るのは難しいです。
専門機器無しで精度がそれなりの点群を取得する際に使用できるものになります。
OSSとして公開していますので、Unityをお持ちの方は是非お試しください。
点群を撮っている最中の画面は以下のようになります。
せっかくUnityで作ってるのでSFチックにしました。Insta映えするかも?
読み方ガイド
盆休みで作成した長い記事です。理論だけ知りたい方は 1〜2 章、手を動かしたい方は 3〜6 章と 13 章、精度の話が見たい方は 8〜10 章へ。11章は同じことをする人へのプレゼントです。
1. 写真からは「大きさ」が決まらない
写真から立体を起こす技術(フォトグラメトリ)には、困ったところがあります。
写真をどれだけ集めても、写ったものの大きさは決まりません。
割と正確な形と長さの比は出せます。
しかし、この壁は何メートルですか?には答えられません。
なので普通、既知の長さのものを一緒に写します。定規、スケールバー、コンベックス、対空標識。現場で三脚立ててるアレ...(トータルステーションといいます)などなど。
本プロジェクトでは写真から大きさまで出します。
1-1. カメラの式
いちばん単純なモデルは針穴写真機です。箱に小さな穴を開けて、奥にスクリーンを置く。ある一点から出た光のうち穴を通った1本だけがその穴を通り、スクリーンに届くので、点が元の点として写ります。
穴を原点にして、カメラの向きを Z 軸、右を X 軸、上を Y 軸とします。スクリーンは穴から距離 $f$の位置にあるとします。
空間の点 $(X, Y, Z)$ から出た光がスクリーンに当たる位置 $(x, y)$ は、相似な三角形で求まります。
$$
\frac{x}{f} = \frac{X}{Z}, \qquad \frac{y}{f} = \frac{Y}{Z}
$$
つまり
$$
x = f \frac{X}{Z}, \qquad y = f \frac{Y}{Z}
$$
奥行き $Z$ で割っている。
つまり遠ければ遠いものほどそれが占めるスクリーンの割合は小さくなりますよね。
画像は画素の集合体ですから、スクリーン(センサー)上の位置 $(x, y)\text{ [mm]}$ を、デジタル画像の左上を基準にした「ピクセル座標」$(u, v)\text{ [px]}$ に変換します。
1mmあたりのピクセル数を $k_x\text{ [px/mm]}$ とすると
$$\text{ピクセルでの位置} = x \times k_x = \left(f \cdot k_x\right) \frac{X}{Z}$$ここで、焦点距離 $f\text{ [mm]}$ $\times$ 1mmあたりのピクセル数をひとまとめにしたものを、ピクセル単位の焦点距離 $f_x\text{ [px]}$ といいます。
Pixel 9aを例に出すと($f_x \approx 1388.6$)であるので、このカメラのレンズで撮影すると、距離 $Z$ に対する横幅 $X$ の比に約1388.6を掛けたピクセル数として写るということになります。カメラ固有の定数と捉えてください。
画像の中心(原点)をずらす数学の計算ではレンズの中心直下を $(0,0)$ としていますが、実際の画像データでは画像の左上が $(0,0)$ です。そのため、画像中心のピクセル座標 $(c_x, c_y)$ を足して位置を補正します。
なので、最終的な式は
$$u = f_x \frac{X}{Z} + c_x, \qquad v = f_y \frac{Y}{Z} + c_y$$
$\frac{X}{Z}$ :遠近感(奥行き $Z$ が深いほど小さくなる)
$f_x$ :カメラの倍率(焦点距離 $\times$ 画素密度)
$c_x$ :画像中心のオフセット(原点を左上へシフト)
1-2. 大きさが消えるのはなぜ ?
一度、空間の点を全て $k$ 倍に拡大して、カメラの位置も同じ $k$ 倍に動かしてみます。
$$
u = f_x \frac{kX}{kZ} + c_x = f_x \frac{X}{Z} + c_x
$$
$k$ が約分されて消えます。
10cmの模型を10cm離れて撮った写真と、10mの実物を10m離れて撮った写真は、画素として完全に同じになります。何枚撮っても、どの角度から撮っても変わりません。
情報が足りないのではなく、原理的に同じ像になってしまうのです。
当たり前の誰もが知る現象ですが、数式で説明される機会が少ないので少し新鮮かと思います。
1-3. 2枚あれば形は出る
大きさは決まりませんが、形は2枚あれば出ます。
同じ点を2つの別々の位置から撮ると、画像上で少しずれた位置に写ります。このずれを視差と呼びます。別のカメラが距離 $B$ だけ離れているとき、距離 $Z$ にある点の視差 $d$(画素)は
$$
d = f_x \frac{B}{Z} \quad \Longrightarrow \quad Z = f_x \frac{B}{d}
$$
視差が分かれば奥行きが出ます。
つまり、基線長$B$離した位置からの画像があれば形を出せるわけです。
どの分野も、既知の長さを1本どこかに入れています。
| 分野 | 既知の長さ |
|---|---|
| 測量 | 基線(別途きちんと測ったもの) |
| 写真測量 | コンベックス、対空標識 |
| ステレオカメラ | 2つのレンズの間隔(工場で固定) |
| LiDAR | 光速 × 往復時間 |
| 今回 | 加速度計の値を積分した移動距離 |
2. 加速度計からカメラの移動距離を計算する
2-1. 中身はばねと錘
加速度計の中身は、ばねに吊るされた小さな錘です
(実際はMEMSという微細加工の櫛歯らしいです)
錘が加速されるとばねが伸びます。伸びた量 $x$ から、運動方程式より、
$$ma = kx \implies a = \frac{k}{m} x$$
ばね定数 $k$ と質量 $m$ は、その部品が持っている物理量なので、工場で校正されています。
静止していても 9.8$m/s^2$ を示す
加速度計が測るのは、重力以外の力による加速度です。机に置いた加速度計は重力に逆らって机から支えられていますよね。この支える力が測られるので、静止時に上向き 9.8 を示します。自由落下させると 0 になります。
2-2. 2回積分すると距離
加速度を時間で積分して速度、もう1回積分して位置が出ます。高校物理ですね。
$$
v(t) = \int a, dt
$$
$$
x(t) = \int v, dt
$$
この $x$ が、先ほどの $B$ です。
$B$ がメートル単位で出るので、
$$Z = \frac{f_x \cdot B}{d}$$
もメートルで出ます。
つまりマーカーも基準尺も置かずにスマホだけで実寸が出るという理屈です。
2-3. なぜカメラと組ませるのか
加速度計は長時間の撮影に向きません。
2回積分するので、誤差が時間の2乗単位で膨らんでしまうからです。
一方でカメラは、形と向きは長時間正確に追えますが、大きさが分かりません。
| 得意 | 不得意 | |
|---|---|---|
| カメラ | 形・向きの測定 | 大きさの測定 |
| 加速度計 | 大きさの測定 | 長時間の測定 |
得意・不得意が嚙み合っていますね。
これらカメラ・加速度計を組み合わせた技術を VIO(Visual-Inertial Odometry)といいます。
Googleが開発したARCoreのモーショントラッキングはこれです。
このプロジェクトは、ARCoreが計算しているカメラ位置をそのまま貰ってきて、写真と一緒に保存することで、計算可能としています。
3. 全体の流れ
スマホ側は姿勢付きの連番JPEGを貯め、
PC側はその姿勢を使って点群を起こし、GPSによる地図座標を付けます。
4. 撮影
4-1. 動画にしなかった理由
動画で撮って、あとからコマに割ればいいのでは?と思いますよね。
私も思いました。
ただ、動画をコマ割りすると、
- フレーム間圧縮で特徴点がつぶされしまう。(動画は前のコマとの差分で容量を減らそうとするので、平坦な部分の情報が落ちてしまいます)
- 1コマずつ姿勢と対応させられない。
という欠点があります。
なのでJPEGを連番で保存し、1枚1枚にカメラの位置と向きの情報を添えています。
4-2. CSVの中身
frames.csv で、1行が1枚に対応しています。
| 列 | 中身 |
|---|---|
local_px/py/pz |
ARCoreの局所座標系でのカメラ位置($m$) |
local_qx/qy/qz/qw |
カメラの向き |
lat, lon, alt_ellipsoid
|
地球座標(緯度・経度・楕円体高) |
acc_h, acc_v, acc_yaw
|
ARCoreが申告する精度 |
fx, fy, cx, cy
|
内部パラメータ。フレームごとに動くので毎回控える |
screen_orientation |
端末の向き |
方針として、計測値だけにしました。
現地で計算処理をさせず、全部持って帰ってPCで考えます。
暑い中何回も計測しに行きたくないですもんね。
私は何回も後悔し、何回も計測にいきました...
4-3. 計測後に後悔したこと
-
XR Origin の Camera Y Offset
既定値1.1176がカメラのYに常時加算されます。立っている人の目の高さを想定した値だと思います。こちらは実測が欲しいので0にします。
しばらく気づきませんでした。罠です(迫真)。 -
解像度ARCore経由のフレームは処理用の解像度で、既定は 640×480 です。
選べたのは 640×480 / 1280×720 / 1920×1080のどれか。
5. 姿勢を渡す
ここではARCoreが出した姿勢を与えます。計算処理が不要で、実測値を取れるからです。
COLMAPから姿勢推定もできるのですが、やはり現場で計算処理をさせたくないので、今回はパスしました。
5-1. 反転を2回かけると回転になる
UnityとCOLMAPでは軸関係が違います。
| Unity / ARCore | COLMAP | |
|---|---|---|
| 世界座標 | 左手系 | 右手系 |
| カメラのY軸 | 上向き | 下向き |
| 姿勢の向き | カメラを世界に置く | 世界の点をカメラから見る |
左手系から右手系にするには軸をどれか反転します。
$$
S = \begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & -1
\end{pmatrix}
$$
カメラのY軸を上から下に変更するのも反転ですね。
$$
C = \begin{pmatrix}
1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & 1
\end{pmatrix}
$$
反転を1回かけたものは回転ではないですね。 鏡映しになってしまいます。
行列式だと $-1$ ですね。
2回かけると $+1$ に戻ります。つまり180度回転させたことになります。
$$
R_{c2w} = S \cdot R_{unity} \cdot C
$$
Unity、ARCore、COLMAPを使用しているとき、片方の反転だけ直そうとすると、鏡映しのデータが完成してしまうということです。
5-2. 縦持ちの90度
もうひとつ回転系で。
AR FoundationがUnityのカメラに与える姿勢は画面の向きに合わせたものです。
一方、画像として取り出せるのはセンサーそのままで横向きに端末を持った状態。
端末を縦に持って撮ると、この2つが90度ずれてしまいます。画像を回さずに保存している以上、姿勢の側で戻すしかありません。視線軸まわりの回転なので、回転行列の列を組み替えます。
| 起きた3次元点 | |
|---|---|
| 補正なし | 2,815点 |
| 補正あり | 30,023点 |
三角測量は、2台のカメラからあそこに点が見えるぞと伸ばした2本の視線の交点を探す作業です。
ところが、カメラの回転が90度ずれていると、視線がぐるっとねじれてしまいます。
中心付近の特徴点は回転軸の近くなので、視線のズレはほぼゼロですが、
画面の端視線の向きが約28度も明当違いな方向を向いてしまいます。
ねじれてしまった2本の視線は、空間上でねじれの位置になるか、カメラの真後ろで交差してしまいます。
COLMAPは、ずれや位置関係を厳しくチェックするため、交わらない視線からできた、でたらめな3D点をすべて捨てます。
6. 点群にする
姿勢を固定した場合と、COLMAPに自力推定させた場合を比べました。
| 姿勢を固定(ARCore) | 自力推定 | |
|---|---|---|
| 登録画像 | 286 / 286 | 245 / 286 |
| 3次元点 | 30,023 | 34,923 |
| 再投影誤差 | 1.80 px | 0.885 px |
| 平均トラック長 | 3.61 | 9.49 |
| 所要時間 | 4.6分 | 46.8分 |
姿勢を与えると全枚数が10倍速く通る代わりに、精度が自力推定に劣っています。
ARCoreの姿勢は向きは合っているが、写真から詰めた姿勢ほどの精度はない」ということになります。
トラック長は、1つの3次元点が何枚に写っていたかを示しています。
長いほど信頼できます。2枚しか写っていない点は誤対応でも成立してしまうので、位置が不安定になります。
6-1. ARのカメラ位置は固定するより、後から合わせる
結論から言うと、スマホ(ARCore)が出したカメラ位置(ARポーズ)は信じ切って固定してはいけないということが分かりました。
COLMAPに自力で位置を解かせ(mapper)、最後にARポーズを使ってスケールを補正する(model_aligner)手法に切り替えたところ、精度が劇的に改善しました。
| モデル | 復元できた画像数 | 点の数 | 平均トラック長 | ずれの小ささ |
|---|---|---|---|---|
| ARポーズ固定 | 2,732 | 702,673 | 2.95 | 1.71 px |
| mapper → model_aligner | 1,365 | 527,252 | 5.56 | 0.830 px |
同じ点を見続けられる回数(トラック長)が2倍近くになり、1点あたりのピクセル誤差は半分以下まで激減しました。
なぜARポーズ固定だと失敗するのか?
一見するとスマホが計算してくれた位置をそのまま使えば速くて正確そうです。
しかし、ARCoreとCOLMAPが要求する精度の次元が全く違っていました。
COLMAP(画像処理)の世界
3m先の物体を撮る場合、わずか1ピクセルのズレが「2.2mm」の違いになります。COLMAPはミリ単位の整合性を求めています。
ARCore(スマホVIO)の世界
手ブレや歩行のズレを抑えるには十分ですが、局所的にも「約4.0cm」の誤差を持っています。
つまり、COLMAPから見ると「ARCoreの指定位置は、18ピクセルもズレた大雑把な位置」だったのです。
このズレた位置を正解として固定してしまったため、COLMAPは「視線が交差しないでたらめな点だ」と判断し、対応点の実に98%をゴミ箱へ捨ててしまっていました。
| ポーズの残差 | |
|---|---|
| 局所150枚の中央値 | 4.0 cm |
| 133m全行程の中央値 | 25.2 cm |
| 同・平均 | 3.26 m |
焦点距離1,388・被写体距離3mだと、1画素が約2.2mm相当です。局所の中央値4.0cmでも、要求精度の約20倍粗い。
だから point_triangulator でポーズを固定すると、対応点の98%が棄却されます。速いけど精度が劣るどころではなく、そもそも足りていなかったわけですね。ARのポーズは点群の実寸を与える基準としては優秀ですが、MVSの入力としては粗すぎる。役割が違うということでした。
「解き直すと実寸が保たれないのでは」という心配は、model_aligner が解決します。ARのカメラ中心を参照位置として渡すと、任意スケールで解かれたモデルを元の実寸・向き・原点へ戻してくれます(3章の図のオレンジの矢印がこれです)。
colmap model_aligner \
--input_path <mapperの出力> \
--output_path <出力> \
--ref_images_path ref_centers.txt \
--ref_is_gps 0 \
--alignment_type custom \
--alignment_max_error 0.10
ref_centers.txt はARポーズのカメラ中心の一覧です。これを通さないと mapper の出力は任意のスケール・向き・原点になり、平面直角座標に変換できません。
6-2. 密ステレオまで通す
ここまでは疎な点群です。特徴点が立った場所にしか点がありません。面として見るには密ステレオ(MVS)が要ります。1枚1枚について深度マップを作り、最後に全部を融合します。
| 項目 | 値 |
|---|---|
| 深度マップを計算した枚数 | 683(1,365枚の1/2) |
| 解像度(長辺) | 1920 px |
| 幾何整合 | 有効 |
| 幾何整合後の有効画素率 | 38.94% |
| 出力 | 6,221,638点 |
疎な点群(12万点)の51倍です。
処理は素直に重いです。RTX 5070 での実績。
| 工程 | 時間 |
|---|---|
| mapper(1,366枚) | 1時間36分 |
| 測光パス(683枚・1920px) | 1時間45分 |
| 幾何整合パス(683枚・1920px) | 約2時間20分 |
| 融合 | 8.9分 |
中間ファイル(深度+法線)は1枚あたり約66MB、683枚で約43GB。ピーク使用量であって保管量ではなく、融合後は捨てられます。
(1) 間引くならモデルごと間引く
patch-match.cfg の参照画像リストだけを間引いて、ソースを __auto__ で全フレームから選ばせる構成は、geom_consistency=true で破綻します。
測光パスはソース「画像」だけで足りるので完走します。ところが幾何整合パスはソース側の「深度マップ」を読むので、無いものを探して落ちる。測光パスに1時間45分費やした後で segfault します。
colmap image_deleter --image_names_path でモデル自体を絞ってください。image_undistorter --image_list_path は書き出す画像を絞るだけでモデルは全画像を保持するので、この用途には使えません。
(2) stereo_fusion の use_cache は既定 0
このときワークスペース全体(深度11GB + 法線32GB)をRAMに読み込もうとします。cache_size は無視されます。WSL2 だとプロセスではなく VMごと落ちるので、dmesg にも stderr にも何も残りません。ログが Loading workspace data... の直後で途切れるだけ。
--StereoFusion.use_cache 1 --StereoFusion.cache_size 6
これを明示すると、683枚・1920px が8.9分で終わりました。
7. 座標を付ける
点群はまだ「起動地点が原点」の局所座標にいます。これを世界地図の座標に移します。
7-1. 楕円体高と標高は違う
GPSが返す高さは楕円体からの高さ(楕円体高)です。一方、地図や標識の「標高」はジオイドという面からの高さ。ジオイドは平均海面を陸まで延長した重力の等しい面で、地下の密度のむらでデコボコしています。
$$
\text{楕円体高} = \text{標高} + \text{ジオイド高}
$$
日本付近のジオイド高は 30〜45m あります。
今回の撮影地は 37.75m でした(国土地理院 GSIGEO2011)。最初は概算の37.0mを定数で使っていて、Zが0.75mずれていました。場所ごとに確かめないと駄目ですね。
7-2. 平行移動で済ませられない理由
局所座標の原点を平面直角座標に変換して、あとは平行移動でいいのでは?と思いますよね。
子午線収差という問題があります。平面直角座標系の「北」と、その地点での本当の北は、原点から東西に離れるほどずれます。
名古屋付近で 0.215度。小さく見えますが、
$$
30,\mathrm{m} \times \tan(0.215°) \approx 11,\mathrm{cm}
$$
30m先で11cm。 無視できません。
なので実装では、点を1つずつ
局所座標 → EUN(東・上・北) → ECEF(地球中心の直交座標)
→ 緯度経度 → 平面直角座標
と厳密に通しています。数千万点でも耐えるよう、リストを作らず1点ずつ流し込む書き方にしました。
なお緯度経度⇔平面直角座標の変換(全19系)は単体で使えるようにしてあります。国土地理院の変換サービスと突き合わせて0.1mm差で一致しました。
Pythonで平面直角座標を扱いたい方はこれだけ持っていってください。
8. 精度をどう測ったか
8-1. 1つの数字にまとめない
この手法の誤差は、性質の違う2つに分かれます。
| 何が決まるか | 誤差源 | |
|---|---|---|
| 相対 | 点群の形と大きさ | 加速度計、カメラ、三角測量 |
| 絶対 | 地球上のどこに、どの向きで置くか | VPS(測位) |
相対はcm級、絶対はm級です。
8-2. 巻尺を撮る
実寸が本当に合っているかを確かめるには、外部の物差しが要ります。
黒い塀に巻尺を 173cm 伸ばして貼り、そのまわりを歩きながら撮りました。
- 各コマで巻尺(暗い壁の上の細長い明るい帯)を見つけ、両端の画素位置を求める
- ARCoreの局所姿勢だけを使って視線を作り、両端をそれぞれ三角測量する
- 2点間の距離を巻尺の読みと比べる
VPSや緯度経度を使わず、 加速度計とカメラで出した長さを巻き尺と突き合わせています。
端点は「明るさが半分に落ちる位置」で決めました。ぼけが左右対称なら物理的な端に一致して、コマが変わっても同じ場所を指すからです。
| 三角測量による長さ | 172.6 cm |
| 巻尺の読み | 173.0 cm |
| 差 | −0.4 cm(−0.22%) |
コマの選び方を変えると 170.7〜172.6cm に散りました。
ばらつき約2cmなので、結論は172 ± 1 cmで、実寸は約1%の精度で正しいということになります。
実際、0.1%級で合っているとは言えないと思います。 どの区分でも173cmより1〜2cm短く出ますが、端点の決め方が内側に寄るためとも説明できて、区別がつきませんでした。
8-3. 危なかった話
最初、巻尺を貼った直後のコマだけで計算しました。結果は 136cm。21%違います。
原因は基線長です。三角測量は、2つの視点がどれだけ離れているかで精度が決まります。奥行き $Z$ の誤差は、視差の誤差 $\sigma_d$ からこう伝わります。
$$
Z = f_x \frac{B}{d} \quad \Longrightarrow \quad \sigma_Z = \frac{Z^2 \sigma_d}{f_x B}
$$
$Z$ の2乗に比例し、基線長 $B$ に反比例する。 遠いものほど急に悪化して、視点が離れていないほど悪化します。
数字を入れます($Z$ = 1.7m、$f_x$ = 1388px、$\sigma_d$ = 10px)。
| 使ったコマ | 基線長 $B$ | 理論上の $\sigma_Z$ | 実際 |
|---|---|---|---|
| 貼った直後(立ち止まって撮影) | 0.115 m | 約 18 cm(10%) | 135.9 cm |
| 歩き回ったコマも含める | 0.723 m | 約 3 cm(1.7%) | 172.6 cm |
長さは奥行きに比例するので、奥行きが10%狂えば長さも10%狂います。
基準物を貼って、その場から撮っただけでは何も測れません。 今回は歩いたあとも巻尺が何度も画面に入っていたので救われました。目安は 基線長 ÷ 対象までの距離 ≧ 0.3。2m先のものなら60cm以上は横に動く必要があります。
現場で対空標識を置いたあと、その場で満足して帰ってはいけない、ということですね。
9. ARCore Geospatial は実際どう動いているか
前節の副産物として、公式には書かれていない挙動がいくつか見えました。
9-1. 新しい測位が入るのは一瞬だけ
acc_yaw(方位精度の自己申告)と、楕円体高から局所座標の上下動を引いた値を並べると、同じ時刻に同じ段差が出ます。
| 時刻 | 高さの下駄 | 方位精度 |
|---|---|---|
| 63〜167秒 | 43.539 m で一定 | 8.66 → 8.84度(じわじわ悪化) |
| 167.0〜167.7秒 | −0.503 m | 8.84 → 4.43度 |
| 168〜378秒 | 43.039 m で一定 | 4.43 → 4.85度(じわじわ悪化) |
読み方はこうです。
- 段差 = VPSが照合し直した瞬間。 ここで位置・方位・高さがまとめて付け替わる
- 段差と段差の間は推測航法。 新しい観測が入らないので、自己申告の精度が毎分0.113度のペースで膨らんでいく
314秒・133mの撮影で、段差は1回だけでした。地球に対する置き方の情報は、実質2回分しか入っていません。
前半と後半で当てはめ直すと、方位が3.17度、位置が1.65m違いました。これが1回の撮影のあいだに「点群全体が回って動いた量」です。自己申告の方位精度4.74度と同じ桁なので、申告値は信用してよさそうです。
実用上の意味: 現地に着いたらすぐ撮らず、歩いて周囲を映してから読む。経験則として言われていることですが、数字で裏が取れました。
9-2. 高さは毎回測っていない
楕円体高から局所座標の上下動を引くと、フレーム間の変化が中央値 1.10mm しかない定数が残りました。
楕円体高は、局所座標の鉛直に定数を足したものです。 その定数はVPSが照合し直したときだけ更新されます。
- 垂直精度の値は「毎フレームの高さの精度」ではなく「足している定数の精度」
- 短時間の高低差は局所座標の精度(cm級)で信用してよい
- 絶対標高は定数まかせなので、別系統で与える必要がある
9-3. その定数はどこから来たのか
「高さは測位ではなく地形の標高データから引いているのでは」と疑っていました。国土地理院の数値標高モデル(DEM、5mメッシュ)と突き合わせます。
DEM標高 中央値 3.15 m
楕円体高 中央値 43.5 m
その差 40.36 m
ジオイド高 37.75 m(GSIGEO2011)
→ 端末の地上高 2.61 m
端末を1.4mの高さで持っていたとすると、ARCoreの高さは +1.21m 高い。
地形データから引いているなら、DEMとはcm級で噛み合うはずです。噛み合わない。「地形データ由来」という疑いは支持されませんでした。
9-4. 仮定を外して測り直す
上の計算は「端末を1.4mで持っていた」という仮定に乗っています。
密な点群ができたので、仮定なしで測れるようになりました。点群には地面そのものが写っているからです。
2m格子に切って、各マスの下位2%を地面とみなし、同じ場所のDEM標高と比べます。
| 点群の地面 − DEM標高 | |
|---|---|
| 中央値 | +1.01 m |
| 25% / 75% | +0.74 / +1.87 m |
| 標準偏差 | 1.35 m |
仮定に頼った +1.21m と、地面から直接測った +1.01m が一致しました。 別の道筋で同じ答えが出たので、この数字は信用してよさそうです。
ばらつきが大きいのは、壁や木しか写っていないマスが「地面」を高く見積もるためです。
一方このずれは、自己申告の垂直精度1.33mの内側に収まります。高さは測位由来で、1m級の偏りを持っていると読むのが素直です。そしてこれは、DEMを使えば机上で補正できます。
10. 数字の見取り図
相対 — 形と大きさ
| 値 | 根拠 | |
|---|---|---|
| 実寸(スケール) | 1%以内 | 巻尺173cm → 172.6cm |
| 再投影誤差 | 1.6〜1.8 px | COLMAP |
絶対 — 置き場所と向き
| 中央値 | 備考 | |
|---|---|---|
| 水平位置 | 2.36 m | 68%信頼半径。3回に1回はこの外 |
| 方位 | 4.74 度 | これが効く |
| 高さ | 約1.0mの偏り | 点群の地面とDEMで実測。DEMで補正可能 |
方位誤差は距離に比例します。
| 基準点からの距離 | 4.74度によるずれ |
|---|---|
| 10 m | 0.83 m |
| 50 m | 4.15 m |
| 100 m | 8.30 m |
| 130 m | 10.78 m |
「形は数%、置き場所は数m」。 これがこの手法の実力です(8-1の図がそのまま見取り図になります)。
用途で言うと、出来形の寸法を見るには足りるが、そのまま座標を信じて杭は打てない、というあたりですね。
11. 踏んだやつ一覧
実際に踏んだものだけ書きます。
| 症状 | 原因 |
|---|---|
| 点がほとんど起きない | 縦持ちの90度ずれ。視線軸まわりの回転補正を忘れている |
| 点群が鏡像になる | 左手系→右手系の反転を片方しかかけていない |
| 実寸が2割ずれる | 基線長不足。立ち止まって撮っただけでは測れない |
| 標高が1〜2mずれる | ジオイド高を概算のままにしている。場所ごとに調べる |
| 点群が2m平行移動する | 基準の緯度経度と基準の局所座標を、別々に中央値で取っていた(同じ1コマから取ること) |
| カメラの高さが1.1m高い | Unity XR Origin の Camera Y Offset の既定値 |
| CSVの1列目が読めない | BOM付き。utf-8-sig で開く |
| 点群がスカスカなのにエラーが出ない | フレームの間引きが2の冪でない(下記) |
| 密ステレオが1時間45分後に落ちる | 参照リストだけ間引いた。image_deleter でモデルごと絞る |
stereo_fusion でWSLごと落ちる |
use_cache が既定0。全ワークスペースをRAMに載せようとする |
| 対応点の98%が棄却される | ARのポーズをMVSに使った。解き直して合わせ直す |
11-1. 間引きは2の冪でないと全損する
これは一番ひどかったので独立させます。
sequential_matcher の quadratic overlap が張るペアは、画像インデックスの差が 1, 2, 4, …, 512 の2の冪だけです。だから等間隔に間引くと、こうなります。
| 間引き | 枚数 | セット内の有効ペア | 判定 |
|---|---|---|---|
| なし | 2,732 | 15,815 (100%) | 可 |
| 1/2 | 1,366 | 6,533 (41.3%) | 可 |
| 1/3 | 911 | 0 (0%) | 不可 |
| 1/4 | 683 | 2,583 (16.3%) | 可 |
1/3で間引くと、ペアが1組も残りません。 そしてエラーも警告も出ず最後まで走りきります。 出てくるのは、点は少ないが一見それらしい点群。
私は最初これで911枚のモデルを作って、平均トラック長2.13という数字を「まあこんなものか」と眺めていました。全損してました。
colmap model_analyzer の平均トラック長で早期に気づけます。2台なら異常だと思ってください。任意の枚数が欲しいなら、間引いた後のセットにマッチャーを掛け直すことです。
12. まだできていないこと
正直に並べておきます。
- レンズ歪みを補正していません。 AR Foundationが係数を返さないので、歪みなしで近似しています
- 撮り方の問題は残っています。 同じ場所を一度しか通らないと、解き直してもトラック長は伸びきりません
- 高さが地形を追っているかは未判定。撮影地の起伏が0.51mしかなくて判定できませんでした。坂のある場所を探してやってみます。
13. 試すには
撮る側
- ARCore対応のAndroid端末(私は Google Pixel 9a)
- Unity 6000.5.6f1 / AR Foundation 6.5.0 / ARCore Extensions 1.54.0
- Geospatial API 用の Google Cloud APIキー(無料枠で十分足ります)
点群にする側
- COLMAP(疎な点群だけならCPUでも動きます。密な点群までやるならCUDA版とGPUが必要です。。私はRTX5070を使用しています。)
- Python 3(追加インストール不要。標準ライブラリだけで動きます)
流れは3章の図の通りで、コマンドにするとこれだけです。
# 1. 撮る(アプリで撮影 → 端末から取り出し)
./tools/pull.ps1 -Latest
# 2. 点群にする(姿勢の変換 → COLMAP)
python tools/frames_to_colmap.py <rec_dir>
colmap feature_extractor ...
colmap sequential_matcher ...
# 3. LASにする(ジオイド高を調べて座標を付ける)
python tools/check_dem.py <rec_dir>
python tools/points_to_las.py <model> <rec_dir> --zone 7 --geoid 37.75
手順の全文はリポジトリの README と各ドキュメント(現地チェックリスト、COLMAP手順、GPU機での密ステレオ)にあります。撮るときのコツは本文の通り、「方位を調整中」のまま撮り始めない、基準物のまわりは歩く、の2つです。
出力は LAS 1.2(平面直角座標+標高)なので、CloudCompare でそのまま開けます。
Unityの場合は私作のOSS(UnityLasImporter)も是非使用してください。
まとめ
- 写真からは原理的に大きさが決まりません。すべてを $k$ 倍しても同じ像になるからです
- どの分野も既知の長さを1本入れています。今回は加速度計がその役をします。ばね定数と質量は工場で校正された絶対量なので、出力の $\mathrm{m/s^2}$ はメートルに直結しています
- 実際に測ったら、実寸は1%以内でした(巻尺173cm → 172.6cm)
- ただし地球上のどこに置くかは数m級です。この2つを混ぜないでください
- ARCoreの緯度経度は連続した測位ではありません。区間内では局所座標に固定の変換を掛けたもので、新しい情報が入るのはVPSが解を組み直す瞬間だけでした
以下、MITライセンスでOSSとして使用可能です。
ツールはPython標準ライブラリだけで動きます。特に平面直角座標の変換(全19系、国土地理院と0.1mm一致)は単体で使えるので、必要な方は持っていってください。
長い記事にお付き合いいただきありがとうございました。
夏季休暇を利用したちょっとした自由研究をしてみました。
私自身、センサーや地理情報の活用方法を再学習することができて楽しかったです。
中学・高校とちょっと大学学部レベルの数学・物理でここまでできてしまうの、
ちょっと面白くないですか?
みなさんも是非点群データを撮影・活用してみてください。
ではでは皆さま、今日も一日、ご安全に!














