1 // Copyright Contributors to the OpenVDB Project 2 // SPDX-License-Identifier: MPL-2.0 3 4 /// @file math/Operators.h 5 6 #ifndef OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED 7 #define OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED 8 9 #include "FiniteDifference.h" 10 #include "Stencils.h" 11 #include "Maps.h" 12 #include "Transform.h" 13 #include <cmath> // for std::sqrt() 14 15 16 namespace openvdb { 17 OPENVDB_USE_VERSION_NAMESPACE 18 namespace OPENVDB_VERSION_NAME { 19 namespace math { 20 21 // Simple tools to help determine when type conversions are needed 22 template<typename Vec3T> struct is_vec3d { static const bool value = false; }; 23 template<> struct is_vec3d<Vec3d> { static const bool value = true; }; 24 25 template<typename T> struct is_double { static const bool value = false; }; 26 template<> struct is_double<double> { static const bool value = true; }; 27 28 29 /// @brief Adapter to associate a map with a world-space operator, 30 /// giving it the same call signature as an index-space operator 31 /// @todo For now, the operator's result type must be specified explicitly, 32 /// but eventually it should be possible, via traits, to derive the result type 33 /// from the operator type. 34 template<typename MapType, typename OpType, typename ResultType> 35 struct MapAdapter { 36 MapAdapter(const MapType& m): map(m) {} 37 38 template<typename AccessorType> 39 inline ResultType 40 result(const AccessorType& grid, const Coord& ijk) { return OpType::result(map, grid, ijk); } 41 42 template<typename StencilType> 43 inline ResultType 44 result(const StencilType& stencil) { return OpType::result(map, stencil); } 45 46 const MapType map; 47 }; 48 49 50 /// Adapter for vector-valued index-space operators to return the vector magnitude 51 template<typename OpType> 52 struct ISOpMagnitude { 53 template<typename AccessorType> 54 static inline double result(const AccessorType& grid, const Coord& ijk) { 55 return double(OpType::result(grid, ijk).length()); 56 } 57 58 template<typename StencilType> 59 static inline double result(const StencilType& stencil) { 60 return double(OpType::result(stencil).length()); 61 } 62 }; 63 64 /// Adapter for vector-valued world-space operators to return the vector magnitude 65 template<typename OpType, typename MapT> 66 struct OpMagnitude { 67 template<typename AccessorType> 68 static inline double result(const MapT& map, const AccessorType& grid, const Coord& ijk) { 69 return double(OpType::result(map, grid, ijk).length()); 70 } 71 72 template<typename StencilType> 73 static inline double result(const MapT& map, const StencilType& stencil) { 74 return double(OpType::result(map, stencil).length()); 75 } 76 }; 77 78 /// @cond OPENVDB_DOCS_INTERNAL 79 80 namespace internal { 81 82 // This additional layer is necessary for Visual C++ to compile. 83 template<typename T> 84 struct ReturnValue { 85 using ValueType = typename T::ValueType; 86 using Vec3Type = math::Vec3<ValueType>; 87 }; 88 89 } // namespace internal 90 91 /// @endcond 92 93 // ---- Operators defined in index space 94 95 96 //@{ 97 /// @brief Gradient operators defined in index space of various orders 98 template<DScheme DiffScheme> 99 struct ISGradient 100 { 101 // random access version 102 template<typename Accessor> static Vec3<typename Accessor::ValueType> 103 result(const Accessor& grid, const Coord& ijk) 104 { 105 using ValueType = typename Accessor::ValueType; 106 using Vec3Type = Vec3<ValueType>; 107 return Vec3Type( D1<DiffScheme>::inX(grid, ijk), 108 D1<DiffScheme>::inY(grid, ijk), 109 D1<DiffScheme>::inZ(grid, ijk) ); 110 } 111 112 // stencil access version 113 template<typename StencilT> static Vec3<typename StencilT::ValueType> 114 result(const StencilT& stencil) 115 { 116 using ValueType = typename StencilT::ValueType; 117 using Vec3Type = Vec3<ValueType>; 118 return Vec3Type( D1<DiffScheme>::inX(stencil), 119 D1<DiffScheme>::inY(stencil), 120 D1<DiffScheme>::inZ(stencil) ); 121 } 122 }; 123 //@} 124 125 /// struct that relates the BiasedGradientScheme to the 126 /// forward and backward difference methods used, as well as to 127 /// the correct stencil type for index space use 128 template<BiasedGradientScheme bgs> 129 struct BIAS_SCHEME { 130 static const DScheme FD = FD_1ST; 131 static const DScheme BD = BD_1ST; 132 133 template<typename GridType, bool IsSafe = true> 134 struct ISStencil { 135 using StencilType = SevenPointStencil<GridType, IsSafe>; 136 }; 137 }; 138 139 template<> struct BIAS_SCHEME<FIRST_BIAS> 140 { 141 static const DScheme FD = FD_1ST; 142 static const DScheme BD = BD_1ST; 143 144 template<typename GridType, bool IsSafe = true> 145 struct ISStencil { 146 using StencilType = SevenPointStencil<GridType, IsSafe>; 147 }; 148 }; 149 150 template<> struct BIAS_SCHEME<SECOND_BIAS> 151 { 152 static const DScheme FD = FD_2ND; 153 static const DScheme BD = BD_2ND; 154 155 template<typename GridType, bool IsSafe = true> 156 struct ISStencil { 157 using StencilType = ThirteenPointStencil<GridType, IsSafe>; 158 }; 159 }; 160 template<> struct BIAS_SCHEME<THIRD_BIAS> 161 { 162 static const DScheme FD = FD_3RD; 163 static const DScheme BD = BD_3RD; 164 165 template<typename GridType, bool IsSafe = true> 166 struct ISStencil { 167 using StencilType = NineteenPointStencil<GridType, IsSafe>; 168 }; 169 }; 170 template<> struct BIAS_SCHEME<WENO5_BIAS> 171 { 172 static const DScheme FD = FD_WENO5; 173 static const DScheme BD = BD_WENO5; 174 175 template<typename GridType, bool IsSafe = true> 176 struct ISStencil { 177 using StencilType = NineteenPointStencil<GridType, IsSafe>; 178 }; 179 }; 180 template<> struct BIAS_SCHEME<HJWENO5_BIAS> 181 { 182 static const DScheme FD = FD_HJWENO5; 183 static const DScheme BD = BD_HJWENO5; 184 185 template<typename GridType, bool IsSafe = true> 186 struct ISStencil { 187 using StencilType = NineteenPointStencil<GridType, IsSafe>; 188 }; 189 }; 190 191 192 //@{ 193 /// @brief Biased Gradient Operators, using upwinding defined by the @c Vec3Bias input 194 195 template<BiasedGradientScheme GradScheme, typename Vec3Bias> 196 struct ISGradientBiased 197 { 198 static const DScheme FD = BIAS_SCHEME<GradScheme>::FD; 199 static const DScheme BD = BIAS_SCHEME<GradScheme>::BD; 200 201 // random access version 202 template<typename Accessor> 203 static Vec3<typename Accessor::ValueType> 204 result(const Accessor& grid, const Coord& ijk, const Vec3Bias& V) 205 { 206 using ValueType = typename Accessor::ValueType; 207 using Vec3Type = Vec3<ValueType>; 208 209 return Vec3Type(V[0]<0 ? D1<FD>::inX(grid,ijk) : D1<BD>::inX(grid,ijk), 210 V[1]<0 ? D1<FD>::inY(grid,ijk) : D1<BD>::inY(grid,ijk), 211 V[2]<0 ? D1<FD>::inZ(grid,ijk) : D1<BD>::inZ(grid,ijk) ); 212 } 213 214 // stencil access version 215 template<typename StencilT> 216 static Vec3<typename StencilT::ValueType> 217 result(const StencilT& stencil, const Vec3Bias& V) 218 { 219 using ValueType = typename StencilT::ValueType; 220 using Vec3Type = Vec3<ValueType>; 221 222 return Vec3Type(V[0]<0 ? D1<FD>::inX(stencil) : D1<BD>::inX(stencil), 223 V[1]<0 ? D1<FD>::inY(stencil) : D1<BD>::inY(stencil), 224 V[2]<0 ? D1<FD>::inZ(stencil) : D1<BD>::inZ(stencil) ); 225 } 226 }; 227 228 229 template<BiasedGradientScheme GradScheme> 230 struct ISGradientNormSqrd 231 { 232 static const DScheme FD = BIAS_SCHEME<GradScheme>::FD; 233 static const DScheme BD = BIAS_SCHEME<GradScheme>::BD; 234 235 236 // random access version 237 template<typename Accessor> 238 static typename Accessor::ValueType 239 result(const Accessor& grid, const Coord& ijk) 240 { 241 using ValueType = typename Accessor::ValueType; 242 using Vec3Type = math::Vec3<ValueType>; 243 244 Vec3Type up = ISGradient<FD>::result(grid, ijk); 245 Vec3Type down = ISGradient<BD>::result(grid, ijk); 246 return math::GodunovsNormSqrd(grid.getValue(ijk)>0, down, up); 247 } 248 249 // stencil access version 250 template<typename StencilT> 251 static typename StencilT::ValueType 252 result(const StencilT& stencil) 253 { 254 using ValueType = typename StencilT::ValueType; 255 using Vec3Type = math::Vec3<ValueType>; 256 257 Vec3Type up = ISGradient<FD>::result(stencil); 258 Vec3Type down = ISGradient<BD>::result(stencil); 259 return math::GodunovsNormSqrd(stencil.template getValue<0, 0, 0>()>0, down, up); 260 } 261 }; 262 263 #ifdef DWA_OPENVDB // for SIMD - note will do the computations in float 264 template<> 265 struct ISGradientNormSqrd<HJWENO5_BIAS> 266 { 267 // random access version 268 template<typename Accessor> 269 static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk) 270 { 271 struct GetValue 272 { 273 const Accessor& acc; 274 GetValue(const Accessor& acc_): acc(acc_) {} 275 // Return the grid value at ijk converted to simd::Float4::value_type (= float). 276 inline simd::Float4::value_type operator()(const Coord& ijk_) { 277 return static_cast<simd::Float4::value_type>(acc.getValue(ijk_)); 278 } 279 } 280 valueAt(grid); 281 282 // SSE optimized 283 const simd::Float4 284 v1(valueAt(ijk.offsetBy(-2, 0, 0)) - valueAt(ijk.offsetBy(-3, 0, 0)), 285 valueAt(ijk.offsetBy( 0,-2, 0)) - valueAt(ijk.offsetBy( 0,-3, 0)), 286 valueAt(ijk.offsetBy( 0, 0,-2)) - valueAt(ijk.offsetBy( 0, 0,-3)), 0), 287 v2(valueAt(ijk.offsetBy(-1, 0, 0)) - valueAt(ijk.offsetBy(-2, 0, 0)), 288 valueAt(ijk.offsetBy( 0,-1, 0)) - valueAt(ijk.offsetBy( 0,-2, 0)), 289 valueAt(ijk.offsetBy( 0, 0,-1)) - valueAt(ijk.offsetBy( 0, 0,-2)), 0), 290 v3(valueAt(ijk ) - valueAt(ijk.offsetBy(-1, 0, 0)), 291 valueAt(ijk ) - valueAt(ijk.offsetBy( 0,-1, 0)), 292 valueAt(ijk ) - valueAt(ijk.offsetBy( 0, 0,-1)), 0), 293 v4(valueAt(ijk.offsetBy( 1, 0, 0)) - valueAt(ijk ), 294 valueAt(ijk.offsetBy( 0, 1, 0)) - valueAt(ijk ), 295 valueAt(ijk.offsetBy( 0, 0, 1)) - valueAt(ijk ), 0), 296 v5(valueAt(ijk.offsetBy( 2, 0, 0)) - valueAt(ijk.offsetBy( 1, 0, 0)), 297 valueAt(ijk.offsetBy( 0, 2, 0)) - valueAt(ijk.offsetBy( 0, 1, 0)), 298 valueAt(ijk.offsetBy( 0, 0, 2)) - valueAt(ijk.offsetBy( 0, 0, 1)), 0), 299 v6(valueAt(ijk.offsetBy( 3, 0, 0)) - valueAt(ijk.offsetBy( 2, 0, 0)), 300 valueAt(ijk.offsetBy( 0, 3, 0)) - valueAt(ijk.offsetBy( 0, 2, 0)), 301 valueAt(ijk.offsetBy( 0, 0, 3)) - valueAt(ijk.offsetBy( 0, 0, 2)), 0), 302 down = math::WENO5(v1, v2, v3, v4, v5), 303 up = math::WENO5(v6, v5, v4, v3, v2); 304 305 return math::GodunovsNormSqrd(grid.getValue(ijk)>0, down, up); 306 } 307 308 // stencil access version 309 template<typename StencilT> 310 static typename StencilT::ValueType result(const StencilT& s) 311 { 312 using F4Val = simd::Float4::value_type; 313 314 // SSE optimized 315 const simd::Float4 316 v1(F4Val(s.template getValue<-2, 0, 0>()) - F4Val(s.template getValue<-3, 0, 0>()), 317 F4Val(s.template getValue< 0,-2, 0>()) - F4Val(s.template getValue< 0,-3, 0>()), 318 F4Val(s.template getValue< 0, 0,-2>()) - F4Val(s.template getValue< 0, 0,-3>()), 0), 319 v2(F4Val(s.template getValue<-1, 0, 0>()) - F4Val(s.template getValue<-2, 0, 0>()), 320 F4Val(s.template getValue< 0,-1, 0>()) - F4Val(s.template getValue< 0,-2, 0>()), 321 F4Val(s.template getValue< 0, 0,-1>()) - F4Val(s.template getValue< 0, 0,-2>()), 0), 322 v3(F4Val(s.template getValue< 0, 0, 0>()) - F4Val(s.template getValue<-1, 0, 0>()), 323 F4Val(s.template getValue< 0, 0, 0>()) - F4Val(s.template getValue< 0,-1, 0>()), 324 F4Val(s.template getValue< 0, 0, 0>()) - F4Val(s.template getValue< 0, 0,-1>()), 0), 325 v4(F4Val(s.template getValue< 1, 0, 0>()) - F4Val(s.template getValue< 0, 0, 0>()), 326 F4Val(s.template getValue< 0, 1, 0>()) - F4Val(s.template getValue< 0, 0, 0>()), 327 F4Val(s.template getValue< 0, 0, 1>()) - F4Val(s.template getValue< 0, 0, 0>()), 0), 328 v5(F4Val(s.template getValue< 2, 0, 0>()) - F4Val(s.template getValue< 1, 0, 0>()), 329 F4Val(s.template getValue< 0, 2, 0>()) - F4Val(s.template getValue< 0, 1, 0>()), 330 F4Val(s.template getValue< 0, 0, 2>()) - F4Val(s.template getValue< 0, 0, 1>()), 0), 331 v6(F4Val(s.template getValue< 3, 0, 0>()) - F4Val(s.template getValue< 2, 0, 0>()), 332 F4Val(s.template getValue< 0, 3, 0>()) - F4Val(s.template getValue< 0, 2, 0>()), 333 F4Val(s.template getValue< 0, 0, 3>()) - F4Val(s.template getValue< 0, 0, 2>()), 0), 334 down = math::WENO5(v1, v2, v3, v4, v5), 335 up = math::WENO5(v6, v5, v4, v3, v2); 336 337 return math::GodunovsNormSqrd(s.template getValue<0, 0, 0>()>0, down, up); 338 } 339 }; 340 #endif //DWA_OPENVDB // for SIMD - note will do the computations in float 341 //@} 342 343 344 //@{ 345 /// @brief Laplacian defined in index space, using various center-difference stencils 346 template<DDScheme DiffScheme> 347 struct ISLaplacian 348 { 349 // random access version 350 template<typename Accessor> 351 static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk); 352 353 // stencil access version 354 template<typename StencilT> 355 static typename StencilT::ValueType result(const StencilT& stencil); 356 }; 357 358 359 template<> 360 struct ISLaplacian<CD_SECOND> 361 { 362 // random access version 363 template<typename Accessor> 364 static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk) 365 { 366 return grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) + 367 grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy(0, -1, 0)) + 368 grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy(0, 0,-1)) 369 - 6*grid.getValue(ijk); 370 } 371 372 // stencil access version 373 template<typename StencilT> 374 static typename StencilT::ValueType result(const StencilT& stencil) 375 { 376 return stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() + 377 stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() + 378 stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() 379 - 6*stencil.template getValue< 0, 0, 0>(); 380 } 381 }; 382 383 template<> 384 struct ISLaplacian<CD_FOURTH> 385 { 386 // random access version 387 template<typename Accessor> 388 static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk) 389 { 390 using ValueT = typename Accessor::ValueType; 391 return static_cast<ValueT>( 392 (-1./12.)*( 393 grid.getValue(ijk.offsetBy(2,0,0)) + grid.getValue(ijk.offsetBy(-2, 0, 0)) + 394 grid.getValue(ijk.offsetBy(0,2,0)) + grid.getValue(ijk.offsetBy( 0,-2, 0)) + 395 grid.getValue(ijk.offsetBy(0,0,2)) + grid.getValue(ijk.offsetBy( 0, 0,-2)) ) 396 + (4./3.)*( 397 grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) + 398 grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy( 0,-1, 0)) + 399 grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy( 0, 0,-1)) ) 400 - 7.5*grid.getValue(ijk)); 401 } 402 403 // stencil access version 404 template<typename StencilT> 405 static typename StencilT::ValueType result(const StencilT& stencil) 406 { 407 using ValueT = typename StencilT::ValueType; 408 return static_cast<ValueT>( 409 (-1./12.)*( 410 stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() + 411 stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() + 412 stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() ) 413 + (4./3.)*( 414 stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() + 415 stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() + 416 stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() ) 417 - 7.5*stencil.template getValue< 0, 0, 0>()); 418 } 419 }; 420 421 template<> 422 struct ISLaplacian<CD_SIXTH> 423 { 424 // random access version 425 template<typename Accessor> 426 static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk) 427 { 428 using ValueT = typename Accessor::ValueType; 429 return static_cast<ValueT>( 430 (1./90.)*( 431 grid.getValue(ijk.offsetBy(3,0,0)) + grid.getValue(ijk.offsetBy(-3, 0, 0)) + 432 grid.getValue(ijk.offsetBy(0,3,0)) + grid.getValue(ijk.offsetBy( 0,-3, 0)) + 433 grid.getValue(ijk.offsetBy(0,0,3)) + grid.getValue(ijk.offsetBy( 0, 0,-3)) ) 434 - (3./20.)*( 435 grid.getValue(ijk.offsetBy(2,0,0)) + grid.getValue(ijk.offsetBy(-2, 0, 0)) + 436 grid.getValue(ijk.offsetBy(0,2,0)) + grid.getValue(ijk.offsetBy( 0,-2, 0)) + 437 grid.getValue(ijk.offsetBy(0,0,2)) + grid.getValue(ijk.offsetBy( 0, 0,-2)) ) 438 + 1.5 *( 439 grid.getValue(ijk.offsetBy(1,0,0)) + grid.getValue(ijk.offsetBy(-1, 0, 0)) + 440 grid.getValue(ijk.offsetBy(0,1,0)) + grid.getValue(ijk.offsetBy( 0,-1, 0)) + 441 grid.getValue(ijk.offsetBy(0,0,1)) + grid.getValue(ijk.offsetBy( 0, 0,-1)) ) 442 - (3*49/18.)*grid.getValue(ijk)); 443 } 444 445 // stencil access version 446 template<typename StencilT> 447 static typename StencilT::ValueType result(const StencilT& stencil) 448 { 449 using ValueT = typename StencilT::ValueType; 450 return static_cast<ValueT>( 451 (1./90.)*( 452 stencil.template getValue< 3, 0, 0>() + stencil.template getValue<-3, 0, 0>() + 453 stencil.template getValue< 0, 3, 0>() + stencil.template getValue< 0,-3, 0>() + 454 stencil.template getValue< 0, 0, 3>() + stencil.template getValue< 0, 0,-3>() ) 455 - (3./20.)*( 456 stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() + 457 stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() + 458 stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() ) 459 + 1.5 *( 460 stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() + 461 stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() + 462 stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() ) 463 - (3*49/18.)*stencil.template getValue< 0, 0, 0>()); 464 } 465 }; 466 //@} 467 468 469 //@{ 470 /// Divergence operator defined in index space using various first derivative schemes 471 template<DScheme DiffScheme> 472 struct ISDivergence 473 { 474 // random access version 475 template<typename Accessor> static typename Accessor::ValueType::value_type 476 result(const Accessor& grid, const Coord& ijk) 477 { 478 return D1Vec<DiffScheme>::inX(grid, ijk, 0) + 479 D1Vec<DiffScheme>::inY(grid, ijk, 1) + 480 D1Vec<DiffScheme>::inZ(grid, ijk, 2); 481 } 482 483 // stencil access version 484 template<typename StencilT> static typename StencilT::ValueType::value_type 485 result(const StencilT& stencil) 486 { 487 return D1Vec<DiffScheme>::inX(stencil, 0) + 488 D1Vec<DiffScheme>::inY(stencil, 1) + 489 D1Vec<DiffScheme>::inZ(stencil, 2); 490 } 491 }; 492 //@} 493 494 495 //@{ 496 /// Curl operator defined in index space using various first derivative schemes 497 template<DScheme DiffScheme> 498 struct ISCurl 499 { 500 // random access version 501 template<typename Accessor> 502 static typename Accessor::ValueType result(const Accessor& grid, const Coord& ijk) 503 { 504 using Vec3Type = typename Accessor::ValueType; 505 return Vec3Type( D1Vec<DiffScheme>::inY(grid, ijk, 2) - //dw/dy - dv/dz 506 D1Vec<DiffScheme>::inZ(grid, ijk, 1), 507 D1Vec<DiffScheme>::inZ(grid, ijk, 0) - //du/dz - dw/dx 508 D1Vec<DiffScheme>::inX(grid, ijk, 2), 509 D1Vec<DiffScheme>::inX(grid, ijk, 1) - //dv/dx - du/dy 510 D1Vec<DiffScheme>::inY(grid, ijk, 0) ); 511 } 512 513 // stencil access version 514 template<typename StencilT> 515 static typename StencilT::ValueType result(const StencilT& stencil) 516 { 517 using Vec3Type = typename StencilT::ValueType; 518 return Vec3Type( D1Vec<DiffScheme>::inY(stencil, 2) - //dw/dy - dv/dz 519 D1Vec<DiffScheme>::inZ(stencil, 1), 520 D1Vec<DiffScheme>::inZ(stencil, 0) - //du/dz - dw/dx 521 D1Vec<DiffScheme>::inX(stencil, 2), 522 D1Vec<DiffScheme>::inX(stencil, 1) - //dv/dx - du/dy 523 D1Vec<DiffScheme>::inY(stencil, 0) ); 524 } 525 }; 526 //@} 527 528 529 //@{ 530 /// Compute the mean curvature in index space 531 template<DDScheme DiffScheme2, DScheme DiffScheme1> 532 struct ISMeanCurvature 533 { 534 /// @brief Random access version 535 /// @return @c true if the gradient is nonzero, in which case the mean curvature 536 /// is returned in two parts, @a alpha and @a beta, where @a alpha is the numerator 537 /// in ∇ · (∇Φ / |∇Φ|) and @a beta is |∇Φ|. 538 template<typename Accessor> 539 static bool result(const Accessor& grid, const Coord& ijk, 540 typename Accessor::ValueType& alpha, 541 typename Accessor::ValueType& beta) 542 { 543 using ValueType = typename Accessor::ValueType; 544 545 const ValueType Dx = D1<DiffScheme1>::inX(grid, ijk); 546 const ValueType Dy = D1<DiffScheme1>::inY(grid, ijk); 547 const ValueType Dz = D1<DiffScheme1>::inZ(grid, ijk); 548 549 const ValueType Dx2 = Dx*Dx; 550 const ValueType Dy2 = Dy*Dy; 551 const ValueType Dz2 = Dz*Dz; 552 const ValueType normGrad = Dx2 + Dy2 + Dz2; 553 if (normGrad <= math::Tolerance<ValueType>::value()) { 554 alpha = beta = 0; 555 return false; 556 } 557 558 const ValueType Dxx = D2<DiffScheme2>::inX(grid, ijk); 559 const ValueType Dyy = D2<DiffScheme2>::inY(grid, ijk); 560 const ValueType Dzz = D2<DiffScheme2>::inZ(grid, ijk); 561 562 const ValueType Dxy = D2<DiffScheme2>::inXandY(grid, ijk); 563 const ValueType Dyz = D2<DiffScheme2>::inYandZ(grid, ijk); 564 const ValueType Dxz = D2<DiffScheme2>::inXandZ(grid, ijk); 565 566 // for return 567 alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz)); 568 beta = ValueType(std::sqrt(double(normGrad))); // * 1/dx 569 return true; 570 } 571 572 /// @brief Stencil access version 573 /// @return @c true if the gradient is nonzero, in which case the mean curvature 574 /// is returned in two parts, @a alpha and @a beta, where @a alpha is the numerator 575 /// in ∇ · (∇Φ / |∇Φ|) and @a beta is |∇Φ|. 576 template<typename StencilT> 577 static bool result(const StencilT& stencil, 578 typename StencilT::ValueType& alpha, 579 typename StencilT::ValueType& beta) 580 { 581 using ValueType = typename StencilT::ValueType; 582 const ValueType Dx = D1<DiffScheme1>::inX(stencil); 583 const ValueType Dy = D1<DiffScheme1>::inY(stencil); 584 const ValueType Dz = D1<DiffScheme1>::inZ(stencil); 585 586 const ValueType Dx2 = Dx*Dx; 587 const ValueType Dy2 = Dy*Dy; 588 const ValueType Dz2 = Dz*Dz; 589 const ValueType normGrad = Dx2 + Dy2 + Dz2; 590 if (normGrad <= math::Tolerance<ValueType>::value()) { 591 alpha = beta = 0; 592 return false; 593 } 594 595 const ValueType Dxx = D2<DiffScheme2>::inX(stencil); 596 const ValueType Dyy = D2<DiffScheme2>::inY(stencil); 597 const ValueType Dzz = D2<DiffScheme2>::inZ(stencil); 598 599 const ValueType Dxy = D2<DiffScheme2>::inXandY(stencil); 600 const ValueType Dyz = D2<DiffScheme2>::inYandZ(stencil); 601 const ValueType Dxz = D2<DiffScheme2>::inXandZ(stencil); 602 603 // for return 604 alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz)); 605 beta = ValueType(std::sqrt(double(normGrad))); // * 1/dx 606 return true; 607 } 608 }; 609 610 //////////////////////////////////////////////////////// 611 612 // --- Operators defined in the Range of a given map 613 614 //@{ 615 /// @brief Center difference gradient operators, defined with respect to 616 /// the range-space of the @c map 617 /// @note This will need to be divided by two in the case of CD_2NDT 618 template<typename MapType, DScheme DiffScheme> 619 struct Gradient 620 { 621 // random access version 622 template<typename Accessor> 623 static typename internal::ReturnValue<Accessor>::Vec3Type 624 result(const MapType& map, const Accessor& grid, const Coord& ijk) 625 { 626 using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type; 627 628 Vec3d iGradient( ISGradient<DiffScheme>::result(grid, ijk) ); 629 return Vec3Type(map.applyIJT(iGradient, ijk.asVec3d())); 630 } 631 632 // stencil access version 633 template<typename StencilT> 634 static typename internal::ReturnValue<StencilT>::Vec3Type 635 result(const MapType& map, const StencilT& stencil) 636 { 637 using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type; 638 639 Vec3d iGradient( ISGradient<DiffScheme>::result(stencil) ); 640 return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d())); 641 } 642 }; 643 644 // Partial template specialization of Gradient 645 // translation, any order 646 template<DScheme DiffScheme> 647 struct Gradient<TranslationMap, DiffScheme> 648 { 649 // random access version 650 template<typename Accessor> 651 static typename internal::ReturnValue<Accessor>::Vec3Type 652 result(const TranslationMap&, const Accessor& grid, const Coord& ijk) 653 { 654 return ISGradient<DiffScheme>::result(grid, ijk); 655 } 656 657 // stencil access version 658 template<typename StencilT> 659 static typename internal::ReturnValue<StencilT>::Vec3Type 660 result(const TranslationMap&, const StencilT& stencil) 661 { 662 return ISGradient<DiffScheme>::result(stencil); 663 } 664 }; 665 666 /// Full template specialization of Gradient 667 /// uniform scale, 2nd order 668 template<> 669 struct Gradient<UniformScaleMap, CD_2ND> 670 { 671 // random access version 672 template<typename Accessor> 673 static typename internal::ReturnValue<Accessor>::Vec3Type 674 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 675 { 676 using ValueType = typename internal::ReturnValue<Accessor>::ValueType; 677 using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type; 678 679 Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) ); 680 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 681 return iGradient * inv2dx; 682 } 683 684 // stencil access version 685 template<typename StencilT> 686 static typename internal::ReturnValue<StencilT>::Vec3Type 687 result(const UniformScaleMap& map, const StencilT& stencil) 688 { 689 using ValueType = typename internal::ReturnValue<StencilT>::ValueType; 690 using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type; 691 692 Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) ); 693 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 694 return iGradient * inv2dx; 695 } 696 }; 697 698 /// Full template specialization of Gradient 699 /// uniform scale translate, 2nd order 700 template<> 701 struct Gradient<UniformScaleTranslateMap, CD_2ND> 702 { 703 // random access version 704 template<typename Accessor> 705 static typename internal::ReturnValue<Accessor>::Vec3Type 706 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 707 { 708 using ValueType = typename internal::ReturnValue<Accessor>::ValueType; 709 using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type; 710 711 Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) ); 712 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 713 return iGradient * inv2dx; 714 } 715 716 // stencil access version 717 template<typename StencilT> 718 static typename internal::ReturnValue<StencilT>::Vec3Type 719 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 720 { 721 using ValueType = typename internal::ReturnValue<StencilT>::ValueType; 722 using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type; 723 724 Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) ); 725 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 726 return iGradient * inv2dx; 727 } 728 }; 729 730 /// Full template specialization of Gradient 731 /// scale, 2nd order 732 template<> 733 struct Gradient<ScaleMap, CD_2ND> 734 { 735 // random access version 736 template<typename Accessor> 737 static typename internal::ReturnValue<Accessor>::Vec3Type 738 result(const ScaleMap& map, const Accessor& grid, const Coord& ijk) 739 { 740 using ValueType = typename internal::ReturnValue<Accessor>::ValueType; 741 using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type; 742 743 Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) ); 744 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 745 const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0]; 746 const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1]; 747 const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2]; 748 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 749 return Vec3Type(ValueType(gradient0), 750 ValueType(gradient1), 751 ValueType(gradient2)); 752 } 753 754 // stencil access version 755 template<typename StencilT> 756 static typename internal::ReturnValue<StencilT>::Vec3Type 757 result(const ScaleMap& map, const StencilT& stencil) 758 { 759 using ValueType = typename internal::ReturnValue<StencilT>::ValueType; 760 using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type; 761 762 Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) ); 763 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 764 const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0]; 765 const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1]; 766 const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2]; 767 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 768 return Vec3Type(ValueType(gradient0), 769 ValueType(gradient1), 770 ValueType(gradient2)); 771 } 772 }; 773 774 /// Full template specialization of Gradient 775 /// scale translate, 2nd order 776 template<> 777 struct Gradient<ScaleTranslateMap, CD_2ND> 778 { 779 // random access version 780 template<typename Accessor> 781 static typename internal::ReturnValue<Accessor>::Vec3Type 782 result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 783 { 784 using ValueType = typename internal::ReturnValue<Accessor>::ValueType; 785 using Vec3Type = typename internal::ReturnValue<Accessor>::Vec3Type; 786 787 Vec3Type iGradient( ISGradient<CD_2NDT>::result(grid, ijk) ); 788 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 789 const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0]; 790 const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1]; 791 const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2]; 792 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 793 return Vec3Type(ValueType(gradient0), 794 ValueType(gradient1), 795 ValueType(gradient2)); 796 } 797 798 // Stencil access version 799 template<typename StencilT> 800 static typename internal::ReturnValue<StencilT>::Vec3Type 801 result(const ScaleTranslateMap& map, const StencilT& stencil) 802 { 803 using ValueType = typename internal::ReturnValue<StencilT>::ValueType; 804 using Vec3Type = typename internal::ReturnValue<StencilT>::Vec3Type; 805 806 Vec3Type iGradient( ISGradient<CD_2NDT>::result(stencil) ); 807 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 808 const auto gradient0 = iGradient[0] * map.getInvTwiceScale()[0]; 809 const auto gradient1 = iGradient[1] * map.getInvTwiceScale()[1]; 810 const auto gradient2 = iGradient[2] * map.getInvTwiceScale()[2]; 811 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 812 return Vec3Type(ValueType(gradient0), 813 ValueType(gradient1), 814 ValueType(gradient2)); 815 } 816 }; 817 //@} 818 819 820 //@{ 821 /// @brief Biased gradient operators, defined with respect to the range-space of the map 822 /// @note This will need to be divided by two in the case of CD_2NDT 823 template<typename MapType, BiasedGradientScheme GradScheme> 824 struct GradientBiased 825 { 826 // random access version 827 template<typename Accessor> static math::Vec3<typename Accessor::ValueType> 828 result(const MapType& map, const Accessor& grid, const Coord& ijk, 829 const Vec3<typename Accessor::ValueType>& V) 830 { 831 using ValueType = typename Accessor::ValueType; 832 using Vec3Type = math::Vec3<ValueType>; 833 834 Vec3d iGradient( ISGradientBiased<GradScheme, Vec3Type>::result(grid, ijk, V) ); 835 return Vec3Type(map.applyIJT(iGradient, ijk.asVec3d())); 836 } 837 838 // stencil access version 839 template<typename StencilT> static math::Vec3<typename StencilT::ValueType> 840 result(const MapType& map, const StencilT& stencil, 841 const Vec3<typename StencilT::ValueType>& V) 842 { 843 using ValueType = typename StencilT::ValueType; 844 using Vec3Type = math::Vec3<ValueType>; 845 846 Vec3d iGradient( ISGradientBiased<GradScheme, Vec3Type>::result(stencil, V) ); 847 return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d())); 848 } 849 }; 850 //@} 851 852 853 //////////////////////////////////////////////////////// 854 855 // Computes |Grad[Phi]| using upwinding 856 template<typename MapType, BiasedGradientScheme GradScheme> 857 struct GradientNormSqrd 858 { 859 static const DScheme FD = BIAS_SCHEME<GradScheme>::FD; 860 static const DScheme BD = BIAS_SCHEME<GradScheme>::BD; 861 862 863 // random access version 864 template<typename Accessor> 865 static typename Accessor::ValueType 866 result(const MapType& map, const Accessor& grid, const Coord& ijk) 867 { 868 using ValueType = typename Accessor::ValueType; 869 using Vec3Type = math::Vec3<ValueType>; 870 871 Vec3Type up = Gradient<MapType, FD>::result(map, grid, ijk); 872 Vec3Type down = Gradient<MapType, BD>::result(map, grid, ijk); 873 return math::GodunovsNormSqrd(grid.getValue(ijk)>0, down, up); 874 } 875 876 // stencil access version 877 template<typename StencilT> 878 static typename StencilT::ValueType 879 result(const MapType& map, const StencilT& stencil) 880 { 881 using ValueType = typename StencilT::ValueType; 882 using Vec3Type = math::Vec3<ValueType>; 883 884 Vec3Type up = Gradient<MapType, FD>::result(map, stencil); 885 Vec3Type down = Gradient<MapType, BD>::result(map, stencil); 886 return math::GodunovsNormSqrd(stencil.template getValue<0, 0, 0>()>0, down, up); 887 } 888 }; 889 890 /// Partial template specialization of GradientNormSqrd 891 template<BiasedGradientScheme GradScheme> 892 struct GradientNormSqrd<UniformScaleMap, GradScheme> 893 { 894 // random access version 895 template<typename Accessor> 896 static typename Accessor::ValueType 897 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 898 { 899 using ValueType = typename Accessor::ValueType; 900 901 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 902 return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk); 903 } 904 905 // stencil access version 906 template<typename StencilT> 907 static typename StencilT::ValueType 908 result(const UniformScaleMap& map, const StencilT& stencil) 909 { 910 using ValueType = typename StencilT::ValueType; 911 912 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 913 return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil); 914 } 915 }; 916 917 /// Partial template specialization of GradientNormSqrd 918 template<BiasedGradientScheme GradScheme> 919 struct GradientNormSqrd<UniformScaleTranslateMap, GradScheme> 920 { 921 // random access version 922 template<typename Accessor> 923 static typename Accessor::ValueType 924 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 925 { 926 using ValueType = typename Accessor::ValueType; 927 928 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 929 return invdxdx * ISGradientNormSqrd<GradScheme>::result(grid, ijk); 930 } 931 932 // stencil access version 933 template<typename StencilT> 934 static typename StencilT::ValueType 935 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 936 { 937 using ValueType = typename StencilT::ValueType; 938 939 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 940 return invdxdx * ISGradientNormSqrd<GradScheme>::result(stencil); 941 } 942 }; 943 944 945 //@{ 946 /// @brief Compute the divergence of a vector-valued grid using differencing 947 /// of various orders, the result defined with respect to the range-space of the map. 948 template<typename MapType, DScheme DiffScheme> 949 struct Divergence 950 { 951 // random access version 952 template<typename Accessor> static typename Accessor::ValueType::value_type 953 result(const MapType& map, const Accessor& grid, const Coord& ijk) 954 { 955 using ValueType = typename Accessor::ValueType::value_type; 956 957 ValueType div(0); 958 for (int i=0; i < 3; i++) { 959 Vec3d vec( D1Vec<DiffScheme>::inX(grid, ijk, i), 960 D1Vec<DiffScheme>::inY(grid, ijk, i), 961 D1Vec<DiffScheme>::inZ(grid, ijk, i) ); 962 div += ValueType(map.applyIJT(vec, ijk.asVec3d())[i]); 963 } 964 return div; 965 } 966 967 // stencil access version 968 template<typename StencilT> static typename StencilT::ValueType::value_type 969 result(const MapType& map, const StencilT& stencil) 970 { 971 using ValueType = typename StencilT::ValueType::value_type; 972 973 ValueType div(0); 974 for (int i=0; i < 3; i++) { 975 Vec3d vec( D1Vec<DiffScheme>::inX(stencil, i), 976 D1Vec<DiffScheme>::inY(stencil, i), 977 D1Vec<DiffScheme>::inZ(stencil, i) ); 978 div += ValueType(map.applyIJT(vec, stencil.getCenterCoord().asVec3d())[i]); 979 } 980 return div; 981 } 982 }; 983 984 /// Partial template specialization of Divergence 985 /// translation, any scheme 986 template<DScheme DiffScheme> 987 struct Divergence<TranslationMap, DiffScheme> 988 { 989 // random access version 990 template<typename Accessor> static typename Accessor::ValueType::value_type 991 result(const TranslationMap&, const Accessor& grid, const Coord& ijk) 992 { 993 using ValueType = typename Accessor::ValueType::value_type; 994 995 ValueType div(0); 996 div =ISDivergence<DiffScheme>::result(grid, ijk); 997 return div; 998 } 999 1000 // stencil access version 1001 template<typename StencilT> static typename StencilT::ValueType::value_type 1002 result(const TranslationMap&, const StencilT& stencil) 1003 { 1004 using ValueType = typename StencilT::ValueType::value_type; 1005 1006 ValueType div(0); 1007 div =ISDivergence<DiffScheme>::result(stencil); 1008 return div; 1009 } 1010 }; 1011 1012 /// Partial template specialization of Divergence 1013 /// uniform scale, any scheme 1014 template<DScheme DiffScheme> 1015 struct Divergence<UniformScaleMap, DiffScheme> 1016 { 1017 // random access version 1018 template<typename Accessor> static typename Accessor::ValueType::value_type 1019 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 1020 { 1021 using ValueType = typename Accessor::ValueType::value_type; 1022 1023 ValueType div(0); 1024 1025 div =ISDivergence<DiffScheme>::result(grid, ijk); 1026 ValueType invdx = ValueType(map.getInvScale()[0]); 1027 return div * invdx; 1028 } 1029 1030 // stencil access version 1031 template<typename StencilT> static typename StencilT::ValueType::value_type 1032 result(const UniformScaleMap& map, const StencilT& stencil) 1033 { 1034 using ValueType = typename StencilT::ValueType::value_type; 1035 1036 ValueType div(0); 1037 1038 div =ISDivergence<DiffScheme>::result(stencil); 1039 ValueType invdx = ValueType(map.getInvScale()[0]); 1040 return div * invdx; 1041 } 1042 }; 1043 1044 /// Partial template specialization of Divergence 1045 /// uniform scale and translation, any scheme 1046 template<DScheme DiffScheme> 1047 struct Divergence<UniformScaleTranslateMap, DiffScheme> 1048 { 1049 // random access version 1050 template<typename Accessor> static typename Accessor::ValueType::value_type 1051 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1052 { 1053 using ValueType = typename Accessor::ValueType::value_type; 1054 1055 ValueType div(0); 1056 1057 div =ISDivergence<DiffScheme>::result(grid, ijk); 1058 ValueType invdx = ValueType(map.getInvScale()[0]); 1059 return div * invdx; 1060 } 1061 1062 // stencil access version 1063 template<typename StencilT> static typename StencilT::ValueType::value_type 1064 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 1065 { 1066 using ValueType = typename StencilT::ValueType::value_type; 1067 1068 ValueType div(0); 1069 1070 div =ISDivergence<DiffScheme>::result(stencil); 1071 ValueType invdx = ValueType(map.getInvScale()[0]); 1072 return div * invdx; 1073 } 1074 }; 1075 1076 /// Full template specialization of Divergence 1077 /// uniform scale 2nd order 1078 template<> 1079 struct Divergence<UniformScaleMap, CD_2ND> 1080 { 1081 // random access version 1082 template<typename Accessor> static typename Accessor::ValueType::value_type 1083 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 1084 { 1085 using ValueType = typename Accessor::ValueType::value_type; 1086 1087 ValueType div(0); 1088 div =ISDivergence<CD_2NDT>::result(grid, ijk); 1089 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 1090 return div * inv2dx; 1091 } 1092 1093 // stencil access version 1094 template<typename StencilT> static typename StencilT::ValueType::value_type 1095 result(const UniformScaleMap& map, const StencilT& stencil) 1096 { 1097 using ValueType = typename StencilT::ValueType::value_type; 1098 1099 ValueType div(0); 1100 div =ISDivergence<CD_2NDT>::result(stencil); 1101 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 1102 return div * inv2dx; 1103 } 1104 }; 1105 1106 /// Full template specialization of Divergence 1107 /// uniform scale translate 2nd order 1108 template<> 1109 struct Divergence<UniformScaleTranslateMap, CD_2ND> 1110 { 1111 // random access version 1112 template<typename Accessor> static typename Accessor::ValueType::value_type 1113 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1114 { 1115 using ValueType = typename Accessor::ValueType::value_type; 1116 1117 ValueType div(0); 1118 1119 div =ISDivergence<CD_2NDT>::result(grid, ijk); 1120 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 1121 return div * inv2dx; 1122 } 1123 1124 // stencil access version 1125 template<typename StencilT> static typename StencilT::ValueType::value_type 1126 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 1127 { 1128 using ValueType = typename StencilT::ValueType::value_type; 1129 1130 ValueType div(0); 1131 1132 div =ISDivergence<CD_2NDT>::result(stencil); 1133 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 1134 return div * inv2dx; 1135 } 1136 }; 1137 1138 /// Partial template specialization of Divergence 1139 /// scale, any scheme 1140 template<DScheme DiffScheme> 1141 struct Divergence<ScaleMap, DiffScheme> 1142 { 1143 // random access version 1144 template<typename Accessor> static typename Accessor::ValueType::value_type 1145 result(const ScaleMap& map, const Accessor& grid, const Coord& ijk) 1146 { 1147 using ValueType = typename Accessor::ValueType::value_type; 1148 1149 ValueType div = ValueType( 1150 D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) + 1151 D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) + 1152 D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2])); 1153 return div; 1154 } 1155 1156 // stencil access version 1157 template<typename StencilT> static typename StencilT::ValueType::value_type 1158 result(const ScaleMap& map, const StencilT& stencil) 1159 { 1160 using ValueType = typename StencilT::ValueType::value_type; 1161 1162 ValueType div(0); 1163 div = ValueType( 1164 D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) + 1165 D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) + 1166 D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) ); 1167 return div; 1168 } 1169 }; 1170 1171 /// Partial template specialization of Divergence 1172 /// scale translate, any scheme 1173 template<DScheme DiffScheme> 1174 struct Divergence<ScaleTranslateMap, DiffScheme> 1175 { 1176 // random access version 1177 template<typename Accessor> static typename Accessor::ValueType::value_type 1178 result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1179 { 1180 using ValueType = typename Accessor::ValueType::value_type; 1181 1182 ValueType div = ValueType( 1183 D1Vec<DiffScheme>::inX(grid, ijk, 0) * (map.getInvScale()[0]) + 1184 D1Vec<DiffScheme>::inY(grid, ijk, 1) * (map.getInvScale()[1]) + 1185 D1Vec<DiffScheme>::inZ(grid, ijk, 2) * (map.getInvScale()[2])); 1186 return div; 1187 } 1188 1189 // stencil access version 1190 template<typename StencilT> static typename StencilT::ValueType::value_type 1191 result(const ScaleTranslateMap& map, const StencilT& stencil) 1192 { 1193 using ValueType = typename StencilT::ValueType::value_type; 1194 1195 ValueType div(0); 1196 div = ValueType( 1197 D1Vec<DiffScheme>::inX(stencil, 0) * (map.getInvScale()[0]) + 1198 D1Vec<DiffScheme>::inY(stencil, 1) * (map.getInvScale()[1]) + 1199 D1Vec<DiffScheme>::inZ(stencil, 2) * (map.getInvScale()[2]) ); 1200 return div; 1201 } 1202 }; 1203 1204 /// Full template specialization Divergence 1205 /// scale 2nd order 1206 template<> 1207 struct Divergence<ScaleMap, CD_2ND> 1208 { 1209 // random access version 1210 template<typename Accessor> static typename Accessor::ValueType::value_type 1211 result(const ScaleMap& map, const Accessor& grid, const Coord& ijk) 1212 { 1213 using ValueType = typename Accessor::ValueType::value_type; 1214 1215 ValueType div = ValueType( 1216 D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) + 1217 D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) + 1218 D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) ); 1219 return div; 1220 } 1221 1222 // stencil access version 1223 template<typename StencilT> static typename StencilT::ValueType::value_type 1224 result(const ScaleMap& map, const StencilT& stencil) 1225 { 1226 using ValueType = typename StencilT::ValueType::value_type; 1227 1228 ValueType div = ValueType( 1229 D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) + 1230 D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) + 1231 D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) ); 1232 return div; 1233 } 1234 }; 1235 1236 /// Full template specialization of Divergence 1237 /// scale and translate, 2nd order 1238 template<> 1239 struct Divergence<ScaleTranslateMap, CD_2ND> 1240 { 1241 // random access version 1242 template<typename Accessor> static typename Accessor::ValueType::value_type 1243 result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1244 { 1245 using ValueType = typename Accessor::ValueType::value_type; 1246 1247 ValueType div = ValueType( 1248 D1Vec<CD_2NDT>::inX(grid, ijk, 0) * (map.getInvTwiceScale()[0]) + 1249 D1Vec<CD_2NDT>::inY(grid, ijk, 1) * (map.getInvTwiceScale()[1]) + 1250 D1Vec<CD_2NDT>::inZ(grid, ijk, 2) * (map.getInvTwiceScale()[2]) ); 1251 return div; 1252 } 1253 1254 // stencil access version 1255 template<typename StencilT> static typename StencilT::ValueType::value_type 1256 result(const ScaleTranslateMap& map, const StencilT& stencil) 1257 { 1258 using ValueType = typename StencilT::ValueType::value_type; 1259 1260 ValueType div = ValueType( 1261 D1Vec<CD_2NDT>::inX(stencil, 0) * (map.getInvTwiceScale()[0]) + 1262 D1Vec<CD_2NDT>::inY(stencil, 1) * (map.getInvTwiceScale()[1]) + 1263 D1Vec<CD_2NDT>::inZ(stencil, 2) * (map.getInvTwiceScale()[2]) ); 1264 return div; 1265 } 1266 }; 1267 //@} 1268 1269 1270 //@{ 1271 /// @brief Compute the curl of a vector-valued grid using differencing 1272 /// of various orders in the space defined by the range of the map. 1273 template<typename MapType, DScheme DiffScheme> 1274 struct Curl 1275 { 1276 // random access version 1277 template<typename Accessor> static typename Accessor::ValueType 1278 result(const MapType& map, const Accessor& grid, const Coord& ijk) 1279 { 1280 using Vec3Type = typename Accessor::ValueType; 1281 Vec3Type mat[3]; 1282 for (int i = 0; i < 3; i++) { 1283 Vec3d vec( 1284 D1Vec<DiffScheme>::inX(grid, ijk, i), 1285 D1Vec<DiffScheme>::inY(grid, ijk, i), 1286 D1Vec<DiffScheme>::inZ(grid, ijk, i)); 1287 // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z) 1288 mat[i] = Vec3Type(map.applyIJT(vec, ijk.asVec3d())); 1289 } 1290 return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3 1291 mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1 1292 mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2 1293 } 1294 1295 // stencil access version 1296 template<typename StencilT> static typename StencilT::ValueType 1297 result(const MapType& map, const StencilT& stencil) 1298 { 1299 using Vec3Type = typename StencilT::ValueType; 1300 Vec3Type mat[3]; 1301 for (int i = 0; i < 3; i++) { 1302 Vec3d vec( 1303 D1Vec<DiffScheme>::inX(stencil, i), 1304 D1Vec<DiffScheme>::inY(stencil, i), 1305 D1Vec<DiffScheme>::inZ(stencil, i)); 1306 // dF_i/dx_j (x_1 = x, x_2 = y, x_3 = z) 1307 mat[i] = Vec3Type(map.applyIJT(vec, stencil.getCenterCoord().asVec3d())); 1308 } 1309 return Vec3Type(mat[2][1] - mat[1][2], // dF_3/dx_2 - dF_2/dx_3 1310 mat[0][2] - mat[2][0], // dF_1/dx_3 - dF_3/dx_1 1311 mat[1][0] - mat[0][1]); // dF_2/dx_1 - dF_1/dx_2 1312 } 1313 }; 1314 1315 /// Partial template specialization of Curl 1316 template<DScheme DiffScheme> 1317 struct Curl<UniformScaleMap, DiffScheme> 1318 { 1319 // random access version 1320 template<typename Accessor> static typename Accessor::ValueType 1321 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 1322 { 1323 using Vec3Type = typename Accessor::ValueType; 1324 using ValueType = typename Vec3Type::value_type; 1325 return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]); 1326 } 1327 1328 // Stencil access version 1329 template<typename StencilT> static typename StencilT::ValueType 1330 result(const UniformScaleMap& map, const StencilT& stencil) 1331 { 1332 using Vec3Type = typename StencilT::ValueType; 1333 using ValueType = typename Vec3Type::value_type; 1334 return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]); 1335 } 1336 }; 1337 1338 /// Partial template specialization of Curl 1339 template<DScheme DiffScheme> 1340 struct Curl<UniformScaleTranslateMap, DiffScheme> 1341 { 1342 // random access version 1343 template<typename Accessor> static typename Accessor::ValueType 1344 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1345 { 1346 using Vec3Type = typename Accessor::ValueType; 1347 using ValueType = typename Vec3Type::value_type; 1348 1349 return ISCurl<DiffScheme>::result(grid, ijk) * ValueType(map.getInvScale()[0]); 1350 } 1351 1352 // stencil access version 1353 template<typename StencilT> static typename StencilT::ValueType 1354 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 1355 { 1356 using Vec3Type = typename StencilT::ValueType; 1357 using ValueType = typename Vec3Type::value_type; 1358 1359 return ISCurl<DiffScheme>::result(stencil) * ValueType(map.getInvScale()[0]); 1360 } 1361 }; 1362 1363 /// Full template specialization of Curl 1364 template<> 1365 struct Curl<UniformScaleMap, CD_2ND> 1366 { 1367 // random access version 1368 template<typename Accessor> static typename Accessor::ValueType 1369 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 1370 { 1371 using Vec3Type = typename Accessor::ValueType; 1372 using ValueType = typename Vec3Type::value_type; 1373 1374 return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]); 1375 } 1376 1377 // stencil access version 1378 template<typename StencilT> static typename StencilT::ValueType 1379 result(const UniformScaleMap& map, const StencilT& stencil) 1380 { 1381 using Vec3Type = typename StencilT::ValueType; 1382 using ValueType = typename Vec3Type::value_type; 1383 1384 return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]); 1385 } 1386 }; 1387 1388 /// Full template specialization of Curl 1389 template<> 1390 struct Curl<UniformScaleTranslateMap, CD_2ND> 1391 { 1392 // random access version 1393 template<typename Accessor> static typename Accessor::ValueType 1394 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1395 { 1396 using Vec3Type = typename Accessor::ValueType; 1397 using ValueType = typename Vec3Type::value_type; 1398 1399 return ISCurl<CD_2NDT>::result(grid, ijk) * ValueType(map.getInvTwiceScale()[0]); 1400 } 1401 1402 // stencil access version 1403 template<typename StencilT> static typename StencilT::ValueType 1404 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 1405 { 1406 using Vec3Type = typename StencilT::ValueType; 1407 using ValueType = typename Vec3Type::value_type; 1408 1409 return ISCurl<CD_2NDT>::result(stencil) * ValueType(map.getInvTwiceScale()[0]); 1410 } 1411 }; 1412 //@} 1413 1414 1415 //@{ 1416 /// @brief Compute the Laplacian at a given location in a grid using finite differencing 1417 /// of various orders. The result is defined in the range of the map. 1418 template<typename MapType, DDScheme DiffScheme> 1419 struct Laplacian 1420 { 1421 // random access version 1422 template<typename Accessor> 1423 static typename Accessor::ValueType result(const MapType& map, 1424 const Accessor& grid, const Coord& ijk) 1425 { 1426 using ValueType = typename Accessor::ValueType; 1427 // all the second derivatives in index space 1428 ValueType iddx = D2<DiffScheme>::inX(grid, ijk); 1429 ValueType iddy = D2<DiffScheme>::inY(grid, ijk); 1430 ValueType iddz = D2<DiffScheme>::inZ(grid, ijk); 1431 1432 ValueType iddxy = D2<DiffScheme>::inXandY(grid, ijk); 1433 ValueType iddyz = D2<DiffScheme>::inYandZ(grid, ijk); 1434 ValueType iddxz = D2<DiffScheme>::inXandZ(grid, ijk); 1435 1436 // second derivatives in index space 1437 Mat3d d2_is(iddx, iddxy, iddxz, 1438 iddxy, iddy, iddyz, 1439 iddxz, iddyz, iddz); 1440 1441 Mat3d d2_rs; // to hold the second derivative matrix in range space 1442 if (is_linear<MapType>::value) { 1443 d2_rs = map.applyIJC(d2_is); 1444 } else { 1445 // compute the first derivatives with 2nd order accuracy. 1446 Vec3d d1_is(static_cast<double>(D1<CD_2ND>::inX(grid, ijk)), 1447 static_cast<double>(D1<CD_2ND>::inY(grid, ijk)), 1448 static_cast<double>(D1<CD_2ND>::inZ(grid, ijk))); 1449 1450 d2_rs = map.applyIJC(d2_is, d1_is, ijk.asVec3d()); 1451 } 1452 1453 // the trace of the second derivative (range space) matrix is laplacian 1454 return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2)); 1455 } 1456 1457 // stencil access version 1458 template<typename StencilT> 1459 static typename StencilT::ValueType result(const MapType& map, const StencilT& stencil) 1460 { 1461 using ValueType = typename StencilT::ValueType; 1462 // all the second derivatives in index space 1463 ValueType iddx = D2<DiffScheme>::inX(stencil); 1464 ValueType iddy = D2<DiffScheme>::inY(stencil); 1465 ValueType iddz = D2<DiffScheme>::inZ(stencil); 1466 1467 ValueType iddxy = D2<DiffScheme>::inXandY(stencil); 1468 ValueType iddyz = D2<DiffScheme>::inYandZ(stencil); 1469 ValueType iddxz = D2<DiffScheme>::inXandZ(stencil); 1470 1471 // second derivatives in index space 1472 Mat3d d2_is(iddx, iddxy, iddxz, 1473 iddxy, iddy, iddyz, 1474 iddxz, iddyz, iddz); 1475 1476 Mat3d d2_rs; // to hold the second derivative matrix in range space 1477 if (is_linear<MapType>::value) { 1478 d2_rs = map.applyIJC(d2_is); 1479 } else { 1480 // compute the first derivatives with 2nd order accuracy. 1481 Vec3d d1_is(D1<CD_2ND>::inX(stencil), 1482 D1<CD_2ND>::inY(stencil), 1483 D1<CD_2ND>::inZ(stencil) ); 1484 1485 d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d()); 1486 } 1487 1488 // the trace of the second derivative (range space) matrix is laplacian 1489 return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2)); 1490 } 1491 }; 1492 1493 1494 template<DDScheme DiffScheme> 1495 struct Laplacian<TranslationMap, DiffScheme> 1496 { 1497 // random access version 1498 template<typename Accessor> 1499 static typename Accessor::ValueType result(const TranslationMap&, 1500 const Accessor& grid, const Coord& ijk) 1501 { 1502 return ISLaplacian<DiffScheme>::result(grid, ijk); 1503 } 1504 1505 // stencil access version 1506 template<typename StencilT> 1507 static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil) 1508 { 1509 return ISLaplacian<DiffScheme>::result(stencil); 1510 } 1511 }; 1512 1513 1514 // The Laplacian is invariant to rotation or reflection. 1515 template<DDScheme DiffScheme> 1516 struct Laplacian<UnitaryMap, DiffScheme> 1517 { 1518 // random access version 1519 template<typename Accessor> 1520 static typename Accessor::ValueType result(const UnitaryMap&, 1521 const Accessor& grid, const Coord& ijk) 1522 { 1523 return ISLaplacian<DiffScheme>::result(grid, ijk); 1524 } 1525 1526 // stencil access version 1527 template<typename StencilT> 1528 static typename StencilT::ValueType result(const UnitaryMap&, const StencilT& stencil) 1529 { 1530 return ISLaplacian<DiffScheme>::result(stencil); 1531 } 1532 }; 1533 1534 1535 template<DDScheme DiffScheme> 1536 struct Laplacian<UniformScaleMap, DiffScheme> 1537 { 1538 // random access version 1539 template<typename Accessor> static typename Accessor::ValueType 1540 result(const UniformScaleMap& map, const Accessor& grid, const Coord& ijk) 1541 { 1542 using ValueType = typename Accessor::ValueType; 1543 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 1544 return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx; 1545 } 1546 1547 // stencil access version 1548 template<typename StencilT> static typename StencilT::ValueType 1549 result(const UniformScaleMap& map, const StencilT& stencil) 1550 { 1551 using ValueType = typename StencilT::ValueType; 1552 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 1553 return ISLaplacian<DiffScheme>::result(stencil) * invdxdx; 1554 } 1555 }; 1556 1557 1558 template<DDScheme DiffScheme> 1559 struct Laplacian<UniformScaleTranslateMap, DiffScheme> 1560 { 1561 // random access version 1562 template<typename Accessor> static typename Accessor::ValueType 1563 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1564 { 1565 using ValueType = typename Accessor::ValueType; 1566 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 1567 return ISLaplacian<DiffScheme>::result(grid, ijk) * invdxdx; 1568 } 1569 1570 // stencil access version 1571 template<typename StencilT> static typename StencilT::ValueType 1572 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 1573 { 1574 using ValueType = typename StencilT::ValueType; 1575 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 1576 return ISLaplacian<DiffScheme>::result(stencil) * invdxdx; 1577 } 1578 }; 1579 1580 1581 template<DDScheme DiffScheme> 1582 struct Laplacian<ScaleMap, DiffScheme> 1583 { 1584 // random access version 1585 template<typename Accessor> static typename Accessor::ValueType 1586 result(const ScaleMap& map, const Accessor& grid, const Coord& ijk) 1587 { 1588 using ValueType = typename Accessor::ValueType; 1589 1590 // compute the second derivatives in index space 1591 ValueType iddx = D2<DiffScheme>::inX(grid, ijk); 1592 ValueType iddy = D2<DiffScheme>::inY(grid, ijk); 1593 ValueType iddz = D2<DiffScheme>::inZ(grid, ijk); 1594 const Vec3d& invScaleSqr = map.getInvScaleSqr(); 1595 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 1596 // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum 1597 const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]; 1598 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 1599 return value; 1600 } 1601 1602 // stencil access version 1603 template<typename StencilT> static typename StencilT::ValueType 1604 result(const ScaleMap& map, const StencilT& stencil) 1605 { 1606 using ValueType = typename StencilT::ValueType; 1607 1608 // compute the second derivatives in index space 1609 ValueType iddx = D2<DiffScheme>::inX(stencil); 1610 ValueType iddy = D2<DiffScheme>::inY(stencil); 1611 ValueType iddz = D2<DiffScheme>::inZ(stencil); 1612 const Vec3d& invScaleSqr = map.getInvScaleSqr(); 1613 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 1614 // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum 1615 const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]; 1616 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 1617 return value; 1618 } 1619 }; 1620 1621 1622 template<DDScheme DiffScheme> 1623 struct Laplacian<ScaleTranslateMap, DiffScheme> 1624 { 1625 // random access version 1626 template<typename Accessor> static typename Accessor::ValueType 1627 result(const ScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 1628 { 1629 using ValueType = typename Accessor::ValueType; 1630 // compute the second derivatives in index space 1631 ValueType iddx = D2<DiffScheme>::inX(grid, ijk); 1632 ValueType iddy = D2<DiffScheme>::inY(grid, ijk); 1633 ValueType iddz = D2<DiffScheme>::inZ(grid, ijk); 1634 const Vec3d& invScaleSqr = map.getInvScaleSqr(); 1635 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 1636 // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum 1637 const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]; 1638 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 1639 return value; 1640 } 1641 1642 // stencil access version 1643 template<typename StencilT> static typename StencilT::ValueType 1644 result(const ScaleTranslateMap& map, const StencilT& stencil) 1645 { 1646 using ValueType = typename StencilT::ValueType; 1647 // compute the second derivatives in index space 1648 ValueType iddx = D2<DiffScheme>::inX(stencil); 1649 ValueType iddy = D2<DiffScheme>::inY(stencil); 1650 ValueType iddz = D2<DiffScheme>::inZ(stencil); 1651 const Vec3d& invScaleSqr = map.getInvScaleSqr(); 1652 OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN 1653 // scale them by the appropriate 1/dx^2, 1/dy^2, 1/dz^2 and sum 1654 const ValueType value = iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]; 1655 OPENVDB_NO_TYPE_CONVERSION_WARNING_END 1656 return value; 1657 } 1658 }; 1659 1660 1661 /// @brief Compute the closest-point transform to a level set. 1662 /// @return the closest point to the surface from which the level set was derived, 1663 /// in the domain space of the map (e.g., voxel space). 1664 template<typename MapType, DScheme DiffScheme> 1665 struct CPT 1666 { 1667 // random access version 1668 template<typename Accessor> static math::Vec3<typename Accessor::ValueType> 1669 result(const MapType& map, const Accessor& grid, const Coord& ijk) 1670 { 1671 using ValueType = typename Accessor::ValueType; 1672 using Vec3Type = Vec3<ValueType>; 1673 1674 // current distance 1675 ValueType d = grid.getValue(ijk); 1676 // compute gradient in physical space where it is a unit normal 1677 // since the grid holds a distance level set. 1678 Vec3d vectorFromSurface(d*Gradient<MapType,DiffScheme>::result(map, grid, ijk)); 1679 if (is_linear<MapType>::value) { 1680 Vec3d result = ijk.asVec3d() - map.applyInverseMap(vectorFromSurface); 1681 return Vec3Type(result); 1682 } else { 1683 Vec3d location = map.applyMap(ijk.asVec3d()); 1684 Vec3d result = map.applyInverseMap(location - vectorFromSurface); 1685 return Vec3Type(result); 1686 } 1687 } 1688 1689 // stencil access version 1690 template<typename StencilT> static math::Vec3<typename StencilT::ValueType> 1691 result(const MapType& map, const StencilT& stencil) 1692 { 1693 using ValueType = typename StencilT::ValueType; 1694 using Vec3Type = Vec3<ValueType>; 1695 1696 // current distance 1697 ValueType d = stencil.template getValue<0, 0, 0>(); 1698 // compute gradient in physical space where it is a unit normal 1699 // since the grid holds a distance level set. 1700 Vec3d vectorFromSurface(d*Gradient<MapType, DiffScheme>::result(map, stencil)); 1701 if (is_linear<MapType>::value) { 1702 Vec3d result = stencil.getCenterCoord().asVec3d() 1703 - map.applyInverseMap(vectorFromSurface); 1704 return Vec3Type(result); 1705 } else { 1706 Vec3d location = map.applyMap(stencil.getCenterCoord().asVec3d()); 1707 Vec3d result = map.applyInverseMap(location - vectorFromSurface); 1708 return Vec3Type(result); 1709 } 1710 } 1711 }; 1712 1713 1714 /// @brief Compute the closest-point transform to a level set. 1715 /// @return the closest point to the surface from which the level set was derived, 1716 /// in the range space of the map (e.g., in world space) 1717 template<typename MapType, DScheme DiffScheme> 1718 struct CPT_RANGE 1719 { 1720 // random access version 1721 template<typename Accessor> static Vec3<typename Accessor::ValueType> 1722 result(const MapType& map, const Accessor& grid, const Coord& ijk) 1723 { 1724 using ValueType = typename Accessor::ValueType; 1725 using Vec3Type = Vec3<ValueType>; 1726 // current distance 1727 ValueType d = grid.getValue(ijk); 1728 // compute gradient in physical space where it is a unit normal 1729 // since the grid holds a distance level set. 1730 Vec3Type vectorFromSurface = 1731 d*Gradient<MapType,DiffScheme>::result(map, grid, ijk); 1732 Vec3d result = map.applyMap(ijk.asVec3d()) - vectorFromSurface; 1733 1734 return Vec3Type(result); 1735 } 1736 1737 // stencil access version 1738 template<typename StencilT> static Vec3<typename StencilT::ValueType> 1739 result(const MapType& map, const StencilT& stencil) 1740 { 1741 using ValueType = typename StencilT::ValueType; 1742 using Vec3Type = Vec3<ValueType>; 1743 // current distance 1744 ValueType d = stencil.template getValue<0, 0, 0>(); 1745 // compute gradient in physical space where it is a unit normal 1746 // since the grid holds a distance level set. 1747 Vec3Type vectorFromSurface = 1748 d*Gradient<MapType, DiffScheme>::result(map, stencil); 1749 Vec3d result = map.applyMap(stencil.getCenterCoord().asVec3d()) - vectorFromSurface; 1750 1751 return Vec3Type(result); 1752 } 1753 }; 1754 1755 1756 /// @brief Compute the mean curvature. 1757 /// @details The mean curvature is returned in two parts, @a alpha and @a beta, 1758 /// where @a alpha is the numerator in ∇ · (∇Φ / |∇Φ|) 1759 /// and @a beta is |∇Φ|. 1760 template<typename MapType, DDScheme DiffScheme2, DScheme DiffScheme1> 1761 struct MeanCurvature 1762 { 1763 /// @brief Random access version 1764 /// @return @c true if the gradient is nonzero, in which case the mean curvature 1765 /// is returned in two parts, @a alpha and @a beta, where @a alpha is the numerator 1766 /// in ∇ · (∇Φ / |∇Φ|) and @a beta is |∇Φ|. 1767 template<typename Accessor> 1768 static bool compute(const MapType& map, const Accessor& grid, const Coord& ijk, 1769 double& alpha, double& beta) 1770 { 1771 using ValueType = typename Accessor::ValueType; 1772 1773 // compute the gradient in index and world space 1774 Vec3d d1_is(static_cast<double>(D1<DiffScheme1>::inX(grid, ijk)), 1775 static_cast<double>(D1<DiffScheme1>::inY(grid, ijk)), 1776 static_cast<double>(D1<DiffScheme1>::inZ(grid, ijk))), d1_ws; 1777 if (is_linear<MapType>::value) {//resolved at compiletime 1778 d1_ws = map.applyIJT(d1_is); 1779 } else { 1780 d1_ws = map.applyIJT(d1_is, ijk.asVec3d()); 1781 } 1782 const double Dx2 = d1_ws(0)*d1_ws(0); 1783 const double Dy2 = d1_ws(1)*d1_ws(1); 1784 const double Dz2 = d1_ws(2)*d1_ws(2); 1785 const double normGrad = Dx2 + Dy2 + Dz2; 1786 if (normGrad <= math::Tolerance<double>::value()) { 1787 alpha = beta = 0; 1788 return false; 1789 } 1790 1791 // all the second derivatives in index space 1792 ValueType iddx = D2<DiffScheme2>::inX(grid, ijk); 1793 ValueType iddy = D2<DiffScheme2>::inY(grid, ijk); 1794 ValueType iddz = D2<DiffScheme2>::inZ(grid, ijk); 1795 1796 ValueType iddxy = D2<DiffScheme2>::inXandY(grid, ijk); 1797 ValueType iddyz = D2<DiffScheme2>::inYandZ(grid, ijk); 1798 ValueType iddxz = D2<DiffScheme2>::inXandZ(grid, ijk); 1799 1800 // second derivatives in index space 1801 Mat3d d2_is(iddx, iddxy, iddxz, 1802 iddxy, iddy, iddyz, 1803 iddxz, iddyz, iddz); 1804 1805 // convert second derivatives to world space 1806 Mat3d d2_ws; 1807 if (is_linear<MapType>::value) {//resolved at compiletime 1808 d2_ws = map.applyIJC(d2_is); 1809 } else { 1810 d2_ws = map.applyIJC(d2_is, d1_is, ijk.asVec3d()); 1811 } 1812 1813 // assemble the nominator and denominator for mean curvature 1814 alpha = (Dx2*(d2_ws(1,1)+d2_ws(2,2))+Dy2*(d2_ws(0,0)+d2_ws(2,2)) 1815 +Dz2*(d2_ws(0,0)+d2_ws(1,1)) 1816 -2*(d1_ws(0)*(d1_ws(1)*d2_ws(0,1)+d1_ws(2)*d2_ws(0,2)) 1817 +d1_ws(1)*d1_ws(2)*d2_ws(1,2))); 1818 beta = std::sqrt(normGrad); // * 1/dx 1819 return true; 1820 } 1821 1822 template<typename Accessor> 1823 static typename Accessor::ValueType result(const MapType& map, 1824 const Accessor& grid, const Coord& ijk) 1825 { 1826 using ValueType = typename Accessor::ValueType; 1827 double alpha, beta; 1828 return compute(map, grid, ijk, alpha, beta) ? 1829 ValueType(alpha/(2. *math::Pow3(beta))) : 0; 1830 } 1831 1832 template<typename Accessor> 1833 static typename Accessor::ValueType normGrad(const MapType& map, 1834 const Accessor& grid, const Coord& ijk) 1835 { 1836 using ValueType = typename Accessor::ValueType; 1837 double alpha, beta; 1838 return compute(map, grid, ijk, alpha, beta) ? 1839 ValueType(alpha/(2. *math::Pow2(beta))) : 0; 1840 } 1841 1842 /// @brief Stencil access version 1843 /// @return @c true if the gradient is nonzero, in which case the mean curvature 1844 /// is returned in two parts, @a alpha and @a beta, where @a alpha is the numerator 1845 /// in ∇ · (∇Φ / |∇Φ|) and @a beta is |∇Φ|. 1846 template<typename StencilT> 1847 static bool compute(const MapType& map, const StencilT& stencil, 1848 double& alpha, double& beta) 1849 { 1850 using ValueType = typename StencilT::ValueType; 1851 1852 // compute the gradient in index and world space 1853 Vec3d d1_is(D1<DiffScheme1>::inX(stencil), 1854 D1<DiffScheme1>::inY(stencil), 1855 D1<DiffScheme1>::inZ(stencil) ), d1_ws; 1856 if (is_linear<MapType>::value) {//resolved at compiletime 1857 d1_ws = map.applyIJT(d1_is); 1858 } else { 1859 d1_ws = map.applyIJT(d1_is, stencil.getCenterCoord().asVec3d()); 1860 } 1861 const double Dx2 = d1_ws(0)*d1_ws(0); 1862 const double Dy2 = d1_ws(1)*d1_ws(1); 1863 const double Dz2 = d1_ws(2)*d1_ws(2); 1864 const double normGrad = Dx2 + Dy2 + Dz2; 1865 if (normGrad <= math::Tolerance<double>::value()) { 1866 alpha = beta = 0; 1867 return false; 1868 } 1869 1870 // all the second derivatives in index space 1871 ValueType iddx = D2<DiffScheme2>::inX(stencil); 1872 ValueType iddy = D2<DiffScheme2>::inY(stencil); 1873 ValueType iddz = D2<DiffScheme2>::inZ(stencil); 1874 1875 ValueType iddxy = D2<DiffScheme2>::inXandY(stencil); 1876 ValueType iddyz = D2<DiffScheme2>::inYandZ(stencil); 1877 ValueType iddxz = D2<DiffScheme2>::inXandZ(stencil); 1878 1879 // second derivatives in index space 1880 Mat3d d2_is(iddx, iddxy, iddxz, 1881 iddxy, iddy, iddyz, 1882 iddxz, iddyz, iddz); 1883 1884 // convert second derivatives to world space 1885 Mat3d d2_ws; 1886 if (is_linear<MapType>::value) {//resolved at compiletime 1887 d2_ws = map.applyIJC(d2_is); 1888 } else { 1889 d2_ws = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d()); 1890 } 1891 1892 // for return 1893 alpha = (Dx2*(d2_ws(1,1)+d2_ws(2,2))+Dy2*(d2_ws(0,0)+d2_ws(2,2)) 1894 +Dz2*(d2_ws(0,0)+d2_ws(1,1)) 1895 -2*(d1_ws(0)*(d1_ws(1)*d2_ws(0,1)+d1_ws(2)*d2_ws(0,2)) 1896 +d1_ws(1)*d1_ws(2)*d2_ws(1,2))); 1897 beta = std::sqrt(normGrad); // * 1/dx 1898 return true; 1899 } 1900 1901 template<typename StencilT> 1902 static typename StencilT::ValueType 1903 result(const MapType& map, const StencilT stencil) 1904 { 1905 using ValueType = typename StencilT::ValueType; 1906 double alpha, beta; 1907 return compute(map, stencil, alpha, beta) ? 1908 ValueType(alpha/(2*math::Pow3(beta))) : 0; 1909 } 1910 1911 template<typename StencilT> 1912 static typename StencilT::ValueType normGrad(const MapType& map, const StencilT stencil) 1913 { 1914 using ValueType = typename StencilT::ValueType; 1915 double alpha, beta; 1916 return compute(map, stencil, alpha, beta) ? 1917 ValueType(alpha/(2*math::Pow2(beta))) : 0; 1918 } 1919 }; 1920 1921 1922 template<DDScheme DiffScheme2, DScheme DiffScheme1> 1923 struct MeanCurvature<TranslationMap, DiffScheme2, DiffScheme1> 1924 { 1925 // random access version 1926 template<typename Accessor> 1927 static typename Accessor::ValueType result(const TranslationMap&, 1928 const Accessor& grid, const Coord& ijk) 1929 { 1930 using ValueType = typename Accessor::ValueType; 1931 1932 ValueType alpha, beta; 1933 return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta) ? 1934 ValueType(alpha /(2*math::Pow3(beta))) : 0; 1935 } 1936 1937 template<typename Accessor> 1938 static typename Accessor::ValueType normGrad(const TranslationMap&, 1939 const Accessor& grid, const Coord& ijk) 1940 { 1941 using ValueType = typename Accessor::ValueType; 1942 1943 ValueType alpha, beta; 1944 return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta) ? 1945 ValueType(alpha/(2*math::Pow2(beta))) : 0; 1946 } 1947 1948 // stencil access version 1949 template<typename StencilT> 1950 static typename StencilT::ValueType result(const TranslationMap&, const StencilT& stencil) 1951 { 1952 using ValueType = typename StencilT::ValueType; 1953 1954 ValueType alpha, beta; 1955 return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta) ? 1956 ValueType(alpha /(2*math::Pow3(beta))) : 0; 1957 } 1958 1959 template<typename StencilT> 1960 static typename StencilT::ValueType normGrad(const TranslationMap&, const StencilT& stencil) 1961 { 1962 using ValueType = typename StencilT::ValueType; 1963 1964 ValueType alpha, beta; 1965 return ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta) ? 1966 ValueType(alpha/(2*math::Pow2(beta))) : 0; 1967 } 1968 }; 1969 1970 1971 template<DDScheme DiffScheme2, DScheme DiffScheme1> 1972 struct MeanCurvature<UniformScaleMap, DiffScheme2, DiffScheme1> 1973 { 1974 // random access version 1975 template<typename Accessor> 1976 static typename Accessor::ValueType result(const UniformScaleMap& map, 1977 const Accessor& grid, const Coord& ijk) 1978 { 1979 using ValueType = typename Accessor::ValueType; 1980 1981 ValueType alpha, beta; 1982 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) { 1983 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 1984 return ValueType(alpha*inv2dx/math::Pow3(beta)); 1985 } 1986 return 0; 1987 } 1988 1989 template<typename Accessor> 1990 static typename Accessor::ValueType normGrad(const UniformScaleMap& map, 1991 const Accessor& grid, const Coord& ijk) 1992 { 1993 using ValueType = typename Accessor::ValueType; 1994 1995 ValueType alpha, beta; 1996 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) { 1997 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 1998 return ValueType(alpha*invdxdx/(2*math::Pow2(beta))); 1999 } 2000 return 0; 2001 } 2002 2003 // stencil access version 2004 template<typename StencilT> 2005 static typename StencilT::ValueType result(const UniformScaleMap& map, const StencilT& stencil) 2006 { 2007 using ValueType = typename StencilT::ValueType; 2008 2009 ValueType alpha, beta; 2010 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) { 2011 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 2012 return ValueType(alpha*inv2dx/math::Pow3(beta)); 2013 } 2014 return 0; 2015 } 2016 2017 template<typename StencilT> 2018 static typename StencilT::ValueType normGrad(const UniformScaleMap& map, const StencilT& stencil) 2019 { 2020 using ValueType = typename StencilT::ValueType; 2021 2022 ValueType alpha, beta; 2023 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) { 2024 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 2025 return ValueType(alpha*invdxdx/(2*math::Pow2(beta))); 2026 } 2027 return 0; 2028 } 2029 }; 2030 2031 2032 template<DDScheme DiffScheme2, DScheme DiffScheme1> 2033 struct MeanCurvature<UniformScaleTranslateMap, DiffScheme2, DiffScheme1> 2034 { 2035 // random access version 2036 template<typename Accessor> static typename Accessor::ValueType 2037 result(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 2038 { 2039 using ValueType = typename Accessor::ValueType; 2040 2041 ValueType alpha, beta; 2042 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) { 2043 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 2044 return ValueType(alpha*inv2dx/math::Pow3(beta)); 2045 } 2046 return 0; 2047 } 2048 2049 template<typename Accessor> static typename Accessor::ValueType 2050 normGrad(const UniformScaleTranslateMap& map, const Accessor& grid, const Coord& ijk) 2051 { 2052 using ValueType = typename Accessor::ValueType; 2053 2054 ValueType alpha, beta; 2055 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(grid, ijk, alpha, beta)) { 2056 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 2057 return ValueType(alpha*invdxdx/(2*math::Pow2(beta))); 2058 } 2059 return 0; 2060 } 2061 2062 // stencil access version 2063 template<typename StencilT> static typename StencilT::ValueType 2064 result(const UniformScaleTranslateMap& map, const StencilT& stencil) 2065 { 2066 using ValueType = typename StencilT::ValueType; 2067 2068 ValueType alpha, beta; 2069 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) { 2070 ValueType inv2dx = ValueType(map.getInvTwiceScale()[0]); 2071 return ValueType(alpha*inv2dx/math::Pow3(beta)); 2072 } 2073 return 0; 2074 } 2075 2076 template<typename StencilT> static typename StencilT::ValueType 2077 normGrad(const UniformScaleTranslateMap& map, const StencilT& stencil) 2078 { 2079 using ValueType = typename StencilT::ValueType; 2080 2081 ValueType alpha, beta; 2082 if (ISMeanCurvature<DiffScheme2, DiffScheme1>::result(stencil, alpha, beta)) { 2083 ValueType invdxdx = ValueType(map.getInvScaleSqr()[0]); 2084 return ValueType(alpha*invdxdx/(2*math::Pow2(beta))); 2085 } 2086 return 0; 2087 } 2088 }; 2089 2090 2091 /// @brief A wrapper that holds a MapBase::ConstPtr and exposes a reduced set 2092 /// of functionality needed by the mathematical operators 2093 /// @details This may be used in some <tt>Map</tt>-templated code, when the overhead of 2094 /// actually resolving the @c Map type is large compared to the map work to be done. 2095 class GenericMap 2096 { 2097 public: 2098 template<typename GridType> 2099 GenericMap(const GridType& g): mMap(g.transform().baseMap()) {} 2100 2101 GenericMap(const Transform& t): mMap(t.baseMap()) {} 2102 GenericMap(MapBase::Ptr map): mMap(ConstPtrCast<const MapBase>(map)) {} 2103 GenericMap(MapBase::ConstPtr map): mMap(map) {} 2104 ~GenericMap() {} 2105 2106 Vec3d applyMap(const Vec3d& in) const { return mMap->applyMap(in); } 2107 Vec3d applyInverseMap(const Vec3d& in) const { return mMap->applyInverseMap(in); } 2108 2109 Vec3d applyIJT(const Vec3d& in) const { return mMap->applyIJT(in); } 2110 Vec3d applyIJT(const Vec3d& in, const Vec3d& pos) const { return mMap->applyIJT(in, pos); } 2111 Mat3d applyIJC(const Mat3d& m) const { return mMap->applyIJC(m); } 2112 Mat3d applyIJC(const Mat3d& m, const Vec3d& v, const Vec3d& pos) const 2113 { return mMap->applyIJC(m,v,pos); } 2114 2115 double determinant() const { return mMap->determinant(); } 2116 double determinant(const Vec3d& in) const { return mMap->determinant(in); } 2117 2118 Vec3d voxelSize() const { return mMap->voxelSize(); } 2119 Vec3d voxelSize(const Vec3d&v) const { return mMap->voxelSize(v); } 2120 2121 private: 2122 MapBase::ConstPtr mMap; 2123 }; 2124 2125 } // end math namespace 2126 } // namespace OPENVDB_VERSION_NAME 2127 } // end openvdb namespace 2128 2129 #endif // OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED 2130