おっぱい山をFOSS4Gでアーカイブする方法

この記事はFOSS4G Advent Calendar 2016の20日目の記事です。

今年のFOSS4G Tokyoの特別セッションは「アーカイブ×FOSS4G」というテーマで、大変刺激を受けました。
ということで私もそれにのっかって、おっぱい山をFOSS4Gでアーカイブする方法を紹介したいと思います。

仙台市上空よりおっぱい山を望む（巨乳化処理済み）

1. おっぱい山を作る

公園の砂場に行き、おっぱい山を作ります。おっさん一人で砂場で遊んでいると通報されますので、自分の子か親戚の子と作業してください。

作業の様子	おっぱい山

2. 写真を撮る

おっぱい山が出来上がったら、UAV(ドローン)を用いた公共測量マニュアル（案）に従って写真を撮ります。ドローンはありませんので、スマホを手に砂場の周りをぐるっと回ってシャッターを押してください。注意点は、写真に位置情報が記録されるよう設定しておくこと、7割以上オーバーラップさせて撮影すること、公園にいるママさん達に不審がられないようにすることです。

撮影の様子。自分の影が写らないように。	撮影したおっぱい山（80枚）

3. OpenDroneMapで3Dモデルとオルソ画像を作る

写真が準備できたら、OpenDroneMapで3Dモデルとオルソ画像を作成します。OpenDroneMapは、ドローンで撮影した空中写真から地形の3Dモデルとオルソ画像を作ることができるオープンソースのソフトウェアです。*1ドローンだけでなく、スマホで撮影した写真からでも作成できます。

3-1. Bash on Windowsの導入

OpenDroneMapをWindowsで使うためにBash on Windowsを導入します。Bash on Windowsとは、Windowsで動作するBashシェルです。これによってOpenDroneMapのようなLinux用のソフトもWindowsでコンパイル＆実行できるようになります。導入方法はここでは説明しないので、適当にググってください。

3-2. OpenDroneMapのコンパイル

Bash on Windowsを起動して、以下のコマンドでGitHubからOpenDroneMapのソースをダウンロードします。

git clone https://github.com/OpenDroneMap/OpenDroneMap.git

~/.bashrcにパスを追加して、sourceコマンドで設定を読み直します。/your/path/の部分はソースのある場所に置き換えてください。

export PYTHONPATH=$PYTHONPATH:/your/path/OpenDroneMap/SuperBuild/install/lib/python2.7/dist-packages
export PYTHONPATH=$PYTHONPATH:/your/path/OpenDroneMap/SuperBuild/src/opensfm
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/your/path/OpenDroneMap/SuperBuild/install/lib

source ~/.bashrc

OpenDroneMapのフォルダに移動してコンパイルします。依存するライブラリなど諸々インストールされるので結構時間がかかります。(1時間ほど)

cd OpenDroneMap
bash configure.sh

3-3. GCP(グラウンドコントロールポイント)の作成

OpenDroneMapを実行する前に、写真の位置あわせをするためのGCPを作成しておきます。スマホの写真には位置情報が記録されているのでGCPがなくてもオルソ画像は出来ますが、精度をあげるために作成します。GCPファイルは、gcp_list.txtという名前で以下のフォーマットで作成します。

coordinate system description
x1 y1 z1 pixelx1 pixely1 imagename1
x2 y2 z2 pixelx2 pixely2 imagename2
x3 y3 z3 pixelx3 pixely3 imagename3

coordinate system descriptionは位置の投影法、x1,y1,z1は位置座標、pixelx1,pixely1は写真のピクセル座標、imagenameは写真ファイル名です。1枚の写真に対して最低5点ぐらいはあると良いようです。具体的には、こんな感じです。

WGS84 UTM 54N
494399.2 4232697.1 5 333 1590 IMG_6785.jpg
494399.2 4232697.1 5 758 1698 IMG_6786.jpg
494399.2 4232697.1 5 877 1942 IMG_6787.jpg

手動でこのファイルを作ろうとすると結構手間ですが、支援ツールもあるようなので、利用してみると良いかもしれません。*2

3-4. OpenDroneMapを実行する

プログラムを実行する準備として、OpenDroneMapフォルダの直下にoppaiyamaフォルダを作成します。さらにoppaiyamaフォルダの中にimagesフォルダを作って、そこにスマホで撮影したおっぱい山の画像を入れておきます。gcp_list.txtはoppaiyamaフォルダの直下に入れます。
その後、以下のコマンドでプログラムを実行します。写真の枚数や設定にもよりますが、80枚で30分ぐらいはかかるので、気長に待ちます。

sh run.sh --project-path oppaiyama --odm_georeferencing-useGcp  --odm_orthophoto-resolution 1000

実行パラメータは、色々ありますが、今回は、gcpを使用して、オルソ画像の解像度を1000pixcel/mと設定しておきます。
パラメータの詳細は、ここを参照してください。https://github.com/OpenDroneMap/OpenDroneMap/wiki/Run-Time-Parameters

3-5. データを確認する

プログラムが終了したら、oppaiyamaフォルダの中に3Dモデルとオルソ画像が出来上がっています。3Dモデルは、odm_texturingフォルダの中に「odm_textured_model_geo.obj」という名前で入っているので、MeshLabなどで読み込んで確認してみましょう。また、オルソ画像は、odm_orthophotoフォルダの中に「odm_orthophoto.tif」という名前で入っているので、QGISに読み込んで確認してみましょう。

MeshLabで3Dモデルを表示	QGISでオルソ画像を表示

4. OpenDroneMapのデータをCesium用に変換する

OpenDroneMapで作成した地形データをCesiumで表示するためデータを変換します。Cesiumとは、地図をブラウザで3D表示できるオープンソースのjavascriptライブラリです。Cesiumの使い方は、このあと説明しますが、とりあえず読み込むためのデータを作成しておきます。

CesiumにOpenDroneMapのデータを読み込む方法は二通りあります。一つは、objファイルをgltf形式の3Dモデルに変換して読み込む方法です。もう一つは、lasファイル（点群データ）をCesium用の標高タイルに変換して読み込む方法です。

4-1.（方法その1） objファイルをgltf形式の3Dモデルに変換する

おっぱい山を3Dモデルとして読み込むために、まずmeshlabでobjファイルをcollada(dae)形式に変換します。
odm_texturingフォルダの中のodm_textured_model_geo.objをMeshLabに読みこんでcollada(dae)形式でエクスポートします。その際、MeshLabでいらない部分を消すなど編集しておくと良いかもしれません。ファイル名は「oppaiyama.dae」としておきます。

次に、collada(dae)形式をgltf形式に変換するためcollada2gltfというソフトをGitHubからダウンロードしてコンパイルします。*3

git clone https://github.com/KhronosGroup/COLLADA2GLTF.git
cd COLLADA2GLTF
cmake .
make
sudo make install

以下のコマンドで、oppaiyama.daeをoppaiyama.gltfに変換します。-eはテクスチャ等も同梱するオプションです。

collada2gltf -e -f oppaiyama.dae -o oppaiyama.gltf

collada2gltfが簡単に利用できるのもBash on Windowsのおかげです。slコマンドも動くよ！

4-2.（方法その2）lasファイルを標高タイルに変換する

おっぱい山を標高タイルとして読み込むために、LastoolsでlasファイルをDEMデータに変換します。まず、Lastoolsを以下からダウンロードます。
http://www.cs.unc.edu/~isenburg/lastools/
Lastoolsの一部のツールは残念ながらフリーウェアではありません。今回使用するlas2dem.exeも商用使用だと有料になりますが、個人的に使用する分にはサイズ制限はありますが使用できるようです。*4ソフトを解凍したら以下のコマンドで変換します。ただし、las2dem.exeはWindows用の実行ファイルなので、Bash on Windowsではなく、別途、PowerShellを起動して実行します。変換するlasファイルは、odm_georeferencingフォルダに入っているodm_georeferenced_model.ply.lasを使用します。

cd Lastools\bin
las2dem.exe -i odm_georeferenced_model.ply.las -o dem.tif -elevation -step 0.001 -utm 54N -nodata 0

オプションの設定は、高さ方向0.001mの解像度で、UTM54の投影法、nodataは0の標高データをdem.tifで出力するとなっています。

次に、DEMデータをCesium用の標高タイル(heightmap-1.0 terrain format)*5に変換するため、Cesium Terrain Builderをコンパイルします。Cesium Terrain Builderのコンパイルには、gdal2.0.0以上が必要なので、以下のコマンドでgdalをアップデートしておきます。(Bash On Windowsでの作業に戻っているので間違えないように)

sudo add-apt-repository ppa:ubuntugis/ppa && sudo apt-get update
sudo apt-get install libgeos-c1v5 gdal-bin

gdalをアップデートできたら、Cesium Terrain Builderのソースをダウンロードしてコンパイルします。

git clone https://github.com/geo-data/cesium-terrain-builder.git
cd cesium-terrain-builder
mkdir build && cd build && cmake ..
sudo make install

ソフトは準備できましたが、今回のデータは対象範囲が狭すぎるためそのままでは上手く変換できません。そのため、Cesium Terrin Builderで変換する前にgdal_translateで巨乳化処理をしておきます。
gdal_translateの-a_ullrオプションで範囲を1000倍に拡大し、-scaleオプションで高さも1000倍にして最低標高を調整しておきます。ついでに-projwinオプションで必要な範囲にクリップします。(nodataを範囲に含むと変換結果がおかしくなる？ので切り取る)

gdal_translate -projwin 494398.737 4232696.642 494400.030 4232694.956 -scale 4.740 5.160 -80 340 -a_nodata -9999 -a_ullr 494398.737 4232696.642 495893.03 4231008.956 dem.tif  dem_scale.tif

Cesium Terrain Builderのctb-tileコマンドで、DEMデータを標高タイルに変換します。標高タイルはterrainフォルダに保存します。

mkdir terrain
ctb-tile -o terrain dem_scale.tif

標高タイルには以下の定義ファイルも必要なので、「layer.json」という名前でterrainフォルダの中に入れておきます。

{
  "tilejson": "2.1.0",
  "format": "heightmap-1.0",
  "version": "1.0.0",
  "scheme": "tms",
  "tiles": ["{z}/{x}/{y}.terrain"]
}

また、0/1/0.terrainと0/0/0.terrainは必須(?)のようなので、0/0/0.terrainを、そのまま0/1/0.terrainとしてコピーしておきます。（今回のデータ範囲は、0/0/0の領域ではないので作成されていないため）

オルソ画像も地図タイルにしておきます。dem.tifと同じように範囲を1000倍に拡大して、png形式の地図タイルに変換しておきます。なお、ctb-tileで作成される地図タイルはtmsなので原点は南になります。

gdal_translate -projwin 494398.737 4232696.642 494400.030 4232694.956 -a_ullr 494398.737 4232696.642 495893.03 4231008.956 odm_orthophoto.tif ortho_scale.tif
mkdir ortho
ctb-tile -p mercator -f png -o ortho ortho_scale.tif

ctb-tileのオプションの詳細は、こちらにあります。https://github.com/geo-data/cesium-terrain-builder

5. Cesiumでおっぱい山を表示する

おっぱい山のデータが出来上がったので、これをCesiumで表示してみます。

5-1. Cesiumを用意する

Cesiumのサイトからライブラリをダウンロードします。
https://cesiumjs.org/

解凍するとファイルとフォルダが色々入っていますが、使用するのはAppsフォルダの中のHelloWorld.htmlとBuildフォルダの中のCesiumフォルダ一式になります。
ひとまず動作確認のために、Bash on Windowsで以下のコマンドを打って、簡易ウェブサーバーを起動します。その後、ブラウザで http://localhost:8000/ にアクセスして「Hello World」のリンクをクリックしてみましょう。Ctrl+Cでサーバーを終了できます。

cd Cesium-1.28
python -m SimpleHTTPServer

5-2. gltfの読み込み

gltfを読み込む場合は、以下のコードをHelloWorld.htmlに追加します。3Dモデルを正しい方向で読み込むには、中心を軸に-90度回転させる必要があります。高さの中心は、今回は見た目で微調整しました。Cesiumの高さの基準は楕円体面（海面0mではない）なので、そのへんの調整も必要です。

    var center = Cesium.Cartesian3.fromDegrees(140.9360018,  38.2422495, -5.5);//地形を読み込む場合は42.5にする
    var heading = 0.0;
    var pitch = 0.0;
    var roll = -Math.PI / 2.0;
    var hpr = new Cesium.HeadingPitchRoll(heading, pitch, roll);
    var transform = Cesium.Transforms.headingPitchRollToFixedFrame(center, hpr);
    var model = viewer.scene.primitives.add(Cesium.Model.fromGltf({
         url:"./oppaiyama.gltf",modelMatrix : transform, scale : 1.0})
    );

5-3. 標高タイルとオルソ画像の読み込み

標高タイルとオルソ画像の地図タイルを読み込む場合は、以下のコードをHelloWorld.htmlに追加します。オルソ画像の地図タイルは、南が原点となっているので{reverseY}と指定します。

var terrainProvider = new Cesium.CesiumTerrainProvider({
    	url : './terrain'
});
viewer.terrainProvider = terrainProvider;

var ortho = new Cesium.UrlTemplateImageryProvider({ 
    url: './ortho/{z}/{x}/{reverseY}.png'
});
viewer.scene.imageryLayers.addImageryProvider(ortho);

5-4. おっぱい山の表示

ついにCesiumで表示できるところまで来ました。5-1.と同様に簡易ウェブサーバーを起動して確認したいところですが、gzip形式で圧縮されている標高タイルは、そのままではSimpleHTTPServerでは配信できません。そこで、元のSimpleHTTPServer.pyをコピーして以下のコードを加えます。こうすると拡張子が.terrainの標高タイルはgzip形式として配信されます。

cp /usr/lib/python2.7/SimpleHTTPServer.py ./Cesium-1.28

self.send_response(200)
self.send_header("Content-type", ctype)
## SimpleHTTPServer.pyの95行目以降に以下の3行加える（send_head関数内）##
base, ext = posixpath.splitext(path)
if ext == ".terrain":
   self.send_header("Content-Encoding", "gzip")
###################################################

なお、通常のウェブサーバーの場合は.htaccessを以下のように記述して設置しておきます。

AddType application/octet-stream .terrain
AddEncoding x-gzip .terrain

SimpleHTTPServer.pyを変更したら簡易サーバーを起動しましょう。5-1.の時と違って、-m オプションは付けずに.pyも含むので間違えないように。

python SimpleHTTPServer.py

http://localhost:8000/Apps/HelloWorld.html にアクセスしてみましょう。（簡易サーバーだと、読み込みに少し時間がかかります。また、時々gltfの読み込みに失敗するので、表示されない場合はShift+リロードしてみます。）

OpenDroneMapで作成した3DモデルをCesiumで表示
http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld1.html



OpenDroneMapで作成した点群データを標高タイルに変換してCesiumで表示
http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld2.html

6. アーカイブをVRで体験

ここまでの作業で、おっぱい山をアーカイブして、Cesiumで表示することができました。Cesiumでは、さらにヴァーチャルリアルティ(VR)表示も可能なので、やってみましょう。以下のようにhtmlを変更すると、画面右下にVRボタンが表示されるようになります。

var viewer = new Cesium.Viewer('cesiumContainer', {
    vrButton : true
});

VRボタンを押すと2画面表示になり、Google Cardboardのようなヘッドセットで見ると臨場感あふれるバーチャルリアリティ体験ができます。

※ VR表示はiPhoneには対応していませんが、PCでVR表示させた画面をリモートデスクトップでスマホに表示すれば可能になります。

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld.html

まとめ

子供と作った砂場のおっぱい山は、夕日とともに崩れてなくなってしまいますが、このようにアーカイブしておけば、思い出と一緒に永遠に残しておくことができます。
プライスレスなこの活動、FOSS4Gがあれば誰でもできます。皆さんもぜひお試しください！それでは、良いクリスマスを。

追伸：
Google Cardboadは、このあと別のおっぱい山の閲覧に有効活用させていただきました。VRすげーよ。ムフフ。

参考リンク

この記事に使用した画像、コード

https://github.com/tmizu23/cesium_oppaiyama

おっぱい山アーカイブ

http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld.html
http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld1.html
http://ecoris.co.jp/map/Cesium-1.27/Apps/HelloWorld2.html

OpenDroneMap関連

https://github.com/OpenDroneMap/OpenDroneMap
https://github.com/OpenDroneMap/OpenDroneMap/wiki/Running-OpenDroneMap
https://github.com/OpenDroneMap/OpenDroneMap/wiki/Run-Time-Parameters
http://lists.osgeo.org/pipermail//opendronemap-users/2016-June/000314.html
https://github.com/wolkstein/OpenDroneMap-GCP_LIST.TXT-generator

データ変換関連

http://www.cs.unc.edu/~isenburg/lastools/
https://www.cs.unc.edu/~isenburg/lastools/download/las2dem_README.txt
https://github.com/KhronosGroup/COLLADA2GLTF
https://github.com/geo-data/cesium-terrain-builder
https://github.com/AnalyticalGraphicsInc/obj2gltf
http://www.sarasafavi.com/installing-gdalogr-on-ubuntu.html
http://blog.mastermaps.com/2014/10/3d-terrains-with-cesium.html
http://blog.mastermaps.com/2016/09/creating-tins-in-saga-gis.html
http://blog.mastermaps.com/2016/09/creating-tin-from-raster-dem.html
https://cesiumjs.org/data-and-assets/terrain/formats/heightmap-1.0.html
http://cesiumjs.org/data-and-assets/terrain/formats/quantized-mesh-1.0.html
http://www.pdal.io/

Cesium関連

https://cesiumjs.org/
http://www.slideshare.net/makinux7/cesium3
http://www.slideshare.net/yuki1225/foss4g-2016-hokkaidocesium
http://www.slideshare.net/Ihizaki/cesium-63857910
http://stackoverflow.com/questions/38709743/cesiums-height-value-is-altitude-or-above-sea-level

そのほか

http://meshlab.sourceforge.net/
https://vr.google.com/intl/ja_jp/cardboard/
http://psgsv2.gsi.go.jp/koukyou/public/uav/
http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q10144883875
http://www.osgeo.jp/events/foss4g-2016/foss4g-2016-tokyo/foss4g-2016-tokyo-coreday

*1:商用だとPhotoScanとかPix4Dなどがありますが、これらはウン十万円もする高価なソフトで、個人ではとても手がでません。そのようなツールと同等なものを使えるのは、ほんとFOSS4Gのおかげです。ただ、まだまだ開発途中で商用のものに及ばない部分は沢山ありますが。

*2:https://github.com/wolkstein/OpenDroneMap-GCP_LIST.TXT-generator

*3:obje形式から直接gltf形式に変換できるobj2gltfというソフトもありますが、テクスチャが表示できない不具合(?)があったので見送りました。

*4:PDALとかOpenDroneMapでDEMに変換できると良いのですが

*5:Cesium用の標高タイルはheightmap-1.0 形式のほかに、quantized-mesh-1.0 形式というものがあります。こちらの形式は標高データをTINにして、それをタイルに切り分ける方法のようですが、まだ使いやすい変換ソフトが無いようなので試せていません