File indexing completed on 2025-01-18 10:10:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #ifndef ROOT_Math_GenVector_Transform3D
0017 #define ROOT_Math_GenVector_Transform3D 1
0018
0019
0020
0021 #include "Math/GenVector/DisplacementVector3D.h"
0022
0023 #include "Math/GenVector/PositionVector3D.h"
0024
0025 #include "Math/GenVector/Rotation3D.h"
0026
0027 #include "Math/GenVector/Translation3D.h"
0028
0029
0030 #include "Math/GenVector/AxisAnglefwd.h"
0031 #include "Math/GenVector/EulerAnglesfwd.h"
0032 #include "Math/GenVector/Quaternionfwd.h"
0033 #include "Math/GenVector/RotationZYXfwd.h"
0034 #include "Math/GenVector/RotationXfwd.h"
0035 #include "Math/GenVector/RotationYfwd.h"
0036 #include "Math/GenVector/RotationZfwd.h"
0037
0038 #include <iostream>
0039 #include <type_traits>
0040 #include <cmath>
0041
0042
0043
0044
0045
0046 namespace ROOT {
0047
0048 namespace Math {
0049
0050 namespace Impl {
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079 template <typename T = double>
0080 class Transform3D {
0081
0082 public:
0083 typedef T Scalar;
0084
0085 typedef DisplacementVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Vector;
0086 typedef PositionVector3D<Cartesian3D<T>, DefaultCoordinateSystemTag> Point;
0087
0088 enum ETransform3DMatrixIndex {
0089 kXX = 0, kXY = 1, kXZ = 2, kDX = 3,
0090 kYX = 4, kYY = 5, kYZ = 6, kDY = 7,
0091 kZX = 8, kZY = 9, kZZ =10, kDZ = 11
0092 };
0093
0094
0095
0096
0097
0098
0099 Transform3D()
0100 {
0101 SetIdentity();
0102 }
0103
0104
0105
0106
0107
0108 template<class IT>
0109 Transform3D(IT begin, IT end)
0110 {
0111 SetComponents(begin,end);
0112 }
0113
0114
0115
0116
0117 Transform3D( const Rotation3D & r, const Vector & v)
0118 {
0119 AssignFrom( r, v );
0120 }
0121
0122
0123
0124 Transform3D(const Rotation3D &r, const Translation3D<T> &t) { AssignFrom(r, t.Vect()); }
0125
0126
0127
0128
0129
0130
0131
0132 template <class ARotation, class CoordSystem, class Tag>
0133 Transform3D( const ARotation & r, const DisplacementVector3D<CoordSystem,Tag> & v)
0134 {
0135 AssignFrom( Rotation3D(r), Vector (v.X(),v.Y(),v.Z()) );
0136 }
0137
0138
0139
0140
0141
0142
0143
0144 template <class ARotation>
0145 Transform3D(const ARotation &r, const Translation3D<T> &t)
0146 {
0147 AssignFrom( Rotation3D(r), t.Vect() );
0148 }
0149
0150
0151 #ifdef OLD_VERSION
0152
0153
0154
0155 Transform3D( const Vector & v, const Rotation3D & r)
0156 {
0157
0158 AssignFrom( r, r(v) );
0159 }
0160 #endif
0161
0162
0163
0164
0165 explicit constexpr Transform3D( const Rotation3D & r) {
0166 AssignFrom(r);
0167 }
0168
0169
0170
0171
0172 explicit constexpr Transform3D( const AxisAngle & r) {
0173 AssignFrom(Rotation3D(r));
0174 }
0175 explicit constexpr Transform3D( const EulerAngles & r) {
0176 AssignFrom(Rotation3D(r));
0177 }
0178 explicit constexpr Transform3D( const Quaternion & r) {
0179 AssignFrom(Rotation3D(r));
0180 }
0181 explicit constexpr Transform3D( const RotationZYX & r) {
0182 AssignFrom(Rotation3D(r));
0183 }
0184
0185
0186
0187 explicit constexpr Transform3D( const RotationX & r) {
0188 AssignFrom(Rotation3D(r));
0189 }
0190 explicit constexpr Transform3D( const RotationY & r) {
0191 AssignFrom(Rotation3D(r));
0192 }
0193 explicit constexpr Transform3D( const RotationZ & r) {
0194 AssignFrom(Rotation3D(r));
0195 }
0196
0197
0198
0199
0200
0201 template<class CoordSystem, class Tag>
0202 explicit constexpr Transform3D( const DisplacementVector3D<CoordSystem,Tag> & v) {
0203 AssignFrom(Vector(v.X(),v.Y(),v.Z()));
0204 }
0205
0206
0207
0208
0209 explicit constexpr Transform3D( const Vector & v) {
0210 AssignFrom(v);
0211 }
0212
0213
0214
0215
0216 explicit constexpr Transform3D(const Translation3D<T> &t) { AssignFrom(t.Vect()); }
0217
0218
0219
0220
0221 #ifdef OLD_VERSION
0222
0223
0224
0225
0226
0227
0228 template <class ARotation, class CoordSystem, class Tag>
0229 Transform3D(const DisplacementVector3D<CoordSystem,Tag> & v , const ARotation & r)
0230 {
0231
0232 Rotation3D r3d(r);
0233 AssignFrom( r3d, r3d( Vector(v.X(),v.Y(),v.Z()) ) );
0234 }
0235 #endif
0236
0237 public:
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250 template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0251 Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
0252 const Point &to2)
0253 {
0254
0255
0256 Vector x1 = (fr1 - fr0).Unit();
0257 Vector y1 = (fr2 - fr0).Unit();
0258 Vector x2 = (to1 - to0).Unit();
0259 Vector y2 = (to2 - to0).Unit();
0260
0261
0262
0263 const T cos1 = x1.Dot(y1);
0264 const T cos2 = x2.Dot(y2);
0265
0266 if (std::fabs(T(1) - cos1) <= T(0.000001) || std::fabs(T(1) - cos2) <= T(0.000001)) {
0267 std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
0268 SetIdentity();
0269 } else {
0270 if (std::fabs(cos1 - cos2) > T(0.000001)) {
0271 std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
0272 }
0273
0274
0275
0276 Vector z1 = (x1.Cross(y1)).Unit();
0277 y1 = z1.Cross(x1);
0278
0279 Vector z2 = (x2.Cross(y2)).Unit();
0280 y2 = z2.Cross(x2);
0281
0282 T x1x = x1.x();
0283 T x1y = x1.y();
0284 T x1z = x1.z();
0285 T y1x = y1.x();
0286 T y1y = y1.y();
0287 T y1z = y1.z();
0288 T z1x = z1.x();
0289 T z1y = z1.y();
0290 T z1z = z1.z();
0291
0292 T x2x = x2.x();
0293 T x2y = x2.y();
0294 T x2z = x2.z();
0295 T y2x = y2.x();
0296 T y2y = y2.y();
0297 T y2z = y2.z();
0298 T z2x = z2.x();
0299 T z2y = z2.y();
0300 T z2z = z2.z();
0301
0302 T detxx = (y1y * z1z - z1y * y1z);
0303 T detxy = -(y1x * z1z - z1x * y1z);
0304 T detxz = (y1x * z1y - z1x * y1y);
0305 T detyx = -(x1y * z1z - z1y * x1z);
0306 T detyy = (x1x * z1z - z1x * x1z);
0307 T detyz = -(x1x * z1y - z1x * x1y);
0308 T detzx = (x1y * y1z - y1y * x1z);
0309 T detzy = -(x1x * y1z - y1x * x1z);
0310 T detzz = (x1x * y1y - y1x * x1y);
0311
0312 T txx = x2x * detxx + y2x * detyx + z2x * detzx;
0313 T txy = x2x * detxy + y2x * detyy + z2x * detzy;
0314 T txz = x2x * detxz + y2x * detyz + z2x * detzz;
0315 T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
0316 T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
0317 T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
0318 T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
0319 T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
0320 T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
0321
0322
0323
0324 T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
0325 T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
0326
0327 SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
0328 dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
0329 }
0330 }
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344 template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0345 Transform3D(const Point &fr0, const Point &fr1, const Point &fr2, const Point &to0, const Point &to1,
0346 const Point &to2)
0347 {
0348
0349
0350 Vector x1 = (fr1 - fr0).Unit();
0351 Vector y1 = (fr2 - fr0).Unit();
0352 Vector x2 = (to1 - to0).Unit();
0353 Vector y2 = (to2 - to0).Unit();
0354
0355
0356
0357 const T cos1 = x1.Dot(y1);
0358 const T cos2 = x2.Dot(y2);
0359
0360 const auto m1 = (abs(T(1) - cos1) <= T(0.000001) || abs(T(1) - cos2) <= T(0.000001));
0361
0362 const auto m2 = (abs(cos1 - cos2) > T(0.000001));
0363 if (any_of(m2)) {
0364 std::cerr << "Transform3D: Warning: angles between axes are not equal" << std::endl;
0365 }
0366
0367
0368
0369 Vector z1 = (x1.Cross(y1)).Unit();
0370 y1 = z1.Cross(x1);
0371
0372 Vector z2 = (x2.Cross(y2)).Unit();
0373 y2 = z2.Cross(x2);
0374
0375 T x1x = x1.x();
0376 T x1y = x1.y();
0377 T x1z = x1.z();
0378 T y1x = y1.x();
0379 T y1y = y1.y();
0380 T y1z = y1.z();
0381 T z1x = z1.x();
0382 T z1y = z1.y();
0383 T z1z = z1.z();
0384
0385 T x2x = x2.x();
0386 T x2y = x2.y();
0387 T x2z = x2.z();
0388 T y2x = y2.x();
0389 T y2y = y2.y();
0390 T y2z = y2.z();
0391 T z2x = z2.x();
0392 T z2y = z2.y();
0393 T z2z = z2.z();
0394
0395 T detxx = (y1y * z1z - z1y * y1z);
0396 T detxy = -(y1x * z1z - z1x * y1z);
0397 T detxz = (y1x * z1y - z1x * y1y);
0398 T detyx = -(x1y * z1z - z1y * x1z);
0399 T detyy = (x1x * z1z - z1x * x1z);
0400 T detyz = -(x1x * z1y - z1x * x1y);
0401 T detzx = (x1y * y1z - y1y * x1z);
0402 T detzy = -(x1x * y1z - y1x * x1z);
0403 T detzz = (x1x * y1y - y1x * x1y);
0404
0405 T txx = x2x * detxx + y2x * detyx + z2x * detzx;
0406 T txy = x2x * detxy + y2x * detyy + z2x * detzy;
0407 T txz = x2x * detxz + y2x * detyz + z2x * detzz;
0408 T tyx = x2y * detxx + y2y * detyx + z2y * detzx;
0409 T tyy = x2y * detxy + y2y * detyy + z2y * detzy;
0410 T tyz = x2y * detxz + y2y * detyz + z2y * detzz;
0411 T tzx = x2z * detxx + y2z * detyx + z2z * detzx;
0412 T tzy = x2z * detxy + y2z * detyy + z2z * detzy;
0413 T tzz = x2z * detxz + y2z * detyz + z2z * detzz;
0414
0415
0416
0417 T dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
0418 T dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
0419
0420 SetComponents(txx, txy, txz, dx2 - txx * dx1 - txy * dy1 - txz * dz1, tyx, tyy, tyz,
0421 dy2 - tyx * dx1 - tyy * dy1 - tyz * dz1, tzx, tzy, tzz, dz2 - tzx * dx1 - tzy * dy1 - tzz * dz1);
0422
0423 if (any_of(m1)) {
0424 std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
0425 SetIdentity(m1);
0426 }
0427 }
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437 template<class ForeignMatrix>
0438 explicit constexpr Transform3D(const ForeignMatrix & m) {
0439 SetComponents(m);
0440 }
0441
0442
0443
0444
0445 Transform3D(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
0446 {
0447 SetComponents (xx, xy, xz, dx, yx, yy, yz, dy, zx, zy, zz, dz);
0448 }
0449
0450
0451
0452
0453
0454
0455
0456
0457 template <class ForeignMatrix>
0458 Transform3D<T> &operator=(const ForeignMatrix &m)
0459 {
0460 SetComponents(m);
0461 return *this;
0462 }
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472 template<class IT>
0473 void SetComponents(IT begin, IT end) {
0474 for (int i = 0; i <12; ++i) {
0475 fM[i] = *begin;
0476 ++begin;
0477 }
0478 (void)end;
0479 assert (end==begin);
0480 }
0481
0482
0483
0484
0485
0486 template<class IT>
0487 void GetComponents(IT begin, IT end) const {
0488 for (int i = 0; i <12; ++i) {
0489 *begin = fM[i];
0490 ++begin;
0491 }
0492 (void)end;
0493 assert (end==begin);
0494 }
0495
0496
0497
0498
0499 template<class IT>
0500 void GetComponents(IT begin) const {
0501 std::copy(fM, fM + 12, begin);
0502 }
0503
0504
0505
0506
0507
0508
0509
0510 template<class ForeignMatrix>
0511 void
0512 SetTransformMatrix (const ForeignMatrix & m) {
0513 fM[kXX]=m(0,0); fM[kXY]=m(0,1); fM[kXZ]=m(0,2); fM[kDX]=m(0,3);
0514 fM[kYX]=m(1,0); fM[kYY]=m(1,1); fM[kYZ]=m(1,2); fM[kDY]=m(1,3);
0515 fM[kZX]=m(2,0); fM[kZY]=m(2,1); fM[kZZ]=m(2,2); fM[kDZ]=m(2,3);
0516 }
0517
0518
0519
0520
0521
0522
0523 template<class ForeignMatrix>
0524 void
0525 GetTransformMatrix (ForeignMatrix & m) const {
0526 m(0,0)=fM[kXX]; m(0,1)=fM[kXY]; m(0,2)=fM[kXZ]; m(0,3)=fM[kDX];
0527 m(1,0)=fM[kYX]; m(1,1)=fM[kYY]; m(1,2)=fM[kYZ]; m(1,3)=fM[kDY];
0528 m(2,0)=fM[kZX]; m(2,1)=fM[kZY]; m(2,2)=fM[kZZ]; m(2,3)=fM[kDZ];
0529 }
0530
0531
0532
0533
0534
0535 void SetComponents(T xx, T xy, T xz, T dx, T yx, T yy, T yz, T dy, T zx, T zy, T zz, T dz)
0536 {
0537 fM[kXX]=xx; fM[kXY]=xy; fM[kXZ]=xz; fM[kDX]=dx;
0538 fM[kYX]=yx; fM[kYY]=yy; fM[kYZ]=yz; fM[kDY]=dy;
0539 fM[kZX]=zx; fM[kZY]=zy; fM[kZZ]=zz; fM[kDZ]=dz;
0540 }
0541
0542
0543
0544
0545 void GetComponents(T &xx, T &xy, T &xz, T &dx, T &yx, T &yy, T &yz, T &dy, T &zx, T &zy, T &zz, T &dz) const
0546 {
0547 xx=fM[kXX]; xy=fM[kXY]; xz=fM[kXZ]; dx=fM[kDX];
0548 yx=fM[kYX]; yy=fM[kYY]; yz=fM[kYZ]; dy=fM[kDY];
0549 zx=fM[kZX]; zy=fM[kZY]; zz=fM[kZZ]; dz=fM[kDZ];
0550 }
0551
0552
0553
0554
0555
0556
0557 template<class AnyRotation, class V>
0558 void GetDecomposition(AnyRotation &r, V &v) const {
0559 GetRotation(r);
0560 GetTranslation(v);
0561 }
0562
0563
0564
0565
0566
0567 void GetDecomposition(Rotation3D &r, Vector &v) const {
0568 GetRotation(r);
0569 GetTranslation(v);
0570 }
0571
0572
0573
0574
0575 Rotation3D Rotation() const {
0576 return Rotation3D( fM[kXX], fM[kXY], fM[kXZ],
0577 fM[kYX], fM[kYY], fM[kYZ],
0578 fM[kZX], fM[kZY], fM[kZZ] );
0579 }
0580
0581
0582
0583
0584 template <class AnyRotation>
0585 AnyRotation Rotation() const {
0586 return AnyRotation(Rotation3D(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]));
0587 }
0588
0589
0590
0591
0592 template <class AnyRotation>
0593 void GetRotation(AnyRotation &r) const {
0594 r = Rotation();
0595 }
0596
0597
0598
0599
0600 Translation3D<T> Translation() const { return Translation3D<T>(fM[kDX], fM[kDY], fM[kDZ]); }
0601
0602
0603
0604
0605
0606 template <class AnyVector>
0607 void GetTranslation(AnyVector &v) const {
0608 v.SetXYZ(fM[kDX], fM[kDY], fM[kDZ]);
0609 }
0610
0611
0612
0613
0614
0615
0616
0617
0618
0619 Point operator() (const Point & p) const {
0620 return Point ( fM[kXX]*p.X() + fM[kXY]*p.Y() + fM[kXZ]*p.Z() + fM[kDX],
0621 fM[kYX]*p.X() + fM[kYY]*p.Y() + fM[kYZ]*p.Z() + fM[kDY],
0622 fM[kZX]*p.X() + fM[kZY]*p.Y() + fM[kZZ]*p.Z() + fM[kDZ] );
0623 }
0624
0625
0626
0627
0628
0629
0630 Vector operator() (const Vector & v) const {
0631 return Vector( fM[kXX]*v.X() + fM[kXY]*v.Y() + fM[kXZ]*v.Z() ,
0632 fM[kYX]*v.X() + fM[kYY]*v.Y() + fM[kYZ]*v.Z() ,
0633 fM[kZX]*v.X() + fM[kZY]*v.Y() + fM[kZZ]*v.Z() );
0634 }
0635
0636
0637
0638
0639
0640 template <class CoordSystem>
0641 PositionVector3D<CoordSystem> operator()(const PositionVector3D<CoordSystem> &p) const
0642 {
0643 return PositionVector3D<CoordSystem>(operator()(Point(p)));
0644 }
0645
0646
0647
0648 template <class CoordSystem>
0649 PositionVector3D<CoordSystem> operator*(const PositionVector3D<CoordSystem> &v) const
0650 {
0651 return operator()(v);
0652 }
0653
0654
0655
0656
0657 template<class CoordSystem >
0658 DisplacementVector3D<CoordSystem> operator() (const DisplacementVector3D <CoordSystem> & v) const {
0659 return DisplacementVector3D<CoordSystem>(operator()(Vector(v)));
0660 }
0661
0662
0663
0664 template <class CoordSystem>
0665 DisplacementVector3D<CoordSystem> operator*(const DisplacementVector3D<CoordSystem> &v) const
0666 {
0667 return operator()(v);
0668 }
0669
0670
0671
0672
0673
0674
0675 Vector ApplyInverse(const Vector &v) const
0676 {
0677 return Vector(fM[kXX] * v.X() + fM[kYX] * v.Y() + fM[kZX] * v.Z(),
0678 fM[kXY] * v.X() + fM[kYY] * v.Y() + fM[kZY] * v.Z(),
0679 fM[kXZ] * v.X() + fM[kYZ] * v.Y() + fM[kZZ] * v.Z());
0680 }
0681
0682
0683
0684
0685
0686
0687
0688 Point ApplyInverse(const Point &p) const
0689 {
0690 Point tmp(p.X() - fM[kDX], p.Y() - fM[kDY], p.Z() - fM[kDZ]);
0691 return Point(fM[kXX] * tmp.X() + fM[kYX] * tmp.Y() + fM[kZX] * tmp.Z(),
0692 fM[kXY] * tmp.X() + fM[kYY] * tmp.Y() + fM[kZY] * tmp.Z(),
0693 fM[kXZ] * tmp.X() + fM[kYZ] * tmp.Y() + fM[kZZ] * tmp.Z());
0694 }
0695
0696
0697
0698
0699
0700
0701 template <class CoordSystem>
0702 PositionVector3D<CoordSystem> ApplyInverse(const PositionVector3D<CoordSystem> &p) const
0703 {
0704 return PositionVector3D<CoordSystem>(ApplyInverse(Point(p)));
0705 }
0706
0707
0708
0709
0710
0711
0712 template <class CoordSystem>
0713 DisplacementVector3D<CoordSystem> ApplyInverse(const DisplacementVector3D<CoordSystem> &p) const
0714 {
0715 return DisplacementVector3D<CoordSystem>(ApplyInverse(Vector(p)));
0716 }
0717
0718
0719
0720
0721 template <class CoordSystem, class Tag1, class Tag2>
0722 void Transform(const PositionVector3D<CoordSystem, Tag1> &p1, PositionVector3D<CoordSystem, Tag2> &p2) const
0723 {
0724 const Point xyzNew = operator()(Point(p1.X(), p1.Y(), p1.Z()));
0725 p2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
0726 }
0727
0728
0729
0730
0731
0732 template <class CoordSystem, class Tag1, class Tag2>
0733 void Transform(const DisplacementVector3D<CoordSystem, Tag1> &v1, DisplacementVector3D<CoordSystem, Tag2> &v2) const
0734 {
0735 const Vector xyzNew = operator()(Vector(v1.X(), v1.Y(), v1.Z()));
0736 v2.SetXYZ( xyzNew.X(), xyzNew.Y(), xyzNew.Z() );
0737 }
0738
0739
0740
0741
0742 template <class CoordSystem >
0743 LorentzVector<CoordSystem> operator() (const LorentzVector<CoordSystem> & q) const {
0744 const Vector xyzNew = operator()(Vector(q.Vect()));
0745 return LorentzVector<CoordSystem>(xyzNew.X(), xyzNew.Y(), xyzNew.Z(), q.E());
0746 }
0747
0748
0749
0750 template <class CoordSystem>
0751 LorentzVector<CoordSystem> operator*(const LorentzVector<CoordSystem> &q) const
0752 {
0753 return operator()(q);
0754 }
0755
0756
0757
0758
0759 template <typename TYPE>
0760 Plane3D<TYPE> operator()(const Plane3D<TYPE> &plane) const
0761 {
0762
0763 const auto n = plane.Normal();
0764
0765
0766 const auto d = plane.HesseDistance();
0767 Point p(-d * n.X(), -d * n.Y(), -d * n.Z());
0768 return Plane3D<TYPE>(operator()(n), operator()(p));
0769 }
0770
0771
0772 template <typename TYPE>
0773 Plane3D<TYPE> operator*(const Plane3D<TYPE> &plane) const
0774 {
0775 return operator()(plane);
0776 }
0777
0778
0779
0780
0781
0782
0783 inline Transform3D<T> &operator*=(const Transform3D<T> &t);
0784
0785
0786
0787
0788 inline Transform3D<T> operator*(const Transform3D<T> &t) const;
0789
0790
0791
0792
0793 template <typename SCALAR = T, typename std::enable_if<std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0794 void Invert()
0795 {
0796
0797
0798
0799
0800
0801
0802 T detxx = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
0803 T detxy = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
0804 T detxz = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
0805 T det = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
0806 if (det == T(0)) {
0807 std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
0808 return;
0809 }
0810 det = T(1) / det;
0811 detxx *= det;
0812 detxy *= det;
0813 detxz *= det;
0814 T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
0815 T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
0816 T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
0817 T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
0818 T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
0819 T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
0820 SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
0821 detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
0822 -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
0823 }
0824
0825
0826
0827
0828 template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0829 void Invert()
0830 {
0831
0832
0833
0834
0835
0836
0837 T detxx = fM[kYY] * fM[kZZ] - fM[kYZ] * fM[kZY];
0838 T detxy = fM[kYX] * fM[kZZ] - fM[kYZ] * fM[kZX];
0839 T detxz = fM[kYX] * fM[kZY] - fM[kYY] * fM[kZX];
0840 T det = fM[kXX] * detxx - fM[kXY] * detxy + fM[kXZ] * detxz;
0841 const auto detZmask = (det == T(0));
0842 if (any_of(detZmask)) {
0843 std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
0844 det(detZmask) = T(1);
0845 }
0846 det = T(1) / det;
0847 detxx *= det;
0848 detxy *= det;
0849 detxz *= det;
0850 T detyx = (fM[kXY] * fM[kZZ] - fM[kXZ] * fM[kZY]) * det;
0851 T detyy = (fM[kXX] * fM[kZZ] - fM[kXZ] * fM[kZX]) * det;
0852 T detyz = (fM[kXX] * fM[kZY] - fM[kXY] * fM[kZX]) * det;
0853 T detzx = (fM[kXY] * fM[kYZ] - fM[kXZ] * fM[kYY]) * det;
0854 T detzy = (fM[kXX] * fM[kYZ] - fM[kXZ] * fM[kYX]) * det;
0855 T detzz = (fM[kXX] * fM[kYY] - fM[kXY] * fM[kYX]) * det;
0856
0857 if (any_of(detZmask)) {
0858 detxx(detZmask) = T(0);
0859 detxy(detZmask) = T(0);
0860 detxz(detZmask) = T(0);
0861 detyx(detZmask) = T(0);
0862 detyy(detZmask) = T(0);
0863 detyz(detZmask) = T(0);
0864 detzx(detZmask) = T(0);
0865 detzy(detZmask) = T(0);
0866 detzz(detZmask) = T(0);
0867 }
0868
0869 SetComponents(detxx, -detyx, detzx, -detxx * fM[kDX] + detyx * fM[kDY] - detzx * fM[kDZ], -detxy, detyy, -detzy,
0870 detxy * fM[kDX] - detyy * fM[kDY] + detzy * fM[kDZ], detxz, -detyz, detzz,
0871 -detxz * fM[kDX] + detyz * fM[kDY] - detzz * fM[kDZ]);
0872 }
0873
0874
0875
0876
0877 Transform3D<T> Inverse() const
0878 {
0879 Transform3D<T> t(*this);
0880 t.Invert();
0881 return t;
0882 }
0883
0884
0885
0886
0887
0888 bool operator==(const Transform3D<T> &rhs) const
0889 {
0890 return (fM[0] == rhs.fM[0] && fM[1] == rhs.fM[1] && fM[2] == rhs.fM[2] && fM[3] == rhs.fM[3] &&
0891 fM[4] == rhs.fM[4] && fM[5] == rhs.fM[5] && fM[6] == rhs.fM[6] && fM[7] == rhs.fM[7] &&
0892 fM[8] == rhs.fM[8] && fM[9] == rhs.fM[9] && fM[10] == rhs.fM[10] && fM[11] == rhs.fM[11]);
0893 }
0894
0895
0896
0897
0898
0899 bool operator!=(const Transform3D<T> &rhs) const { return !operator==(rhs); }
0900
0901 protected:
0902
0903
0904
0905
0906 void AssignFrom(const Rotation3D &r, const Vector &v)
0907 {
0908
0909
0910 T rotData[9];
0911 r.GetComponents(rotData, rotData + 9);
0912
0913 for (int i = 0; i < 3; ++i) fM[i] = rotData[i];
0914
0915 for (int i = 0; i < 3; ++i) fM[kYX + i] = rotData[3 + i];
0916
0917 for (int i = 0; i < 3; ++i) fM[kZX + i] = rotData[6 + i];
0918
0919
0920 T vecData[3];
0921 v.GetCoordinates(vecData, vecData + 3);
0922 fM[kDX] = vecData[0];
0923 fM[kDY] = vecData[1];
0924 fM[kDZ] = vecData[2];
0925 }
0926
0927
0928
0929
0930 void AssignFrom(const Rotation3D &r)
0931 {
0932
0933 T rotData[9];
0934 r.GetComponents(rotData, rotData + 9);
0935 for (int i = 0; i < 3; ++i) {
0936 for (int j = 0; j < 3; ++j) fM[4 * i + j] = rotData[3 * i + j];
0937
0938 fM[4 * i + 3] = T(0);
0939 }
0940 }
0941
0942
0943
0944
0945 void AssignFrom(const Vector &v)
0946 {
0947
0948 fM[kXX] = T(1);
0949 fM[kXY] = T(0);
0950 fM[kXZ] = T(0);
0951 fM[kDX] = v.X();
0952 fM[kYX] = T(0);
0953 fM[kYY] = T(1);
0954 fM[kYZ] = T(0);
0955 fM[kDY] = v.Y();
0956 fM[kZX] = T(0);
0957 fM[kZY] = T(0);
0958 fM[kZZ] = T(1);
0959 fM[kDZ] = v.Z();
0960 }
0961
0962
0963
0964
0965 void SetIdentity()
0966 {
0967
0968 fM[kXX] = T(1);
0969 fM[kXY] = T(0);
0970 fM[kXZ] = T(0);
0971 fM[kDX] = T(0);
0972 fM[kYX] = T(0);
0973 fM[kYY] = T(1);
0974 fM[kYZ] = T(0);
0975 fM[kDY] = T(0);
0976 fM[kZX] = T(0);
0977 fM[kZY] = T(0);
0978 fM[kZZ] = T(1);
0979 fM[kDZ] = T(0);
0980 }
0981
0982
0983
0984
0985
0986 template <typename SCALAR = T, typename std::enable_if<!std::is_arithmetic<SCALAR>::value>::type * = nullptr>
0987 void SetIdentity(const typename SCALAR::mask_type m)
0988 {
0989
0990 fM[kXX](m) = T(1);
0991 fM[kXY](m) = T(0);
0992 fM[kXZ](m) = T(0);
0993 fM[kDX](m) = T(0);
0994 fM[kYX](m) = T(0);
0995 fM[kYY](m) = T(1);
0996 fM[kYZ](m) = T(0);
0997 fM[kDY](m) = T(0);
0998 fM[kZX](m) = T(0);
0999 fM[kZY](m) = T(0);
1000 fM[kZZ](m) = T(1);
1001 fM[kDZ](m) = T(0);
1002 }
1003
1004 private:
1005 T fM[12];
1006 };
1007
1008
1009
1010
1011
1012
1013 template <class T>
1014 inline Transform3D<T> &Transform3D<T>::operator*=(const Transform3D<T> &t)
1015 {
1016
1017
1018 SetComponents(fM[kXX]*t.fM[kXX]+fM[kXY]*t.fM[kYX]+fM[kXZ]*t.fM[kZX],
1019 fM[kXX]*t.fM[kXY]+fM[kXY]*t.fM[kYY]+fM[kXZ]*t.fM[kZY],
1020 fM[kXX]*t.fM[kXZ]+fM[kXY]*t.fM[kYZ]+fM[kXZ]*t.fM[kZZ],
1021 fM[kXX]*t.fM[kDX]+fM[kXY]*t.fM[kDY]+fM[kXZ]*t.fM[kDZ]+fM[kDX],
1022
1023 fM[kYX]*t.fM[kXX]+fM[kYY]*t.fM[kYX]+fM[kYZ]*t.fM[kZX],
1024 fM[kYX]*t.fM[kXY]+fM[kYY]*t.fM[kYY]+fM[kYZ]*t.fM[kZY],
1025 fM[kYX]*t.fM[kXZ]+fM[kYY]*t.fM[kYZ]+fM[kYZ]*t.fM[kZZ],
1026 fM[kYX]*t.fM[kDX]+fM[kYY]*t.fM[kDY]+fM[kYZ]*t.fM[kDZ]+fM[kDY],
1027
1028 fM[kZX]*t.fM[kXX]+fM[kZY]*t.fM[kYX]+fM[kZZ]*t.fM[kZX],
1029 fM[kZX]*t.fM[kXY]+fM[kZY]*t.fM[kYY]+fM[kZZ]*t.fM[kZY],
1030 fM[kZX]*t.fM[kXZ]+fM[kZY]*t.fM[kYZ]+fM[kZZ]*t.fM[kZZ],
1031 fM[kZX]*t.fM[kDX]+fM[kZY]*t.fM[kDY]+fM[kZZ]*t.fM[kDZ]+fM[kDZ]);
1032
1033 return *this;
1034 }
1035
1036 template <class T>
1037 inline Transform3D<T> Transform3D<T>::operator*(const Transform3D<T> &t) const
1038 {
1039
1040
1041 return Transform3D<T>(fM[kXX] * t.fM[kXX] + fM[kXY] * t.fM[kYX] + fM[kXZ] * t.fM[kZX],
1042 fM[kXX] * t.fM[kXY] + fM[kXY] * t.fM[kYY] + fM[kXZ] * t.fM[kZY],
1043 fM[kXX] * t.fM[kXZ] + fM[kXY] * t.fM[kYZ] + fM[kXZ] * t.fM[kZZ],
1044 fM[kXX] * t.fM[kDX] + fM[kXY] * t.fM[kDY] + fM[kXZ] * t.fM[kDZ] + fM[kDX],
1045
1046 fM[kYX] * t.fM[kXX] + fM[kYY] * t.fM[kYX] + fM[kYZ] * t.fM[kZX],
1047 fM[kYX] * t.fM[kXY] + fM[kYY] * t.fM[kYY] + fM[kYZ] * t.fM[kZY],
1048 fM[kYX] * t.fM[kXZ] + fM[kYY] * t.fM[kYZ] + fM[kYZ] * t.fM[kZZ],
1049 fM[kYX] * t.fM[kDX] + fM[kYY] * t.fM[kDY] + fM[kYZ] * t.fM[kDZ] + fM[kDY],
1050
1051 fM[kZX] * t.fM[kXX] + fM[kZY] * t.fM[kYX] + fM[kZZ] * t.fM[kZX],
1052 fM[kZX] * t.fM[kXY] + fM[kZY] * t.fM[kYY] + fM[kZZ] * t.fM[kZY],
1053 fM[kZX] * t.fM[kXZ] + fM[kZY] * t.fM[kYZ] + fM[kZZ] * t.fM[kZZ],
1054 fM[kZX] * t.fM[kDX] + fM[kZY] * t.fM[kDY] + fM[kZZ] * t.fM[kDZ] + fM[kDZ]);
1055 }
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070 template <class T>
1071 inline Transform3D<T> operator*(const Rotation3D &r, const Translation3D<T> &t)
1072 {
1073 return Transform3D<T>(r, r(t.Vect()));
1074 }
1075 template <class T>
1076 inline Transform3D<T> operator*(const RotationX &r, const Translation3D<T> &t)
1077 {
1078 Rotation3D r3(r);
1079 return Transform3D<T>(r3, r3(t.Vect()));
1080 }
1081 template <class T>
1082 inline Transform3D<T> operator*(const RotationY &r, const Translation3D<T> &t)
1083 {
1084 Rotation3D r3(r);
1085 return Transform3D<T>(r3, r3(t.Vect()));
1086 }
1087 template <class T>
1088 inline Transform3D<T> operator*(const RotationZ &r, const Translation3D<T> &t)
1089 {
1090 Rotation3D r3(r);
1091 return Transform3D<T>(r3, r3(t.Vect()));
1092 }
1093 template <class T>
1094 inline Transform3D<T> operator*(const RotationZYX &r, const Translation3D<T> &t)
1095 {
1096 Rotation3D r3(r);
1097 return Transform3D<T>(r3, r3(t.Vect()));
1098 }
1099 template <class T>
1100 inline Transform3D<T> operator*(const AxisAngle &r, const Translation3D<T> &t)
1101 {
1102 Rotation3D r3(r);
1103 return Transform3D<T>(r3, r3(t.Vect()));
1104 }
1105 template <class T>
1106 inline Transform3D<T> operator*(const EulerAngles &r, const Translation3D<T> &t)
1107 {
1108 Rotation3D r3(r);
1109 return Transform3D<T>(r3, r3(t.Vect()));
1110 }
1111 template <class T>
1112 inline Transform3D<T> operator*(const Quaternion &r, const Translation3D<T> &t)
1113 {
1114 Rotation3D r3(r);
1115 return Transform3D<T>(r3, r3(t.Vect()));
1116 }
1117
1118
1119
1120
1121
1122
1123
1124 template <class T>
1125 inline Transform3D<T> operator*(const Translation3D<T> &t, const Rotation3D &r)
1126 {
1127 return Transform3D<T>(r, t.Vect());
1128 }
1129 template <class T>
1130 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationX &r)
1131 {
1132 return Transform3D<T>(Rotation3D(r), t.Vect());
1133 }
1134 template <class T>
1135 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationY &r)
1136 {
1137 return Transform3D<T>(Rotation3D(r), t.Vect());
1138 }
1139 template <class T>
1140 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZ &r)
1141 {
1142 return Transform3D<T>(Rotation3D(r), t.Vect());
1143 }
1144 template <class T>
1145 inline Transform3D<T> operator*(const Translation3D<T> &t, const RotationZYX &r)
1146 {
1147 return Transform3D<T>(Rotation3D(r), t.Vect());
1148 }
1149 template <class T>
1150 inline Transform3D<T> operator*(const Translation3D<T> &t, const EulerAngles &r)
1151 {
1152 return Transform3D<T>(Rotation3D(r), t.Vect());
1153 }
1154 template <class T>
1155 inline Transform3D<T> operator*(const Translation3D<T> &t, const Quaternion &r)
1156 {
1157 return Transform3D<T>(Rotation3D(r), t.Vect());
1158 }
1159 template <class T>
1160 inline Transform3D<T> operator*(const Translation3D<T> &t, const AxisAngle &r)
1161 {
1162 return Transform3D<T>(Rotation3D(r), t.Vect());
1163 }
1164
1165
1166
1167
1168
1169
1170
1171 template <class T>
1172 inline Transform3D<T> operator*(const Transform3D<T> &t, const Translation3D<T> &d)
1173 {
1174 Rotation3D r = t.Rotation();
1175 return Transform3D<T>(r, r(d.Vect()) + t.Translation().Vect());
1176 }
1177
1178
1179
1180
1181
1182 template <class T>
1183 inline Transform3D<T> operator*(const Translation3D<T> &d, const Transform3D<T> &t)
1184 {
1185 return Transform3D<T>(t.Rotation(), t.Translation().Vect() + d.Vect());
1186 }
1187
1188
1189
1190
1191
1192
1193
1194
1195 template <class T>
1196 inline Transform3D<T> operator*(const Transform3D<T> &t, const Rotation3D &r)
1197 {
1198 return Transform3D<T>(t.Rotation() * r, t.Translation());
1199 }
1200 template <class T>
1201 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationX &r)
1202 {
1203 return Transform3D<T>(t.Rotation() * r, t.Translation());
1204 }
1205 template <class T>
1206 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationY &r)
1207 {
1208 return Transform3D<T>(t.Rotation() * r, t.Translation());
1209 }
1210 template <class T>
1211 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZ &r)
1212 {
1213 return Transform3D<T>(t.Rotation() * r, t.Translation());
1214 }
1215 template <class T>
1216 inline Transform3D<T> operator*(const Transform3D<T> &t, const RotationZYX &r)
1217 {
1218 return Transform3D<T>(t.Rotation() * r, t.Translation());
1219 }
1220 template <class T>
1221 inline Transform3D<T> operator*(const Transform3D<T> &t, const EulerAngles &r)
1222 {
1223 return Transform3D<T>(t.Rotation() * r, t.Translation());
1224 }
1225 template <class T>
1226 inline Transform3D<T> operator*(const Transform3D<T> &t, const AxisAngle &r)
1227 {
1228 return Transform3D<T>(t.Rotation() * r, t.Translation());
1229 }
1230 template <class T>
1231 inline Transform3D<T> operator*(const Transform3D<T> &t, const Quaternion &r)
1232 {
1233 return Transform3D<T>(t.Rotation() * r, t.Translation());
1234 }
1235
1236
1237
1238
1239
1240
1241
1242 template <class T>
1243 inline Transform3D<T> operator*(const Rotation3D &r, const Transform3D<T> &t)
1244 {
1245 return Transform3D<T>(r * t.Rotation(), r * t.Translation().Vect());
1246 }
1247 template <class T>
1248 inline Transform3D<T> operator*(const RotationX &r, const Transform3D<T> &t)
1249 {
1250 Rotation3D r3d(r);
1251 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1252 }
1253 template <class T>
1254 inline Transform3D<T> operator*(const RotationY &r, const Transform3D<T> &t)
1255 {
1256 Rotation3D r3d(r);
1257 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1258 }
1259 template <class T>
1260 inline Transform3D<T> operator*(const RotationZ &r, const Transform3D<T> &t)
1261 {
1262 Rotation3D r3d(r);
1263 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1264 }
1265 template <class T>
1266 inline Transform3D<T> operator*(const RotationZYX &r, const Transform3D<T> &t)
1267 {
1268 Rotation3D r3d(r);
1269 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1270 }
1271 template <class T>
1272 inline Transform3D<T> operator*(const EulerAngles &r, const Transform3D<T> &t)
1273 {
1274 Rotation3D r3d(r);
1275 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1276 }
1277 template <class T>
1278 inline Transform3D<T> operator*(const AxisAngle &r, const Transform3D<T> &t)
1279 {
1280 Rotation3D r3d(r);
1281 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1282 }
1283 template <class T>
1284 inline Transform3D<T> operator*(const Quaternion &r, const Transform3D<T> &t)
1285 {
1286 Rotation3D r3d(r);
1287 return Transform3D<T>(r3d * t.Rotation(), r3d * t.Translation().Vect());
1288 }
1289
1290
1291
1292
1293
1294
1295
1296
1297 template <class T>
1298 std::ostream &operator<<(std::ostream &os, const Transform3D<T> &t)
1299 {
1300
1301
1302
1303 T m[12];
1304 t.GetComponents(m, m + 12);
1305 os << "\n" << m[0] << " " << m[1] << " " << m[2] << " " << m[3];
1306 os << "\n" << m[4] << " " << m[5] << " " << m[6] << " " << m[7];
1307 os << "\n" << m[8] << " " << m[9] << " " << m[10] << " " << m[11] << "\n";
1308 return os;
1309 }
1310
1311 }
1312
1313
1314 typedef Impl::Transform3D<double> Transform3D;
1315 typedef Impl::Transform3D<float> Transform3DF;
1316
1317 }
1318
1319 }
1320
1321
1322 #endif