地震被害の VR 映像生成のための三次元点群の平滑化
Smoothing of Point Cloud for Generating VR Video Image of Earthquake
Damage○山下 拓三*1
Takuzo Yamashita＊1*1 防災科学技術研究所地震減災実験研究部門

主任研究員

博士(工学)

Chief Researcher, Earthquake Disaster Mitigation Research Division, National Research Institute for Earth Science and Disaster
Resilience, Dr. Eng.
キーワード：三次元点群; バーチャルリアリティ; 最適化問題;
Keywords: Point cloud; virtual reality; optimization problem.

1. はじめに

なお，,は画像座標，下付き添字 i は点群のインデッ

防災科学技術研究所では，Ｅ－ディフェンス実験を活用

クスで，上付きの は計測データを表している．以降，x，した地震被害の仮想現実（VR）体験システム 1)の開発を進y，z 成分について同じ手順のため，x 成分を例に説明する．

めている．著者は，2015 年実施の 10 層 RC 建物実験 2)に3 次元点群の計測ノイズに起因する凹凸を低減するために，おいて，空間センサー（Kinect v2）を用いて三次元点群の計測値からの誤差が許容できる範囲

時刻歴データを取得し，ヘッドマウントディスプレイて，画像座標に対して 3 次元座標の各成分が滑らかに分布（HMD）で表示する高い臨場感のある室内映像の作成に取するように最適化問題を構成している．

で点群を動かしり組んでいる．リアルタイムで高品質の映像を HMD で表

添字 i は画像座標の 1 次元のインデックスであるが，画示するためには，三次元点群からメッシュを作成し，その像座標の縦方向のインデックス，横方向のインデックスをメッシュを粗化してデータを削減する方法が考えられる．

それぞれ ， として，必要に応じて 1 対 1 対応となる 2 種

しかし，計測ノイズの影響により得られた三次元点群から類のインデックスを使い分けて表記する．

メッシュを生成すると凹凸のある面が生成されて，十分な続いて，次式により偏微分項を差分化する．

映像の品質を確保できない．また，そのようなメッシュを

𝜕 2 𝑥(𝝃)≈ 𝑥(𝝃𝛼+1,𝛽 ) − 2𝑥(𝝃𝛼,𝛽 ) + 𝑥(𝝃𝛼−1,𝛽 )𝜕𝜉 2 |𝝃=𝝃𝛼,𝛽粗化しても品質は改善されない．
本報告では，3 次元点群を平滑化する最適化手法を提案≔ 𝑥𝛼+1,𝛽 − 2𝑥𝛼,𝛽 + 𝑥𝛼−1,𝛽 =: 𝑎𝛼,𝛽(4)し，10 層 RC 実験で取得した 3 次元点群に適用する．平滑化された点群に対して，メッシュを生成して，そのメッシュを粗化し，得られたメッシュを可視化して提案手法の有𝜕 2 𝑥(𝝃)≈ 𝑥(𝝃𝛼,𝛽+1 ) − 2𝑥(𝝃𝛼,𝛽 ) + 𝑥(𝝃𝛼,𝛽−1 )𝜕𝜂 2 |𝝃=𝝃𝛼,𝛽≔ 𝑥𝛼,𝛽+1 − 2𝑥𝛼,𝛽 + 𝑥𝛼,𝛽−1 =: 𝑏𝛼,𝛽効性を検証する．

(5)

なお，以降の式を見やすくするために係数 ， を導入している．差分化された目的関数は次式で表される．

3. 3 次元点群の平滑化手法
3.1. 最適化問題

𝐹 =3 次元座標の各成分 x，y，z に対して，次式の最適化問

𝑎+ 𝑏∑{| 𝛼,𝛽 | | 𝛼,𝛽 |}𝛼,𝛽(6)

題を構成する．
min𝒙

𝜕 2 𝑥(𝝃)𝜕 2 𝑥(𝝃)+𝑠. 𝑡. |𝑥(𝝃𝑖 ) − 𝑥̂(𝝃𝑖 )| ≤ 𝜖 𝑥 (𝝃𝑖 ) (1)∑ | 𝜕𝜉 2 || 𝜕𝜂 2 |𝝃=𝝃
𝝃=𝝃𝑖𝑖

3.2. 目的関数の近似
目的関数は原点で微分できない絶対値の関数を含んで𝑖いることから，解析的に設計変数を求めることができない．
𝜕 2 𝑦(𝝃)𝜕 2 𝑦(𝝃)min+𝑠. 𝑡. |𝑦(𝝃𝑖 ) − 𝑦̂(𝝃𝑖 )| ≤ 𝜖 𝑦 (𝝃𝑖 ) (2)|𝒚 ∑ | 𝜕𝜉 2 |𝜕𝜂 2 |𝝃=𝝃𝝃=𝝃
𝑖𝑖𝑖

𝜕 2 𝑧(𝝃)𝜕 2 𝑧(𝝃)min+𝑠. 𝑡. |𝑧(𝝃𝑖 ) − 𝑧̂(𝝃𝑖 )| ≤ 𝜖 𝑧 (𝝃𝑖 ) (3)|𝒛 ∑ | 𝜕𝜉 2 |𝜕𝜂 2 |𝝃=𝝃𝑖
𝝃=𝝃𝑖𝑖
日本建築学会情報システム技術委員会
第 42 回情報･システム･利用･技術シンポジウム論文集， 20-23 ， 2019 年 12 月，東京
Proceedings of the 42nd Symposium on Computer Technology of Information ，Systems and Applications ， AIJ ， 20-23 ， Dec. ， 2019 ， Tokyo

そこで絶対値の関数を次式のように 2 次関数で近似する．
2𝑎𝛼,𝛽𝑐𝛼,𝛽|𝑎𝛼,𝛽 | ≈ 2𝑐 + 2 for 𝑐𝛼,𝛽 > 0(7)𝛼,𝛽|𝑏𝛼,𝛽 | ≈𝑏2𝛼,𝛽2𝑑𝛼,𝛽+𝑑𝛼,𝛽2for 𝑑𝛼,𝛽 > 0(8)－ 20 －

絶対値の関数を 2 次関数で近似するため，反復計算によりと設定している．

近似最適化問題を解くことになる．そこで，n 回目の反復
3.4. 最適化問題の求解での値を上添字(n)で表すと，この時，
(𝑛)

(𝑛)

(𝑛)

(𝑛)

(𝑛)

𝑐𝛼,𝛽 ≔ |𝑎𝛼,𝛽 | = |𝑥𝛼+1,𝛽 − 2𝑥𝛼,𝛽 + 𝑥𝛼−1,𝛽 |

(9)

(𝑛)

(10)

(𝑛)

(𝑛)

(𝑛)

(𝑛)

𝑑𝛼,𝛽 ≔ |𝑏𝛼,𝛽 | = |𝑥𝛼,𝛽+1 − 2𝑥𝛼,𝛽 + 𝑥𝛼,𝛽−1 |とすると，絶対値の関数と近似に用いた 2 次関数の関係は
図 1 のようになる．

(𝑛+1)
目的関数𝐹 (𝑛+1) = ∑𝛼,𝛽 𝐹𝛼,𝛽は，式(11)より，𝒙(𝑛+1) についての 2 次形式の項と定数の項で表される．そのことから，(𝑛 + 1)回目の反復で解くべき最適化問題は以下の 2 次計画問題となる．
𝑇min {𝒙(𝑛+1) } 𝑷 (𝑛) 𝒙(𝑛+1) s. t 𝑮𝒙(𝑛+1) ≤ 𝒉(𝑛)
𝒙(𝑛)

(15)

(𝑛)

ここに，𝑷 (𝑛) は𝑐𝛼,𝛽 ，𝑑𝛼,𝛽 から計算でき，また，𝑮 = [𝑰−𝑰 ]𝑇𝒉(𝑛) = [𝒉(𝑛)+𝒉(𝑛)− ]𝑇(16)(𝑛)である．ただし，𝒉+ , 𝒉(𝑛)− の i 成分は次式で表される．
(𝑛)

(0)

ℎ+,𝑖 = 𝑥𝑖 + 𝜖 𝑥(𝑛) (𝝃𝑖 ),(𝑛)(0)ℎ−,𝑖 = 𝑥𝑖 − 𝜖 𝑥(𝑛) (𝝃𝑖 )(17)

2 次計画問題の求解には，Python の最適化ライブラリ
CVXOPT3)の関数を利用する．
Figure 1. Approximation of an absolute value function using aquadratic function.

4. 適用事例
本節では，10 層 RC 建物の加振実験の室内において空間以上より，n+1 回目の反復での点( , )についての目的関

(𝑛+1)

𝐹𝛼,𝛽=

(𝑛+1) 2{𝑎𝛼,𝛽 }(𝑛)(𝑛)+𝑐𝛼,𝛽22𝑐𝛼,𝛽+(𝑛+1) 2{𝑏𝛼,𝛽 }(𝑛)2𝑑𝛼,𝛽フレームのデータを検証用のデータとして用いる．検証用(𝑛)+𝑑𝛼,𝛽のデータは，RGB 画像データ（解像度：1920×1080），お2よび，深度画像データ（解像度：512×424）である．深度1⎡ 1(𝑛+1)(𝑛+1)(𝑛+1) 2(𝑛)= ⎢ (𝑛) (𝑥𝛼+1,𝛽 − 2𝑥𝛼,𝛽 + 𝑥𝛼−1,𝛽 ) + 𝑐𝛼,𝛽⎢2 𝑐⎣ 𝛼,𝛽(11)

(𝑛+1)
(𝑛+1)
(𝑛+1) 2(𝑛) ⎤+ (𝑛) (𝑥𝛼,𝛽+1 − 2𝑥𝛼,𝛽 + 𝑥𝛼,𝛽+1 ) + 𝑑𝛼,𝛽 ⎥⎥𝑑𝛼,𝛽⎦1,は，点群が動ける範囲を指定するものであり，対象とする点の周囲の点群の値をもとに次式により求める．
(12)

ここに，𝑛𝑠 は計算に用いる周辺の点群の個数を指定する変2の時 25 個の点群の値が用いられる．また，⎧⎪(𝑥𝑘,𝑙 − 𝑥𝛼,𝛽 )2⎪
𝑓𝑘,𝑙 ≔ ⎨⎪0⎪⎩

217,088 である．
作成したメッシュはオープンソースの可視化ソフトである MeshLab_64bit_fp v2016.12（以下，MeshLab）を用いて可チャとして用いた．生成されたメッシュでは境界部で品質の悪いメッシュが見られたため MeshLab の機能を用いて，メッシュの面積が Case1 では 5mm2 以下，Case2 では 3mm2以下の場合，メッシュのエッジの長さが 80mm 以上の場合，√√ 1 𝛼+𝑛𝑠 𝛽+𝑛𝑠𝑥𝑓𝜖 (𝝃𝛼,𝛽 ) = √𝑁 ∑ ∑ 𝑘,𝑙⎷ 𝑘=𝛼−𝑛𝑠 𝑙=𝛽−𝑛𝑠数で，データの解像度より，解くべき最適化問題の次元は視化した．可視化に際しては RGB 画像をメッシュのテクス3.3. 制約条件式の閾値の計算
制約条件式での閾値センサー（Kinect v2）で撮影された家具が揺れる前のフレーム(Case1)と家具が揺れている際のフレーム(Case2)の 2数の近似値は次式で表される．

𝑖𝑓
𝑖𝑓|𝑑𝑘,𝑙 − 𝑑𝛼,𝛽 |< 𝑡𝑜𝑙1𝑑𝛼,𝛽|𝑑𝑘,𝑙 − 𝑑𝛼,𝛽 |≥ 𝑡𝑜𝑙1𝑑𝛼,𝛽𝑑𝑘,𝑙 ≔ √𝑥2𝑘,𝑙 + 𝑦2𝑘,𝑙 + 𝑧2𝑘,𝑙または，接続関係にあるメッシュ群のメッシュの数が 200 個以下の場合は該当するメッシュを削除している．
メッシュの粗化には，MeshLab の Quadratic Edge CollapseDecimation(with texture)の機能を用いた．設定したパラメー
タは削除後の目標メッシュ数（Target number of faces）のみであり，その値を 30,000 に設定した．

(13)

作成したメッシュは，平滑化無し，粗化無しのメッシュ（non-smoothing fine mesh），平滑化無し，粗化有りのメッシ(14)である．式(12)は標準偏差に類似した式であるが，式(10)にュ（non-smoothing coarse mesh），平滑化有り，粗化無しのメッシュ（smoothing fine mesh），平滑化有り，粗化有りのメッシュ（smoothing coarse mesh）の 4 パターンである．

示すように点(𝑘, 𝑙)が対象の点(𝛼, 𝛽)と 3 次元空間上で距離表 1 にこれら 4 パターンのメッシュの節点とメッシュのが離れている場合は，その値は加算されない．ここで，𝑡𝑜𝑙1総数を示す．粗化についていずれのケースでも，元のメッはユーザーが指定する値であり，4 節の適用事例では 0.05

シュの総数の 1/10 以下である設定パラメータとほぼ同じ

報告 H09－ 21 －Table 1. Total number of vertices and faces.

も凹凸があることが確認できる．一方，平滑化することに(a) Case1

より凹凸が解消されていることが確認できる．粗化した場
合，家具の平らな面では平滑化したメッシュのほうがよりNon-smoothing

SmoothingFine

CoarseFine

Coarse

粗くメッシュ分割ができていることがわかる．図 5 にVertices

163,385

17,447

176,432

19,312

Case2 での 4 パターンの本棚の 3 次元画像の比較を示す．

Faces

313,047

29,999336,775

30,000加振時の家具が動いている場合でも，Case1 と同様の結果が得られている．Case2 のほうが平滑化の有無による見た(b) Case2目の違いはそれほど大きくないがこれは，空間センサーとNon-smoothing

SmoothingFine

CoarseFine

CoarseVertices

158,747

17,714

179,507

19,828Faces

303,380

29,999340,872

30,000のメッシュ数に削減できている．Case1 と Case2 での室内の RGB 画像と 3 次元画像をそれぞれ図 2，図 3 に示す．
ここで 3 次元画像として smoothing fine mesh の結果を示している．オクルージョンにより後ろの家具の映像は生成できていないが，実験では，複数台の空間センサーで時刻同期をして撮影していることから，これらのデータを統合した映像の作成が今後の課題となる．
図 4 に Case1 での 4 パターンの本棚の 3 次元画像の比較を示す．平滑化することにより，家具のエッジが滑らかに表現されて映像の見た目が良くなっていることがわかる．
また，粗化した場合，元の映像と見た目に大きな違いがな
く，平滑化したメッシュでは見た目が損なわれずデータ量を 1/10 以下に低減することができている．メッシュをみると平滑化しない場合，元のメッシュは，家具の平らな面でFig 2. RGB image (upper) and 3D image (lower) for Case1.

報告 H09

家具の距離が近く計測ノイズが小さいためと解釈できる．
5. おわりに本研究では，計測ノイズを有する三次元点群を平滑化する最適化手法を提案した．10 層 RC 建物のＥ－ディフェン
ス実験で取得された室内の 3 次元点群データに提案手法を適用して点群の平滑化を行い，その点群データからテクスチャ付きメッシュを作成して 3 次元映像を作成した．その結果，計測ノイズに由来する面の凹凸を解消した品質の高いメッシュを作成でき 3 次元映像の品質向上を確認した．
更に，そのメッシュをメッシュ数が 1/10 程度に削減されるように疎化しても映像の品質が変わらないことを確認した．
本手法は，三次元点群データを 3 次元画像計測に利用するためのノイズ除去手法としても応用の可能性がある．今後の研究として，振動台実験において本手法をノイズ除去に用いて空間センサーによる 3 次元変位計測を試み，モーションキャプチャ等の他の変位計測手法との比較により，応用の可能性を検証する予定である．

Fig 3. RGB image (upper) and 3D image(lower) for Case2.

－ 22 －

(a)

(b)

(c)

(d)

Figure 4. Shelves without meshes (upper pictures) and with meshes (lower pictures) for Case1. (a) nonsmoothing fine mesh; (b) non-smoothing coarse mesh; (c) smoothing fine mesh; (d) smoothing coarse mesh.

(a)

(b)

(c)

(d)

Figure 5. Shelves without meshes (upper pictures) and with meshes (lower pictures) for Case2. (a) nonsmoothing fine mesh; (b) non-smoothing coarse mesh; (c) smoothing fine mesh; (d) smoothing coarse mesh.
謝辞
本研究は JSPS 科研費 JP19H02286 の助成を受けたものです．
ここに記して謝意を表します．

［参考文献］
1） Yamashita, T. et al.: Development of a virtual reality experiencesystem for interior damage due to an earthquake: Utilizing E-報告 H09

Defense shake table test, {¥it Journal of disaster research}, Vol.12,
No.5, pp.882-890, 2017.
2） 土佐内優介，他：大型震動台による 10 階建て鉄筋コンクリート造建物の三次元振動実験（2015） ‐基礎すべり建物の実験方法と応答‐，日本建築学会構造系論文集，83 巻，750号，pp.1139-1149，2018.
3） http://cvxopt.org/, 2019 年 10 月 2 日アクセス．

－ 23 －