| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // Copyright Contributors to the OpenVDB Project | ||
| 2 | // SPDX-License-Identifier: MPL-2.0 | ||
| 3 | /// | ||
| 4 | /// @file RayIntersector.h | ||
| 5 | /// | ||
| 6 | /// @author Ken Museth | ||
| 7 | /// | ||
| 8 | /// @brief Accelerated intersection of a ray with a narrow-band level | ||
| 9 | /// set or a generic (e.g. density) volume. This will of course be | ||
| 10 | /// useful for respectively surface and volume rendering. | ||
| 11 | /// | ||
| 12 | /// @details This file defines two main classes, | ||
| 13 | /// LevelSetRayIntersector and VolumeRayIntersector, as well as the | ||
| 14 | /// three support classes LevelSetHDDA, VolumeHDDA and LinearSearchImpl. | ||
| 15 | /// The LevelSetRayIntersector is templated on the LinearSearchImpl class | ||
| 16 | /// and calls instances of the LevelSetHDDA class. The reason to split | ||
| 17 | /// level set ray intersection into three classes is twofold. First | ||
| 18 | /// LevelSetRayIntersector defines the public API for client code and | ||
| 19 | /// LinearSearchImpl defines the actual algorithm used for the | ||
| 20 | /// ray level-set intersection. In other words this design will allow | ||
| 21 | /// for the public API to be fixed while the intersection algorithm | ||
| 22 | /// can change without resolving to (slow) virtual methods. Second, | ||
| 23 | /// LevelSetHDDA, which implements a hierarchical Differential Digital | ||
| 24 | /// Analyzer, relies on partial template specialization, so it has to | ||
| 25 | /// be a standalone class (as opposed to a member class of | ||
| 26 | /// LevelSetRayIntersector). The VolumeRayIntersector is conceptually | ||
| 27 | /// much simpler than the LevelSetRayIntersector, and hence it only | ||
| 28 | /// depends on VolumeHDDA that implements the hierarchical | ||
| 29 | /// Differential Digital Analyzer. | ||
| 30 | |||
| 31 | |||
| 32 | #ifndef OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED | ||
| 33 | #define OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED | ||
| 34 | |||
| 35 | #include <openvdb/math/DDA.h> | ||
| 36 | #include <openvdb/math/Math.h> | ||
| 37 | #include <openvdb/math/Ray.h> | ||
| 38 | #include <openvdb/math/Stencils.h> | ||
| 39 | #include <openvdb/Grid.h> | ||
| 40 | #include <openvdb/Types.h> | ||
| 41 | #include "Morphology.h" | ||
| 42 | #include <iostream> | ||
| 43 | #include <type_traits> | ||
| 44 | |||
| 45 | |||
| 46 | namespace openvdb { | ||
| 47 | OPENVDB_USE_VERSION_NAMESPACE | ||
| 48 | namespace OPENVDB_VERSION_NAME { | ||
| 49 | namespace tools { | ||
| 50 | |||
| 51 | // Helper class that implements the actual search of the zero-crossing | ||
| 52 | // of the level set along the direction of a ray. This particular | ||
| 53 | // implementation uses iterative linear search. | ||
| 54 | template<typename GridT, int Iterations = 0, typename RealT = double> | ||
| 55 | class LinearSearchImpl; | ||
| 56 | |||
| 57 | |||
| 58 | ///////////////////////////////////// LevelSetRayIntersector ///////////////////////////////////// | ||
| 59 | |||
| 60 | |||
| 61 | /// @brief This class provides the public API for intersecting a ray | ||
| 62 | /// with a narrow-band level set. | ||
| 63 | /// | ||
| 64 | /// @details It wraps a SearchImplT with a simple public API and | ||
| 65 | /// performs the actual hierarchical tree node and voxel traversal. | ||
| 66 | /// | ||
| 67 | /// @warning Use the (default) copy-constructor to make sure each | ||
| 68 | /// computational thread has their own instance of this class. This is | ||
| 69 | /// important since the SearchImplT contains a ValueAccessor that is | ||
| 70 | /// not thread-safe. However copying is very efficient. | ||
| 71 | /// | ||
| 72 | /// @see tools/RayTracer.h for examples of intended usage. | ||
| 73 | /// | ||
| 74 | /// @todo Add TrilinearSearchImpl, as an alternative to LinearSearchImpl, | ||
| 75 | /// that performs analytical 3D trilinear intersection tests, i.e., solves | ||
| 76 | /// cubic equations. This is slower but also more accurate than the 1D | ||
| 77 | /// linear interpolation in LinearSearchImpl. | ||
| 78 | template<typename GridT, | ||
| 79 | typename SearchImplT = LinearSearchImpl<GridT>, | ||
| 80 | int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL, | ||
| 81 | typename RayT = math::Ray<Real> > | ||
| 82 |
9/18✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
|
9 | class LevelSetRayIntersector |
| 83 | { | ||
| 84 | public: | ||
| 85 | using GridType = GridT; | ||
| 86 | using RayType = RayT; | ||
| 87 | using RealType = typename RayT::RealType; | ||
| 88 | using Vec3Type = typename RayT::Vec3T; | ||
| 89 | using ValueT = typename GridT::ValueType; | ||
| 90 | using TreeT = typename GridT::TreeType; | ||
| 91 | |||
| 92 | static_assert(NodeLevel >= -1 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range"); | ||
| 93 | static_assert(std::is_floating_point<ValueT>::value, | ||
| 94 | "level set grids must have scalar, floating-point value types"); | ||
| 95 | |||
| 96 | /// @brief Constructor | ||
| 97 | /// @param grid level set grid to intersect rays against. | ||
| 98 | /// @param isoValue optional iso-value for the ray-intersection. | ||
| 99 | 18 | LevelSetRayIntersector(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>()) | |
| 100 | 18 | : mTester(grid, isoValue) | |
| 101 | { | ||
| 102 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
|
18 | if (!grid.hasUniformVoxels() ) { |
| 103 | ✗ | OPENVDB_THROW(RuntimeError, | |
| 104 | "LevelSetRayIntersector only supports uniform voxels!"); | ||
| 105 | } | ||
| 106 |
2/4✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
|
18 | if (grid.getGridClass() != GRID_LEVEL_SET) { |
| 107 | ✗ | OPENVDB_THROW(RuntimeError, | |
| 108 | "LevelSetRayIntersector only supports level sets!" | ||
| 109 | "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)"); | ||
| 110 | } | ||
| 111 | 18 | } | |
| 112 | |||
| 113 | /// @brief Return the iso-value used for ray-intersections | ||
| 114 | const ValueT& getIsoValue() const { return mTester.getIsoValue(); } | ||
| 115 | |||
| 116 | /// @brief Return @c true if the index-space ray intersects the level set. | ||
| 117 | /// @param iRay ray represented in index space. | ||
| 118 | bool intersectsIS(const RayType& iRay) const | ||
| 119 | { | ||
| 120 | if (!mTester.setIndexRay(iRay)) return false;//missed bbox | ||
| 121 | return math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester); | ||
| 122 | } | ||
| 123 | |||
| 124 | /// @brief Return @c true if the index-space ray intersects the level set | ||
| 125 | /// @param iRay ray represented in index space. | ||
| 126 | /// @param iTime if an intersection was found it is assigned the time of the | ||
| 127 | /// intersection along the index ray. | ||
| 128 | bool intersectsIS(const RayType& iRay, RealType &iTime) const | ||
| 129 | { | ||
| 130 | if (!mTester.setIndexRay(iRay)) return false;//missed bbox | ||
| 131 | iTime = mTester.getIndexTime(); | ||
| 132 | return math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester); | ||
| 133 | } | ||
| 134 | |||
| 135 | /// @brief Return @c true if the index-space ray intersects the level set. | ||
| 136 | /// @param iRay ray represented in index space. | ||
| 137 | /// @param xyz if an intersection was found it is assigned the | ||
| 138 | /// intersection point in index space, otherwise it is unchanged. | ||
| 139 | bool intersectsIS(const RayType& iRay, Vec3Type& xyz) const | ||
| 140 | { | ||
| 141 | if (!mTester.setIndexRay(iRay)) return false;//missed bbox | ||
| 142 | if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set | ||
| 143 | mTester.getIndexPos(xyz); | ||
| 144 | return true; | ||
| 145 | } | ||
| 146 | |||
| 147 | /// @brief Return @c true if the index-space ray intersects the level set. | ||
| 148 | /// @param iRay ray represented in index space. | ||
| 149 | /// @param xyz if an intersection was found it is assigned the | ||
| 150 | /// intersection point in index space, otherwise it is unchanged. | ||
| 151 | /// @param iTime if an intersection was found it is assigned the time of the | ||
| 152 | /// intersection along the index ray. | ||
| 153 | bool intersectsIS(const RayType& iRay, Vec3Type& xyz, RealType &iTime) const | ||
| 154 | { | ||
| 155 | if (!mTester.setIndexRay(iRay)) return false;//missed bbox | ||
| 156 | if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set | ||
| 157 | mTester.getIndexPos(xyz); | ||
| 158 | iTime = mTester.getIndexTime(); | ||
| 159 | return true; | ||
| 160 | } | ||
| 161 | |||
| 162 | /// @brief Return @c true if the world-space ray intersects the level set. | ||
| 163 | /// @param wRay ray represented in world space. | ||
| 164 | bool intersectsWS(const RayType& wRay) const | ||
| 165 | { | ||
| 166 | if (!mTester.setWorldRay(wRay)) return false;//missed bbox | ||
| 167 | return math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester); | ||
| 168 | } | ||
| 169 | |||
| 170 | /// @brief Return @c true if the world-space ray intersects the level set. | ||
| 171 | /// @param wRay ray represented in world space. | ||
| 172 | /// @param wTime if an intersection was found it is assigned the time of the | ||
| 173 | /// intersection along the world ray. | ||
| 174 | bool intersectsWS(const RayType& wRay, RealType &wTime) const | ||
| 175 | { | ||
| 176 | if (!mTester.setWorldRay(wRay)) return false;//missed bbox | ||
| 177 | wTime = mTester.getWorldTime(); | ||
| 178 | return math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester); | ||
| 179 | } | ||
| 180 | |||
| 181 | /// @brief Return @c true if the world-space ray intersects the level set. | ||
| 182 | /// @param wRay ray represented in world space. | ||
| 183 | /// @param world if an intersection was found it is assigned the | ||
| 184 | /// intersection point in world space, otherwise it is unchanged | ||
| 185 | bool intersectsWS(const RayType& wRay, Vec3Type& world) const | ||
| 186 | { | ||
| 187 | if (!mTester.setWorldRay(wRay)) return false;//missed bbox | ||
| 188 | if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set | ||
| 189 | mTester.getWorldPos(world); | ||
| 190 | return true; | ||
| 191 | } | ||
| 192 | |||
| 193 | /// @brief Return @c true if the world-space ray intersects the level set. | ||
| 194 | /// @param wRay ray represented in world space. | ||
| 195 | /// @param world if an intersection was found it is assigned the | ||
| 196 | /// intersection point in world space, otherwise it is unchanged. | ||
| 197 | /// @param wTime if an intersection was found it is assigned the time of the | ||
| 198 | /// intersection along the world ray. | ||
| 199 | 2097168 | bool intersectsWS(const RayType& wRay, Vec3Type& world, RealType &wTime) const | |
| 200 | { | ||
| 201 |
2/2✓ Branch 1 taken 266264 times.
✓ Branch 2 taken 782320 times.
|
2097168 | if (!mTester.setWorldRay(wRay)) return false;//missed bbox |
| 202 |
2/2✓ Branch 1 taken 205869 times.
✓ Branch 2 taken 60395 times.
|
532528 | if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set |
| 203 | 411738 | mTester.getWorldPos(world); | |
| 204 | 411738 | wTime = mTester.getWorldTime(); | |
| 205 | 411738 | return true; | |
| 206 | } | ||
| 207 | |||
| 208 | /// @brief Return @c true if the world-space ray intersects the level set. | ||
| 209 | /// @param wRay ray represented in world space. | ||
| 210 | /// @param world if an intersection was found it is assigned the | ||
| 211 | /// intersection point in world space, otherwise it is unchanged. | ||
| 212 | /// @param normal if an intersection was found it is assigned the normal | ||
| 213 | /// of the level set surface in world space, otherwise it is unchanged. | ||
| 214 | ✗ | bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal) const | |
| 215 | { | ||
| 216 | ✗ | if (!mTester.setWorldRay(wRay)) return false;//missed bbox | |
| 217 | ✗ | if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set | |
| 218 | ✗ | mTester.getWorldPosAndNml(world, normal); | |
| 219 | ✗ | return true; | |
| 220 | } | ||
| 221 | |||
| 222 | /// @brief Return @c true if the world-space ray intersects the level set. | ||
| 223 | /// @param wRay ray represented in world space. | ||
| 224 | /// @param world if an intersection was found it is assigned the | ||
| 225 | /// intersection point in world space, otherwise it is unchanged. | ||
| 226 | /// @param normal if an intersection was found it is assigned the normal | ||
| 227 | /// of the level set surface in world space, otherwise it is unchanged. | ||
| 228 | /// @param wTime if an intersection was found it is assigned the time of the | ||
| 229 | /// intersection along the world ray. | ||
| 230 | bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal, RealType &wTime) const | ||
| 231 | { | ||
| 232 | if (!mTester.setWorldRay(wRay)) return false;//missed bbox | ||
| 233 | if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set | ||
| 234 | mTester.getWorldPosAndNml(world, normal); | ||
| 235 | wTime = mTester.getWorldTime(); | ||
| 236 | return true; | ||
| 237 | } | ||
| 238 | |||
| 239 | private: | ||
| 240 | |||
| 241 | mutable SearchImplT mTester; | ||
| 242 | |||
| 243 | };// LevelSetRayIntersector | ||
| 244 | |||
| 245 | |||
| 246 | ////////////////////////////////////// VolumeRayIntersector ////////////////////////////////////// | ||
| 247 | |||
| 248 | |||
| 249 | /// @brief This class provides the public API for intersecting a ray | ||
| 250 | /// with a generic (e.g. density) volume. | ||
| 251 | /// @details Internally it performs the actual hierarchical tree node traversal. | ||
| 252 | /// @warning Use the (default) copy-constructor to make sure each | ||
| 253 | /// computational thread has their own instance of this class. This is | ||
| 254 | /// important since it contains a ValueAccessor that is | ||
| 255 | /// not thread-safe and a CoordBBox of the active voxels that should | ||
| 256 | /// not be re-computed for each thread. However copying is very efficient. | ||
| 257 | /// @par Example: | ||
| 258 | /// @code | ||
| 259 | /// // Create an instance for the master thread | ||
| 260 | /// VolumeRayIntersector inter(grid); | ||
| 261 | /// // For each additional thread use the copy constructor. This | ||
| 262 | /// // amortizes the overhead of computing the bbox of the active voxels! | ||
| 263 | /// VolumeRayIntersector inter2(inter); | ||
| 264 | /// // Before each ray-traversal set the index ray. | ||
| 265 | /// iter.setIndexRay(ray); | ||
| 266 | /// // or world ray | ||
| 267 | /// iter.setWorldRay(ray); | ||
| 268 | /// // Now you can begin the ray-marching using consecutive calls to VolumeRayIntersector::march | ||
| 269 | /// double t0=0, t1=0;// note the entry and exit times are with respect to the INDEX ray | ||
| 270 | /// while ( inter.march(t0, t1) ) { | ||
| 271 | /// // perform line-integration between t0 and t1 | ||
| 272 | /// }} | ||
| 273 | /// @endcode | ||
| 274 | template<typename GridT, | ||
| 275 | int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL, | ||
| 276 | typename RayT = math::Ray<Real> > | ||
| 277 | class VolumeRayIntersector | ||
| 278 | { | ||
| 279 | public: | ||
| 280 | using GridType = GridT; | ||
| 281 | using RayType = RayT; | ||
| 282 | using RealType = typename RayT::RealType; | ||
| 283 | using RootType = typename GridT::TreeType::RootNodeType; | ||
| 284 | using TreeT = tree::Tree<typename RootType::template ValueConverter<bool>::Type>; | ||
| 285 | |||
| 286 | static_assert(NodeLevel >= 0 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range"); | ||
| 287 | |||
| 288 | /// @brief Grid constructor | ||
| 289 | /// @param grid Generic grid to intersect rays against. | ||
| 290 | /// @param dilationCount The number of voxel dilations performed | ||
| 291 | /// on (a boolean copy of) the input grid. This allows the | ||
| 292 | /// intersector to account for the size of interpolation kernels | ||
| 293 | /// in client code. | ||
| 294 | /// @throw RuntimeError if the voxels of the grid are not uniform | ||
| 295 | /// or the grid is empty. | ||
| 296 | 13 | VolumeRayIntersector(const GridT& grid, int dilationCount = 0) | |
| 297 | : mIsMaster(true) | ||
| 298 |
2/4✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
|
13 | , mTree(new TreeT(grid.tree(), false, TopologyCopy())) |
| 299 | , mGrid(&grid) | ||
| 300 |
1/2✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
|
26 | , mAccessor(*mTree) |
| 301 | { | ||
| 302 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
|
13 | if (!grid.hasUniformVoxels() ) { |
| 303 | ✗ | OPENVDB_THROW(RuntimeError, | |
| 304 | "VolumeRayIntersector only supports uniform voxels!"); | ||
| 305 | } | ||
| 306 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
|
13 | if ( grid.empty() ) { |
| 307 | ✗ | OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids"); | |
| 308 | } | ||
| 309 | |||
| 310 | // Dilate active voxels to better account for the size of interpolation kernels | ||
| 311 | 13 | tools::dilateActiveValues(*mTree, dilationCount, tools::NN_FACE, tools::IGNORE_TILES); | |
| 312 | |||
| 313 | 13 | mTree->root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false); | |
| 314 | |||
| 315 | mBBox.max().offset(1);//padding so the bbox of a node becomes (origin,origin + node_dim) | ||
| 316 | 13 | } | |
| 317 | |||
| 318 | /// @brief Grid and BBox constructor | ||
| 319 | /// @param grid Generic grid to intersect rays against. | ||
| 320 | /// @param bbox The axis-aligned bounding-box in the index space of the grid. | ||
| 321 | /// @warning It is assumed that bbox = (min, min + dim) where min denotes | ||
| 322 | /// to the smallest grid coordinates and dim are the integer length of the bbox. | ||
| 323 | /// @throw RuntimeError if the voxels of the grid are not uniform | ||
| 324 | /// or the grid is empty. | ||
| 325 | VolumeRayIntersector(const GridT& grid, const math::CoordBBox& bbox) | ||
| 326 | : mIsMaster(true) | ||
| 327 | , mTree(new TreeT(grid.tree(), false, TopologyCopy())) | ||
| 328 | , mGrid(&grid) | ||
| 329 | , mAccessor(*mTree) | ||
| 330 | , mBBox(bbox) | ||
| 331 | { | ||
| 332 | if (!grid.hasUniformVoxels() ) { | ||
| 333 | OPENVDB_THROW(RuntimeError, | ||
| 334 | "VolumeRayIntersector only supports uniform voxels!"); | ||
| 335 | } | ||
| 336 | if ( grid.empty() ) { | ||
| 337 | OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids"); | ||
| 338 | } | ||
| 339 | } | ||
| 340 | |||
| 341 | /// @brief Shallow copy constructor | ||
| 342 | /// @warning This copy constructor creates shallow copies of data | ||
| 343 | /// members of the instance passed as the argument. For | ||
| 344 | /// performance reasons we are not using shared pointers (their | ||
| 345 | /// mutex-lock impairs multi-threading). | ||
| 346 | ✗ | VolumeRayIntersector(const VolumeRayIntersector& other) | |
| 347 | : mIsMaster(false) | ||
| 348 | ✗ | , mTree(other.mTree)//shallow copy | |
| 349 | ✗ | , mGrid(other.mGrid)//shallow copy | |
| 350 | , mAccessor(*mTree)//initialize new (vs deep copy) | ||
| 351 | , mRay(other.mRay)//deep copy | ||
| 352 | ✗ | , mTmax(other.mTmax)//deep copy | |
| 353 | ✗ | , mBBox(other.mBBox)//deep copy | |
| 354 | { | ||
| 355 | } | ||
| 356 | |||
| 357 | /// @brief Destructor | ||
| 358 |
28/108✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 1 times.
✗ Branch 27 not taken.
✓ Branch 30 taken 1 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 1 times.
✗ Branch 33 not taken.
✓ Branch 36 taken 1 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 1 times.
✗ Branch 39 not taken.
✓ Branch 42 taken 1 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 1 times.
✗ Branch 45 not taken.
✓ Branch 48 taken 1 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 1 times.
✗ Branch 51 not taken.
✓ Branch 54 taken 1 times.
✗ Branch 55 not taken.
✓ Branch 56 taken 1 times.
✗ Branch 57 not taken.
✓ Branch 60 taken 1 times.
✗ Branch 61 not taken.
✓ Branch 62 taken 1 times.
✗ Branch 63 not taken.
✓ Branch 65 taken 1 times.
✗ Branch 66 not taken.
✓ Branch 67 taken 1 times.
✗ Branch 68 not taken.
✓ Branch 69 taken 1 times.
✗ Branch 70 not taken.
✓ Branch 72 taken 1 times.
✗ Branch 73 not taken.
✓ Branch 74 taken 1 times.
✗ Branch 75 not taken.
✓ Branch 76 taken 1 times.
✗ Branch 77 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
✗ Branch 98 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 120 not taken.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 132 not taken.
✗ Branch 133 not taken.
✗ Branch 135 not taken.
✗ Branch 136 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 140 not taken.
✗ Branch 141 not taken.
✗ Branch 142 not taken.
✗ Branch 143 not taken.
|
13 | ~VolumeRayIntersector() { if (mIsMaster) delete mTree; } |
| 359 | |||
| 360 | /// @brief Return @c false if the index ray misses the bbox of the grid. | ||
| 361 | /// @param iRay Ray represented in index space. | ||
| 362 | /// @warning Call this method (or setWorldRay) before the ray | ||
| 363 | /// traversal starts and use the return value to decide if further | ||
| 364 | /// marching is required. | ||
| 365 | 13 | inline bool setIndexRay(const RayT& iRay) | |
| 366 | { | ||
| 367 | 13 | mRay = iRay; | |
| 368 | 13 | const bool hit = mRay.clip(mBBox); | |
| 369 |
1/2✓ Branch 0 taken 13 times.
✗ Branch 1 not taken.
|
13 | if (hit) mTmax = mRay.t1(); |
| 370 | 13 | return hit; | |
| 371 | } | ||
| 372 | |||
| 373 | /// @brief Return @c false if the world ray misses the bbox of the grid. | ||
| 374 | /// @param wRay Ray represented in world space. | ||
| 375 | /// @warning Call this method (or setIndexRay) before the ray | ||
| 376 | /// traversal starts and use the return value to decide if further | ||
| 377 | /// marching is required. | ||
| 378 | /// @details Since hit times are computed with respect to the ray | ||
| 379 | /// represented in index space of the current grid, it is | ||
| 380 | /// recommended that either the client code uses getIndexPos to | ||
| 381 | /// compute index position from hit times or alternatively keeps | ||
| 382 | /// an instance of the index ray and instead uses setIndexRay to | ||
| 383 | /// initialize the ray. | ||
| 384 | ✗ | inline bool setWorldRay(const RayT& wRay) | |
| 385 | { | ||
| 386 | ✗ | return this->setIndexRay(wRay.worldToIndex(*mGrid)); | |
| 387 | } | ||
| 388 | |||
| 389 | 22 | inline typename RayT::TimeSpan march() | |
| 390 | { | ||
| 391 | 22 | const typename RayT::TimeSpan t = mHDDA.march(mRay, mAccessor); | |
| 392 |
3/4✓ Branch 0 taken 13 times.
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
|
22 | if (t.t1>0) mRay.setTimes(t.t1 + math::Delta<RealType>::value(), mTmax); |
| 393 | 22 | return t; | |
| 394 | } | ||
| 395 | |||
| 396 | /// @brief Return @c true if the ray intersects active values, | ||
| 397 | /// i.e. either active voxels or tiles. Only when a hit is | ||
| 398 | /// detected are t0 and t1 updated with the corresponding entry | ||
| 399 | /// and exit times along the INDEX ray! | ||
| 400 | /// @note Note that t0 and t1 are only resolved at the node level | ||
| 401 | /// (e.g. a LeafNode with active voxels) as opposed to the individual | ||
| 402 | /// active voxels. | ||
| 403 | /// @param t0 If the return value > 0 this is the time of the | ||
| 404 | /// first hit of an active tile or leaf. | ||
| 405 | /// @param t1 If the return value > t0 this is the time of the | ||
| 406 | /// first hit (> t0) of an inactive tile or exit point of the | ||
| 407 | /// BBOX for the leaf nodes. | ||
| 408 | /// @warning t0 and t1 are computed with respect to the ray represented in | ||
| 409 | /// index space of the current grid, not world space! | ||
| 410 | inline bool march(RealType& t0, RealType& t1) | ||
| 411 | { | ||
| 412 | 22 | const typename RayT::TimeSpan t = this->march(); | |
| 413 | t.get(t0, t1); | ||
| 414 | return t.valid(); | ||
| 415 | } | ||
| 416 | |||
| 417 | /// @brief Generates a list of hits along the ray. | ||
| 418 | /// | ||
| 419 | /// @param list List of hits represented as time spans. | ||
| 420 | /// | ||
| 421 | /// @note ListType is a list of RayType::TimeSpan and is required to | ||
| 422 | /// have the two methods: clear() and push_back(). Thus, it could | ||
| 423 | /// be std::vector<typename RayType::TimeSpan> or | ||
| 424 | /// std::deque<typename RayType::TimeSpan>. | ||
| 425 | template <typename ListType> | ||
| 426 | inline void hits(ListType& list) | ||
| 427 | { | ||
| 428 |
2/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
2 | mHDDA.hits(mRay, mAccessor, list); |
| 429 | 2 | } | |
| 430 | |||
| 431 | /// @brief Return the floating-point index position along the | ||
| 432 | /// current index ray at the specified time. | ||
| 433 | inline Vec3R getIndexPos(RealType time) const { return mRay(time); } | ||
| 434 | |||
| 435 | /// @brief Return the floating-point world position along the | ||
| 436 | /// current index ray at the specified time. | ||
| 437 | ✗ | inline Vec3R getWorldPos(RealType time) const { return mGrid->indexToWorld(mRay(time)); } | |
| 438 | |||
| 439 | inline RealType getWorldTime(RealType time) const | ||
| 440 | { | ||
| 441 | return time*mGrid->transform().baseMap()->applyJacobian(mRay.dir()).length(); | ||
| 442 | } | ||
| 443 | |||
| 444 | /// @brief Return a const reference to the input grid. | ||
| 445 | ✗ | const GridT& grid() const { return *mGrid; } | |
| 446 | |||
| 447 | /// @brief Return a const reference to the (potentially dilated) | ||
| 448 | /// bool tree used to accelerate the ray marching. | ||
| 449 | const TreeT& tree() const { return *mTree; } | ||
| 450 | |||
| 451 | /// @brief Return a const reference to the BBOX of the grid | ||
| 452 | const math::CoordBBox& bbox() const { return mBBox; } | ||
| 453 | |||
| 454 | /// @brief Print bbox, statistics, memory usage and other information. | ||
| 455 | /// @param os a stream to which to write textual information | ||
| 456 | /// @param verboseLevel 1: print bbox only; 2: include boolean tree | ||
| 457 | /// statistics; 3: include memory usage | ||
| 458 | ✗ | void print(std::ostream& os = std::cout, int verboseLevel = 1) | |
| 459 | { | ||
| 460 | ✗ | if (verboseLevel>0) { | |
| 461 | ✗ | os << "BBox: " << mBBox << std::endl; | |
| 462 | ✗ | if (verboseLevel==2) { | |
| 463 | ✗ | mTree->print(os, 1); | |
| 464 | ✗ | } else if (verboseLevel>2) { | |
| 465 | ✗ | mTree->print(os, 2); | |
| 466 | } | ||
| 467 | } | ||
| 468 | } | ||
| 469 | |||
| 470 | private: | ||
| 471 | using AccessorT = typename tree::ValueAccessor<const TreeT,/*IsSafe=*/false>; | ||
| 472 | |||
| 473 | const bool mIsMaster; | ||
| 474 | TreeT* mTree; | ||
| 475 | const GridT* mGrid; | ||
| 476 | AccessorT mAccessor; | ||
| 477 | RayT mRay; | ||
| 478 | RealType mTmax; | ||
| 479 | math::CoordBBox mBBox; | ||
| 480 | math::VolumeHDDA<TreeT, RayType, NodeLevel> mHDDA; | ||
| 481 | |||
| 482 | };// VolumeRayIntersector | ||
| 483 | |||
| 484 | |||
| 485 | //////////////////////////////////////// LinearSearchImpl //////////////////////////////////////// | ||
| 486 | |||
| 487 | |||
| 488 | /// @brief Implements linear iterative search for an iso-value of | ||
| 489 | /// the level set along the direction of the ray. | ||
| 490 | /// | ||
| 491 | /// @note Since this class is used internally in | ||
| 492 | /// LevelSetRayIntersector (define above) and LevelSetHDDA (defined below) | ||
| 493 | /// client code should never interact directly with its API. This also | ||
| 494 | /// explains why we are not concerned with the fact that several of | ||
| 495 | /// its methods are unsafe to call unless roots were already detected. | ||
| 496 | /// | ||
| 497 | /// @details It is approximate due to the limited number of iterations | ||
| 498 | /// which can can be defined with a template parameter. However the default value | ||
| 499 | /// has proven surprisingly accurate and fast. In fact more iterations | ||
| 500 | /// are not guaranteed to give significantly better results. | ||
| 501 | /// | ||
| 502 | /// @warning Since the root-searching algorithm is approximate | ||
| 503 | /// (first-order) it is possible to miss intersections if the | ||
| 504 | /// iso-value is too close to the inside or outside of the narrow | ||
| 505 | /// band (typically a distance less than a voxel unit). | ||
| 506 | /// | ||
| 507 | /// @warning Since this class internally stores a ValueAccessor it is NOT thread-safe, | ||
| 508 | /// so make sure to give each thread its own instance. This of course also means that | ||
| 509 | /// the cost of allocating an instance should (if possible) be amortized over | ||
| 510 | /// as many ray intersections as possible. | ||
| 511 | template<typename GridT, int Iterations, typename RealT> | ||
| 512 | ✗ | class LinearSearchImpl | |
| 513 | { | ||
| 514 | public: | ||
| 515 | using RayT = math::Ray<RealT>; | ||
| 516 | using VecT = math::Vec3<RealT>; | ||
| 517 | using ValueT = typename GridT::ValueType; | ||
| 518 | using AccessorT = typename GridT::ConstAccessor; | ||
| 519 | using StencilT = math::BoxStencil<GridT>; | ||
| 520 | |||
| 521 | /// @brief Constructor from a grid. | ||
| 522 | /// @throw RunTimeError if the grid is empty. | ||
| 523 | /// @throw ValueError if the isoValue is not inside the narrow-band. | ||
| 524 | 18 | LinearSearchImpl(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>()) | |
| 525 | : mStencil(grid), | ||
| 526 | mIsoValue(isoValue), | ||
| 527 | 18 | mMinValue(isoValue - ValueT(2 * grid.voxelSize()[0])), | |
| 528 |
2/4✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
|
54 | mMaxValue(isoValue + ValueT(2 * grid.voxelSize()[0])) |
| 529 | { | ||
| 530 |
2/4✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
|
18 | if ( grid.empty() ) { |
| 531 | ✗ | OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids"); | |
| 532 | } | ||
| 533 |
2/4✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
|
18 | if (mIsoValue<= -grid.background() || |
| 534 | mIsoValue>= grid.background() ){ | ||
| 535 | ✗ | OPENVDB_THROW(ValueError, "The iso-value must be inside the narrow-band!"); | |
| 536 | } | ||
| 537 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
18 | grid.tree().root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false); |
| 538 | 18 | } | |
| 539 | |||
| 540 | /// @brief Return the iso-value used for ray-intersections | ||
| 541 | const ValueT& getIsoValue() const { return mIsoValue; } | ||
| 542 | |||
| 543 | /// @brief Return @c false if the ray misses the bbox of the grid. | ||
| 544 | /// @param iRay Ray represented in index space. | ||
| 545 | /// @warning Call this method before the ray traversal starts. | ||
| 546 | inline bool setIndexRay(const RayT& iRay) | ||
| 547 | { | ||
| 548 | mRay = iRay; | ||
| 549 | return mRay.clip(mBBox);//did it hit the bbox | ||
| 550 | } | ||
| 551 | |||
| 552 | /// @brief Return @c false if the ray misses the bbox of the grid. | ||
| 553 | /// @param wRay Ray represented in world space. | ||
| 554 | /// @warning Call this method before the ray traversal starts. | ||
| 555 | 2097168 | inline bool setWorldRay(const RayT& wRay) | |
| 556 | { | ||
| 557 | 2097168 | mRay = wRay.worldToIndex(mStencil.grid()); | |
| 558 | 2097168 | return mRay.clip(mBBox);//did it hit the bbox | |
| 559 | } | ||
| 560 | |||
| 561 | /// @brief Get the intersection point in index space. | ||
| 562 | /// @param xyz The position in index space of the intersection. | ||
| 563 | ✗ | inline void getIndexPos(VecT& xyz) const { xyz = mRay(mTime); } | |
| 564 | |||
| 565 | /// @brief Get the intersection point in world space. | ||
| 566 | /// @param xyz The position in world space of the intersection. | ||
| 567 | 411738 | inline void getWorldPos(VecT& xyz) const { xyz = mStencil.grid().indexToWorld(mRay(mTime)); } | |
| 568 | |||
| 569 | /// @brief Get the intersection point and normal in world space | ||
| 570 | /// @param xyz The position in world space of the intersection. | ||
| 571 | /// @param nml The surface normal in world space of the intersection. | ||
| 572 | ✗ | inline void getWorldPosAndNml(VecT& xyz, VecT& nml) | |
| 573 | { | ||
| 574 | this->getIndexPos(xyz); | ||
| 575 | ✗ | mStencil.moveTo(xyz); | |
| 576 | ✗ | nml = mStencil.gradient(xyz); | |
| 577 | ✗ | nml.normalize(); | |
| 578 | ✗ | xyz = mStencil.grid().indexToWorld(xyz); | |
| 579 | } | ||
| 580 | |||
| 581 | /// @brief Return the time of intersection along the index ray. | ||
| 582 | inline RealT getIndexTime() const { return mTime; } | ||
| 583 | |||
| 584 | /// @brief Return the time of intersection along the world ray. | ||
| 585 | 411738 | inline RealT getWorldTime() const | |
| 586 | { | ||
| 587 |
2/4✓ Branch 2 taken 205869 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 205869 times.
✗ Branch 5 not taken.
|
823476 | return mTime*mStencil.grid().transform().baseMap()->applyJacobian(mRay.dir()).length(); |
| 588 | } | ||
| 589 | |||
| 590 | private: | ||
| 591 | |||
| 592 | /// @brief Initiate the local voxel intersection test. | ||
| 593 | /// @warning Make sure to call this method before the local voxel intersection test. | ||
| 594 | inline void init(RealT t0) | ||
| 595 | { | ||
| 596 | 537330 | mT[0] = t0; | |
| 597 | 537330 | mV[0] = static_cast<ValueT>(this->interpValue(t0)); | |
| 598 | 537330 | } | |
| 599 | |||
| 600 | inline void setRange(RealT t0, RealT t1) { mRay.setTimes(t0, t1); } | ||
| 601 | |||
| 602 | /// @brief Return a const reference to the ray. | ||
| 603 | 1585724 | inline const RayT& ray() const { return mRay; } | |
| 604 | |||
| 605 | /// @brief Return true if a node of the specified type exists at ijk. | ||
| 606 | template <typename NodeT> | ||
| 607 | inline bool hasNode(const Coord& ijk) | ||
| 608 | { | ||
| 609 | 7670073 | return mStencil.accessor().template probeConstNode<NodeT>(ijk) != nullptr; | |
| 610 | } | ||
| 611 | |||
| 612 | /// @brief Return @c true if an intersection is detected. | ||
| 613 | /// @param ijk Grid coordinate of the node origin or voxel being tested. | ||
| 614 | /// @param time Time along the index ray being tested. | ||
| 615 | /// @warning Only if an intersection is detected is it safe to | ||
| 616 | /// call getIndexPos, getWorldPos and getWorldPosAndNml! | ||
| 617 | 7285620 | inline bool operator()(const Coord& ijk, RealT time) | |
| 618 | { | ||
| 619 | ValueT V; | ||
| 620 | 7285620 | if (mStencil.accessor().probeValue(ijk, V) &&//within narrow band | |
| 621 |
6/6✓ Branch 0 taken 985768 times.
✓ Branch 1 taken 2657042 times.
✓ Branch 2 taken 985767 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 985764 times.
✓ Branch 5 taken 3 times.
|
7285620 | V>mMinValue && V<mMaxValue) {// and close to iso-value? |
| 622 | 1971528 | mT[1] = time; | |
| 623 |
2/2✓ Branch 1 taken 205869 times.
✓ Branch 2 taken 779895 times.
|
1971528 | mV[1] = static_cast<ValueT>(this->interpValue(time)); |
| 624 |
2/2✓ Branch 0 taken 205869 times.
✓ Branch 1 taken 779895 times.
|
1971528 | if (math::ZeroCrossing(mV[0], mV[1])) { |
| 625 | 411738 | mTime = this->interpTime(); | |
| 626 | OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN | ||
| 627 |
2/2✓ Branch 0 taken 411722 times.
✓ Branch 1 taken 205861 times.
|
1235166 | for (int n=0; Iterations>0 && n<Iterations; ++n) {//resolved at compile-time |
| 628 |
2/2✓ Branch 1 taken 207730 times.
✓ Branch 2 taken 203992 times.
|
823444 | V = static_cast<ValueT>(this->interpValue(mTime)); |
| 629 |
2/2✓ Branch 0 taken 207730 times.
✓ Branch 1 taken 203992 times.
|
823444 | const int m = math::ZeroCrossing(mV[0], V) ? 1 : 0; |
| 630 | 823444 | mV[m] = V; | |
| 631 | 823444 | mT[m] = mTime; | |
| 632 | 823444 | mTime = this->interpTime(); | |
| 633 | } | ||
| 634 | OPENVDB_NO_UNREACHABLE_CODE_WARNING_END | ||
| 635 | 16 | return true; | |
| 636 | } | ||
| 637 | 1559790 | mT[0] = mT[1]; | |
| 638 | 1559790 | mV[0] = mV[1]; | |
| 639 | } | ||
| 640 | return false; | ||
| 641 | } | ||
| 642 | |||
| 643 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 617591 times.
|
1235182 | inline RealT interpTime() |
| 644 | { | ||
| 645 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 617591 times.
|
1235182 | assert( math::isApproxLarger(mT[1], mT[0], RealT(1e-6) ) ); |
| 646 | 1235182 | return mT[0]+(mT[1]-mT[0])*mV[0]/(mV[0]-mV[1]); | |
| 647 | } | ||
| 648 | |||
| 649 | 1934816 | inline RealT interpValue(RealT time) | |
| 650 | { | ||
| 651 | 1934816 | const VecT pos = mRay(time); | |
| 652 | 1934816 | mStencil.moveTo(pos); | |
| 653 | 1934816 | return mStencil.interpolation(pos) - mIsoValue; | |
| 654 | } | ||
| 655 | |||
| 656 | template<typename, int> friend struct math::LevelSetHDDA; | ||
| 657 | |||
| 658 | RayT mRay; | ||
| 659 | StencilT mStencil; | ||
| 660 | RealT mTime;//time of intersection | ||
| 661 | ValueT mV[2]; | ||
| 662 | RealT mT[2]; | ||
| 663 | const ValueT mIsoValue, mMinValue, mMaxValue; | ||
| 664 | math::CoordBBox mBBox; | ||
| 665 | };// LinearSearchImpl | ||
| 666 | |||
| 667 | } // namespace tools | ||
| 668 | } // namespace OPENVDB_VERSION_NAME | ||
| 669 | } // namespace openvdb | ||
| 670 | |||
| 671 | #endif // OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED | ||
| 672 |