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