From 3885a4c6f9f5d2962ad6327efc724ed3193ced29 Mon Sep 17 00:00:00 2001 From: Wenzel Jakob Date: Tue, 23 Oct 2012 12:00:52 -0400 Subject: [PATCH] new branch for motion blur and spacetime kd-trees --- include/mitsuba/render/gkdtree.h | 4 +- include/mitsuba/render/sahkdtree3.h | 3 +- include/mitsuba/render/sahkdtree4.h | 304 +++++++++++++++++++ src/shapes/SConscript | 1 + src/shapes/pointcache.cpp | 435 ++++++++++++++++++++++++++++ 5 files changed, 745 insertions(+), 2 deletions(-) create mode 100644 include/mitsuba/render/sahkdtree4.h create mode 100644 src/shapes/pointcache.cpp diff --git a/include/mitsuba/render/gkdtree.h b/include/mitsuba/render/gkdtree.h index c0509d97..4f19ec4e 100644 --- a/include/mitsuba/render/gkdtree.h +++ b/include/mitsuba/render/gkdtree.h @@ -1321,7 +1321,9 @@ protected: return a.axis < b.axis; if (a.pos != b.pos) return a.pos < b.pos; - return a.type < b.type; + if (a.type != b.type) + return a.type < b.type; + return a.index < b.index; } }; diff --git a/include/mitsuba/render/sahkdtree3.h b/include/mitsuba/render/sahkdtree3.h index 9db4338d..8053d121 100644 --- a/include/mitsuba/render/sahkdtree3.h +++ b/include/mitsuba/render/sahkdtree3.h @@ -63,7 +63,8 @@ public: * Given a split on axis \a axis that produces children having extents * \a leftWidth and \a rightWidth along \a axis, compute the probability * of traversing the left and right child during a typical query - * operation. + * operation. In the case of the surface area heuristic, this is simply + * the ratio of surface areas. */ inline std::pair operator()(int axis, Float leftWidth, Float rightWidth) const { return std::pair( diff --git a/include/mitsuba/render/sahkdtree4.h b/include/mitsuba/render/sahkdtree4.h new file mode 100644 index 00000000..6900bd80 --- /dev/null +++ b/include/mitsuba/render/sahkdtree4.h @@ -0,0 +1,304 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2012 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#if !defined(__SAH_KDTREE4_H) +#define __SAH_KDTREE4_H + +#include + +MTS_NAMESPACE_BEGIN + +typedef TAABB AABB4; + +/** + * \brief Implements the four-dimensional surface area heuristic for use + * by the \ref GenericKDTree construction algorithm. + */ +class SurfaceAreaHeuristic4 { +public: + /** + * \brief Initialize the surface area heuristic with the bounds of + * a parent node + * + * Precomputes some information so that traversal probabilities + * of potential split planes can be evaluated efficiently + */ + inline SurfaceAreaHeuristic4(const AABB4 &aabb) : m_aabb(aabb) { + const Vector4 extents(aabb.getExtents()); + const Float temp = 1.0f / (extents.x * extents.y + + extents.y*extents.z + extents.x*extents.z); + + m_temp0 = Vector4( + extents.y * extents.z * temp, + extents.x * extents.z * temp, + extents.x * extents.y * temp, + 0.0f); + + m_temp1 = Vector4( + (extents.y + extents.z) * temp, + (extents.x + extents.z) * temp, + (extents.x + extents.y) * temp, + 1.0f / extents.w); + } + + /** + * Given a split on axis \a axis that produces children having extents + * \a leftWidth and \a rightWidth along \a axis, compute the probability + * of traversing the left and right child during a typical query + * operation. + */ + inline std::pair operator()(int axis, Float leftWidth, Float rightWidth) const { + if (axis == 3 && m_temp1.w == std::numeric_limits::infinity()) { + return std::pair( + std::numeric_limits::infinity(), + std::numeric_limits::infinity() + ); + } + + return std::pair( + m_temp0[axis] + m_temp1[axis] * leftWidth, + m_temp0[axis] + m_temp1[axis] * rightWidth); + } + + /** + * Compute the underlying quantity used by the tree construction + * heuristic. This is used to compute the final cost of a kd-tree. + */ + inline static Float getQuantity(const AABB4 &aabb) { + const Vector4 extents(aabb.getExtents()); + Float result = 2 * (extents[0] * extents[1] + extents[1] * extents[2] + + extents[2] * extents[0]); + if (extents[3] != 0) + result *= extents[3]; + return result; + } +private: + Vector4 m_temp0, m_temp1; + AABB4 m_aabb; +}; + +/** + * This class specializes \ref GenericKDTree to a four-dimensional + * tree to be used for spacetime ray tracing. One additional function call + * must be implemented by subclasses: + * + * /// Check whether a primitive is intersected by the given ray. + * /// Some temporary space is supplied, which can be used to cache + * /// information about the intersection + * bool intersect(const Ray &ray, IndexType idx, + * Float mint, Float maxt, Float &t, void *tmp); + * + * This class implements an epsilon-free version of the optimized ray + * traversal algorithm (TA^B_{rec}), which is explained in Vlastimil + * Havran's PhD thesis "Heuristic Ray Shooting Algorithms". + * + * \author Wenzel Jakob + */ +template + class SAHKDTree4D : public GenericKDTree { +public: + typedef typename KDTreeBase::SizeType SizeType; + typedef typename KDTreeBase::IndexType IndexType; + typedef typename KDTreeBase::KDNode KDNode; + +protected: + void buildInternal() { + SizeType primCount = this->cast()->getPrimitiveCount(); + KDLog(EInfo, "Constructing a 4D SAH kd-tree (%i primitives) ..", primCount); + GenericKDTree::buildInternal(); + } + + /** + * \brief Hashed mailbox implementation + */ + struct HashedMailbox { + inline HashedMailbox() { + memset(entries, 0xFF, sizeof(IndexType)*MTS_KD_MAILBOX_SIZE); + } + + inline void put(IndexType primIndex) { + entries[primIndex & MTS_KD_MAILBOX_MASK] = primIndex; + } + + inline bool contains(IndexType primIndex) const { + return entries[primIndex & MTS_KD_MAILBOX_MASK] == primIndex; + } + + IndexType entries[MTS_KD_MAILBOX_SIZE]; + }; + + /// Ray traversal stack entry for Havran-style incoherent ray tracing + struct KDStackEntryHavran { + /* Pointer to the far child */ + const KDNode * __restrict node; + /* Distance traveled along the ray (entry or exit) */ + Float t; + /* Previous stack item */ + uint32_t prev; + /* Associated point */ + Point p; + }; + + /** + * \brief Ray tracing kd-tree traversal loop (Havran variant) + * + * This is generally the most robust and fastest traversal routine + * of the methods implemented in this class. + */ + template FINLINE + bool rayIntersectHavran(const Ray &ray, Float mint, Float maxt, + Float &t, void *temp) const { + KDStackEntryHavran stack[MTS_KD_MAXDEPTH]; + #if 0 + static const int prevAxisTable[] = { 2, 0, 1 }; + static const int nextAxisTable[] = { 1, 2, 0 }; + #endif + + #if defined(MTS_KD_MAILBOX_ENABLED) + HashedMailbox mailbox; + #endif + + /* Set up the entry point */ + uint32_t enPt = 0; + stack[enPt].t = mint; + stack[enPt].p = ray(mint); + + /* Set up the exit point */ + uint32_t exPt = 1; + stack[exPt].t = maxt; + stack[exPt].p = ray(maxt); + stack[exPt].node = NULL; + + bool foundIntersection = false; + const KDNode * __restrict currNode = this->m_nodes; + while (currNode != NULL) { + while (EXPECT_TAKEN(!currNode->isLeaf())) { + const Float splitVal = (Float) currNode->getSplit(); + const int axis = currNode->getAxis(); + const KDNode * __restrict farChild; + + if (axis == 3) { + if (ray.time <= splitVal) + currNode = currNode->getLeft(); + else + currNode = currNode->getRight(); + continue; + } else if (stack[enPt].p[axis] <= splitVal) { + if (stack[exPt].p[axis] <= splitVal) { + /* Cases N1, N2, N3, P5, Z2 and Z3 (see thesis) */ + currNode = currNode->getLeft(); + continue; + } + + /* Typo in Havran's thesis: + (it specifies "stack[exPt].p == splitVal", which + is clearly incorrect) */ + if (stack[enPt].p[axis] == splitVal) { + /* Case Z1 */ + currNode = currNode->getRight(); + continue; + } + + /* Case N4 */ + currNode = currNode->getLeft(); + farChild = currNode + 1; // getRight() + } else { /* stack[enPt].p[axis] > splitVal */ + if (splitVal < stack[exPt].p[axis]) { + /* Cases P1, P2, P3 and N5 */ + currNode = currNode->getRight(); + continue; + } + /* Case P4 */ + farChild = currNode->getLeft(); + currNode = farChild + 1; // getRight() + } + + /* Cases P4 and N4 -- calculate the distance to the split plane */ + Float distToSplit = (splitVal - ray.o[axis]) * ray.dRcp[axis]; + + /* Set up a new exit point */ + const uint32_t tmp = exPt++; + if (exPt == enPt) /* Do not overwrite the entry point */ + ++exPt; + + KDAssert(exPt < MTS_KD_MAXDEPTH); + stack[exPt].prev = tmp; + stack[exPt].t = distToSplit; + stack[exPt].node = farChild; + + #if 1 + /* Intrestingly, this appears to be faster than the + original code with the prevAxis & nextAxis table */ + stack[exPt].p = ray(distToSplit); + stack[exPt].p[axis] = splitVal; + #else + const int nextAxis = nextAxisTable[axis]; + const int prevAxis = prevAxisTable[axis]; + stack[exPt].p[axis] = splitVal; + stack[exPt].p[nextAxis] = ray.o[nextAxis] + + distToSplit*ray.d[nextAxis]; + stack[exPt].p[prevAxis] = ray.o[prevAxis] + + distToSplit*ray.d[prevAxis]; + #endif + + } + + /* Reached a leaf node */ + for (IndexType entry=currNode->getPrimStart(), + last = currNode->getPrimEnd(); entry != last; entry++) { + const IndexType primIdx = this->m_indices[entry]; + + #if defined(MTS_KD_MAILBOX_ENABLED) + if (mailbox.contains(primIdx)) + continue; + #endif + + bool result; + if (!shadowRay) + result = this->cast()->intersect(ray, primIdx, mint, maxt, t, temp); + else + result = this->cast()->intersect(ray, primIdx, mint, maxt); + + if (result) { + if (shadowRay) + return true; + maxt = t; + foundIntersection = true; + } + + #if defined(MTS_KD_MAILBOX_ENABLED) + mailbox.put(primIdx); + #endif + } + + if (stack[exPt].t > maxt) + break; + + /* Pop from the stack and advance to the next node on the interval */ + enPt = exPt; + currNode = stack[exPt].node; + exPt = stack[enPt].prev; + } + + return foundIntersection; + } +}; + +MTS_NAMESPACE_END + +#endif /* __SAH_KDTREE4_H */ diff --git a/src/shapes/SConscript b/src/shapes/SConscript index 6abc53ce..2561dc83 100644 --- a/src/shapes/SConscript +++ b/src/shapes/SConscript @@ -12,5 +12,6 @@ plugins += env.SharedLibrary('hair', ['hair.cpp']) plugins += env.SharedLibrary('shapegroup', ['shapegroup.cpp']) plugins += env.SharedLibrary('instance', ['instance.cpp']) plugins += env.SharedLibrary('animatedinstance', ['animatedinstance.cpp']) +plugins += env.SharedLibrary('pointcache', ['pointcache.cpp']) Export('plugins') diff --git a/src/shapes/pointcache.cpp b/src/shapes/pointcache.cpp new file mode 100644 index 00000000..f208690d --- /dev/null +++ b/src/shapes/pointcache.cpp @@ -0,0 +1,435 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2012 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include +#include +#include +#include +#include +#include +#include + +#define SHAPE_PER_SEGMENT 1 +#define NO_CLIPPING_SUPPORT 1 + +MTS_NAMESPACE_BEGIN + +class SpaceTimeKDTree : public SAHKDTree4D { + friend class GenericKDTree; + friend class SAHKDTree4D; +public: + /// Temporarily holds some intersection information + struct IntersectionCache { + Point p[3]; + Float u, v; + }; + + SpaceTimeKDTree(const std::vector &frameTimes, std::vector &positions, + Triangle *triangles, size_t vertexCount, size_t triangleCount) + : m_frameTimes(frameTimes), m_positions(positions), m_triangles(triangles), + m_vertexCount(vertexCount), m_triangleCount(triangleCount) { + + Log(EInfo, "Total amount of vertex data: %s", + memString(vertexCount*frameTimes.size()*sizeof(float)*3).c_str()); + + //setClip(false); + //setExactPrimitiveThreshold(10); + buildInternal(); + + /* Collect some statistics */ + std::stack stack; + + stack.push(m_nodes); + size_t spatialSplits = 0, timeSplits = 0; + while (!stack.empty()) { + const KDNode *node = stack.top(); + stack.pop(); + if (!node->isLeaf()) { + if (node->getAxis() == 3) { + timeSplits++; + } else { + spatialSplits++; + } + stack.push((const KDNode *) node->getLeft()); + stack.push((const KDNode *) node->getRight()); + } + } + + KDLog(EInfo, "Spacetime kd-tree statistics"); + KDLog(EInfo, " Time interval = [%f, %f]" , m_tightAABB.min.w, m_tightAABB.max.w); + KDLog(EInfo, " Spatial splits = " SIZE_T_FMT, spatialSplits); + KDLog(EInfo, " Time splits = " SIZE_T_FMT, timeSplits); + KDLog(EInfo, ""); + + m_spatialAABB = AABB( + Point(m_aabb.min.x, m_aabb.min.y, m_aabb.min.z), + Point(m_aabb.max.x, m_aabb.max.y, m_aabb.max.z) + ); + } + + /// Return one of the points stored in the point cache + inline Point getPoint(uint32_t frame, uint32_t index) const { + float *ptr = m_positions[frame] + index*3; +#if defined(__LITTLE_ENDIAN__) + return Point( + (Float) endianness_swap(ptr[0]), + (Float) endianness_swap(ptr[1]), + (Float) endianness_swap(ptr[2])); +#else + return Point((Float) ptr[0], (Float) ptr[1], (Float) ptr[2]); +#endif + } + + // ======================================================================== + // Implementation of functions required by the parent class + // ======================================================================== + + /// Return the total number of primitives that are organized in the tree + inline SizeType getPrimitiveCount() const { +#ifdef SHAPE_PER_SEGMENT + return m_triangleCount * (m_frameTimes.size() - 1); +#else + return m_triangleCount; +#endif + } + + /// Return the 4D extents for one of the primitives contained in the tree + AABB4 getAABB(IndexType index) const { +#ifdef SHAPE_PER_SEGMENT + int frameIdx = index / m_triangleCount; + int triangleIdx = index % m_triangleCount; + const Triangle &tri = m_triangles[triangleIdx]; + + AABB aabb; + for (int i=0; i<3; ++i) { + aabb.expandBy(getPoint(frameIdx, tri.idx[i])); + aabb.expandBy(getPoint(frameIdx+1, tri.idx[i])); + } + + return AABB4( + Point4(aabb.min.x, aabb.min.y, aabb.min.z, m_frameTimes[frameIdx]), + Point4(aabb.max.x, aabb.max.y, aabb.max.z, m_frameTimes[frameIdx+1]) + ); +#else + AABB aabb; + const Triangle &tri = m_triangles[index]; + for (size_t i=0; i 1) + return false; + + /* Compute interpolated positions */ + Point p[3]; + for (int i=0; i<3; ++i) + p[i] = (1 - alpha) * getPoint(frameIdx, tri.idx[i]) + + alpha * getPoint(frameIdx+1, tri.idx[i]); + + Float tempU, tempV, tempT; + if (!Triangle::rayIntersect(p[0], p[1], p[2], ray, tempU, tempV, tempT)) + return false; + if (tempT < mint || tempT > maxt) + return false; + + if (tmp != NULL) { + IntersectionCache *cache = + static_cast(tmp); + t = tempT; + memcpy(cache->p, p, sizeof(Point)*3); + cache->u = tempU; + cache->v = tempV; + } + return true; + } + + /// Cast a shadow ray against a specific triangle + inline bool intersect(const Ray &ray, IndexType idx, + Float mint, Float maxt) const { + Float tempT; + /* No optimized version for shadow rays yet */ + return intersect(ray, idx, mint, maxt, tempT, NULL); + } + + // ======================================================================== + // Miscellaneous + // ======================================================================== + + /// Intersect a ray with all primitives stored in the kd-tree + inline bool rayIntersect(const Ray &ray, Float _mint, Float _maxt, + Float &t, void *temp) const { + Float tempT = std::numeric_limits::infinity(); + Float mint, maxt; + + if (m_spatialAABB.rayIntersect(ray, mint, maxt)) { + if (_mint > mint) mint = _mint; + if (_maxt < maxt) maxt = _maxt; + + if (EXPECT_TAKEN(maxt > mint && ray.time >= m_aabb.min.w && ray.time <= m_aabb.max.w)) { + if (rayIntersectHavran(ray, mint, maxt, tempT, temp)) { + t = tempT; + return true; + } + } + } + return false; + } + + /** + * \brief Intersect a ray with all primitives stored in the kd-tree + * (Visiblity query version) + */ + inline bool rayIntersect(const Ray &ray, Float _mint, Float _maxt) const { + Float tempT = std::numeric_limits::infinity(); + Float mint, maxt; + + if (m_spatialAABB.rayIntersect(ray, mint, maxt)) { + if (_mint > mint) mint = _mint; + if (_maxt < maxt) maxt = _maxt; + + if (EXPECT_TAKEN(maxt > mint && ray.time >= m_aabb.min.w && ray.time <= m_aabb.max.w)) + if (rayIntersectHavran(ray, mint, maxt, tempT, NULL)) + return true; + } + return false; + } + + inline const Triangle *getTriangles() const { + return m_triangles; + } + + /// Return an AABB with the spatial extents + inline const AABB &getSpatialAABB() const { + return m_spatialAABB; + } + + MTS_DECLARE_CLASS() +protected: + std::vector m_frameTimes; + std::vector m_positions; + Triangle *m_triangles; + size_t m_vertexCount; + size_t m_triangleCount; + AABB m_spatialAABB; +}; + +class PointCache : public Shape { +public: + PointCache(const Properties &props) : Shape(props) { + FileResolver *fResolver = Thread::getThread()->getFileResolver(); + fs::path path = fResolver->resolve(props.getString("filename")); + if (path.extension() != ".mdd") + Log(EError, "Point cache files must have the extension \".mdd\""); + + m_mmap = new MemoryMappedFile(path); + + ref mStream = new MemoryStream((uint8_t *) m_mmap->getData(), + m_mmap->getSize()); + mStream->setByteOrder(Stream::EBigEndian); + + uint32_t frameCount = mStream->readUInt(); + m_vertexCount = mStream->readUInt(); + + Log(EInfo, "Point cache has %i frames and %i vertices", frameCount, m_vertexCount); + + for (uint32_t i=0; ireadSingle()); + + for (uint32_t i=0; i(mStream->getCurrentData())); + mStream->skip(m_vertexCount * 3 * sizeof(float)); + } + Assert(mStream->getPos() == mStream->getSize()); + } + + PointCache(Stream *stream, InstanceManager *manager) + : Shape(stream, manager) { + /// TBD + } + + void serialize(Stream *stream, InstanceManager *manager) const { + Shape::serialize(stream, manager); + /// TBD + } + + void configure() { + Shape::configure(); + + if (m_mesh == NULL) + Log(EError, "A nested triangle mesh is required so that " + "connectivity information can be extracted!"); + if (m_mesh->getVertexCount() != m_vertexCount) + Log(EError, "Point cache and nested geometry have mismatched " + "numbers of vertices!"); + + m_kdtree = new SpaceTimeKDTree(m_frameTimes, m_positions, m_mesh->getTriangles(), + m_vertexCount, m_mesh->getTriangleCount()); + m_aabb = m_kdtree->getSpatialAABB(); + } + + bool rayIntersect(const Ray &ray, Float mint, + Float maxt, Float &t, void *temp) const { + return m_kdtree->rayIntersect(ray, mint, maxt, t, temp); + } + + bool rayIntersect(const Ray &ray, Float mint, Float maxt) const { + return m_kdtree->rayIntersect(ray, mint, maxt); + } + + void fillIntersectionRecord(const Ray &ray, + const void *temp, Intersection &its) const { + const SpaceTimeKDTree::IntersectionCache *cache + = reinterpret_cast(temp); + + const Vector b(1 - cache->u - cache->v, cache->u, cache->v); + const Point p0 = cache->p[0]; + const Point p1 = cache->p[1]; + const Point p2 = cache->p[2]; + + Normal faceNormal(cross(p1-p0, p2-p0)); + Float length = faceNormal.length(); + if (!faceNormal.isZero()) + faceNormal /= length; + + /* Just the basic attributes for now and geometric normals */ + its.p = ray(its.t); + its.geoFrame = Frame(faceNormal); + its.shFrame = its.geoFrame; + its.wi = its.toLocal(-ray.d); + its.shape = this; + its.hasUVPartials = false; + its.time = ray.time; + } + + AABB getAABB() const { + return m_kdtree->getSpatialAABB(); + } + + Float getSurfaceArea() const { + Log(EError, "PointCache::getSurfaceArea(): Not implemented."); + return -1; + } + + size_t getPrimitiveCount() const { + return m_mesh->getTriangleCount(); + } + + size_t getEffectivePrimitiveCount() const { + return m_mesh->getTriangleCount(); + } + + void addChild(const std::string &name, ConfigurableObject *child) { + const Class *cClass = child->getClass(); + if (cClass->derivesFrom(TriMesh::m_theClass)) { + Assert(m_mesh == NULL); + m_mesh = static_cast(child); + } else { + Shape::addChild(name, child); + } + } + + std::string toString() const { + std::ostringstream oss; + oss << "PointCache[" << endl + << "]"; + return oss.str(); + } + + MTS_DECLARE_CLASS() +private: + ref m_mmap; + ref m_kdtree; + std::vector m_frameTimes; + std::vector m_positions; + ref m_mesh; + uint32_t m_vertexCount; + AABB m_aabb; +}; + +MTS_IMPLEMENT_CLASS(SpaceTimeKDTree, false, KDTreeBase) +MTS_IMPLEMENT_CLASS_S(PointCache, false, Shape) +MTS_EXPORT_PLUGIN(PointCache, "Point cache"); +MTS_NAMESPACE_END