Mitsuba Renderer  0.5.0
sahkdtree4.h
Go to the documentation of this file.
1 /*
2  This file is part of Mitsuba, a physically based rendering system.
3
4  Copyright (c) 2007-2014 by Wenzel Jakob and others.
5
6  Mitsuba is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License Version 3
9
10  Mitsuba is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #if !defined(__SAH_KDTREE4_H)
20 #define __SAH_KDTREE4_H
21
23
25
27
28 /**
29  * \brief Implements the four-dimensional surface area heuristic for use
30  * by the \ref GenericKDTree construction algorithm.
31  */
33 public:
34  /**
35  * \brief Initialize the surface area heuristic with the bounds of
36  * a parent node
37  *
38  * Precomputes some information so that traversal probabilities
39  * of potential split planes can be evaluated efficiently
40  */
41  inline SurfaceAreaHeuristic4(const AABB4 &aabb) : m_aabb(aabb) {
42  const Vector4 extents(aabb.getExtents());
43  const Float temp = 1.0f / (extents.x * extents.y
44  + extents.y*extents.z + extents.x*extents.z);
45
46  m_temp0 = Vector4(
47  extents.y * extents.z * temp,
48  extents.x * extents.z * temp,
49  extents.x * extents.y * temp,
50  0.0f);
51
52  m_temp1 = Vector4(
53  (extents.y + extents.z) * temp,
54  (extents.x + extents.z) * temp,
55  (extents.x + extents.y) * temp,
56  1.0f / extents.w);
57  }
58
59  /**
60  * Given a split on axis \a axis that produces children having extents
61  * \a leftWidth and \a rightWidth along \a axis, compute the probability
62  * of traversing the left and right child during a typical query
63  * operation.
64  */
65  inline std::pair<Float, Float> operator()(int axis, Float leftWidth, Float rightWidth) const {
66  if (axis == 3 && m_temp1.w == std::numeric_limits<Float>::infinity()) {
67  return std::pair<Float, Float>(
68  std::numeric_limits<Float>::infinity(),
69  std::numeric_limits<Float>::infinity()
70  );
71  }
72
73  return std::pair<Float, Float>(
74  m_temp0[axis] + m_temp1[axis] * leftWidth,
75  m_temp0[axis] + m_temp1[axis] * rightWidth);
76  }
77
78  /**
79  * Compute the underlying quantity used by the tree construction
80  * heuristic. This is used to compute the final cost of a kd-tree.
81  */
82  inline static Float getQuantity(const AABB4 &aabb) {
83  const Vector4 extents(aabb.getExtents());
84  Float result = 2 * (extents[0] * extents[1] + extents[1] * extents[2]
85  + extents[2] * extents[0]);
86  if (extents[3] != 0)
87  result *= extents[3];
88  return result;
89  }
90 private:
91  Vector4 m_temp0, m_temp1;
92  AABB4 m_aabb;
93 };
94
95 /**
96  * This class specializes \ref GenericKDTree to a four-dimensional
97  * tree to be used for spacetime ray tracing. One additional function call
98  * must be implemented by subclasses:
99  *
100  * /// Check whether a primitive is intersected by the given ray.
101  * /// Some temporary space is supplied, which can be used to cache
102  * /// information about the intersection
103  * bool intersect(const Ray &ray, IndexType idx,
104  * Float mint, Float maxt, Float &t, void *tmp);
105  *
106  * This class implements an epsilon-free version of the optimized ray
107  * traversal algorithm (TA^B_{rec}), which is explained in Vlastimil
108  * Havran's PhD thesis "Heuristic Ray Shooting Algorithms".
109  *
110  * \author Wenzel Jakob
111  */
112 template <typename Derived>
113  class SAHKDTree4D : public GenericKDTree<AABB4, SurfaceAreaHeuristic4, Derived> {
114 public:
118
119 protected:
120  void buildInternal() {
121  SizeType primCount = this->cast()->getPrimitiveCount();
122  KDLog(EInfo, "Constructing a 4D SAH kd-tree (%i primitives) ..", primCount);
124  }
125
126  /**
127  * \brief Hashed mailbox implementation
128  */
129  struct HashedMailbox {
130  inline HashedMailbox() {
131  memset(entries, 0xFF, sizeof(IndexType)*MTS_KD_MAILBOX_SIZE);
132  }
133
134  inline void put(IndexType primIndex) {
135  entries[primIndex & MTS_KD_MAILBOX_MASK] = primIndex;
136  }
137
138  inline bool contains(IndexType primIndex) const {
139  return entries[primIndex & MTS_KD_MAILBOX_MASK] == primIndex;
140  }
141
143  };
144
145  /// Ray traversal stack entry for Havran-style incoherent ray tracing
147  /* Pointer to the far child */
148  const KDNode * __restrict node;
149  /* Distance traveled along the ray (entry or exit) */
151  /* Previous stack item */
153  /* Associated point */
155  };
156
157  /**
158  * \brief Ray tracing kd-tree traversal loop (Havran variant)
159  *
160  * This is generally the most robust and fastest traversal routine
161  * of the methods implemented in this class.
162  */
164  bool rayIntersectHavran(const Ray &ray, Float mint, Float maxt,
165  Float &t, void *temp) const {
167  #if 0
168  static const int prevAxisTable[] = { 2, 0, 1 };
169  static const int nextAxisTable[] = { 1, 2, 0 };
170  #endif
171
172  #if defined(MTS_KD_MAILBOX_ENABLED)
173  HashedMailbox mailbox;
174  #endif
175
176  /* Set up the entry point */
177  uint32_t enPt = 0;
178  stack[enPt].t = mint;
179  stack[enPt].p = ray(mint);
180
181  /* Set up the exit point */
182  uint32_t exPt = 1;
183  stack[exPt].t = maxt;
184  stack[exPt].p = ray(maxt);
185  stack[exPt].node = NULL;
186
187  bool foundIntersection = false;
188  const KDNode * __restrict currNode = this->m_nodes;
189  while (currNode != NULL) {
190  while (EXPECT_TAKEN(!currNode->isLeaf())) {
191  const Float splitVal = (Float) currNode->getSplit();
192  const int axis = currNode->getAxis();
193  const KDNode * __restrict farChild;
194
195  if (axis == 3) {
196  if (ray.time <= splitVal)
197  currNode = currNode->getLeft();
198  else
199  currNode = currNode->getRight();
200  continue;
201  } else if (stack[enPt].p[axis] <= splitVal) {
202  if (stack[exPt].p[axis] <= splitVal) {
203  /* Cases N1, N2, N3, P5, Z2 and Z3 (see thesis) */
204  currNode = currNode->getLeft();
205  continue;
206  }
207
208  /* Typo in Havran's thesis:
209  (it specifies "stack[exPt].p == splitVal", which
210  is clearly incorrect) */
211  if (stack[enPt].p[axis] == splitVal) {
212  /* Case Z1 */
213  currNode = currNode->getRight();
214  continue;
215  }
216
217  /* Case N4 */
218  currNode = currNode->getLeft();
219  farChild = currNode + 1; // getRight()
220  } else { /* stack[enPt].p[axis] > splitVal */
221  if (splitVal < stack[exPt].p[axis]) {
222  /* Cases P1, P2, P3 and N5 */
223  currNode = currNode->getRight();
224  continue;
225  }
226  /* Case P4 */
227  farChild = currNode->getLeft();
228  currNode = farChild + 1; // getRight()
229  }
230
231  /* Cases P4 and N4 -- calculate the distance to the split plane */
232  Float distToSplit = (splitVal - ray.o[axis]) * ray.dRcp[axis];
233
234  /* Set up a new exit point */
235  const uint32_t tmp = exPt++;
236  if (exPt == enPt) /* Do not overwrite the entry point */
237  ++exPt;
238
239  KDAssert(exPt < MTS_KD_MAXDEPTH);
240  stack[exPt].prev = tmp;
241  stack[exPt].t = distToSplit;
242  stack[exPt].node = farChild;
243
244  #if 1
245  /* Intrestingly, this appears to be faster than the
246  original code with the prevAxis & nextAxis table */
247  stack[exPt].p = ray(distToSplit);
248  stack[exPt].p[axis] = splitVal;
249  #else
250  const int nextAxis = nextAxisTable[axis];
251  const int prevAxis = prevAxisTable[axis];
252  stack[exPt].p[axis] = splitVal;
253  stack[exPt].p[nextAxis] = ray.o[nextAxis]
254  + distToSplit*ray.d[nextAxis];
255  stack[exPt].p[prevAxis] = ray.o[prevAxis]
256  + distToSplit*ray.d[prevAxis];
257  #endif
258
259  }
260
261  /* Reached a leaf node */
262  for (IndexType entry=currNode->getPrimStart(),
263  last = currNode->getPrimEnd(); entry != last; entry++) {
264  const IndexType primIdx = this->m_indices[entry];
265
266  #if defined(MTS_KD_MAILBOX_ENABLED)
267  if (mailbox.contains(primIdx))
268  continue;
269  #endif
270
271  bool result;
273  result = this->cast()->intersect(ray, primIdx, mint, maxt, t, temp);
274  else
275  result = this->cast()->intersect(ray, primIdx, mint, maxt);
276
277  if (result) {
279  return true;
280  maxt = t;
281  foundIntersection = true;
282  }
283
284  #if defined(MTS_KD_MAILBOX_ENABLED)
285  mailbox.put(primIdx);
286  #endif
287  }
288
289  if (stack[exPt].t > maxt)
290  break;
291
292  /* Pop from the stack and advance to the next node on the interval */
293  enPt = exPt;
294  currNode = stack[exPt].node;
295  exPt = stack[enPt].prev;
296  }
297
298  return foundIntersection;
299  }
300 };
301
303
304 #endif /* __SAH_KDTREE4_H */
void put(IndexType primIndex)
Definition: sahkdtree4.h:134
FINLINE bool rayIntersectHavran(const Ray &ray, Float mint, Float maxt, Float &t, void *temp) const
Ray tracing kd-tree traversal loop (Havran variant)
Definition: sahkdtree4.h:164
#define KDLog(level, fmt,...)
Definition: gkdtree.h:642
KDTreeBase< AABB4 >::KDNode KDNode
Definition: sahkdtree4.h:117
std::pair< Float, Float > operator()(int axis, Float leftWidth, Float rightWidth) const
Definition: sahkdtree4.h:65
bool contains(IndexType primIndex) const
Definition: sahkdtree4.h:138
const KDNode *__restrict node
Definition: sahkdtree4.h:148
KDTreeBase< AABB4 >::SizeType SizeType
Definition: sahkdtree4.h:115
static Float getQuantity(const AABB4 &aabb)
Definition: sahkdtree4.h:82
#define KDAssert(expr)
Definition: gkdtree.h:56
void buildInternal()
Definition: sahkdtree4.h:120
Base class of all kd-trees.
Definition: gkdtree.h:441
Float t
Definition: sahkdtree4.h:150
TAABB< Point4 > AABB4
Definition: fwd.h:157
#define MTS_NAMESPACE_BEGIN
Definition: platform.h:137
#define MTS_KD_MAXDEPTH
Activate lots of extra checks.
Definition: gkdtree.h:37
More relevant debug / information message.
Definition: formatter.h:31
TVector4< Float > Vector4
Definition: fwd.h:122
HashedMailbox()
Definition: sahkdtree4.h:130
Optimized KD-tree acceleration data structure for n-dimensional (n&lt;=4) shapes and various queries on ...
Definition: gkdtree.h:706
Definition: fwd.h:65
Point p
Definition: sahkdtree4.h:154
VectorType getExtents() const
Calculate the bounding box extents.
Definition: aabb.h:292
Implements the four-dimensional surface area heuristic for use by the GenericKDTree construction algo...
Definition: sahkdtree4.h:32
#define MTS_KD_MAILBOX_SIZE
Definition: sahkdtree3.h:31
Hashed mailbox implementation.
Definition: sahkdtree4.h:129
Definition: fwd.h:100
Ray traversal stack entry for Havran-style incoherent ray tracing.
Definition: sahkdtree4.h:146