元伊勢大饗石のLiDAR測量 LiDAR measurement of Oo-Ae-ishi, a big table rock for a banquet, near the former Grand Shrine of Ise

はじめに

 このpostingでは,このテーマに限定する。iPhone 12 Proによって大饗石のLiDAR測量結果は得ている。さらに,幾つかの手順で外観だけでなく,数値データも整理する必要がある。そのアプリとして最も利用されているのが,CloudCompareである。Windows版のオープンソースであり,ドーネーションが期待されている。ここではインストールされていることを前提とする。メニューの一部が日本語化されているバージョンもあるが翻訳が不完全で,日本語化されているものはお勧めできない。
 ぼくはmouseのWindows 10とmac内のWindows 10 on Parallels Desktopにインストール (2022年3月実施)している。このウェブページでのコンテンツ利用には,mac内が便利ではあるが,多々支障はある。なお,mac版も提供されているが以前使って問題があって,ぼくのmacからは削除している。以下を実行する上で,次のページを読み込む必要性がある。

 iPhone 12 Pro撮影の3Dスキャン画像の座標を捉える
 裸岩露頭のiPhone 12 Proを使った点群撮影
 二つの3Dスキャンマップを繋ぐ

1  iPhone 12 ProからCloudCompareにファイル転送する

 iPhone 12 Proからファイル出力の際には,USDZ以外,Metascan Proの契約をしていないとできない。唯一有効なPC処理アプリCloudcompareは,残念ながら,USDZに対応していない。大饗石のLiDAR撮影結果は,すでにMeshファイルとしてはFBX(121.5MB),Point CloudファイルとしてはPLY(81.8MB)とXYZ(zip 89.9MB)を出力し,macに取り込んでいる。XYZ出力はzipファイルとして出力されているので,解凍する必要がある。この順でCloudCompareに取り込んだ結果が次の図1である。
 左端のツールのうち,虫眼鏡直下のSet top viewが図1のコンソールに現れている。このコンソールの右下隅には,3D軸が見え,これは右手座標系でZ軸は手前方向に延びている。

 図1の中央の赤線で囲んだ範囲が大饗石で,この上に何か白いキノコのようなものが見える。図2ではこの図1を少し回転したもので,赤丸で囲んだこの白いキノコのようなものは大饗石よりも高い位置にあることがわかる。スキャンニングの際に,大饗石よりも高い位置の「障害物」が捉えられているのである。大饗石だけを残す形で切り取ることになる。

図1 CloudCompareでのPLYファイル表示
図2 ゴミは大饗石より上方に

以上,Dec. 31, 2022記。

2 必要な部分を切り抜く

 このLiDAR測量の際の意図をはっきりとは,もう覚えていない! 大饗石のみを詳細にスキャンしたつもりではあるが,大饗石周辺もかなり捉えられている。大饗石のみを切り出そうと考えていたが,まずはより広い範囲をまずは抜き出したいと思う。そこで,topビューで大饗石とその周辺を切り出すのではなくて,地上の木本以外の地形を切り出すべく,z軸に直交する平面で切り取りたいと思う。
 CloudCompareで不要なポイント群または面群の削除 にはその手法を整理しているので,参照して欲しい。この参照を前提に,以下,記述する。

1 まずはこの例では,set lift sideして,ツール左端の‖をクリックして,Segmentation [PAUSED] にして,フロートや邪魔になりそうなスギの頭などがsegmentしやすいように回転などして,ツール左端のトグルスウィッチ‖をクリックして,segmentationを実行してゆく(図3→図4)。segmentationを実行すると,Segmentation [PAUSED]にもどるので,位置決めをして,また,トグルスウィッチ‖をクリックして,segmentationを実行してゆく。

図3 polygonでsegmentationの範囲を選ぶ
図4 Segment Outを実行後

 一応,終了しよう。ツールの✅️(confirm segmentation Enter)した。その直後が図5で,remainingとsegmentedがともに表示されている。前者の✓を外したのが図6であり,作業結果が見える。set right side view,表示である。ほぼ側面図にあたる。中央のトップが水平に近い岩塊が凝灰岩(fall)からなる大饗石で,左手の岩塊はデイサイトの大岩塊である。いずれも前期中新世の溶岩または凝灰岩層から,崩落したものである。
 保存作業であるが,remainingファイルはこれを右クリックで選んでdeleteした。

図5 confirm直後
図6 remainingを非表示に

 作業結果のクラウドを保存するには,segmented Fileを選んで,保存することになるが。File > saveの形で,ファイル名に_segmentedを追加して,元ファイルのフォルダに保存すべく実行して,確認のために,このtxtファイルを読み込んで表示したのが図7であるが,残念ながらRGB表示ができない。
 図8のように,保存する際に下段にファイル名とファイルの種類を指定する必要性がある。デフォールトでは元々のファイル名がそのまま使用され,ファイルの種類もtxtファイルになってしまう。図7はファイル名には_segmentedを追加したがtxtファイルで保存してしまった。図8ではファイルの種類をbinつまり,ClouCompare entities,にして保存したのである。binとは,binary拡張子で,バイナリー形式のデータであることを表している。txtファイル形式ではCloudCompareの機能を発揮できないのである。

図7 txtファイルで保存した場合
図8 binファイルで保存した場合

 それゆえ,処理のそれぞれの画期で,ファイル名に作業を付加したbinファイルを保存してゆく必要がある。図6の形で完成とも考えることができる。大饗石だけよりも,立地がよく理解できる。とはいえ,大饗石の3D情報を得たいと考えての調査であったので,この大饗石だけを刳り抜いて見たい。図9は大饗石のズームインである。

図9 大饗石ズームイン

3 大饗石だけの刳り貫き

大饗石3Dクラウドを確保したい。図10はset top viewにして表示したもので,この形でsegmentedしたが不明瞭な部分が多く,図11〜13のようにパースで見て,segmentedしていった。図14は新たにbinファイルとして保存したものである。binファイル名は,metascan_20221121-1440大饗石周辺のロックフォール_onlyOoaeishi.binというように名付けることができるが,この下位のクラウドの名称は自動的に,metascan_20221121-1440大饗石周辺のロックフォール – Cloud.segmented.segmented,となる。

図10 Set top view
図11 Segmented 3D
図12 Segmented 3D
図13 Segmented 3D
図14 Only Ooaeishiのbinファイル

以上,Jan. 1, 2022記。

4 大饗石の体積を求める1: まずはデフォルトで

 3Dファイルにはメッシュmeshと点群point cloudがある。画像的な表現では前者がより優れているが,右座標系など幾何学的な分析には後者がすぐれているようにぼくは感じている。体積計算に関しては,公式マニュアルは全く役に立たない。
 では,http://www.cloudcompare.org/doc/wiki/index.php/Main_Page,はどうだろうか。2.5Dで検索すると,2.5D Volume のページに行きつく。このページを理解した上での使用例として,下記を挙げる。CloudCompareでの体積計算 その2 (その1はこの「その2」の予告だけ)
(引用1)2.5D Volumeの説明に従って,大饗石の体積計算に挑む。このアプリは次の機能がある。
This tool can be used to compute the volume between a 2.5D cloud and an arbitrary plane (constant height) or between two 2.5D clouds.
 二つのクラウド間の扱う例としては,例えば工事前Beforeと工事後Afterで,土量のプラスマイナスを求めることができる。一つだけのクラウドで体積を求めるというのは,この大饗石の体積を求めたい場合である。図6のように,大饗石のトップはほぼ水平で下部もまあ水平と見做すことができるので,下部の最も低い部分の高度を使えばいいと思う。(引用1)の「an arbitrary plane (constant height) 」をその高度とすれば良いだろうということになる。

1 フォルダと対応ファイルに✓して,対応ファイルを選ぶ(ハイライト)。そして,Tools > Volume > Compute 2.5D volume 。
 図15はleft side viewである。奥行きの軸がXで,縦長の軸がY,そして垂直軸はZとなっている。左のペーンのPropertiesのBox dimensionsを見ると,X: 5.39993, Y: 7.43598, Z: 2.96101,となっている。iPhone 12 ProのLiDAR測量は海抜高度を反映していないので,上記の例とはかなり異なり,数字が大きくない。大饗石は地中に多少埋まっているであろうが,凝灰岩の層理を反映しており,後掲するVolume caluculationのパネルでのGround(空間)/Before(時間)のGroundと考えて良いだろう。

図15 left side view

 この大饗石の座標値が未だ理解できないので,Point list picking(6点)を実施してみた。図16がそれである。この画像表示はpoint間の位置関係がわかるように岩の上部を少し手前に回転している。Point list pickingの表では,座標値を一目では3行しか見ることができないので,右上の表のツールのフロッピーアイコンをクリックして保存して,次のように,そのテキストを,参照した。

図16 point list picking 6点

Point #0,-0.939450204372,2.67226147652,-0.758838474751
Point #1,-0.485350400209,3.48180985451,0.743179500103
Point #2,-0.7682762146,0.985117554665,2.12557053566
Point #3,-2.89230465889,-3.08031964302,1.72712016106
Point #4,-2.35597205162,-3.85543513298,1.05731666088
Point #5,-3.39330339432,-1.8211826086,-0.585319936275

 この6点の座標値からXYZ軸を理解することができる。今,必要なのはGround値で,Point #0のZ値と見做して問題ないだろう。Ground値を-0.758 (m) とする。ミリメートル単位である。図17ではパラメータを入力している。図18は,図17の赤いベルトのUpdateを実行した結果である。図18下部のResultsが一気には見えないので,”Copy to clipboard”ボタンを押して,Google Chromeを使って,macのmailにgmailした。
 その結果を図19の左手のテキスト群に示している。図18の右手のrelative heightベルトの左手の模様は,大饗石の高度分布を示している。図19は大饗石のtop viewであるが,これと図18の模様と対応しているのである。

図17 パラメータ値の入力
図18 体積計算の結果

Volume: 70.765
Surface: 36.000
———————-

Added volume: (+)70.765
Removed volume: (-)0.000

———————-

Matching cells: 75.0%
Non-matching cells:
   ground = 25.0%
    ceil = 0.0%

Average neighbors per cell: 7.1 / 8.0

図19 top viewとrelative height

以上,Jan. 2, 2022記。

5 大饗石の体積を求める2: step (メッシュサイズ)= 0.008 m

 とにかく,デフォルトでやってみた。図19で言うと,長軸方向(上辺が斜面の下手で下辺が上手),これの直交方向が短軸方向,そして高さは図16で推定できる。大雑把にいうと,6.5 m x 3.6 m x 2.7 m = 63 ㎥ 。Volume: 70.765 ㎥ は悪い値ではない。ただ,Matching cells: 75.0%という値は適切ではない。
 さて,この体積計算アプリ(図17)で,ユーザーが設定できるパラメータは実質,”Grid”では”step”のみ,Ceil/Afterとしては指定したクラウドの”Empty cells”だけである。stepに関する説明がマニュアルにないが,図19から見ると,stepは,まさにgrid(方眼)を構成する正方形squareの一辺の長さである。step = 1.000000としたので,size = 6 x 8となったのだが,適切な英語表現に替えると,”step”はside-size of a square,”size”はnumber of squaresである。point cloudとmeshとの対比からすると,混乱するが,stepはmesh sizeで,sizeはmesh numbersである。
 体積を求めるための地物に対して,メッシュサイズ(step)を小さくするほどメッシュ数は大きくなってゆく。この図17の場合,size: 6 x 8とある。この設定での計算結果が図19に見えるが,横軸(左右)X 6,縦軸(上下)Y 8のグリッド(白いドット)が見える。かなり粗っぽく計算されたことがわかる。
 図20と21は,step = 0.008(自動的に,size = 676 x 930)で計算した結果である。サイズは1000 x 1000以内に収まるようにというアドバイスがある利用サイトにあったので,その意味の根拠はわからないが,そういう制限のもとで,step = 0.008が決まった。
 計算結果をクリップボードにコピーして,このウィンドウズからメールでペーストして,macに送った結果を次に示す。
Volume: 31.317 Surface: 14.883
———————-
Added volume: (+)31.317
Removed volume: (-)0.000
———————-
Matching cells: 37.0% Non-matching cells: ground = 63.0% ceil = 0.0%
Average neighbors per cell: 5.6 / 8.0

 体積は第4章に比べてほぼ半分になってしまったマッチングセル率は37%で,非マッチングセル率はceilではゼロだが,groundでは63%に及んでいる。これはGroundを斜面下縁(図16に続くPoint#0のZ値=-0.758)にしたことに依っているのであろう。この体積 31㎥ を提示されると,なるほど,と思ってしまう。計算のアルゴリズムは公開されていないが,メッシュ数676 x 930 = 628,680で,メッシュサイズ = 0.008 m だから,(676 x 0.008) x (930x 0.008) = 5.408 m x 7.44 m = 40.235 ㎡ になる。図21の矩形域の面積がこれにあたる。それぞれのメッシュでクラウドの深度幅を掛けて,全体の和を求めると,体積となる。粗っぽく目算すると地物はメッシュ枠の半分ほどであり,深度幅が1.5mほどだすると,体積30 ㎥ 余りとなるのである。

図20 step 0.008, size 676 x 930
図21 同 DB TreeとProperties

 さて,図22と23は,この地物の空っぽの裏が見えている。この体積計算過程では,LiDAR測量結果が正しく把握されている。色分けはZ値を示すもので,もちろんトップ面のほとんどが赤色である。Z = -0.7588からの比高によって色付けされているのである。オーバーハングの場なども正確に表現されている。

図22 3Dパース1
図23 3Dパース2

 CloudCompare octree8分木空間分割を最適化が参考になる)を表示する。DB Treeには,体積計算を実行すると,Height difference raster (0.008)のサブディレクトリにOctree(8分木空間分割)が作成される。これに✓して,PropertiesのOctree欄でDisplay levelを5にした場合が,図24の表示である。多数の平行6面体が見えるだろう。このレベルだと,Current level情報に,Cell size : 0.232482 m, Cell count: 1,822, Filled Volume: 22.8938 ㎥ が見える。

図24 Octree Display level =5, Cell size : 0.232482 m, Cell count: 1,822, Filled Volume: 22.8938㎤
図25 CloudCompare octree 掲載図

 Cloudcompareでは,この図25のウサギのような3Dの体積を求めることができる。

6 大饗石の体積を求める3: 点群欠損部の解消

 作業過程でこの地物にさらにフロートなどが残っていた。図22にその滓を示している。滓を排除したあと,自動的に作成されるファイル名は,segmented.segmented.setmentedと三つ並ぶ。これが現段階での地物の最新クラウドである。LiDARスキャンした際の不手際で,点群データの欠損が見られる。図21や図23の赤丸の部分である。図21の赤丸部分はトップ面に属し,図23の残りの二つの赤丸部分は側面に属している。

図22 滓が残っていた
図23 点群データの欠損の例

 この作業は,フィールドで確認作業をすれば不要だったかも知れない。1日目にこのLiDAR測量をしたので,当日のホテルで確認をすれば良かったとは思う。フィールドではiPhone 12 Proの画面での作業になるので見落とす可能性もある。図20左のパネル内に赤丸で示した赤字のメッセージは,クラウドに欠損があり解決すべきとしている。この欠損解決するため,iwamaによる iPhone LiDAR×CloudCompareで体積算出しようぜ!を参考にした。
 iwamaのサイトでは点群データをメッシュデータに変換して,点群の欠損部をなくした後に,このメッシュデータを点群データに戻している。点群データをメッシュデータに変換するだけで,何故,元の点群の欠損を消すことができるのか,それはメッシュデータは点群の欠損部を否応無くメッシュでカバーしてしまうからのようである。
 Hajime Saitoさんの メッシュのデータ構造 によれば(不要と思われる部分は削除),「メッシュは突き詰めれば点の集まりで、点と点を結んで線を作り、線と線を結んで面を作る。点は頂点、線はエッジ、面はフェイスと呼ばれる。」,とある。
 この文で,「点の集まり」のところが点群データなのだが,点群データからメッシュデータに変換する意義などが示された,点群データの3Dモデル化!点群をCADデータ・メッシュデータ・サーフェスデータに変換 ,が参考になる。「点群データは、3Dスキャンによって取得した段階の場合、単純に点の集合の状態でしかありません。そのためCADデータとして編集を行う際には、点群データを「点」から「面」の状態に変換していく必要があります。」,とある。点から面の状態に変換する場合,メッシュ(ポリゴン)データとサーフェス(ジオメトリ)データがあり,ここではメッシュデータに限定する。その説明の一部を次に引用する。

メッシュデータ

メッシュデータ(ポリゴンデータ)の場合は、点群のそれぞれの点を頂点としてとらえ、辺と面で接続した状態の3Dモデルの形状データを表していきます。ただ、要するに点同士をつなぎ合わせた状態の形状にあたるため、カクカクとした状態のぎこちない見た目になるのが特徴です。

 このメッシュデータの説明から,点群データの欠落部が埋められるのは容易に理解できる。
 それでは,iPhone LiDAR×CloudCompareで体積算出しようぜ!のプロセスを大饗石の点群データに適用してゆく。表現は替えている。

手順1 DB Treeから該当の点群データを選択し、Edit > Nomals > Computeで,特にパラメータには触らず,OKを選択すると,点群データに「表裏」ができる(図24)。
手順2 点群データをメッシュデータに変換するため、Edit > Mesh > Delaunay2.5D(XYPlane)で,特にパラメータには触らず,OKを選択する。ここで,図25のように,Keep old normals ? のアラートが現れるので,Yesを選択する。

図24 手順1の結果
図25 手順2でのアラート

 観察すると,大饗石のすべての表面には,図26のようにメッシュが貼られている。

手順3 メッシュデータを再度点群データに変換するため、ツールバーのPoints Sampling on mesh(図26に見える上段のツールバー群の左端(選択中なのでブルーになっている))を実行する。Points Number(デフォルト)を選択し,OK。この実行の結果,図27に施した赤丸のように, ⋯⋯⋯⋯⋯.mesh.sampled という点群データファイルが生まれているが,変換前のメッシュデータも合わせて表示されているので,点群データが見えない。図28では,変換前のファイルの前にある□で✓を外しているので,点群データだけが見えている。

図26 欠損部分もメッシュになっている
図27 点群データになった
図28 メッシュデータの変換から生まれた新たな点群データ

 ここまでの作業結果を確認したい。点群の欠損を果たして埋めることができたのか。図29はメッシュデータであるが,図23に見られる欠損は埋められている。
 図30を見ると,この変換過程に入る前のもともとの点群データと比べて,トップ面の点密度が小さくなったような印象を受ける。Propertiesで見ると,Points 1,000, 623であり,滓を取った点群データ(図23)Points 615,237と比べると,むしろ増えており,問題無いのだろう。

省略手順 iwamaのサイトでは,点群密度を間引く作業が用意されているが当方は不要なので省略する。

図29 メッシュデータでの点群データ
図30 新たな点群データの天井部分

同省略手順 iwamaのサイトでは,地物のベースの部分を切り取って,ground側のクラウドとする作業が実施されている。二つのクラウドの差を取るCloudCompareの手法を実現するためのなかなか良い工夫だと思うが,その必要性があるのか疑問なので,一応,この手順も省略したい。
 後に実行しようとしたが,segmentation 機能と

手順4 「5 大饗石の体積を求める2」と全く同じように計算した結果を図31に示す。Volume: 53.826,Surface 23.716,Matching cells: 59.2%である。Empty Cellsはinterpolateに替えたが,「5 大饗石の体積を求める2」でもinterpolateに替えたが改善は見られず図20の赤字でのエラーメッセージが出たのであるが,今回はこのメッセージは出ていない。

図31 体積計算

 図31を画面一杯に表示することもできる。Resultsも一気に見ることができる。Octreeの詳細も見える。図31と32で計算結果が異なるが,図31の出力の後,処理を修正したことが反映しているようだ。図33はHeight difference rasterの分布を見ている。図22〜24の赤色の3D分布と比べて,圧倒的に事実を反映していることがわかる。この章で実行した「点群ファイル→メッシュファイル→点群ファイル」のプロセスは体積計算には欠かせないものだろうと思う。

図32 図31の拡大版
図33 Heigt difference raster (0.08)のみ表示

 本章の作業結果と,前章の作業結果では,かなりの大きな違いがあり,本章の作業は体積計算に欠かせないものだと思う。

7 大饗石の体積を求める4: メッシュデータの面積と体積の計算

 メッシュデータの場合には,Edit > Mesh > Measure Surface,そして,Measure volume,のコマンドを使って,計算が可能で,結果は下部のConsoleに表示される(マニュアルp. 65)この場合の結果は,
時刻 [Mesh Surface] Mesh ファイル名: S = 1141.79 (square units)
時刻 [Mesh Surface] Average triangle surface: 0.000928008 (square units)

時刻 [Mesh Volume] Mesh ファイル名: V = 35.2555 (cube units)
時刻 [Mesh Volume] The above volume might be invalid (mesh has holes)

面積は何故か,非現実的な値となっている。1141.79 = 0.000928008 x 1230367で求めることができる。1230367は,このメッシュファイルのProperties > Mesh > Faces = 1,230,367に対応している。非常に詳細な凹凸を認識しての表面積であり,誤りでは無いが,一般に認知されている表面積ではない。一般に認知されているのはいわば表面のデコボコを無視して平滑な面で構成される地物の表面積なのだろう。

 体積も同様の計算過程を経たとしても点群ファイルとの差は大きくは無いはずで,V = 35 ㎥は適当と考えて良いのであるが,なんと,点群ファイルの欠損部分埋めた筈のメッシュファイルにもホールがあるという。というわけで,このメッセージ故に信頼できないとなるのである。

おわりに

 エラーを回避できた計算結果は, 第6章での53.8㎥ のみであった。この結果を受け容れることになる。ただ,CloudCompareの体積計算機能が適切なのか,簡単に計測が可能な直方体を通じて知る必要がある。次のポスティングに引き継ぎたいと思う。

身近な直方体の体積をiPhone 12 Pro+ LiDAR写真測量で求める

以上,Jan. 7, 2023記。