写真測量：MeshLabによる点群ファイルのジオリファレンシング

MicMacやWebODMで何らかの理由でジオリファレンシングできなかったり、COLMAPやRegard3Dのようにジオリファレンシングできない3Dモデリングソフトを使ったりしてplyファイルだけが得られた場合を考えます。このplyファイルをMeshLabでジオリファレンシングし、GISで表示可能なcsvファイルにします。DEMはGIS上で作れます。

標定点（GCP）データの準備
Reference sceneツールを使う
変換結果の確定（Freeze）
変換結果のエクスポート
必要ならオフセットを元に戻す
QGISに読み込んでみる

標定点（GCP）データの準備

下図のMeshLabには、Rigard3Dで作成した密な点群を読み込んであります。

MeshLabでのジオリファレンシングには最低４点の標定点が必要です。実際にどう測量するかは別に譲るとして、ここでは平面直角座標系で４点準備しました。（緯度経度はだめです。）

なお、MeshLabの標定点座標は大きな桁数を扱えないようですので、(-80000, -215000)をオフセットしたgcp2.txtファイルを使います。

標定点座標とオフセットした標定点座標 — gcp.txt: 標定点座標、gcp2.txt: オフセットした標定点座標

Reference sceneツールを使う

Reference sceneツールを開き、「Load Reference Points From File」ボタンでgcp2.txtファイルを開きます。すると下図右側のように、標定点座標が読み込まれます。あとは、plyファイルの座標を設定したい行番号をクリックして選択し、「Pick current point on MOVING」ボタンでMeshLab画面上の対応する点をクリックしていけば良いのです。

下図は標定点１を画面上でダブルクリックしてplyの座標値を読み込んだところです。

すだくん

やっぱり標定点を画像じゃなく点群上で探すのは難しいよ。大きな対象物を選ばないと大変だ。GUIで標定点を探せるWebODMはやっぱりいいね。

４点すべて選択し、「Calculate Rototranslation」ボタンを押すと、下図のようにError値が表示されます。また、右下の画面には座標変換パラメータが表示されています。この誤差でよければ、「Apply」ボタンを押します。

かわのさん

誤差は12cm～38cmだけど、見にくい密な点群の上で標定点を捜したにしては上出来じゃない？

すだくん

それそれ。この元画像は100m上空から撮影したものだけど、手持ちカメラで川などを撮影したものなら誤差はもっと小さくできるよね。

変換結果の確定（Freeze）

ここまでの作業で、MeshLab上の表示の上では座標変換が実行され、例えば距離を測定すると実距離が出てきますが、実際にはデータそのものは変更されていません。ウインドウ右側のレイヤパネルを見ると、レイヤ名にアスタリスク（＊）がついていることで未確定であることがわかります。変換結果をエクスポートするには、まず変換結果を確定する必要があります。

レイヤパネルでレイヤを選択し、右ボタンメニューでMatrix: Freeze Current Matrixを選びます。出てきたダイアログでApplyボタンをクリックすると結果が確定され、レイヤのアスタリスクも消えます。

変換結果のエクスポート

「File」 – 「Export Mesh As…」メニューで、「XYZ Point Cloud (with or without normal) (*.xyz)」を選び、ファイル名を決めます。すると右側のダイアログが表示されるので、Normalのチェックを外してOKします。以上で、ジオリファレンスがされた点の座標が、下図の例ではpointcloud.xyzというファイルとして出力されます。このファイルは、単にx,y,z値がスペースで区切られて１行に１点記録されたテキストファイルです。

以上でMeshLabでの作業は終わりです。

必要ならオフセットを元に戻す

xyzファイルは単純なテキストファイルですので、Excelなどで簡単にオフセットを加算して平面直角座標系に戻すことができるのですが、点の数が多すぎてExcelの行数に収まらないことも起こります。Windows上でできるだけ簡単にオフセットを加算したCSVファイルに変換するために、VBScriptを作成してみました。

以下のコードをメモ帳にコピペしてxoffset, yoffsetの値を修正し、convert.vbsというファイル名で保存し、コマンドプロンプトで以下の冒頭コメントのように実行します。InFile.xyzと書いてあるところが、今回の例ではpointcloud.xyzということになります。

'VBScriptでxyzファイルにオフセットを加算したcsvを作成
' あらかじめxoffset, yoffsetをこのスクリプトconvert.vbsに書き込んでおく
' CMD.EXEを起動してxyzファイルのあるフォルダに移動し、
' cscript convert.vbs < InFile.xyz > OutFile.csv
' でオフセットが加算される。"<" ">"は標準入出力へのリダイレクション

Option Explicit
Const xoffset = -80000
Const yoffset = -215000

On Error Resume Next
Dim objStdIn
Dim objStdOut
Dim intExitCode

Set objStdIn = Wscript.StdIn

Set objStdOut = Wscript.StdOut

Dim line
Do While objStdIn.AtEndOfStream = false
  line = Split(objStdIn.ReadLine)
  objStdOut.WriteLine  (line(0)+xoffset)&","&(line(1)+yoffset)&","&line(2)
Loop

Set objStdIn   = Nothing
Set objStdOut  = Nothing

If Err.Number <> 0 Then
    intExitCode = 1
Else
    intExitCode = 0
End If

Wscript.Quit(intExitCode)

ファイルに出さずに画面に出すと以下のようになります。Copyright行が出ていますので、変換結果を利用するときにはエディタで削除するかQGIS側の読み込み設定で読み飛ばすかしてください。