Photon Engine 2.0.0-beta
A physically based renderer.
Loading...
Searching...
No Matches
TIndexedKdtree.ipp
Go to the documentation of this file.
1#pragma once
2
3#include "Math/Algorithm/IndexedKdtree/TIndexedKdtree.h"
4#include "Math/Algorithm/acceleration_structure_basics.h"
5
6#include <limits>
7#include <array>
8#include <cmath>
9#include <algorithm>
10#include <optional>
11
12namespace ph::math
13{
14
15template<
16 typename IndexToItem,
17 typename ItemToAABB,
18 typename Index>
19inline TIndexedKdtree<IndexToItem, ItemToAABB, Index>
20::TIndexedKdtree(
21 const std::size_t numItems,
22 IndexToItem indexToItem,
23 ItemToAABB itemToAABB,
24 IndexedKdtreeParams params) :
25
26 m_numItems (numItems),
27 m_indexToItem(std::move(indexToItem)),
28 m_rootAABB (),
29 m_nodeBuffer (),
30 m_numNodes (0),
31 m_itemIndices()
32{
33 if(numItems > 0)
34 {
35 build(std::move(itemToAABB), params);
36 }
37}
38
39template<
40 typename IndexToItem,
41 typename ItemToAABB,
42 typename Index>
43template<
44 typename TesterFunc>
45inline auto TIndexedKdtree<IndexToItem, ItemToAABB, Index>
46::nearestTraversal(const TLineSegment<real>& segment, TesterFunc&& intersectionTester) const
47-> bool
48{
49 static_assert(CItemSegmentIntersectionTester<TesterFunc, Item>);
50
51 PH_ASSERT_GT(m_numNodes, 0);
52
53 struct NodeState
54 {
55 const Node* node;
56 real minT;
57 real maxT;
58 };
59
60 constexpr bool TEST_WITH_ITEM_IDX = CItemSegmentIntersectionTesterWithIndex<TesterFunc, Item>;
61 constexpr int MAX_STACK_HEIGHT = 64;
62
63 real minT, maxT;
64 if(!m_rootAABB.isIntersectingVolume(segment, &minT, &maxT))
65 {
66 return false;
67 }
68 TLineSegment<real> intersectSegment(segment.getOrigin(), segment.getDir(), minT, maxT);
69
70 const Vector3R rcpRayDir(segment.getDir().rcp());
71
72 std::array<NodeState, MAX_STACK_HEIGHT> nodeStack;
73 int stackHeight = 0;
74 bool hasHit = false;
75 const Node* currentNode = &(m_nodeBuffer[0]);
76 while(true)
77 {
78 if(!currentNode->isLeaf())
79 {
80 const auto splitAxis = currentNode->getSplitAxis();
81 const real splitPlaneT = (currentNode->getSplitPos() - segment.getOrigin()[splitAxis]) * rcpRayDir[splitAxis];
82
83 const Node* nearHitNode;
84 const Node* farHitNode;
85 if(
86 (segment.getOrigin()[splitAxis] < currentNode->getSplitPos()) ||
87 (segment.getOrigin()[splitAxis] == currentNode->getSplitPos() &&
88 segment.getDir()[splitAxis] <= 0))
89 {
90 nearHitNode = currentNode + 1;
91 farHitNode = &(m_nodeBuffer[currentNode->getPositiveChildIndex()]);
92 }
93 else
94 {
95 nearHitNode = &(m_nodeBuffer[currentNode->getPositiveChildIndex()]);
96 farHitNode = currentNode + 1;
97 }
98
99 // Case I: Only near node is hit.
100 // (split plane is beyond t-max or behind ray origin)
101 if(splitPlaneT > maxT || splitPlaneT < 0)
102 {
103 currentNode = nearHitNode;
104 }
105 // Case II: Only far node is hit.
106 // (split plane is between ray origin and t-min)
107 else if(splitPlaneT < minT)
108 {
109 currentNode = farHitNode;
110 }
111 // Case III: Both near and far nodes are hit.
112 // (split plane is within [t-min, t-max])
113 else
114 {
115 PH_ASSERT_LT(stackHeight, MAX_STACK_HEIGHT);
116
117 nodeStack[stackHeight].node = farHitNode;
118 nodeStack[stackHeight].minT = splitPlaneT;
119 nodeStack[stackHeight].maxT = maxT;
120 ++stackHeight;
121
122 currentNode = nearHitNode;
123 maxT = splitPlaneT;
124 }
125 }
126 // current node is leaf
127 else
128 {
129 const std::size_t numItems = currentNode->numItems();
130
131 if(numItems == 1)
132 {
133 const auto itemIndex = currentNode->getSingleItemDirectIndex();
134 const Item& item = getItem(itemIndex);
135
136 std::optional<real> hitT;
137 if constexpr(TEST_WITH_ITEM_IDX)
138 {
139 hitT = intersectionTester(item, intersectSegment, itemIndex);
140 }
141 else
142 {
143 hitT = intersectionTester(item, intersectSegment);
144 }
145
146 if(hitT)
147 {
148 intersectSegment.setMaxT(*hitT);
149 hasHit = true;
150 }
151 }
152 else
153 {
154 for(std::size_t i = 0; i < numItems; ++i)
155 {
156 const Index itemIndex = m_itemIndices[currentNode->getIndexBufferOffset() + i];
157 const Item& item = getItem(itemIndex);
158
159 std::optional<real> hitT;
160 if constexpr(TEST_WITH_ITEM_IDX)
161 {
162 hitT = intersectionTester(item, intersectSegment, itemIndex);
163 }
164 else
165 {
166 hitT = intersectionTester(item, intersectSegment);
167 }
168
169 if(hitT)
170 {
171 intersectSegment.setMaxT(*hitT);
172 hasHit = true;
173 }
174 }
175 }
176
177 if(stackHeight > 0)
178 {
179 --stackHeight;
180 currentNode = nodeStack[stackHeight].node;
181 minT = nodeStack[stackHeight].minT;
182 maxT = nodeStack[stackHeight].maxT;
183
184 // Early-out if the test segment cannot reach the next node
185 if(intersectSegment.getMaxT() < minT)
186 {
187 break;
188 }
189 }
190 else
191 {
192 break;
193 }
194 }// end is leaf node
195 }// end infinite loop
196
197 return hasHit;
198}
199
200template<
201 typename IndexToItem,
202 typename ItemToAABB,
203 typename Index>
204inline auto TIndexedKdtree<IndexToItem, ItemToAABB, Index>
205::getAABB() const
206-> AABB3D
207{
208 PH_ASSERT(!isEmpty());
209
210 return m_rootAABB;
211}
212
213template<
214 typename IndexToItem,
215 typename ItemToAABB,
216 typename Index>
217inline auto TIndexedKdtree<IndexToItem, ItemToAABB, Index>
218::isEmpty() const
219-> bool
220{
221 return m_numItems == 0;
222}
223
224template<
225 typename IndexToItem,
226 typename ItemToAABB,
227 typename Index>
228inline auto TIndexedKdtree<IndexToItem, ItemToAABB, Index>
229::getItem(const std::size_t idx) const
230-> Item
231{
232 PH_ASSERT_LT(idx, m_numItems);
233
234 return m_indexToItem(lossless_cast<Index>(idx));
235}
236
237template<
238 typename IndexToItem,
239 typename ItemToAABB,
240 typename Index>
241inline void TIndexedKdtree<IndexToItem, ItemToAABB, Index>
242::build(ItemToAABB itemToAABB, IndexedKdtreeParams params)
243{
244 PH_ASSERT_GT(m_numItems, 0);
245
246 const auto maxNodeDepth = static_cast<std::size_t>(8 + 1.3 * std::log2(m_numItems) + 0.5);
247
248 // Cache item AABB and calculate root AABB
249 std::vector<AABB3D> itemAABBs;
250 m_rootAABB = itemToAABB(getItem(0));
251 for(std::size_t i = 0; i < m_numItems; ++i)
252 {
253 const AABB3D aabb = itemToAABB(getItem(i));
254
255 itemAABBs.push_back(aabb);
256 m_rootAABB.unionWith(aabb);
257 }
258
259 std::unique_ptr<Index[]> negativeItemIndicesCache(new Index[m_numItems]);
260 std::unique_ptr<Index[]> positiveItemIndicesCache(new Index[m_numItems * maxNodeDepth]);
261
262 PH_ASSERT_LE(m_numItems - 1, std::numeric_limits<Index>::max());
263 for(std::size_t i = 0; i < m_numItems; ++i)
264 {
265 negativeItemIndicesCache[i] = lossless_cast<Index>(i);
266 }
267
268 std::array<std::unique_ptr<ItemEndpoint[]>, 3> endPointsCache;
269 for(auto&& cache : endPointsCache)
270 {
271 cache = std::unique_ptr<ItemEndpoint[]>(new ItemEndpoint[m_numItems * 2]);
272 }
273
274 buildNodeRecursive(
275 0,
276 m_rootAABB,
277 negativeItemIndicesCache.get(),
278 m_numItems,
279 0,
280 0,
281 itemAABBs,
282 maxNodeDepth,
283 params,
284 negativeItemIndicesCache.get(),
285 positiveItemIndicesCache.get(),
286 endPointsCache);
287}
288
289template<
290 typename IndexToItem,
291 typename ItemToAABB,
292 typename Index>
293inline void TIndexedKdtree<IndexToItem, ItemToAABB, Index>
294::buildNodeRecursive(
295 std::size_t nodeIndex,
296 const AABB3D& nodeAABB,
297 const Index* nodeItemIndices,
298 std::size_t numNodeItems,
299 std::size_t currentNodeDepth,
300 std::size_t currentBadRefines,
301 const std::vector<AABB3D>& itemAABBs,
302 const std::size_t maxNodeDepth,
303 IndexedKdtreeParams params,
304 Index* negativeItemIndicesCache,
305 Index* positiveItemIndicesCache,
306 std::array<std::unique_ptr<ItemEndpoint[]>, 3>& endpointsCache)
307{
308 ++m_numNodes;
309 if(m_numNodes > m_nodeBuffer.size())
310 {
311 m_nodeBuffer.resize(m_numNodes * 2);
312 }
313 PH_ASSERT_LT(nodeIndex, m_nodeBuffer.size());
314
315 if(currentNodeDepth == maxNodeDepth || numNodeItems <= params.getMaxNodeItems())
316 {
317 m_nodeBuffer[nodeIndex] = Node::makeLeaf({nodeItemIndices, numNodeItems}, m_itemIndices);
318 return;
319 }
320
321 const real noSplitCost = params.getInteractCost() * static_cast<real>(numNodeItems);
322 const real rcpNodeSurfaceArea = 1.0_r / nodeAABB.getSurfaceArea();
323 const Vector3R nodeExtents = nodeAABB.getExtents();
324
325 real bestSplitCost = std::numeric_limits<real>::max();
326 int bestAxis = -1;
327 std::size_t bestEndpointIndex = std::numeric_limits<std::size_t>::max();
328 int axis = static_cast<int>(nodeExtents.maxDimension());
329 int numSplitTrials = 0;
330 while(bestAxis == -1 && numSplitTrials < 3)
331 {
332 for(std::size_t i = 0; i < numNodeItems; ++i)
333 {
334 const Index itemIndex = nodeItemIndices[i];
335 const AABB3D& itemAABB = itemAABBs[itemIndex];
336 endpointsCache[axis][2 * i] = ItemEndpoint{itemAABB.getMinVertex()[axis], itemIndex, EEndpoint::MIN};
337 endpointsCache[axis][2 * i + 1] = ItemEndpoint{itemAABB.getMaxVertex()[axis], itemIndex, EEndpoint::MAX};
338 }
339
340 std::sort(&(endpointsCache[axis][0]), &(endpointsCache[axis][2 * numNodeItems]),
341 [](const ItemEndpoint& a, const ItemEndpoint& b) -> bool
342 {
343 return a.position != b.position ? a.position < b.position
344 : a.type < b.type;
345 });
346
347 std::size_t numNegativeItems = 0;
348 std::size_t numPositiveItems = numNodeItems;
349 for(std::size_t e = 0; e < 2 * numNodeItems; ++e)
350 {
351 if(endpointsCache[axis][e].type == EEndpoint::MAX)
352 {
353 --numPositiveItems;
354 }
355
356 // check if the split point is a reasonable one (within node AABB)
357 real endpoint = endpointsCache[axis][e].position;
358 if(endpoint > nodeAABB.getMinVertex()[axis] &&
359 endpoint < nodeAABB.getMaxVertex()[axis])
360 {
361 Vector3R endpointMinVertex = nodeAABB.getMinVertex();
362 Vector3R endpointMaxVertex = nodeAABB.getMaxVertex();
363 endpointMinVertex[axis] = endpoint;
364 endpointMaxVertex[axis] = endpoint;
365
366 const real probNegative = AABB3D(nodeAABB.getMinVertex(), endpointMaxVertex).getSurfaceArea() * rcpNodeSurfaceArea;
367 const real probPositive = AABB3D(endpointMinVertex, nodeAABB.getMaxVertex()).getSurfaceArea() * rcpNodeSurfaceArea;
368 const real emptyBonus = (numNegativeItems == 0 || numPositiveItems == 0) ? params.getEmptyBonus() : 0.0_r;
369 const real currentSplitCost = params.getTraversalCost() + (1.0_r - emptyBonus) * params.getInteractCost() *
370 (probNegative * static_cast<real>(numNegativeItems) + probPositive * static_cast<real>(numPositiveItems));
371
372 if(currentSplitCost < bestSplitCost)
373 {
374 bestSplitCost = currentSplitCost;
375 bestAxis = axis;
376 bestEndpointIndex = e;
377 }
378 }
379
380 if(endpointsCache[axis][e].type == EEndpoint::MIN)
381 {
382 ++numNegativeItems;
383 }
384 }
385
386 ++numSplitTrials;
387 axis = (axis + 1) % 3;
388 }
389
390 std::size_t newNumBadRefines = currentBadRefines;
391 if(bestSplitCost > noSplitCost)
392 {
393 ++newNumBadRefines;
394 }
395
396 if((bestSplitCost > 4 * noSplitCost && numNodeItems < 16) ||
397 bestAxis == -1 ||
398 newNumBadRefines == 3)
399 {
400 m_nodeBuffer[nodeIndex] = Node::makeLeaf({nodeItemIndices, numNodeItems}, m_itemIndices);
401 return;
402 }
403
404 std::size_t numNegativeItems = 0;
405 for(std::size_t e = 0; e < bestEndpointIndex; ++e)
406 {
407 if(endpointsCache[bestAxis][e].type == EEndpoint::MIN)
408 {
409 negativeItemIndicesCache[numNegativeItems++] = endpointsCache[bestAxis][e].index;
410 }
411 }
412 std::size_t numPositiveItems = 0;
413 for(std::size_t e = bestEndpointIndex + 1; e < 2 * numNodeItems; ++e)
414 {
415 if(endpointsCache[bestAxis][e].type == EEndpoint::MAX)
416 {
417 positiveItemIndicesCache[numPositiveItems++] = endpointsCache[bestAxis][e].index;
418 }
419 }
420 PH_ASSERT(numNegativeItems <= numNodeItems && numPositiveItems <= numNodeItems);
421
422 const real bestSplitPos = endpointsCache[bestAxis][bestEndpointIndex].position;
423
424 Vector3R splitPosMinVertex = nodeAABB.getMinVertex();
425 Vector3R splitPosMaxVertex = nodeAABB.getMaxVertex();
426 splitPosMinVertex[bestAxis] = bestSplitPos;
427 splitPosMaxVertex[bestAxis] = bestSplitPos;
428 const AABB3D negativeNodeAABB(nodeAABB.getMinVertex(), splitPosMaxVertex);
429 const AABB3D positiveNodeAABB(splitPosMinVertex, nodeAABB.getMaxVertex());
430
431 buildNodeRecursive(
432 nodeIndex + 1,
433 negativeNodeAABB,
434 negativeItemIndicesCache,
435 numNegativeItems,
436 currentNodeDepth + 1,
437 newNumBadRefines,
438 itemAABBs,
439 maxNodeDepth,
440 params,
441 negativeItemIndicesCache,
442 positiveItemIndicesCache + numNodeItems,
443 endpointsCache);
444
445 const std::size_t positiveChildIndex = m_numNodes;
446 m_nodeBuffer[nodeIndex] = Node::makeInner(bestSplitPos, bestAxis, positiveChildIndex);
447
448 buildNodeRecursive(
449 positiveChildIndex,
450 positiveNodeAABB,
451 positiveItemIndicesCache,
452 numPositiveItems,
453 currentNodeDepth + 1,
454 newNumBadRefines,
455 itemAABBs,
456 maxNodeDepth,
457 params,
458 negativeItemIndicesCache,
459 positiveItemIndicesCache + numNodeItems,
460 endpointsCache);
461}
462
463}// end namespace ph::math
acceleration_structure_basics.h
ph::math::IndexedKdtreeParams
Definition IndexedKdtreeParams.h:9
ph::math::TAABB3D::getSurfaceArea
T getSurfaceArea() const
Get the surface area of the bound.
Definition TAABB3D.ipp:177
ph::math::TAABB3D::getMinVertex
const TVector3< T > & getMinVertex() const
Get the corner vertex of the minimum (—) octant.
Definition TAABB3D.ipp:146
ph::math::TIndexedKdtree
Definition TIndexedKdtree.h:28
ph::math::TIndexedKdtree::Item
decltype(std::declval< IndexToItem >()(std::declval< Index >())) Item
Definition TIndexedKdtree.h:31
ph::math::TIndexedKdtreeNode
An indexed kD-tree node with compacted memory layout.
Definition TIndexedKdtreeNode.h:23
ph::math::TLineSegment
Represents a line segment in space.
Definition TLineSegment.h:25
ph::math::CItemSegmentIntersectionTester
Definition acceleration_structure_basics.h:38
ph::math::CItemSegmentIntersectionTesterWithIndex
Definition acceleration_structure_basics.h:28
ph::math
Math functions and utilities.
Definition TransformInfo.h:10
ph::math::AABB3D
TAABB3D< real > AABB3D
Definition TAABB3D.h:21
ph::math::Vector3R
TVector3< real > Vector3R
Definition math_fwd.h:52
std
Definition TAABB2D.h:96