| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // Copyright Contributors to the OpenVDB Project | ||
| 2 | // SPDX-License-Identifier: MPL-2.0 | ||
| 3 | |||
| 4 | /// @author Ken Museth | ||
| 5 | /// | ||
| 6 | /// @file tools/LevelSetMorph.h | ||
| 7 | /// | ||
| 8 | /// @brief Shape morphology of level sets. Morphing from a source | ||
| 9 | /// narrow-band level sets to a target narrow-band level set. | ||
| 10 | |||
| 11 | #ifndef OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED | ||
| 12 | #define OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED | ||
| 13 | |||
| 14 | #include "LevelSetTracker.h" | ||
| 15 | #include "Interpolation.h" // for BoxSampler, etc. | ||
| 16 | #include <openvdb/math/FiniteDifference.h> | ||
| 17 | #include <functional> | ||
| 18 | #include <limits> | ||
| 19 | |||
| 20 | |||
| 21 | namespace openvdb { | ||
| 22 | OPENVDB_USE_VERSION_NAMESPACE | ||
| 23 | namespace OPENVDB_VERSION_NAME { | ||
| 24 | namespace tools { | ||
| 25 | |||
| 26 | /// @brief Shape morphology of level sets. Morphing from a source | ||
| 27 | /// narrow-band level sets to a target narrow-band level set. | ||
| 28 | /// | ||
| 29 | /// @details | ||
| 30 | /// The @c InterruptType template argument below refers to any class | ||
| 31 | /// with the following interface: | ||
| 32 | /// @code | ||
| 33 | /// class Interrupter { | ||
| 34 | /// ... | ||
| 35 | /// public: | ||
| 36 | /// void start(const char* name = nullptr) // called when computations begin | ||
| 37 | /// void end() // called when computations end | ||
| 38 | /// bool wasInterrupted(int percent=-1) // return true to break computation | ||
| 39 | /// }; | ||
| 40 | /// @endcode | ||
| 41 | /// | ||
| 42 | /// @note If no template argument is provided for this InterruptType, | ||
| 43 | /// the util::NullInterrupter is used, which implies that all interrupter | ||
| 44 | /// calls are no-ops (i.e., they incur no computational overhead). | ||
| 45 | template<typename GridT, | ||
| 46 | typename InterruptT = util::NullInterrupter> | ||
| 47 | class LevelSetMorphing | ||
| 48 | { | ||
| 49 | public: | ||
| 50 | using GridType = GridT; | ||
| 51 | using TreeType = typename GridT::TreeType; | ||
| 52 | using TrackerT = LevelSetTracker<GridT, InterruptT>; | ||
| 53 | using LeafRange = typename TrackerT::LeafRange; | ||
| 54 | using LeafType = typename TrackerT::LeafType; | ||
| 55 | using BufferType = typename TrackerT::BufferType; | ||
| 56 | using ValueType = typename TrackerT::ValueType; | ||
| 57 | |||
| 58 | /// Main constructor | ||
| 59 | 1 | LevelSetMorphing(GridT& sourceGrid, const GridT& targetGrid, InterruptT* interrupt = nullptr) | |
| 60 | : mTracker(sourceGrid, interrupt) | ||
| 61 | , mTarget(&targetGrid) | ||
| 62 | , mMask(nullptr) | ||
| 63 | , mSpatialScheme(math::HJWENO5_BIAS) | ||
| 64 | , mTemporalScheme(math::TVD_RK2) | ||
| 65 | , mMinMask(0) | ||
| 66 | , mDeltaMask(1) | ||
| 67 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | , mInvertMask(false) |
| 68 | { | ||
| 69 | } | ||
| 70 | |||
| 71 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | virtual ~LevelSetMorphing() {} |
| 72 | |||
| 73 | /// Redefine the target level set | ||
| 74 | ✗ | void setTarget(const GridT& targetGrid) { mTarget = &targetGrid; } | |
| 75 | |||
| 76 | /// Define the alpha mask | ||
| 77 | ✗ | void setAlphaMask(const GridT& maskGrid) { mMask = &maskGrid; } | |
| 78 | |||
| 79 | /// Return the spatial finite-difference scheme | ||
| 80 | ✗ | math::BiasedGradientScheme getSpatialScheme() const { return mSpatialScheme; } | |
| 81 | /// Set the spatial finite-difference scheme | ||
| 82 | ✗ | void setSpatialScheme(math::BiasedGradientScheme scheme) { mSpatialScheme = scheme; } | |
| 83 | |||
| 84 | /// Return the temporal integration scheme | ||
| 85 | ✗ | math::TemporalIntegrationScheme getTemporalScheme() const { return mTemporalScheme; } | |
| 86 | /// Set the temporal integration scheme | ||
| 87 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | void setTemporalScheme(math::TemporalIntegrationScheme scheme) { mTemporalScheme = scheme; } |
| 88 | |||
| 89 | /// Return the spatial finite-difference scheme | ||
| 90 | ✗ | math::BiasedGradientScheme getTrackerSpatialScheme() const | |
| 91 | { | ||
| 92 | ✗ | return mTracker.getSpatialScheme(); | |
| 93 | } | ||
| 94 | /// Set the spatial finite-difference scheme | ||
| 95 | ✗ | void setTrackerSpatialScheme(math::BiasedGradientScheme scheme) | |
| 96 | { | ||
| 97 | mTracker.setSpatialScheme(scheme); | ||
| 98 | } | ||
| 99 | /// Return the temporal integration scheme | ||
| 100 | ✗ | math::TemporalIntegrationScheme getTrackerTemporalScheme() const | |
| 101 | { | ||
| 102 | ✗ | return mTracker.getTemporalScheme(); | |
| 103 | } | ||
| 104 | /// Set the temporal integration scheme | ||
| 105 | ✗ | void setTrackerTemporalScheme(math::TemporalIntegrationScheme scheme) | |
| 106 | { | ||
| 107 | mTracker.setTemporalScheme(scheme); | ||
| 108 | } | ||
| 109 | /// Return the number of normalizations performed per track or normalize call. | ||
| 110 | ✗ | int getNormCount() const { return mTracker.getNormCount(); } | |
| 111 | /// Set the number of normalizations performed per track or normalize call. | ||
| 112 | ✗ | void setNormCount(int n) { mTracker.setNormCount(n); } | |
| 113 | |||
| 114 | /// Return the grain size used for multithreading | ||
| 115 | ✗ | int getGrainSize() const { return mTracker.getGrainSize(); } | |
| 116 | /// @brief Set the grain size used for multithreading. | ||
| 117 | /// @note A grain size of 0 or less disables multithreading! | ||
| 118 | ✗ | void setGrainSize(int grainsize) { mTracker.setGrainSize(grainsize); } | |
| 119 | |||
| 120 | /// @brief Return the minimum value of the mask to be used for the | ||
| 121 | /// derivation of a smooth alpha value. | ||
| 122 | ✗ | ValueType minMask() const { return mMinMask; } | |
| 123 | |||
| 124 | /// @brief Return the maximum value of the mask to be used for the | ||
| 125 | /// derivation of a smooth alpha value. | ||
| 126 | ✗ | ValueType maxMask() const { return mDeltaMask + mMinMask; } | |
| 127 | |||
| 128 | /// @brief Define the range for the (optional) scalar mask. | ||
| 129 | /// @param min Minimum value of the range. | ||
| 130 | /// @param max Maximum value of the range. | ||
| 131 | /// @details Mask values outside the range maps to alpha values of | ||
| 132 | /// respectfully zero and one, and values inside the range maps | ||
| 133 | /// smoothly to 0->1 (unless of course the mask is inverted). | ||
| 134 | /// @throw ValueError if @a min is not smaller than @a max. | ||
| 135 | ✗ | void setMaskRange(ValueType min, ValueType max) | |
| 136 | { | ||
| 137 | ✗ | if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)"); | |
| 138 | ✗ | mMinMask = min; | |
| 139 | ✗ | mDeltaMask = max-min; | |
| 140 | } | ||
| 141 | |||
| 142 | /// @brief Return true if the mask is inverted, i.e. min->max in the | ||
| 143 | /// original mask maps to 1->0 in the inverted alpha mask. | ||
| 144 | ✗ | bool isMaskInverted() const { return mInvertMask; } | |
| 145 | /// @brief Invert the optional mask, i.e. min->max in the original | ||
| 146 | /// mask maps to 1->0 in the inverted alpha mask. | ||
| 147 | ✗ | void invertMask(bool invert=true) { mInvertMask = invert; } | |
| 148 | |||
| 149 | /// @brief Advect the level set from its current time, @a time0, to its | ||
| 150 | /// final time, @a time1. If @a time0 > @a time1, perform backward advection. | ||
| 151 | /// | ||
| 152 | /// @return the number of CFL iterations used to advect from @a time0 to @a time1 | ||
| 153 | size_t advect(ValueType time0, ValueType time1); | ||
| 154 | |||
| 155 | private: | ||
| 156 | |||
| 157 | // disallow copy construction and copy by assignment! | ||
| 158 | LevelSetMorphing(const LevelSetMorphing&);// not implemented | ||
| 159 | LevelSetMorphing& operator=(const LevelSetMorphing&);// not implemented | ||
| 160 | |||
| 161 | template<math::BiasedGradientScheme SpatialScheme> | ||
| 162 | size_t advect1(ValueType time0, ValueType time1); | ||
| 163 | |||
| 164 | template<math::BiasedGradientScheme SpatialScheme, | ||
| 165 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 166 | size_t advect2(ValueType time0, ValueType time1); | ||
| 167 | |||
| 168 | template<math::BiasedGradientScheme SpatialScheme, | ||
| 169 | math::TemporalIntegrationScheme TemporalScheme, | ||
| 170 | typename MapType> | ||
| 171 | size_t advect3(ValueType time0, ValueType time1); | ||
| 172 | |||
| 173 | TrackerT mTracker; | ||
| 174 | const GridT *mTarget, *mMask; | ||
| 175 | math::BiasedGradientScheme mSpatialScheme; | ||
| 176 | math::TemporalIntegrationScheme mTemporalScheme; | ||
| 177 | ValueType mMinMask, mDeltaMask; | ||
| 178 | bool mInvertMask; | ||
| 179 | |||
| 180 | // This templated private class implements all the level set magic. | ||
| 181 | template<typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 182 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 183 | struct Morph | ||
| 184 | { | ||
| 185 | /// Main constructor | ||
| 186 | Morph(LevelSetMorphing<GridT, InterruptT>& parent); | ||
| 187 | /// Shallow copy constructor called by tbb::parallel_for() threads | ||
| 188 | Morph(const Morph& other); | ||
| 189 | /// Shallow copy constructor called by tbb::parallel_reduce() threads | ||
| 190 | Morph(Morph& other, tbb::split); | ||
| 191 | /// destructor | ||
| 192 |
2/288✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✓ Branch 60 taken 32 times.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
✗ Branch 120 not taken.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 132 not taken.
✗ Branch 133 not taken.
✗ Branch 134 not taken.
✗ Branch 135 not taken.
✗ Branch 136 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 139 not taken.
✗ Branch 140 not taken.
✗ Branch 141 not taken.
✗ Branch 142 not taken.
✗ Branch 143 not taken.
✗ Branch 144 not taken.
✗ Branch 145 not taken.
✗ Branch 146 not taken.
✗ Branch 147 not taken.
✗ Branch 148 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✗ Branch 151 not taken.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 154 not taken.
✗ Branch 155 not taken.
✗ Branch 156 not taken.
✗ Branch 157 not taken.
✗ Branch 158 not taken.
✗ Branch 159 not taken.
✗ Branch 160 not taken.
✗ Branch 161 not taken.
✗ Branch 162 not taken.
✗ Branch 163 not taken.
✗ Branch 164 not taken.
✗ Branch 165 not taken.
✗ Branch 166 not taken.
✗ Branch 167 not taken.
✗ Branch 168 not taken.
✗ Branch 169 not taken.
✗ Branch 170 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 176 not taken.
✗ Branch 177 not taken.
✗ Branch 178 not taken.
✗ Branch 179 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 182 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✗ Branch 185 not taken.
✗ Branch 186 not taken.
✗ Branch 187 not taken.
✗ Branch 188 not taken.
✗ Branch 189 not taken.
✗ Branch 190 not taken.
✗ Branch 191 not taken.
✗ Branch 192 not taken.
✗ Branch 193 not taken.
✗ Branch 194 not taken.
✗ Branch 195 not taken.
✗ Branch 196 not taken.
✗ Branch 197 not taken.
✗ Branch 198 not taken.
✗ Branch 199 not taken.
✗ Branch 200 not taken.
✗ Branch 201 not taken.
✗ Branch 202 not taken.
✗ Branch 203 not taken.
✗ Branch 204 not taken.
✗ Branch 205 not taken.
✗ Branch 206 not taken.
✗ Branch 207 not taken.
✗ Branch 208 not taken.
✗ Branch 209 not taken.
✗ Branch 210 not taken.
✗ Branch 211 not taken.
✗ Branch 212 not taken.
✗ Branch 213 not taken.
✗ Branch 214 not taken.
✗ Branch 215 not taken.
✗ Branch 216 not taken.
✗ Branch 217 not taken.
✗ Branch 218 not taken.
✗ Branch 219 not taken.
✗ Branch 220 not taken.
✗ Branch 221 not taken.
✗ Branch 222 not taken.
✗ Branch 223 not taken.
✗ Branch 224 not taken.
✗ Branch 225 not taken.
✗ Branch 226 not taken.
✗ Branch 227 not taken.
✗ Branch 228 not taken.
✗ Branch 229 not taken.
✗ Branch 230 not taken.
✗ Branch 231 not taken.
✗ Branch 232 not taken.
✗ Branch 233 not taken.
✗ Branch 234 not taken.
✗ Branch 235 not taken.
✗ Branch 236 not taken.
✗ Branch 237 not taken.
✗ Branch 238 not taken.
✗ Branch 239 not taken.
✗ Branch 240 not taken.
✗ Branch 241 not taken.
✗ Branch 242 not taken.
✗ Branch 243 not taken.
✗ Branch 244 not taken.
✗ Branch 245 not taken.
✗ Branch 246 not taken.
✗ Branch 247 not taken.
✗ Branch 248 not taken.
✗ Branch 249 not taken.
✗ Branch 250 not taken.
✗ Branch 251 not taken.
✗ Branch 252 not taken.
✗ Branch 253 not taken.
✗ Branch 254 not taken.
✗ Branch 255 not taken.
✗ Branch 256 not taken.
✗ Branch 257 not taken.
✗ Branch 258 not taken.
✗ Branch 259 not taken.
✗ Branch 260 not taken.
✗ Branch 261 not taken.
✗ Branch 262 not taken.
✗ Branch 263 not taken.
✗ Branch 264 not taken.
✗ Branch 265 not taken.
✗ Branch 266 not taken.
✗ Branch 267 not taken.
✗ Branch 268 not taken.
✗ Branch 269 not taken.
✗ Branch 270 not taken.
✗ Branch 271 not taken.
✗ Branch 272 not taken.
✗ Branch 273 not taken.
✗ Branch 274 not taken.
✗ Branch 275 not taken.
✗ Branch 276 not taken.
✗ Branch 277 not taken.
✗ Branch 278 not taken.
✗ Branch 279 not taken.
✗ Branch 280 not taken.
✗ Branch 281 not taken.
✗ Branch 282 not taken.
✗ Branch 283 not taken.
✗ Branch 284 not taken.
✗ Branch 285 not taken.
✗ Branch 286 not taken.
✗ Branch 287 not taken.
|
160 | virtual ~Morph() {} |
| 193 | /// Advect the level set from its current time, time0, to its final time, time1. | ||
| 194 | /// @return number of CFL iterations | ||
| 195 | size_t advect(ValueType time0, ValueType time1); | ||
| 196 | /// Used internally by tbb::parallel_for() | ||
| 197 |
1/2✓ Branch 0 taken 4866 times.
✗ Branch 1 not taken.
|
9732 | void operator()(const LeafRange& r) const |
| 198 | { | ||
| 199 |
1/2✓ Branch 0 taken 4866 times.
✗ Branch 1 not taken.
|
9732 | if (mTask) mTask(const_cast<Morph*>(this), r); |
| 200 | ✗ | else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly"); | |
| 201 | } | ||
| 202 | /// Used internally by tbb::parallel_reduce() | ||
| 203 |
1/2✓ Branch 0 taken 1622 times.
✗ Branch 1 not taken.
|
3244 | void operator()(const LeafRange& r) |
| 204 | { | ||
| 205 |
1/2✓ Branch 0 taken 1622 times.
✗ Branch 1 not taken.
|
3244 | if (mTask) mTask(this, r); |
| 206 | ✗ | else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly"); | |
| 207 | } | ||
| 208 | /// This is only called by tbb::parallel_reduce() threads | ||
| 209 |
2/96✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 35 times.
✓ Branch 7 taken 29 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
|
99 | void join(const Morph& other) { mMaxAbsS = math::Max(mMaxAbsS, other.mMaxAbsS); } |
| 210 | |||
| 211 | /// Enum to define the type of multithreading | ||
| 212 | enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE }; // for internal use | ||
| 213 | // method calling tbb | ||
| 214 | void cook(ThreadingMode mode, size_t swapBuffer = 0); | ||
| 215 | |||
| 216 | /// Sample field and return the CFT time step | ||
| 217 | typename GridT::ValueType sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer); | ||
| 218 | void sampleXformedSpeed(const LeafRange& r, Index speedBuffer); | ||
| 219 | void sampleAlignedSpeed(const LeafRange& r, Index speedBuffer); | ||
| 220 | |||
| 221 | // Convex combination of Phi and a forward Euler advection steps: | ||
| 222 | // Phi(result) = alpha * Phi(phi) + (1-alpha) * (Phi(0) - dt * Speed(speed)*|Grad[Phi(0)]|); | ||
| 223 | template <int Nominator, int Denominator> | ||
| 224 | void euler(const LeafRange&, ValueType, Index, Index, Index); | ||
| 225 | 3244 | inline void euler01(const LeafRange& r, ValueType t, Index s) {this->euler<0,1>(r,t,0,1,s);} | |
| 226 | ✗ | inline void euler12(const LeafRange& r, ValueType t) {this->euler<1,2>(r, t, 1, 1, 2);} | |
| 227 | 3244 | inline void euler34(const LeafRange& r, ValueType t) {this->euler<3,4>(r, t, 1, 2, 3);} | |
| 228 | 3244 | inline void euler13(const LeafRange& r, ValueType t) {this->euler<1,3>(r, t, 1, 2, 3);} | |
| 229 | |||
| 230 | using FuncType = typename std::function<void (Morph*, const LeafRange&)>; | ||
| 231 | LevelSetMorphing* mParent; | ||
| 232 | ValueType mMinAbsS, mMaxAbsS; | ||
| 233 | const MapT* mMap; | ||
| 234 | FuncType mTask; | ||
| 235 | }; // end of private Morph struct | ||
| 236 | |||
| 237 | };//end of LevelSetMorphing | ||
| 238 | |||
| 239 | template<typename GridT, typename InterruptT> | ||
| 240 | inline size_t | ||
| 241 | ✗ | LevelSetMorphing<GridT, InterruptT>::advect(ValueType time0, ValueType time1) | |
| 242 | { | ||
| 243 | ✗ | switch (mSpatialScheme) { | |
| 244 | ✗ | case math::FIRST_BIAS: | |
| 245 | ✗ | return this->advect1<math::FIRST_BIAS >(time0, time1); | |
| 246 | //case math::SECOND_BIAS: | ||
| 247 | //return this->advect1<math::SECOND_BIAS >(time0, time1); | ||
| 248 | //case math::THIRD_BIAS: | ||
| 249 | //return this->advect1<math::THIRD_BIAS >(time0, time1); | ||
| 250 | //case math::WENO5_BIAS: | ||
| 251 | //return this->advect1<math::WENO5_BIAS >(time0, time1); | ||
| 252 | ✗ | case math::HJWENO5_BIAS: | |
| 253 | ✗ | return this->advect1<math::HJWENO5_BIAS>(time0, time1); | |
| 254 | case math::SECOND_BIAS: | ||
| 255 | case math::THIRD_BIAS: | ||
| 256 | case math::WENO5_BIAS: | ||
| 257 | case math::UNKNOWN_BIAS: | ||
| 258 | default: | ||
| 259 | ✗ | OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!"); | |
| 260 | } | ||
| 261 | return 0; | ||
| 262 | } | ||
| 263 | |||
| 264 | template<typename GridT, typename InterruptT> | ||
| 265 | template<math::BiasedGradientScheme SpatialScheme> | ||
| 266 | inline size_t | ||
| 267 | 64 | LevelSetMorphing<GridT, InterruptT>::advect1(ValueType time0, ValueType time1) | |
| 268 | { | ||
| 269 |
1/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
|
64 | switch (mTemporalScheme) { |
| 270 | ✗ | case math::TVD_RK1: | |
| 271 | ✗ | return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1); | |
| 272 | ✗ | case math::TVD_RK2: | |
| 273 | ✗ | return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1); | |
| 274 | 64 | case math::TVD_RK3: | |
| 275 | 64 | return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1); | |
| 276 | case math::UNKNOWN_TIS: | ||
| 277 | default: | ||
| 278 | ✗ | OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!"); | |
| 279 | } | ||
| 280 | return 0; | ||
| 281 | } | ||
| 282 | |||
| 283 | template<typename GridT, typename InterruptT> | ||
| 284 | template<math::BiasedGradientScheme SpatialScheme, | ||
| 285 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 286 | inline size_t | ||
| 287 | 64 | LevelSetMorphing<GridT, InterruptT>::advect2(ValueType time0, ValueType time1) | |
| 288 | { | ||
| 289 | const math::Transform& trans = mTracker.grid().transform(); | ||
| 290 |
1/2✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
|
64 | if (trans.mapType() == math::UniformScaleMap::mapType()) { |
| 291 | 64 | return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1); | |
| 292 | ✗ | } else if (trans.mapType() == math::UniformScaleTranslateMap::mapType()) { | |
| 293 | return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>( | ||
| 294 | ✗ | time0, time1); | |
| 295 | ✗ | } else if (trans.mapType() == math::UnitaryMap::mapType()) { | |
| 296 | ✗ | return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1); | |
| 297 | ✗ | } else if (trans.mapType() == math::TranslationMap::mapType()) { | |
| 298 | ✗ | return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1); | |
| 299 | } else { | ||
| 300 | ✗ | OPENVDB_THROW(ValueError, "MapType not supported!"); | |
| 301 | } | ||
| 302 | return 0; | ||
| 303 | } | ||
| 304 | |||
| 305 | template<typename GridT, typename InterruptT> | ||
| 306 | template<math::BiasedGradientScheme SpatialScheme, | ||
| 307 | math::TemporalIntegrationScheme TemporalScheme, | ||
| 308 | typename MapT> | ||
| 309 | inline size_t | ||
| 310 | 64 | LevelSetMorphing<GridT, InterruptT>::advect3(ValueType time0, ValueType time1) | |
| 311 | { | ||
| 312 | 64 | Morph<MapT, SpatialScheme, TemporalScheme> tmp(*this); | |
| 313 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
128 | return tmp.advect(time0, time1); |
| 314 | } | ||
| 315 | |||
| 316 | |||
| 317 | /////////////////////////////////////////////////////////////////////// | ||
| 318 | |||
| 319 | template<typename GridT, typename InterruptT> | ||
| 320 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 321 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 322 | inline | ||
| 323 | 64 | LevelSetMorphing<GridT, InterruptT>:: | |
| 324 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 325 | Morph(LevelSetMorphing<GridT, InterruptT>& parent) | ||
| 326 | : mParent(&parent) | ||
| 327 | , mMinAbsS(ValueType(1e-6)) | ||
| 328 | , mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()) | ||
| 329 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
64 | , mTask(nullptr) |
| 330 | { | ||
| 331 | } | ||
| 332 | |||
| 333 | template<typename GridT, typename InterruptT> | ||
| 334 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 335 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 336 | inline | ||
| 337 | 1947 | LevelSetMorphing<GridT, InterruptT>:: | |
| 338 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 339 | Morph(const Morph& other) | ||
| 340 | 1947 | : mParent(other.mParent) | |
| 341 | 1947 | , mMinAbsS(other.mMinAbsS) | |
| 342 | 1947 | , mMaxAbsS(other.mMaxAbsS) | |
| 343 | 1947 | , mMap(other.mMap) | |
| 344 | 1947 | , mTask(other.mTask) | |
| 345 | { | ||
| 346 | } | ||
| 347 | |||
| 348 | template<typename GridT, typename InterruptT> | ||
| 349 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 350 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 351 | inline | ||
| 352 | 64 | LevelSetMorphing<GridT, InterruptT>:: | |
| 353 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 354 | Morph(Morph& other, tbb::split) | ||
| 355 | 64 | : mParent(other.mParent) | |
| 356 | 64 | , mMinAbsS(other.mMinAbsS) | |
| 357 | 64 | , mMaxAbsS(other.mMaxAbsS) | |
| 358 | 64 | , mMap(other.mMap) | |
| 359 | 64 | , mTask(other.mTask) | |
| 360 | { | ||
| 361 | } | ||
| 362 | |||
| 363 | template<typename GridT, typename InterruptT> | ||
| 364 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 365 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 366 | inline size_t | ||
| 367 | 64 | LevelSetMorphing<GridT, InterruptT>:: | |
| 368 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 369 | advect(ValueType time0, ValueType time1) | ||
| 370 | { | ||
| 371 | namespace ph = std::placeholders; | ||
| 372 | |||
| 373 | // Make sure we have enough temporal auxiliary buffers for the time | ||
| 374 | // integration AS WELL AS an extra buffer with the speed function! | ||
| 375 | static const Index auxBuffers = 1 + (TemporalScheme == math::TVD_RK3 ? 2 : 1); | ||
| 376 | size_t countCFL = 0; | ||
| 377 |
3/4✓ Branch 0 taken 32 times.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
|
128 | while (time0 < time1 && mParent->mTracker.checkInterrupter()) { |
| 378 | 64 | mParent->mTracker.leafs().rebuildAuxBuffers(auxBuffers); | |
| 379 | |||
| 380 | 64 | const ValueType dt = this->sampleSpeed(time0, time1, auxBuffers); | |
| 381 |
1/2✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
|
64 | if ( math::isZero(dt) ) break;//V is essentially zero so terminate |
| 382 | |||
| 383 | OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN //switch is resolved at compile-time | ||
| 384 | switch(TemporalScheme) { | ||
| 385 | ✗ | case math::TVD_RK1: | |
| 386 | // Perform one explicit Euler step: t1 = t0 + dt | ||
| 387 | // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]| | ||
| 388 | ✗ | mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2); | |
| 389 | |||
| 390 | // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1) | ||
| 391 | ✗ | this->cook(PARALLEL_FOR, 1); | |
| 392 | break; | ||
| 393 | ✗ | case math::TVD_RK2: | |
| 394 | // Perform one explicit Euler step: t1 = t0 + dt | ||
| 395 | // Phi_t1(1) = Phi_t0(0) - dt * Speed(2) * |Grad[Phi(0)]| | ||
| 396 | ✗ | mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/2); | |
| 397 | |||
| 398 | // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1) | ||
| 399 | ✗ | this->cook(PARALLEL_FOR, 1); | |
| 400 | |||
| 401 | // Convex combine explict Euler step: t2 = t0 + dt | ||
| 402 | // Phi_t2(1) = 1/2 * Phi_t0(1) + 1/2 * (Phi_t1(0) - dt * Speed(2) * |Grad[Phi(0)]|) | ||
| 403 | ✗ | mTask = std::bind(&Morph::euler12, ph::_1, ph::_2, dt); | |
| 404 | |||
| 405 | // Cook and swap buffer 0 and 1 such that Phi_t2(0) and Phi_t1(1) | ||
| 406 | ✗ | this->cook(PARALLEL_FOR, 1); | |
| 407 | break; | ||
| 408 | 64 | case math::TVD_RK3: | |
| 409 | // Perform one explicit Euler step: t1 = t0 + dt | ||
| 410 | // Phi_t1(1) = Phi_t0(0) - dt * Speed(3) * |Grad[Phi(0)]| | ||
| 411 | 64 | mTask = std::bind(&Morph::euler01, ph::_1, ph::_2, dt, /*speed*/3); | |
| 412 | |||
| 413 | // Cook and swap buffer 0 and 1 such that Phi_t1(0) and Phi_t0(1) | ||
| 414 | 64 | this->cook(PARALLEL_FOR, 1); | |
| 415 | |||
| 416 | // Convex combine explict Euler step: t2 = t0 + dt/2 | ||
| 417 | // Phi_t2(2) = 3/4 * Phi_t0(1) + 1/4 * (Phi_t1(0) - dt * Speed(3) * |Grad[Phi(0)]|) | ||
| 418 | 64 | mTask = std::bind(&Morph::euler34, ph::_1, ph::_2, dt); | |
| 419 | |||
| 420 | // Cook and swap buffer 0 and 2 such that Phi_t2(0) and Phi_t1(2) | ||
| 421 | 64 | this->cook(PARALLEL_FOR, 2); | |
| 422 | |||
| 423 | // Convex combine explict Euler step: t3 = t0 + dt | ||
| 424 | // Phi_t3(2) = 1/3 * Phi_t0(1) + 2/3 * (Phi_t2(0) - dt * Speed(3) * |Grad[Phi(0)]|) | ||
| 425 | 64 | mTask = std::bind(&Morph::euler13, ph::_1, ph::_2, dt); | |
| 426 | |||
| 427 | // Cook and swap buffer 0 and 2 such that Phi_t3(0) and Phi_t2(2) | ||
| 428 | 64 | this->cook(PARALLEL_FOR, 2); | |
| 429 | break; | ||
| 430 | case math::UNKNOWN_TIS: | ||
| 431 | default: | ||
| 432 | OPENVDB_THROW(ValueError, "Temporal integration scheme not supported!"); | ||
| 433 | }//end of compile-time resolved switch | ||
| 434 | OPENVDB_NO_UNREACHABLE_CODE_WARNING_END | ||
| 435 | |||
| 436 | 64 | time0 += dt; | |
| 437 | 64 | ++countCFL; | |
| 438 | 64 | mParent->mTracker.leafs().removeAuxBuffers(); | |
| 439 | |||
| 440 | // Track the narrow band | ||
| 441 | 64 | mParent->mTracker.track(); | |
| 442 | }//end wile-loop over time | ||
| 443 | |||
| 444 | 64 | return countCFL;//number of CLF propagation steps | |
| 445 | } | ||
| 446 | |||
| 447 | template<typename GridT, typename InterruptT> | ||
| 448 | template<typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 449 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 450 | inline typename GridT::ValueType | ||
| 451 | 64 | LevelSetMorphing<GridT, InterruptT>:: | |
| 452 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 453 | sampleSpeed(ValueType time0, ValueType time1, Index speedBuffer) | ||
| 454 | { | ||
| 455 | namespace ph = std::placeholders; | ||
| 456 | |||
| 457 | 64 | mMaxAbsS = mMinAbsS; | |
| 458 |
1/2✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
|
64 | const size_t leafCount = mParent->mTracker.leafs().leafCount(); |
| 459 |
2/4✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
|
64 | if (leafCount==0 || time0 >= time1) return ValueType(0); |
| 460 | |||
| 461 | const math::Transform& xform = mParent->mTracker.grid().transform(); | ||
| 462 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
64 | if (mParent->mTarget->transform() == xform && |
| 463 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
64 | (mParent->mMask == nullptr || mParent->mMask->transform() == xform)) { |
| 464 | 64 | mTask = std::bind(&Morph::sampleAlignedSpeed, ph::_1, ph::_2, speedBuffer); | |
| 465 | } else { | ||
| 466 | ✗ | mTask = std::bind(&Morph::sampleXformedSpeed, ph::_1, ph::_2, speedBuffer); | |
| 467 | } | ||
| 468 | 64 | this->cook(PARALLEL_REDUCE); | |
| 469 |
1/2✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
|
64 | if (math::isApproxEqual(mMinAbsS, mMaxAbsS)) return ValueType(0);//speed is essentially zero |
| 470 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 31 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
64 | static const ValueType CFL = (TemporalScheme == math::TVD_RK1 ? ValueType(0.3) : |
| 471 | TemporalScheme == math::TVD_RK2 ? ValueType(0.9) : | ||
| 472 | ValueType(1.0))/math::Sqrt(ValueType(3.0)); | ||
| 473 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
64 | const ValueType dt = math::Abs(time1 - time0), dx = mParent->mTracker.voxelSize(); |
| 474 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
64 | return math::Min(dt, ValueType(CFL*dx/mMaxAbsS)); |
| 475 | } | ||
| 476 | |||
| 477 | template<typename GridT, typename InterruptT> | ||
| 478 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 479 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 480 | inline void | ||
| 481 | ✗ | LevelSetMorphing<GridT, InterruptT>:: | |
| 482 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 483 | sampleXformedSpeed(const LeafRange& range, Index speedBuffer) | ||
| 484 | { | ||
| 485 | using VoxelIterT = typename LeafType::ValueOnCIter; | ||
| 486 | using SamplerT = tools::GridSampler<typename GridT::ConstAccessor, tools::BoxSampler>; | ||
| 487 | |||
| 488 | ✗ | const MapT& map = *mMap; | |
| 489 | ✗ | mParent->mTracker.checkInterrupter(); | |
| 490 | |||
| 491 | ✗ | typename GridT::ConstAccessor targetAcc = mParent->mTarget->getAccessor(); | |
| 492 | ✗ | SamplerT target(targetAcc, mParent->mTarget->transform()); | |
| 493 | ✗ | if (mParent->mMask == nullptr) { | |
| 494 | ✗ | for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) { | |
| 495 | ✗ | ValueType* speed = leafIter.buffer(speedBuffer).data(); | |
| 496 | bool isZero = true; | ||
| 497 | ✗ | for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { | |
| 498 | ✗ | ValueType& s = speed[voxelIter.pos()]; | |
| 499 | ✗ | s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d())); | |
| 500 | if (!math::isApproxZero(s)) isZero = false; | ||
| 501 | ✗ | mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s)); | |
| 502 | } | ||
| 503 | ✗ | if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel | |
| 504 | } | ||
| 505 | } else { | ||
| 506 | ✗ | const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask); | |
| 507 | const bool invMask = mParent->isMaskInverted(); | ||
| 508 | typename GridT::ConstAccessor maskAcc = mParent->mMask->getAccessor(); | ||
| 509 | ✗ | SamplerT mask(maskAcc, mParent->mMask->transform()); | |
| 510 | ✗ | for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) { | |
| 511 | ✗ | ValueType* speed = leafIter.buffer(speedBuffer).data(); | |
| 512 | bool isZero = true; | ||
| 513 | ✗ | for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { | |
| 514 | ✗ | const Vec3R xyz = map.applyMap(voxelIter.getCoord().asVec3d());//world space | |
| 515 | ✗ | const ValueType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm); | |
| 516 | ✗ | ValueType& s = speed[voxelIter.pos()]; | |
| 517 | ✗ | s -= target.wsSample(xyz); | |
| 518 | ✗ | s *= invMask ? 1 - a : a; | |
| 519 | if (!math::isApproxZero(s)) isZero = false; | ||
| 520 | ✗ | mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s)); | |
| 521 | } | ||
| 522 | ✗ | if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel | |
| 523 | } | ||
| 524 | } | ||
| 525 | } | ||
| 526 | |||
| 527 | template<typename GridT, typename InterruptT> | ||
| 528 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 529 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 530 | inline void | ||
| 531 | 3244 | LevelSetMorphing<GridT, InterruptT>:: | |
| 532 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 533 | sampleAlignedSpeed(const LeafRange& range, Index speedBuffer) | ||
| 534 | { | ||
| 535 | using VoxelIterT = typename LeafType::ValueOnCIter; | ||
| 536 | |||
| 537 | 3244 | mParent->mTracker.checkInterrupter(); | |
| 538 | |||
| 539 | 3244 | typename GridT::ConstAccessor target = mParent->mTarget->getAccessor(); | |
| 540 | |||
| 541 |
1/2✓ Branch 0 taken 1622 times.
✗ Branch 1 not taken.
|
3244 | if (mParent->mMask == nullptr) { |
| 542 |
2/2✓ Branch 1 taken 1622 times.
✓ Branch 2 taken 1622 times.
|
6488 | for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) { |
| 543 |
1/2✓ Branch 2 taken 1622 times.
✗ Branch 3 not taken.
|
3244 | ValueType* speed = leafIter.buffer(speedBuffer).data(); |
| 544 | bool isZero = true; | ||
| 545 |
2/2✓ Branch 0 taken 221022 times.
✓ Branch 1 taken 1622 times.
|
445288 | for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { |
| 546 | 442044 | ValueType& s = speed[voxelIter.pos()]; | |
| 547 |
4/8✓ Branch 1 taken 221022 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 221022 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 98930 times.
✓ Branch 7 taken 122092 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
442044 | s -= target.getValue(voxelIter.getCoord()); |
| 548 | if (!math::isApproxZero(s)) isZero = false; | ||
| 549 |
2/2✓ Branch 0 taken 1512 times.
✓ Branch 1 taken 219510 times.
|
445068 | mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s)); |
| 550 | } | ||
| 551 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1622 times.
|
3244 | if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel |
| 552 | } | ||
| 553 | } else { | ||
| 554 | ✗ | const ValueType min = mParent->mMinMask, invNorm = 1.0f/(mParent->mDeltaMask); | |
| 555 | const bool invMask = mParent->isMaskInverted(); | ||
| 556 | typename GridT::ConstAccessor mask = mParent->mMask->getAccessor(); | ||
| 557 | ✗ | for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) { | |
| 558 | ✗ | ValueType* speed = leafIter.buffer(speedBuffer).data(); | |
| 559 | bool isZero = true; | ||
| 560 | ✗ | for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { | |
| 561 | ✗ | const Coord ijk = voxelIter.getCoord();//index space | |
| 562 | ✗ | const ValueType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm); | |
| 563 | ✗ | ValueType& s = speed[voxelIter.pos()]; | |
| 564 | ✗ | s -= target.getValue(ijk); | |
| 565 | ✗ | s *= invMask ? 1 - a : a; | |
| 566 | if (!math::isApproxZero(s)) isZero = false; | ||
| 567 | ✗ | mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s)); | |
| 568 | } | ||
| 569 | ✗ | if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel | |
| 570 | } | ||
| 571 | } | ||
| 572 | } | ||
| 573 | |||
| 574 | template<typename GridT, typename InterruptT> | ||
| 575 | template <typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 576 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 577 | inline void | ||
| 578 | 256 | LevelSetMorphing<GridT, InterruptT>:: | |
| 579 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 580 | cook(ThreadingMode mode, size_t swapBuffer) | ||
| 581 | { | ||
| 582 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 128 times.
|
256 | mParent->mTracker.startInterrupter("Morphing level set"); |
| 583 | |||
| 584 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 128 times.
|
256 | const int grainSize = mParent->mTracker.getGrainSize(); |
| 585 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 128 times.
|
256 | const LeafRange range = mParent->mTracker.leafs().leafRange(grainSize); |
| 586 | |||
| 587 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 128 times.
|
256 | if (mParent->mTracker.getGrainSize()==0) { |
| 588 | ✗ | (*this)(range); | |
| 589 |
2/2✓ Branch 0 taken 96 times.
✓ Branch 1 taken 32 times.
|
256 | } else if (mode == PARALLEL_FOR) { |
| 590 | 192 | tbb::parallel_for(range, *this); | |
| 591 |
1/2✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
|
64 | } else if (mode == PARALLEL_REDUCE) { |
| 592 | 64 | tbb::parallel_reduce(range, *this); | |
| 593 | } else { | ||
| 594 | ✗ | OPENVDB_THROW(ValueError, "expected threading mode " << int(PARALLEL_FOR) | |
| 595 | << " or " << int(PARALLEL_REDUCE) << ", got " << int(mode)); | ||
| 596 | } | ||
| 597 | |||
| 598 | 256 | mParent->mTracker.leafs().swapLeafBuffer(swapBuffer, grainSize == 0); | |
| 599 | |||
| 600 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 128 times.
|
256 | mParent->mTracker.endInterrupter(); |
| 601 | } | ||
| 602 | |||
| 603 | template<typename GridT, typename InterruptT> | ||
| 604 | template<typename MapT, math::BiasedGradientScheme SpatialScheme, | ||
| 605 | math::TemporalIntegrationScheme TemporalScheme> | ||
| 606 | template <int Nominator, int Denominator> | ||
| 607 | inline void | ||
| 608 | 9732 | LevelSetMorphing<GridT,InterruptT>:: | |
| 609 | Morph<MapT, SpatialScheme, TemporalScheme>:: | ||
| 610 | euler(const LeafRange& range, ValueType dt, | ||
| 611 | Index phiBuffer, Index resultBuffer, Index speedBuffer) | ||
| 612 | { | ||
| 613 | using SchemeT = math::BIAS_SCHEME<SpatialScheme>; | ||
| 614 | using StencilT = typename SchemeT::template ISStencil<GridType>::StencilType; | ||
| 615 | using VoxelIterT = typename LeafType::ValueOnCIter; | ||
| 616 | using NumGrad = math::GradientNormSqrd<MapT, SpatialScheme>; | ||
| 617 | |||
| 618 | static const ValueType Alpha = ValueType(Nominator)/ValueType(Denominator); | ||
| 619 | static const ValueType Beta = ValueType(1) - Alpha; | ||
| 620 | |||
| 621 | 9732 | mParent->mTracker.checkInterrupter(); | |
| 622 | 9732 | const MapT& map = *mMap; | |
| 623 | 9732 | StencilT stencil(mParent->mTracker.grid()); | |
| 624 | |||
| 625 |
2/2✓ Branch 1 taken 4866 times.
✓ Branch 2 taken 4866 times.
|
19464 | for (typename LeafRange::Iterator leafIter = range.begin(); leafIter; ++leafIter) { |
| 626 |
1/2✓ Branch 2 taken 4866 times.
✗ Branch 3 not taken.
|
9732 | const ValueType* speed = leafIter.buffer(speedBuffer).data(); |
| 627 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4866 times.
|
9732 | if (math::isExactlyEqual(speed[0], std::numeric_limits<ValueType>::max())) continue; |
| 628 |
1/2✓ Branch 2 taken 4866 times.
✗ Branch 3 not taken.
|
9732 | const ValueType* phi = leafIter.buffer(phiBuffer).data(); |
| 629 |
1/2✓ Branch 2 taken 4866 times.
✗ Branch 3 not taken.
|
9732 | ValueType* result = leafIter.buffer(resultBuffer).data(); |
| 630 |
2/2✓ Branch 0 taken 663066 times.
✓ Branch 1 taken 4866 times.
|
1335864 | for (VoxelIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { |
| 631 | const Index n = voxelIter.pos(); | ||
| 632 |
2/2✓ Branch 0 taken 296790 times.
✓ Branch 1 taken 366276 times.
|
1326144 | if (math::isApproxZero(speed[n])) continue; |
| 633 |
1/2✓ Branch 1 taken 663060 times.
✗ Branch 2 not taken.
|
1326120 | stencil.moveTo(voxelIter); |
| 634 | 1326120 | const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil); | |
| 635 | 1326120 | result[n] = Nominator ? Alpha * phi[n] + Beta * v : v; | |
| 636 | }//loop over active voxels in the leaf of the mask | ||
| 637 | }//loop over leafs of the level set | ||
| 638 | } | ||
| 639 | |||
| 640 | |||
| 641 | //////////////////////////////////////// | ||
| 642 | |||
| 643 | |||
| 644 | // Explicit Template Instantiation | ||
| 645 | |||
| 646 | #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION | ||
| 647 | |||
| 648 | #ifdef OPENVDB_INSTANTIATE_LEVELSETMORPH | ||
| 649 | #include <openvdb/util/ExplicitInstantiation.h> | ||
| 650 | #endif | ||
| 651 | |||
| 652 | OPENVDB_INSTANTIATE_CLASS LevelSetMorphing<FloatGrid, util::NullInterrupter>; | ||
| 653 | OPENVDB_INSTANTIATE_CLASS LevelSetMorphing<DoubleGrid, util::NullInterrupter>; | ||
| 654 | |||
| 655 | #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION | ||
| 656 | |||
| 657 | |||
| 658 | } // namespace tools | ||
| 659 | } // namespace OPENVDB_VERSION_NAME | ||
| 660 | } // namespace openvdb | ||
| 661 | |||
| 662 | #endif // OPENVDB_TOOLS_LEVEL_SET_MORPH_HAS_BEEN_INCLUDED | ||
| 663 |