33 #ifndef OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
34 #define OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED
47 template<
typename Vec3T>
struct is_vec3d {
static const bool value =
false; };
50 template<
typename T>
struct is_double {
static const bool value =
false; };
51 template<>
struct is_double<double> {
static const bool value =
true; };
69 template<DScheme DiffScheme>
74 result(
const Accessor& grid,
const Coord& ijk)
76 typedef typename Accessor::ValueType ValueType;
85 result(
const StencilT& stencil)
87 typedef typename StencilT::ValueType ValueType;
99 template<BiasedGradientScheme bgs>
105 template<
typename Gr
idType>
116 template<
typename Gr
idType>
127 template<
typename Gr
idType>
137 template<
typename Gr
idType>
147 template<
typename Gr
idType>
157 template<
typename Gr
idType>
167 template<BiasedGradientScheme GradScheme,
typename Vec3Bias>
176 result(
const Accessor& grid,
const Coord& ijk,
const Vec3Bias& V)
178 typedef typename Accessor::ValueType ValueType;
188 result(
const StencilT& stencil,
const Vec3Bias& V)
190 typedef typename StencilT::ValueType ValueType;
200 template<BiasedGradientScheme GradScheme>
209 template<
typename Accessor>
210 static typename Accessor::ValueType
211 result(
const Accessor& grid,
const Coord& ijk)
213 typedef typename Accessor::ValueType ValueType;
222 template<
typename StencilT>
223 static typename StencilT::ValueType
224 result(
const StencilT& stencil)
226 typedef typename StencilT::ValueType ValueType;
235 #ifdef DWA_OPENVDB // for SIMD - note will do the computations in float
241 template<
typename Accessor>
242 static typename Accessor::ValueType
243 result(
const Accessor& grid,
const Coord& ijk)
245 typedef typename Accessor::ValueType ValueType;
250 v1(grid.getValue(ijk.
offsetBy(-2, 0, 0)) - grid.getValue(ijk.
offsetBy(-3, 0, 0)),
251 grid.getValue(ijk.
offsetBy( 0,-2, 0)) - grid.getValue(ijk.
offsetBy( 0,-3, 0)),
252 grid.getValue(ijk.
offsetBy( 0, 0,-2)) - grid.getValue(ijk.
offsetBy( 0, 0,-3)), 0),
253 v2(grid.getValue(ijk.
offsetBy(-1, 0, 0)) - grid.getValue(ijk.
offsetBy(-2, 0, 0)),
254 grid.getValue(ijk.
offsetBy( 0,-1, 0)) - grid.getValue(ijk.
offsetBy( 0,-2, 0)),
255 grid.getValue(ijk.
offsetBy( 0, 0,-1)) - grid.getValue(ijk.
offsetBy( 0, 0,-2)), 0),
256 v3(grid.getValue(ijk ) - grid.getValue(ijk.
offsetBy(-1, 0, 0)),
257 grid.getValue(ijk ) - grid.getValue(ijk.
offsetBy( 0,-1, 0)),
258 grid.getValue(ijk ) - grid.getValue(ijk.
offsetBy( 0, 0,-1)), 0),
259 v4(grid.getValue(ijk.
offsetBy( 1, 0, 0)) - grid.getValue(ijk ),
260 grid.getValue(ijk.
offsetBy( 0, 1, 0)) - grid.getValue(ijk ),
261 grid.getValue(ijk.
offsetBy( 0, 0, 1)) - grid.getValue(ijk ), 0),
262 v5(grid.getValue(ijk.
offsetBy( 2, 0, 0)) - grid.getValue(ijk.
offsetBy( 1, 0, 0)),
263 grid.getValue(ijk.
offsetBy( 0, 2, 0)) - grid.getValue(ijk.
offsetBy( 0, 1, 0)),
264 grid.getValue(ijk.
offsetBy( 0, 0, 2)) - grid.getValue(ijk.
offsetBy( 0, 0, 1)), 0),
265 v6(grid.getValue(ijk.
offsetBy( 3, 0, 0)) - grid.getValue(ijk.
offsetBy( 2, 0, 0)),
266 grid.getValue(ijk.
offsetBy( 0, 3, 0)) - grid.getValue(ijk.
offsetBy( 0, 2, 0)),
267 grid.getValue(ijk.
offsetBy( 0, 0, 3)) - grid.getValue(ijk.
offsetBy( 0, 0, 2)), 0),
275 template<
typename StencilT>
276 static typename StencilT::ValueType
277 result(
const StencilT& s)
279 typedef typename StencilT::ValueType ValueType;
284 v1(s.template getValue<-2, 0, 0>() - s.template getValue<-3, 0, 0>(),
285 s.template getValue< 0,-2, 0>() - s.template getValue< 0,-3, 0>(),
286 s.template getValue< 0, 0,-2>() - s.template getValue< 0, 0,-3>(), 0),
287 v2(s.template getValue<-1, 0, 0>() - s.template getValue<-2, 0, 0>(),
288 s.template getValue< 0,-1, 0>() - s.template getValue< 0,-2, 0>(),
289 s.template getValue< 0, 0,-1>() - s.template getValue< 0, 0,-2>(), 0),
290 v3(s.template getValue< 0, 0, 0>() - s.template getValue<-1, 0, 0>(),
291 s.template getValue< 0, 0, 0>() - s.template getValue< 0,-1, 0>(),
292 s.template getValue< 0, 0, 0>() - s.template getValue< 0, 0,-1>(), 0),
293 v4(s.template getValue< 1, 0, 0>() - s.template getValue< 0, 0, 0>(),
294 s.template getValue< 0, 1, 0>() - s.template getValue< 0, 0, 0>(),
295 s.template getValue< 0, 0, 1>() - s.template getValue< 0, 0, 0>(), 0),
296 v5(s.template getValue< 2, 0, 0>() - s.template getValue< 1, 0, 0>(),
297 s.template getValue< 0, 2, 0>() - s.template getValue< 0, 1, 0>(),
298 s.template getValue< 0, 0, 2>() - s.template getValue< 0, 0, 1>(), 0),
299 v6(s.template getValue< 3, 0, 0>() - s.template getValue< 2, 0, 0>(),
300 s.template getValue< 0, 3, 0>() - s.template getValue< 0, 2, 0>(),
301 s.template getValue< 0, 0, 3>() - s.template getValue< 0, 0, 2>(), 0),
308 #endif //DWA_OPENVDB // for SIMD - note will do the computations in float
314 template<DDScheme DiffScheme>
318 template<
typename Accessor>
319 static typename Accessor::ValueType result(
const Accessor& grid,
const Coord& ijk);
322 template<
typename StencilT>
323 static typename StencilT::ValueType result(
const StencilT& stencil);
331 template<
typename Accessor>
332 static typename Accessor::ValueType result(
const Accessor& grid,
const Coord& ijk)
334 return grid.getValue(ijk.
offsetBy(1,0,0)) + grid.getValue(ijk.
offsetBy(-1, 0, 0)) +
335 grid.getValue(ijk.
offsetBy(0,1,0)) + grid.getValue(ijk.
offsetBy(0, -1, 0)) +
337 - 6*grid.getValue(ijk);
341 template<
typename StencilT>
342 static typename StencilT::ValueType result(
const StencilT& stencil)
344 return stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
345 stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
346 stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>()
347 - 6*stencil.template getValue< 0, 0, 0>();
356 template<
typename Accessor>
357 static typename Accessor::ValueType result(
const Accessor& grid,
const Coord& ijk)
360 grid.getValue(ijk.
offsetBy(2,0,0)) + grid.getValue(ijk.
offsetBy(-2, 0, 0)) +
361 grid.getValue(ijk.
offsetBy(0,2,0)) + grid.getValue(ijk.
offsetBy( 0,-2, 0)) +
362 grid.getValue(ijk.
offsetBy(0,0,2)) + grid.getValue(ijk.
offsetBy( 0, 0,-2)) )
364 grid.getValue(ijk.
offsetBy(1,0,0)) + grid.getValue(ijk.
offsetBy(-1, 0, 0)) +
365 grid.getValue(ijk.
offsetBy(0,1,0)) + grid.getValue(ijk.
offsetBy( 0,-1, 0)) +
366 grid.getValue(ijk.
offsetBy(0,0,1)) + grid.getValue(ijk.
offsetBy( 0, 0,-1)) )
367 - 7.5*grid.getValue(ijk);
371 template<
typename StencilT>
372 static typename StencilT::ValueType result(
const StencilT& stencil)
375 stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() +
376 stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() +
377 stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() )
379 stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
380 stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
381 stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() )
382 - 7.5*stencil.template getValue< 0, 0, 0>();
390 template<
typename Accessor>
391 static typename Accessor::ValueType result(
const Accessor& grid,
const Coord& ijk)
394 grid.getValue(ijk.
offsetBy(3,0,0)) + grid.getValue(ijk.
offsetBy(-3, 0, 0)) +
395 grid.getValue(ijk.
offsetBy(0,3,0)) + grid.getValue(ijk.
offsetBy( 0,-3, 0)) +
396 grid.getValue(ijk.
offsetBy(0,0,3)) + grid.getValue(ijk.
offsetBy( 0, 0,-3)) )
398 grid.getValue(ijk.
offsetBy(2,0,0)) + grid.getValue(ijk.
offsetBy(-2, 0, 0)) +
399 grid.getValue(ijk.
offsetBy(0,2,0)) + grid.getValue(ijk.
offsetBy( 0,-2, 0)) +
400 grid.getValue(ijk.
offsetBy(0,0,2)) + grid.getValue(ijk.
offsetBy( 0, 0,-2)) )
402 grid.getValue(ijk.
offsetBy(1,0,0)) + grid.getValue(ijk.
offsetBy(-1, 0, 0)) +
403 grid.getValue(ijk.
offsetBy(0,1,0)) + grid.getValue(ijk.
offsetBy( 0,-1, 0)) +
404 grid.getValue(ijk.
offsetBy(0,0,1)) + grid.getValue(ijk.
offsetBy( 0, 0,-1)) )
405 - (3*49/18.)*grid.getValue(ijk);
409 template<
typename StencilT>
410 static typename StencilT::ValueType result(
const StencilT& stencil)
413 stencil.template getValue< 3, 0, 0>() + stencil.template getValue<-3, 0, 0>() +
414 stencil.template getValue< 0, 3, 0>() + stencil.template getValue< 0,-3, 0>() +
415 stencil.template getValue< 0, 0, 3>() + stencil.template getValue< 0, 0,-3>() )
417 stencil.template getValue< 2, 0, 0>() + stencil.template getValue<-2, 0, 0>() +
418 stencil.template getValue< 0, 2, 0>() + stencil.template getValue< 0,-2, 0>() +
419 stencil.template getValue< 0, 0, 2>() + stencil.template getValue< 0, 0,-2>() )
421 stencil.template getValue< 1, 0, 0>() + stencil.template getValue<-1, 0, 0>() +
422 stencil.template getValue< 0, 1, 0>() + stencil.template getValue< 0,-1, 0>() +
423 stencil.template getValue< 0, 0, 1>() + stencil.template getValue< 0, 0,-1>() )
424 - (3*49/18.)*stencil.template getValue< 0, 0, 0>();
432 template<DScheme DiffScheme>
436 template<
typename Accessor>
static typename Accessor::ValueType::value_type
437 result(
const Accessor& grid,
const Coord& ijk)
445 template<
typename StencilT>
static typename StencilT::ValueType::value_type
446 result(
const StencilT& stencil)
458 template<DScheme DiffScheme>
462 template<
typename Accessor>
463 static typename Accessor::ValueType result(
const Accessor& grid,
const Coord& ijk)
465 typedef typename Accessor::ValueType Vec3Type;
475 template<
typename StencilT>
476 static typename StencilT::ValueType result(
const StencilT& stencil)
478 typedef typename StencilT::ValueType Vec3Type;
492 template<DDScheme DiffScheme2, DScheme DiffScheme1>
496 template<
typename Accessor>
497 static void result(
const Accessor& grid,
const Coord& ijk,
498 typename Accessor::ValueType& alpha,
499 typename Accessor::ValueType& beta)
501 typedef typename Accessor::ValueType ValueType;
507 ValueType Dx2 = Dx*Dx;
508 ValueType Dy2 = Dy*Dy;
509 ValueType Dz2 = Dz*Dz;
520 alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
521 beta = ValueType(std::sqrt(
double(Dx2 + Dy2 + Dz2)));
525 template<
typename StencilT>
526 static void result(
const StencilT& stencil,
527 typename StencilT::ValueType& alpha,
528 typename StencilT::ValueType& beta)
530 typedef typename StencilT::ValueType ValueType;
535 ValueType Dx2 = Dx*Dx;
536 ValueType Dy2 = Dy*Dy;
537 ValueType Dz2 = Dz*Dz;
548 alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
549 beta = ValueType(std::sqrt(
double(Dx2 + Dy2 + Dz2)));
561 template<
typename MapType, DScheme DiffScheme>
565 template<
typename Accessor>
567 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk)
572 return Vec3Type(map.applyIJT(iGradient, ijk.
asVec3d()));
576 template<
typename StencilT>
578 result(
const MapType& map,
const StencilT& stencil)
583 return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d()));
589 template<DScheme DiffScheme>
593 template<
typename Accessor>
601 template<
typename StencilT>
615 template<
typename Accessor>
624 return iGradient * inv2dx;
628 template<
typename StencilT>
637 return iGradient * inv2dx;
647 template<
typename Accessor>
656 return iGradient * inv2dx;
660 template<
typename StencilT>
669 return iGradient * inv2dx;
679 template<
typename Accessor>
693 template<
typename StencilT>
695 result(
const ScaleMap& map,
const StencilT& stencil)
713 template<
typename Accessor>
727 template<
typename StencilT>
746 template<
typename MapType, BiasedGradientScheme GradScheme>
751 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk,
754 typedef typename Accessor::ValueType ValueType;
758 return Vec3Type(map.applyIJT(iGradient, ijk.
asVec3d()));
763 result(
const MapType& map,
const StencilT& stencil,
766 typedef typename StencilT::ValueType ValueType;
770 return Vec3Type(map.applyIJT(iGradient, stencil.getCenterCoord().asVec3d()));
777 template<
typename MapType, BiasedGradientScheme GradScheme>
786 template<
typename Accessor>
787 static typename Accessor::ValueType
788 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk)
790 typedef typename Accessor::ValueType ValueType;
799 template<
typename StencilT>
800 static typename StencilT::ValueType
801 result(
const MapType& map,
const StencilT& stencil)
803 typedef typename StencilT::ValueType ValueType;
814 template<BiasedGradientScheme GradScheme>
819 template<
typename Accessor>
820 static typename Accessor::ValueType
823 typedef typename Accessor::ValueType ValueType;
830 template<
typename StencilT>
831 static typename StencilT::ValueType
834 typedef typename StencilT::ValueType ValueType;
841 template<BiasedGradientScheme GradScheme>
846 template<
typename Accessor>
847 static typename Accessor::ValueType
850 typedef typename Accessor::ValueType ValueType;
857 template<
typename StencilT>
858 static typename StencilT::ValueType
861 typedef typename StencilT::ValueType ValueType;
872 template<
typename MapType, DScheme DiffScheme>
876 template<
typename Accessor>
static typename Accessor::ValueType::value_type
877 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk)
879 typedef typename Accessor::ValueType Vec3Type;
880 typedef typename Accessor::ValueType::value_type ValueType;
883 for (
int i=0; i < 3; i++) {
887 div += ValueType(map.applyIJT(vec, ijk.
asVec3d())[i]);
893 template<
typename StencilT>
static typename StencilT::ValueType::value_type
894 result(
const MapType& map,
const StencilT& stencil)
896 typedef typename StencilT::ValueType Vec3Type;
897 typedef typename StencilT::ValueType::value_type ValueType;
900 for (
int i=0; i < 3; i++) {
904 div += ValueType(map.applyIJT(vec, stencil.getCenterCoord().asVec3d())[i]);
912 template<DScheme DiffScheme>
916 template<
typename Accessor>
static typename Accessor::ValueType::value_type
919 typedef typename Accessor::ValueType Vec3Type;
920 typedef typename Accessor::ValueType::value_type ValueType;
928 template<
typename StencilT>
static typename StencilT::ValueType::value_type
931 typedef typename StencilT::ValueType Vec3Type;
932 typedef typename StencilT::ValueType::value_type ValueType;
942 template<DScheme DiffScheme>
946 template<
typename Accessor>
static typename Accessor::ValueType::value_type
949 typedef typename Accessor::ValueType Vec3Type;
950 typedef typename Accessor::ValueType::value_type ValueType;
960 template<
typename StencilT>
static typename StencilT::ValueType::value_type
963 typedef typename StencilT::ValueType Vec3Type;
964 typedef typename StencilT::ValueType::value_type ValueType;
976 template<DScheme DiffScheme>
980 template<
typename Accessor>
static typename Accessor::ValueType::value_type
983 typedef typename Accessor::ValueType Vec3Type;
984 typedef typename Accessor::ValueType::value_type ValueType;
994 template<
typename StencilT>
static typename StencilT::ValueType::value_type
997 typedef typename StencilT::ValueType Vec3Type;
998 typedef typename StencilT::ValueType::value_type ValueType;
1003 ValueType invdx = ValueType(map.
getInvScale()[0]);
1014 template<
typename Accessor>
static typename Accessor::ValueType::value_type
1017 typedef typename Accessor::ValueType Vec3Type;
1018 typedef typename Accessor::ValueType::value_type ValueType;
1023 return div * inv2dx;
1027 template<
typename StencilT>
static typename StencilT::ValueType::value_type
1030 typedef typename StencilT::ValueType Vec3Type;
1031 typedef typename StencilT::ValueType::value_type ValueType;
1036 return div * inv2dx;
1046 template<
typename Accessor>
static typename Accessor::ValueType::value_type
1049 typedef typename Accessor::ValueType Vec3Type;
1050 typedef typename Accessor::ValueType::value_type ValueType;
1056 return div * inv2dx;
1060 template<
typename StencilT>
static typename StencilT::ValueType::value_type
1063 typedef typename StencilT::ValueType Vec3Type;
1064 typedef typename StencilT::ValueType::value_type ValueType;
1070 return div * inv2dx;
1076 template<DScheme DiffScheme>
1080 template<
typename Accessor>
static typename Accessor::ValueType::value_type
1083 typedef typename Accessor::ValueType Vec3Type;
1084 typedef typename Accessor::ValueType::value_type ValueType;
1086 ValueType div = ValueType(
1094 template<
typename StencilT>
static typename StencilT::ValueType::value_type
1097 typedef typename StencilT::ValueType Vec3Type;
1098 typedef typename StencilT::ValueType::value_type ValueType;
1111 template<DScheme DiffScheme>
1115 template<
typename Accessor>
static typename Accessor::ValueType::value_type
1118 typedef typename Accessor::ValueType Vec3Type;
1119 typedef typename Accessor::ValueType::value_type ValueType;
1121 ValueType div = ValueType(
1129 template<
typename StencilT>
static typename StencilT::ValueType::value_type
1132 typedef typename StencilT::ValueType Vec3Type;
1133 typedef typename StencilT::ValueType::value_type ValueType;
1150 template<
typename Accessor>
static typename Accessor::ValueType::value_type
1153 typedef typename Accessor::ValueType Vec3Type;
1154 typedef typename Accessor::ValueType::value_type ValueType;
1156 ValueType div = ValueType(
1164 template<
typename StencilT>
static typename StencilT::ValueType::value_type
1167 typedef typename StencilT::ValueType Vec3Type;
1168 typedef typename StencilT::ValueType::value_type ValueType;
1170 ValueType div = ValueType(
1184 template<
typename Accessor>
static typename Accessor::ValueType::value_type
1187 typedef typename Accessor::ValueType Vec3Type;
1188 typedef typename Accessor::ValueType::value_type ValueType;
1190 ValueType div = ValueType(
1198 template<
typename StencilT>
static typename StencilT::ValueType::value_type
1201 typedef typename StencilT::ValueType Vec3Type;
1202 typedef typename StencilT::ValueType::value_type ValueType;
1204 ValueType div = ValueType(
1217 template<
typename MapType, DScheme DiffScheme>
1221 template<
typename Accessor>
static typename Accessor::ValueType
1222 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk)
1224 typedef typename Accessor::ValueType Vec3Type;
1226 for (
int i = 0; i < 3; i++)
1233 mat[i] = Vec3Type(map.applyIJT(vec, ijk.
asVec3d()));
1235 return Vec3Type(mat[2][1] - mat[1][2],
1236 mat[0][2] - mat[2][0],
1237 mat[1][0] - mat[0][1]);
1241 template<
typename StencilT>
static typename StencilT::ValueType
1242 result(
const MapType& map,
const StencilT& stencil)
1244 typedef typename StencilT::ValueType Vec3Type;
1246 for (
int i = 0; i < 3; i++)
1253 mat[i] = Vec3Type(map.applyIJT(vec, stencil.getCenterCoord().asVec3d()));
1255 return Vec3Type(mat[2][1] - mat[1][2],
1256 mat[0][2] - mat[2][0],
1257 mat[1][0] - mat[0][1]);
1262 template<DScheme DiffScheme>
1266 template<
typename Accessor>
static typename Accessor::ValueType
1269 typedef typename Accessor::ValueType Vec3Type;
1270 typedef typename Vec3Type::value_type ValueType;
1275 template<
typename StencilT>
static typename StencilT::ValueType
1278 typedef typename StencilT::ValueType Vec3Type;
1279 typedef typename Vec3Type::value_type ValueType;
1285 template<DScheme DiffScheme>
1289 template<
typename Accessor>
static typename Accessor::ValueType
1293 typedef typename Accessor::ValueType Vec3Type;
1294 typedef typename Vec3Type::value_type ValueType;
1300 template<
typename StencilT>
static typename StencilT::ValueType
1304 typedef typename StencilT::ValueType Vec3Type;
1305 typedef typename Vec3Type::value_type ValueType;
1315 template<
typename Accessor>
static typename Accessor::ValueType
1318 typedef typename Accessor::ValueType Vec3Type;
1319 typedef typename Vec3Type::value_type ValueType;
1325 template<
typename StencilT>
static typename StencilT::ValueType
1328 typedef typename StencilT::ValueType Vec3Type;
1329 typedef typename Vec3Type::value_type ValueType;
1339 template<
typename Accessor>
static typename Accessor::ValueType
1342 typedef typename Accessor::ValueType Vec3Type;
1343 typedef typename Vec3Type::value_type ValueType;
1349 template<
typename StencilT>
static typename StencilT::ValueType
1352 typedef typename StencilT::ValueType Vec3Type;
1353 typedef typename Vec3Type::value_type ValueType;
1364 template<
typename MapType, DDScheme DiffScheme>
1368 template<
typename Accessor>
1369 static typename Accessor::ValueType result(
const MapType& map,
1370 const Accessor& grid,
const Coord& ijk)
1372 typedef typename Accessor::ValueType ValueType;
1383 Mat3d d2_is(iddx, iddxy, iddxz,
1385 iddxz, iddyz, iddz);
1389 d2_rs = map.applyIJC(d2_is);
1396 d2_rs = map.applyIJC(d2_is, d1_is, ijk.
asVec3d());
1400 return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1404 template<
typename StencilT>
1405 static typename StencilT::ValueType result(
const MapType& map,
const StencilT& stencil)
1407 typedef typename StencilT::ValueType ValueType;
1418 Mat3d d2_is(iddx, iddxy, iddxz,
1420 iddxz, iddyz, iddz);
1424 d2_rs = map.applyIJC(d2_is);
1431 d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1435 return ValueType(d2_rs(0,0) + d2_rs(1,1) + d2_rs(2,2));
1440 template<DDScheme DiffScheme>
1444 template<
typename Accessor>
1446 const Accessor& grid,
const Coord& ijk)
1452 template<
typename StencilT>
1453 static typename StencilT::ValueType result(
const TranslationMap&,
const StencilT& stencil)
1461 template<DDScheme DiffScheme>
1465 template<
typename Accessor>
1467 const Accessor& grid,
const Coord& ijk)
1473 template<
typename StencilT>
1474 static typename StencilT::ValueType result(
const UnitaryMap&,
const StencilT& stencil)
1481 template<DDScheme DiffScheme>
1485 template<
typename Accessor>
static typename Accessor::ValueType
1488 typedef typename Accessor::ValueType ValueType;
1494 template<
typename StencilT>
static typename StencilT::ValueType
1497 typedef typename StencilT::ValueType ValueType;
1504 template<DDScheme DiffScheme>
1508 template<
typename Accessor>
static typename Accessor::ValueType
1511 typedef typename Accessor::ValueType ValueType;
1517 template<
typename StencilT>
static typename StencilT::ValueType
1520 typedef typename StencilT::ValueType ValueType;
1527 template<DDScheme DiffScheme>
1531 template<
typename Accessor>
static typename Accessor::ValueType
1534 typedef typename Accessor::ValueType ValueType;
1542 return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1546 template<
typename StencilT>
static typename StencilT::ValueType
1549 typedef typename StencilT::ValueType ValueType;
1557 return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1562 template<DDScheme DiffScheme>
1566 template<
typename Accessor>
static typename Accessor::ValueType
1569 typedef typename Accessor::ValueType ValueType;
1576 return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1580 template<
typename StencilT>
static typename StencilT::ValueType
1583 typedef typename StencilT::ValueType ValueType;
1590 return ValueType(iddx * invScaleSqr[0] + iddy * invScaleSqr[1] + iddz * invScaleSqr[2]);
1598 template<
typename MapType, DScheme DiffScheme>
1603 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk)
1605 typedef typename Accessor::ValueType ValueType;
1609 ValueType d = grid.getValue(ijk);
1614 Vec3d result = ijk.
asVec3d() - map.applyInverseMap(vectorFromSurface);
1615 return Vec3Type(result);
1618 Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1619 return Vec3Type(result);
1625 result(
const MapType& map,
const StencilT& stencil)
1627 typedef typename StencilT::ValueType ValueType;
1631 ValueType d = stencil.template getValue< 0, 0, 0>();
1636 Vec3d result = stencil.getCenterCoord().asVec3d() - map.applyInverseMap(vectorFromSurface);
1637 return Vec3Type(result);
1639 Vec3d location = map.applyMap(stencil.getCenterCoord().asVec3d());
1640 Vec3d result = map.applyInverseMap(location - vectorFromSurface);
1641 return Vec3Type(result);
1650 template<
typename MapType, DScheme DiffScheme>
1655 result(
const MapType& map,
const Accessor& grid,
const Coord& ijk)
1657 typedef typename Accessor::ValueType ValueType;
1660 ValueType d = grid.getValue(ijk);
1663 Vec3Type vectorFromSurface =
1665 Vec3d result = map.applyMap(ijk.
asVec3d()) - vectorFromSurface;
1667 return Vec3Type(result);
1672 result(
const MapType& map,
const StencilT& stencil)
1674 typedef typename StencilT::ValueType ValueType;
1677 ValueType d = stencil.template getValue< 0, 0, 0>();
1680 Vec3Type vectorFromSurface =
1682 Vec3d result = map.applyMap(stencil.getCenterCoord().asVec3d()) - vectorFromSurface;
1684 return Vec3Type(result);
1692 template<
typename MapType, DDScheme DiffScheme2, DScheme DiffScheme1>
1696 template<
typename Accessor>
1697 static void compute(
const MapType& map,
const Accessor& grid,
const Coord& ijk,
1698 double& alpha,
double& beta)
1700 typedef typename Accessor::ValueType ValueType;
1718 Mat3d d2_is(iddx, iddxy, iddxz,
1720 iddxz, iddyz, iddz);
1726 d2_rs = map.applyIJC(d2_is);
1727 d1_rs = map.applyIJT(d1_is);
1729 d2_rs = map.applyIJC(d2_is, d1_is, ijk.
asVec3d());
1730 d1_rs = map.applyIJT(d1_is, ijk.
asVec3d());
1734 double Dx2 = d1_rs(0)*d1_rs(0);
1735 double Dy2 = d1_rs(1)*d1_rs(1);
1736 double Dz2 = d1_rs(2)*d1_rs(2);
1739 alpha = (Dx2*(d2_rs(1,1)+d2_rs(2,2))+Dy2*(d2_rs(0,0)+d2_rs(2,2))
1740 +Dz2*(d2_rs(0,0)+d2_rs(1,1))
1741 -2*(d1_rs(0)*(d1_rs(1)*d2_rs(0,1)+d1_rs(2)*d2_rs(0,2))
1742 +d1_rs(1)*d1_rs(2)*d2_rs(1,2)));
1743 beta = std::sqrt(Dx2 + Dy2 + Dz2);
1746 template<
typename Accessor>
1747 static typename Accessor::ValueType result(
const MapType& map,
1748 const Accessor& grid,
const Coord& ijk)
1750 typedef typename Accessor::ValueType ValueType;
1752 compute(map, grid, ijk, alpha, beta);
1754 return ValueType(alpha/(2. *
math::Pow3(beta)));
1757 template<
typename Accessor>
1758 static typename Accessor::ValueType normGrad(
const MapType& map,
1759 const Accessor& grid,
const Coord& ijk)
1761 typedef typename Accessor::ValueType ValueType;
1763 compute(map, grid, ijk, alpha, beta);
1765 return ValueType(alpha/(2. *
math::Pow2(beta)));
1769 template<
typename StencilT>
1770 static void compute(
const MapType& map,
const StencilT& stencil,
double& alpha,
double& beta)
1772 typedef typename StencilT::ValueType ValueType;
1790 Mat3d d2_is(iddx, iddxy, iddxz,
1792 iddxz, iddyz, iddz);
1798 d2_rs = map.applyIJC(d2_is);
1799 d1_rs = map.applyIJT(d1_is);
1801 d2_rs = map.applyIJC(d2_is, d1_is, stencil.getCenterCoord().asVec3d());
1802 d1_rs = map.applyIJT(d1_is, stencil.getCenterCoord().asVec3d());
1806 double Dx2 = d1_rs(0)*d1_rs(0);
1807 double Dy2 = d1_rs(1)*d1_rs(1);
1808 double Dz2 = d1_rs(2)*d1_rs(2);
1811 alpha = (Dx2*(d2_rs(1,1)+d2_rs(2,2))+Dy2*(d2_rs(0,0)+d2_rs(2,2))
1812 +Dz2*(d2_rs(0,0)+d2_rs(1,1))
1813 -2*(d1_rs(0)*(d1_rs(1)*d2_rs(0,1)+d1_rs(2)*d2_rs(0,2))
1814 +d1_rs(1)*d1_rs(2)*d2_rs(1,2)));
1815 beta = std::sqrt(Dx2 + Dy2 + Dz2);
1818 template<
typename StencilT>
1819 static typename StencilT::ValueType
1820 result(
const MapType& map,
const StencilT stencil)
1822 typedef typename StencilT::ValueType ValueType;
1824 compute(map, stencil, alpha, beta);
1825 return ValueType(alpha/(2*
math::Pow3(beta)));
1828 template<
typename StencilT>
1829 static typename StencilT::ValueType normGrad(
const MapType& map,
const StencilT stencil)
1831 typedef typename StencilT::ValueType ValueType;
1833 compute(map, stencil, alpha, beta);
1835 return ValueType(alpha/(2*
math::Pow2(beta)));
1840 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1844 template<
typename Accessor>
1846 const Accessor& grid,
const Coord& ijk)
1848 typedef typename Accessor::ValueType ValueType;
1850 ValueType alpha, beta;
1853 return ValueType(alpha /(2*
math::Pow3(beta)));
1856 template<
typename Accessor>
1858 const Accessor& grid,
const Coord& ijk)
1860 typedef typename Accessor::ValueType ValueType;
1862 ValueType alpha, beta;
1865 return ValueType(alpha/(2*
math::Pow2(beta)));
1869 template<
typename StencilT>
1870 static typename StencilT::ValueType result(
const TranslationMap&,
const StencilT& stencil)
1872 typedef typename StencilT::ValueType ValueType;
1874 ValueType alpha, beta;
1877 return ValueType(alpha /(2*
math::Pow3(beta)));
1880 template<
typename StencilT>
1881 static typename StencilT::ValueType normGrad(
const TranslationMap&,
const StencilT& stencil)
1883 typedef typename StencilT::ValueType ValueType;
1885 ValueType alpha, beta;
1888 return ValueType(alpha/(2*
math::Pow2(beta)));
1893 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1897 template<
typename Accessor>
1899 const Accessor& grid,
const Coord& ijk)
1901 typedef typename Accessor::ValueType ValueType;
1903 ValueType alpha, beta;
1907 return ValueType(alpha*inv2dx/
math::Pow3(beta));
1910 template<
typename Accessor>
1912 const Accessor& grid,
const Coord& ijk)
1914 typedef typename Accessor::ValueType ValueType;
1916 ValueType alpha, beta;
1920 return ValueType(alpha*invdxdx/(2*
math::Pow2(beta)));
1924 template<
typename StencilT>
1925 static typename StencilT::ValueType result(
const UniformScaleMap& map,
const StencilT& stencil)
1927 typedef typename StencilT::ValueType ValueType;
1929 ValueType alpha, beta;
1933 return ValueType(alpha*inv2dx/
math::Pow3(beta));
1936 template<
typename StencilT>
1937 static typename StencilT::ValueType normGrad(
const UniformScaleMap& map,
const StencilT& stencil)
1939 typedef typename StencilT::ValueType ValueType;
1941 ValueType alpha, beta;
1945 return ValueType(alpha*invdxdx/(2*
math::Pow2(beta)));
1950 template<DDScheme DiffScheme2, DScheme DiffScheme1>
1954 template<
typename Accessor>
static typename Accessor::ValueType
1957 typedef typename Accessor::ValueType ValueType;
1959 ValueType alpha, beta;
1963 return ValueType(alpha*inv2dx/
math::Pow3(beta));
1966 template<
typename Accessor>
static typename Accessor::ValueType
1969 typedef typename Accessor::ValueType ValueType;
1971 ValueType alpha, beta;
1975 return ValueType(alpha*invdxdx/(2*
math::Pow2(beta)));
1979 template<
typename StencilT>
static typename StencilT::ValueType
1982 typedef typename StencilT::ValueType ValueType;
1984 ValueType alpha, beta;
1988 return ValueType(alpha*inv2dx/
math::Pow3(beta));
1991 template<
typename StencilT>
static typename StencilT::ValueType
1994 typedef typename StencilT::ValueType ValueType;
1996 ValueType alpha, beta;
2000 return ValueType(alpha*invdxdx/(2*
math::Pow2(beta)));
2013 template<
typename Gr
idType>
2028 {
return mMap->applyIJC(m,v,pos); }
2045 #endif // OPENVDB_MATH_OPERATORS_HAS_BEEN_INCLUDED