Surface Area Heuristic (SAH) kd-tree of triangles - flat cells
<p>I've implemented a SAH kd-tree based upon the paper <a href="http://dcgi.felk.cvut.cz/home/havran/ARTICLES/ingo06rtKdtree.pdf" rel="noreferrer" title="On building fast kd-Trees for Ray Tracing, and on doing that in O&#40;N log N&#41;">On building fast kd-Trees for Ray Tracing, and on doing that in O(N log N)</a> by Wald and Havran. Note I haven't done their splicing and merging suggestion right at the end to speed up the building of the tree, just the SAH part.</p>

<p>I'm testing the algorithm with an axis-aligned cube where each face is split into two triangles, so I have <code>N = 12</code> triangles in total. The bottom-left vertex (i.e. the one nearest the axes) is actually at the origin.</p>

<p><a href="https://i.stack.imgur.com/hPozx.jpg" rel="noreferrer"><img src="https://i.stack.imgur.com/hPozx.jpg" alt="cube of 5 units with each face split into 2 triangles"></a></p>

<pre><code>Face    Triangle
----------------
Front:  0, 1
Left:   6, 7
Right:  2, 3
Top:    4, 5
Bottom: 10, 11
Back:   8, 9
</code></pre>

<p>Assuming node traversal cost <code>Ct = 1.0</code> and intersection cost <code>Ci = 10.0</code>. We first find the cost of not splitting which is <code>Cns = N * Ci = 12 * 10.0 = 120.0</code>.</p>

<p>Now we go through each axis in turn and do an incremental sweep through the candidate split planes to see if the cost of splitting is cheaper. The first split plane is <code>p = &lt;x,0&gt;</code>. We have <code>Nl = 0</code>, <code>Np = 2</code> and <code>Nr = 10</code> (that is, the number of triangles on the left, in the plane, and to the right of the plane). The two triangles in the plane are number 6 and 7 (the left face). All others are to the right.</p>

<p>The <code>SAH(p,V,Nl,Nr,Np)</code> function is now executed. This takes the split plane, the voxel <code>V</code> to be split, and the numbers of left/right/plane triangles. It computes the probability of hitting the left (flat) voxel as <code>Pl = SA(Vl)/SA(V) = 50/150 = 1/3</code>, the right probability as <code>Pr = SA(Vr)/SA(V) = 150/150 = 1</code>. Now it runs the cost function twice; first with the planar triangles on the left, then with the planar triangles on the right to get <code>Cl</code> and <code>Cr</code> respectively.</p>

<p>The cost function <code>C(Pl,Pr,Nl,Nr)</code> returns <code>bias * (Ct + Ci(Pl * Nl + Pr * Nr))</code></p>

<h2><code>Cl</code>: cost with planar triangles on the left (<code>Nl = 2</code>, <code>Nr = 10</code>)</h2>

<p><code>bias = 1</code> we aren't biasing as neither left nor right voxel is empty.<br/>
<code>Cl = 1 * (1 + 10(1/3 * 2 + 1 * 10)) = 107.666</code></p>

<h2><code>Cr</code>: cost with planar triangles on the right (<code>Nl = 0</code>, <code>Nr = 12</code>)</h2>

<p><code>bias = 0.8</code> empty cell bias comes into play.<br/>
<code>Cr = 0.8 * (1 + 10(1/3 * 0 + 1 * 12)) = 96.8</code></p>

<p>The algorithm determines that <code>Cr = 96.8</code> is better than splitting the two triangles off into a flat cell <code>Cl = 107.666</code> and also better than not splitting the voxel at all <code>Cns = 120.0</code>. No other candidate splits are found to be cheaper. We therefore split into an empty left child, and a <strong>right child containing all the triangles</strong>. When we recurse into the right child to continue the tree building, we perform exactly the same steps as above. It's only because of a max depth termination criterion that this doesn't cause a stack overflow. The resultant tree is very deep.</p>

<p>The paper claims to have thought about this sort of thing:</p>

<blockquote>
  <p>During clipping, special care has to be taken to correctly handle
  special cases like “flat” (i.e., zero-volume) cells, or cases where
  numerical inaccuracies may occur (e.g., for cells that are very thin
  compared to the size of the triangle). For example, we must make sure
  not to “clip away” triangles lying in a flat cell. Note that such
  cases are not rare exceptions, but are in fact encouraged by the SAH,
  as they often produce minimal expected cost.</p>
  
  <p>This local approximation can easily get stuck in a local minimum: As
  the local greedy SAH overestimates CV (p), it might stop subdivision
  even if the correct cost would have indicated further subdivision. In
  particular, the local approximation can lead to premature termination
  for voxels that require splitting off flat cells on the sides: many
  scenes (in particular, architectural ones) contain geometry in the
  form of axis-aligned boxes (a light fixture, a table leg or table top,
  . . . ), in which case the sides have to be “shaved off” until the
  empty interior is exposed. For wrongly chosen parameters, or when
  using cost functions different from the ones we use (in particular,
  ones in which a constant cost is added to the leaf estimate), the
  recursion can terminated prematurely. Though this pre-mature exit
  could also be avoided in a hardcoded way—e.g., only performing the
  automatic termination test for non-flat cells—we propose to follow our
  formulas, in which case no premature exit will happen.</p>
</blockquote>

<p>Am I doing anything wrong? How can I correct this?</p>