フォトンマップ（kd木）を実装した。
動作確認として、カメラから画面中心に向かってレイを飛ばし、物体とぶつかった地点pからいちばん近いフォトンTOP10をkd木を使って検索し、pから各フォトンの位置に線を引くようにしてみた。
線の色はフォトンの色に対応。
いつもの如くリアルタイムでグリグリ動くけど、画面録画がキビしいので静止画だけ。
前回と比べ見栄え的には後退したけれど、実装的には大きな一歩となった。

これでレイトレースして放射輝度推定を行うための材料はほぼ揃ったかたちだが、前回も書いた通り、先にベイクの仕組みを考えたい。
その結果によってレイトレース処理を実装するか、もしくはいきなりベイク処理を実装するか決めようと思う。

kd木はJensen本を参考に実装した。
ただ、巻末に載ってたC++のソースはほとんど見ずに自分でいちから作成した。

あのソース、書き方がかなり古臭い（書かれた時代的に当たり前なんだけど、さらに敢えてCっぽく書いてある気がする）し、配列の最初のインデクスに1を使ってる（0番目は多分無視されてる）のが気に入らなかった。
STLを使えばかなりシンプルに書けるんじゃないかとも思ったし。

でまぁやってみたらかなり苦労したw

苦労したポイントを挙げると以下のような感じ。

1. kd木を二分ヒープ的1次元配列でどう表現するか
2. ソートされたポインタの配列を元にした、実体の配列のソート
3. n近傍探索のためのヒープの扱い（これは単なるうっかりミス）

kd木を二分ヒープ的1次元配列でどう表現するか

以下のような構造体を用いれば普通のkd木の構築処理は難しくない。

struct Node {
  Photon *photon;
  unsigned char plane;
  Node *left_child;
  Node *right_child;
};

しかし、例えばフォトン数が1000万個あるとき、左右の子を指すためのポインタだけで160MBもメモリを消費してしまう。（64bit環境の場合。8バイトx2x1000万 = 160MB）
さらに言うと、自分の場合すでにPhotonの実体は配列で確保されてるのでさらに*photon分のメモリも無駄に保持し続けなければならない。
（plane は、左右の子をどの平面で分割するかを表すidのようなもの。yz平面なら0, xz平面なら1, xy平面なら2とかにする。）

実はこの左右の子を指すポインタ（とPhotonの実体を指すポインタ）分のメモリ消費を0に出来る凄い方法がある。
kd木を二分ヒープの条件（葉ノードは隙間なく左詰めされなければいけない）を満たすように構築してあげれば、kd木を1次元配列で表現することができ、そうすることによって、あるノードのインデクス i から、その左右の子ノードのインデクスを i * 2 + 1, i * 2 + 2 という計算だけで求められるようになるのである。
(plane の情報はPhotonのデータ構造の中に移してしまう。)

こうやって書くと簡単なように思えるけど、kd木の構築ではある中央値を元に再帰的に分割していくわけで、この中央値を単純に (minIndex + maxIndex) / 2 としても二分ヒープの条件を満たすものは出来上がらない。
具体的には以下の図のように、全体的にバランスされた木にはなるけど葉ノードに隙間が出来てしまう。

# 10個の要素を単純なメディアンで再帰的に分割した場合の図（xが隙間、数値はこの木を1次元配列化したときのインデックス）
               0


       1               2

   3       4       5       6
 x   7   x   8   x   9

# 正しくは以下のとおり（最下段のみ表示）
 7   8   9

こうなるとこれを1次元配列化しても隙間付近で i * 2 + 1, i * 2 + 2 が成り立たなくなる。
これを無理やり1次元配列化して扱うと、初期化されない要素が出てきたりしてしまう。
（上の例だと、7で初期化されるべき所が隙間になってるし、8で初期化されるべきところが7で初期化されている。9で初期化されるべき所も隙間になってるし、そもそも8と9は配列外になる。）

これを回避するために、葉ノードのアンバランスさを先読みしてメディアンを決めていく方法を採用した。
で、その方法を言葉で説明するのは大変なんだけど、実装を考えるときに書いた擬似コードがあるのでそれを載せとく。
（入力は配列のどの部分を処理対象にするかを表すmaxIndexとminIndexのみ。）

// 二分ヒープにするためのメディアンの求め方
	size = maxIndex – minIndex + 1
	// 末端の葉ノード数を求める
	nodeCount = 1
	while nodeCount * 2 + 1 < size
		nodeCount = nodeCount * 2 + 1
	leafCount = size – nodeCount
	// 理想状態における葉ノード数を求める
	idealLeafCount = (nodeCount * 2 + 1) – nodeCount
	// 理想状態から見た場合の左右の葉ノード数
	if leafCount > idealLeafCount / 2
		leftLeafCount = idealLeafCount / 2
		rightLeafCount = leafCount – leftLeafCount
	else
		leftLeafCount = leafCount
		rightLeafCount = 0
	// メディアン（インデックス）の算出
	medianIndex = minIndex + (maxIndex – minIndex + (leftLeafCount – rightLeafCount)) / 2

感覚的に説明すれば、「葉ノードを左詰めするために、左側の子ノードに多めに割り振る。そのためにちょっと右よりにメディアンを選択する。」という感じ。

これで二分ヒープの条件を満たすkd木が1次元配列上に構築できた。

ソートされたポインタの配列を元にした、実体の配列のソート

これは大した話じゃないんだけど、高速化のために一時的にPhotonのポインタの配列を用意して、そのポインタ配列に対してkd木の構築処理を行うようにしたもんだから、後でポインタ配列の並びを実体の配列に反映させなきゃいけなくて、その方法の実装に苦労した。
まぁ省メモリを考えなければポインタ配列を保持したままでもいいんだけど、今回はそういうわけには行かなかったので。

最初考えたのが、どの実体がどのポインタから参照されてるかを表す、インデックスの配列（逆参照配列）を一時的に作って、それを参照しながらstd::sortで実体のソートを行う方法。

しかしstd::sortのカスタムプレディケート（lambda使用）の中では、処理対象のインデックスを知る方法がない。
そこで考えたのが、配列の先頭アドレスとプレディケートに引数として渡されてくるオブジェクトとのポインタ演算でインデックスを求める方法。

しかし実際にやってみると、2つの引数のうち片方は変なアドレスが入ってくる。
std::sort内部の一時変数のアドレスだろうなと思ったけど正しいところは調べてない。

そこでstd::sortを諦めて、qsort_bを試してみた。
どうしてもインプレースでプレディケートを書きたいのだw

しかしブロックの中では、その外のオブジェクトへの書き込みができなくて撃沈。
（実体の並び替えだけでなくインデックスの配列も同時に並び替えないといけないが、書き込みできないから並び替えもできない。）

結局インプレースでプレディケートを書くのは諦めて、独自にクイックソートを行う関数を実装した。
基本的すぎて自分で書くなんて事はめったになくて、実装時のうっかりミスで微妙な問題が発生してデバッグに手間取ったというオチ。
うまい具合にqsortを使えば、クイックソート本体は実装せずに済んだかも。

n近傍探索のためのヒープの扱い

上にも書いた通り、これもうっかりミス。
ただし、クイックソートの話のような前フリは無くて、最初から自信満々で実装したものにうっかりミスが紛れ込んでて、ちょっとデバッグに苦労したという話。

ここでいうヒープはkd木の二分ヒープとは別物、というか構造は一緒だけど、こちらのヒープはkd木によるn個の検索結果を保持するためのもの。
一番遠いフォトンが常にルートに来るようにしとくことで、それより近いフォトンが見つかったときに即破棄出来るようにするために使っている。

ヒープの要素を並び替える再帰関数の中のif文の条件が間違ってて、ヒープの配列境界を飛び越えて書き込みをしたせいで、タイミングによって多種多様な現象が引き起こされるという、メモリ破壊の典型をやらかした。

間違ってた部分に絶対の自信を持ってたのでデバッグに時間が掛かったが、ふと思い立って Xcode のスキーマの設定で Enable Guard Edges をオンにして実行したら間違ってるところ（の近く）をバッチリ教えてくれた。
これ開発中は常にオンにしとく勢いでいいなぁ。