点群分割を実装してみる（点群表示機能を自作してみる その６・10億点モデル）

[「点群表示機能を自作してみる」](https://qiita.com/yumimatoba14/items/a89b01145b3af4ce9451)のその６です。10億点のモデルを描画してみました。
[その４](https://qiita.com/yumimatoba14/items/61ea891b303f2063d263#%E8%A9%A6%E4%BD%9C%E6%96%B9%E9%87%9D%E3%81%BE%E3%81%A8%E3%82%81) で案1とした、点群データを場所によって分割する機能を実装します。 木構造の構築はさぼりました。

# 出来上がり図

[GitHub](https://github.com/yumimatoba14/DirectXMfc.git) の以下のリビジョンを参照してください。 

> Revision: 3497790442d7a5842530c2e6eb9339d27a24ba36
> Message:
> modified test code to make huge data.

実行すると10億点のデータ(15GB超のファイル)を作り始めるので注意してください。 ローカルに出力ファイルを作るほか、作業用の一時ファイルも計30GB位、環境変数 TMP の下あたりに作られます。これもご注意ください。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3124939/d6c36094-a7e3-bac5-d118-6f1bf7cdbc9a.png)

一枚100万点の板を奥行方向に1m間隔で100枚、右手方向に10枚並べて合計10億点です。
データを作り終わり、表示してしまえば割とサクサク動くと思います。
ちょっとデータ作成部分が遅い他は、当初狙った程度のものはできたように思います。

2023/5/4 追記
100億点を試そうとしてモデルファイル作成ロジックの不具合を見つけ修正しました。
> Revision: b713d1ddd8d616dfbc5bd3525ec3510c769d2914
> Message:
> fixed some bugs
> - avoided an exception which had been thrown unnecessarily.
> - modified to suppress redrawing under modal dialog to avoid making files repeatedly.

100億点でも一応自分が意図した通りに動いてました。(SSD の容量が足りず HDD 上にファイルを構築しましたが。) ですがプログレッシブ表示のイマイチさが目立つようになってくるようです。  100億点の話はまたいずれ。

# 説明
## 概要
ファイルフォーマット的には前回作ったデータを複数個分を一つにまとめた、という以上には新しいことはありません。 (今回のデータでは1000個強あります。) 但し一つの点群を複数に分割するところに今回の実装の目玉があります。

コードをご覧になる場合は MultiPointListSampleModel2 がテストモデルクラス、構築コードは MultiPointListSampleModel2::PrepareFile() で、実体は D3DPointBlockListBuilder クラスになります。

描画に関しては2種類の描画データの簡略化が実装されています。これらは領域分割を実装したおかげで実現可能となったものです。
* 視界外の点群の描画を避けるコードを入れています。
* 視点から遠い点群を間引く処理が実装されています。 (遠くにある点群は小さく描画されるだけなので、間引いてよいことは直観的に分かると思います。)

併せて領域分割した点群のうち、近くのモノを優先して描画するようにしています。 これも直感的に理解可能かと思います。

## 領域分割
点群の場所による分割のコードは D3DPointBlockListBuilder::Build1Level() にあります。
この関数は格子の順にそこに入る点を並べたファイルを出力します。

1. 最初に点群の AABB を求める。
1. 1格子100万点程度になるように AABB の分割数を決定する。(但し100万点以下になる保証はなし。適当。)
1. 全点を操作し、各格子に点が何点入るかを数えます。 各格子の点の数が決まれば、出力ファイル中のどこの範囲に点を配置すればよいかが決まります。(同一格子中の点の順番は不定。）
1. 格子毎に点を並べた出力ファイルを作成します。

## LOD データの作成
領域分割の出力ファイルを入力として D3DExclusiveLodPointListCreator クラスで LOD データを構築します。 入力ファイルの途中から入力点群が始まることを除けば本質的な変更はありません。

## 領域外の点群の描画除外
PointListSampleModel.cpp の CalcBoxDistanceInProjection() に実装があります。
処理の方針は以下の通りです。
* AABB が視野外か否かを判定します。
* AABB の8個の頂点を投影後の座標値に変換して再度 AABB を構築し、可視範囲との干渉チェックをします。
    * スクリーンの X, Y については -1~1 の範囲が可視範囲です。
    * 奥行方向 Z については視点より手前のものは不可視です。奥行的に視点の先の部分を持つものは可視と判定します。 詳細は透視投影写像について理解することが必要です。

## 各サブ点群の描画精度
点の描画サイズを3D空間での長さで指定している場合、点を間引くと点の向こう側が透けて見えてしまいます。 ですが今回はスクリーン上で最低1ドットの大きさになるように描画しています。(即ち遠い点はスクリーン上でのサイズ指定になっている。)この場合、1ドット分の長さより近くにある点は重なりますので、点を間引いても描画結果が劣化しない(少なくとも向こう側が透けて見えることがない)ことが期待されます。

今回のコードでは AABB の一番近いZ値(0の場合あり)を基準に、必要な描画精度を求めます。

# その他の詳細
## 透視投影写像についてのメモ
透視投影については [XMMatrixPerspectiveFovRH()](https://learn.microsoft.com/ja-jp/windows/win32/api/directxmath/nf-directxmath-xmmatrixperspectivefovrh) が担っています。この関数の仕様を理解しないと、今回の AABB が領域外かの判定ができません。
[その１](https://qiita.com/yumimatoba14/items/e2eb77f82cac7a161c30#projectionmatrix--xmmatrixperspectivefovrh)でも簡単に触れていますが、もう少し詳しくメモします。

以下のように記号を定義します。
```math
P_{ph} = P_v M_p
```
$M_p$ : 透視投影を表す 4x4 行列。
$P_{ph} = (x_{ph}, y_{ph}, z_{ph}, w_{ph})$ : スクリーン座標系(同次座標系)の点
$P_v = (x_v, y_v, z_v, 1)$ : ビュー座標系の点

$P_v$ はビュー座標系の点ということで、以下の前提条件があります。
* 視点は原点に一致。
* 右手座標系を使用しているため、視線方向は$-z_{v}$ 方向。

ビュー座標系は向きと位置を調整したのみで、普通の直交座標系です。 それに対してスクリーン座標系は同次座標系になります。 同次座標系というのはx,y,zに加え4次元目の w(weight)を持つ座標系です。 同次座標系上の点である $P_{ph}$ を通常の3次元座標系 $P_p$ に変換するには以下のようにします。
```math
P_p = (x_p, y_p, z_p) = (x_{ph}/w_{ph}, y_{ph}/w_{ph}, z_{ph}/w_{ph})
```
同次座標系は四角錐形の視体積を表現するために使われているのだと思います。(ビュー座標系で視点から離れた遠くにある2点ほど、スクリーン座標系に変換した時の2点間の距離は小さくなる。)

もう少し詳細をメモします。以下はデバッガで XMMatrixPerspectiveFovRH() の値を見て解釈した内容です。 (どこかに定義がありそうなものですが、見つけられなかったので。)

* $x_p$ については -1～1 が可視範囲。 ($y_p$ も同様。)
* $z_p$ については、 $z_v= -near$ の時に 0、 $z_v = -far$ の時に 1 になるように変換される。
    * 従って $z_p$ については値が小さい方が手間に描画されることになります。 これが D3D11_DEPTH_STENCIL_DESC の DepthFunc で D3D11_COMPARISON_LESS が使われる理由かと思います。
    * 0～1 の値に定義されるのは Viewport の MinDepth、MaxDepth の値をそのように指定しているからかもしれません。(未確認。)
    * $z_{ph}$ は具体的には $z_{ph} = (z0 - z) * (-z1) / (z0-z1)$ です。(但し $z0 = -near$, $z1 = -far$。)

* $w_{ph} = -z_v$ となる。
    * 可視範囲では $z_v < 0$ ですから、 $w_{ph} > 0$ ということになります。 


同次座標系で考える際に w が0や負になることがあるのはおかしいと思いますが、今回は AABB の可視性を考えたいので考慮する必要があります。 ややアドホックですが以下のように可視性判定をしました。
* AABB が視点より $+z_v$ 側にある場合は不可視確定。 $-z_v$ 側にある場合は $w_{ph} > 0$ となるので普通に扱える。
* AABB が視点を包含する場合は、距離0 として可視とみなす。
* AABB が $z_v = 0$ をまたぐ場合、 $z_v=0$ 近傍では $x_p$, $y_p$ は必ず不可視範囲に飛んでいく。 可視性は $z_v << 0$ の個所で $x_p$, $y_p$ の値が可視範囲かどうかによる。
    * 計算上は、$z_v > 0$ となる点については $w_{ph}$ を $0 < w_{ph} << 1$ のようにして、無理やり有効な同次座標系上の点に $P_{ph}$ を移動させて判定しました。

# 参考資料
* [GetTempPath()](https://learn.microsoft.com/ja-jp/windows/win32/api/fileapi/nf-fileapi-gettemppatha) Microsoft
    * 一時ファイルの作成に使用しています。環境変数、TMP、TEMP などが参照されます。