ROOT   6.08/07 Reference Guide
Transform3D.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: W. Brown, M. Fischler, L. Moneta 2005
3
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10
11 // implementation file for class Transform3D
12 //
13 // Created by: Lorenzo Moneta October 27 2005
14 //
15 //
17
19 #include "Math/GenVector/Plane3D.h"
20
21 #include <cmath>
22 #include <algorithm>
23
24
25
26
27 namespace ROOT {
28
29 namespace Math {
30
31
34
35
36 // ========== Constructors and Assignment =====================
37
38
39
40 // construct from two ref frames
41 Transform3D::Transform3D(const XYZPoint & fr0, const XYZPoint & fr1, const XYZPoint & fr2,
42  const XYZPoint & to0, const XYZPoint & to1, const XYZPoint & to2 )
43 {
44  // takes impl. from CLHEP ( E.Chernyaev). To be checked
45
46  XYZVector x1,y1,z1, x2,y2,z2;
47  x1 = (fr1 - fr0).Unit();
48  y1 = (fr2 - fr0).Unit();
49  x2 = (to1 - to0).Unit();
50  y2 = (to2 - to0).Unit();
51
52  // C H E C K A N G L E S
53
54  double cos1, cos2;
55  cos1 = x1.Dot(y1);
56  cos2 = x2.Dot(y2);
57
58  if (std::fabs(1.0-cos1) <= 0.000001 || std::fabs(1.0-cos2) <= 0.000001) {
59  std::cerr << "Transform3D: Error : zero angle between axes" << std::endl;
60  SetIdentity();
61  } else {
62  if (std::fabs(cos1-cos2) > 0.000001) {
63  std::cerr << "Transform3D: Warning: angles between axes are not equal"
64  << std::endl;
65  }
66
67  // F I N D R O T A T I O N M A T R I X
68
69  z1 = (x1.Cross(y1)).Unit();
70  y1 = z1.Cross(x1);
71
72  z2 = (x2.Cross(y2)).Unit();
73  y2 = z2.Cross(x2);
74
75  double x1x = x1.x(); double x1y = x1.y(); double x1z = x1.z();
76  double y1x = y1.x(); double y1y = y1.y(); double y1z = y1.z();
77  double z1x = z1.x(); double z1y = z1.y(); double z1z = z1.z();
78
79  double x2x = x2.x(); double x2y = x2.y(); double x2z = x2.z();
80  double y2x = y2.x(); double y2y = y2.y(); double y2z = y2.z();
81  double z2x = z2.x(); double z2y = z2.y(); double z2z = z2.z();
82
83  double detxx = (y1y *z1z - z1y *y1z );
84  double detxy = -(y1x *z1z - z1x *y1z );
85  double detxz = (y1x *z1y - z1x *y1y );
86  double detyx = -(x1y *z1z - z1y *x1z );
87  double detyy = (x1x *z1z - z1x *x1z );
88  double detyz = -(x1x *z1y - z1x *x1y );
89  double detzx = (x1y *y1z - y1y *x1z );
90  double detzy = -(x1x *y1z - y1x *x1z );
91  double detzz = (x1x *y1y - y1x *x1y );
92
93  double txx = x2x *detxx + y2x *detyx + z2x *detzx;
94  double txy = x2x *detxy + y2x *detyy + z2x *detzy;
95  double txz = x2x *detxz + y2x *detyz + z2x *detzz;
96  double tyx = x2y *detxx + y2y *detyx + z2y *detzx;
97  double tyy = x2y *detxy + y2y *detyy + z2y *detzy;
98  double tyz = x2y *detxz + y2y *detyz + z2y *detzz;
99  double tzx = x2z *detxx + y2z *detyx + z2z *detzx;
100  double tzy = x2z *detxy + y2z *detyy + z2z *detzy;
101  double tzz = x2z *detxz + y2z *detyz + z2z *detzz;
102
103  // S E T T R A N S F O R M A T I O N
104
105  double dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
106  double dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
107
108  SetComponents(txx, txy, txz, dx2-txx*dx1-txy*dy1-txz*dz1,
109  tyx, tyy, tyz, dy2-tyx*dx1-tyy*dy1-tyz*dz1,
110  tzx, tzy, tzz, dz2-tzx*dx1-tzy*dy1-tzz*dz1);
111  }
112 }
113
114
115 // inversion (from CLHEP)
117 {
118  //
119  // Name: Transform3D::inverse Date: 24.09.96
120  // Author: E.Chernyaev (IHEP/Protvino) Revised:
121  //
122  // Function: Find inverse affine transformation.
123
124  double detxx = fM[kYY]*fM[kZZ] - fM[kYZ]*fM[kZY];
125  double detxy = fM[kYX]*fM[kZZ] - fM[kYZ]*fM[kZX];
126  double detxz = fM[kYX]*fM[kZY] - fM[kYY]*fM[kZX];
127  double det = fM[kXX]*detxx - fM[kXY]*detxy + fM[kXZ]*detxz;
128  if (det == 0) {
129  std::cerr << "Transform3D::inverse error: zero determinant" << std::endl;
130  return;
131  }
132  det = 1./det; detxx *= det; detxy *= det; detxz *= det;
133  double detyx = (fM[kXY]*fM[kZZ] - fM[kXZ]*fM[kZY] )*det;
134  double detyy = (fM[kXX]*fM[kZZ] - fM[kXZ]*fM[kZX] )*det;
135  double detyz = (fM[kXX]*fM[kZY] - fM[kXY]*fM[kZX] )*det;
136  double detzx = (fM[kXY]*fM[kYZ] - fM[kXZ]*fM[kYY] )*det;
137  double detzy = (fM[kXX]*fM[kYZ] - fM[kXZ]*fM[kYX] )*det;
138  double detzz = (fM[kXX]*fM[kYY] - fM[kXY]*fM[kYX] )*det;
140  (detxx, -detyx, detzx, -detxx*fM[kDX]+detyx*fM[kDY]-detzx*fM[kDZ],
141  -detxy, detyy, -detzy, detxy*fM[kDX]-detyy*fM[kDY]+detzy*fM[kDZ],
142  detxz, -detyz, detzz, -detxz*fM[kDX]+detyz*fM[kDY]-detzz*fM[kDZ]);
143 }
144
145
146
147
149 {
150  //set identity ( identity rotation and zero translation)
151  fM[kXX] = 1.0; fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kDX] = 0.0;
152  fM[kYX] = 0.0; fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kDY] = 0.0;
153  fM[kZX] = 0.0; fM[kZY] = 0.0; fM[kZZ] = 1.0; fM[kDZ] = 0.0;
154 }
155
156
157 void Transform3D::AssignFrom (const Rotation3D & r, const XYZVector & v)
158 {
159  // assignment from rotation + translation
160
161  double rotData[9];
162  r.GetComponents(rotData, rotData +9);
163  // first raw
164  for (int i = 0; i < 3; ++i)
165  fM[i] = rotData[i];
166  // second raw
167  for (int i = 0; i < 3; ++i)
168  fM[kYX+i] = rotData[3+i];
169  // third raw
170  for (int i = 0; i < 3; ++i)
171  fM[kZX+i] = rotData[6+i];
172
173  // translation data
174  double vecData[3];
175  v.GetCoordinates(vecData, vecData+3);
176  fM[kDX] = vecData[0];
177  fM[kDY] = vecData[1];
178  fM[kDZ] = vecData[2];
179 }
180
181
183 {
184  // assign from only a rotation (null translation)
185  double rotData[9];
186  r.GetComponents(rotData, rotData +9);
187  for (int i = 0; i < 3; ++i) {
188  for (int j = 0; j < 3; ++j)
189  fM[4*i + j] = rotData[3*i+j];
190  // empty vector data
191  fM[4*i + 3] = 0;
192  }
193 }
194
195 void Transform3D::AssignFrom(const XYZVector & v)
196 {
197  // assign from a translation only (identity rotations)
198  fM[kXX] = 1.0; fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kDX] = v.X();
199  fM[kYX] = 0.0; fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kDY] = v.Y();
200  fM[kZX] = 0.0; fM[kZY] = 0.0; fM[kZZ] = 1.0; fM[kDZ] = v.Z();
201 }
202
204 {
205  // transformations on a 3D plane
206  XYZVector n = plane.Normal();
207  // take a point on the plane. Use origin projection on the plane
208  // ( -ad, -bd, -cd) if (a**2 + b**2 + c**2 ) = 1
209  double d = plane.HesseDistance();
210  XYZPoint p( - d * n.X() , - d *n.Y(), -d *n.Z() );
211  return Plane3D ( operator() (n), operator() (p) );
212 }
213
214 std::ostream & operator<< (std::ostream & os, const Transform3D & t)
215 {
216  // TODO - this will need changing for machine-readable issues
217  // and even the human readable form needs formatting improvements
218
219  double m[12];
220  t.GetComponents(m, m+12);
221  os << "\n" << m[0] << " " << m[1] << " " << m[2] << " " << m[3] ;
222  os << "\n" << m[4] << " " << m[5] << " " << m[6] << " " << m[7] ;
223  os << "\n" << m[8] << " " << m[9] << " " << m[10]<< " " << m[11] << "\n";
224  return os;
225 }
226
227 } // end namespace Math
228 } // end namespace ROOT
Class describing a geometrical plane in 3 dimensions.
Definition: Plane3D.h:47
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
void GetComponents(ForeignVector &v1, ForeignVector &v2, ForeignVector &v3) const
Get components into three vectors which will be the (orthonormal) columns of the rotation matrix...
Definition: Rotation3D.h:249
PositionVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > Point
Definition: Transform3D.h:91
PositionVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZPoint
3D Point based on the cartesian coordinates x,y,z in double precision
Definition: Point3Dfwd.h:33
void SetIdentity()
Set identity transformation (identity rotation , zero translation)
void GetComponents(IT begin, IT end) const
Get the 12 matrix components into data specified by an iterator begin and another to the end of the d...
Definition: Transform3D.h:331
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > Vector
Definition: Transform3D.h:90
std::ostream & operator<<(std::ostream &os, const AxisAngle &a)
Stream Output and Input.
Definition: AxisAngle.cxx:91
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
static const double x2[5]
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
Vector Normal() const
Return normal vector to the plane as Cartesian DisplacementVector.
Definition: Plane3D.h:166
Transform3D()
Default constructor (identy rotation) + zero translation.
Definition: Transform3D.h:105
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
void AssignFrom(const Rotation3D &r, const Vector &v)
make transformation from first a rotation then a translation
DisplacementVector3D Cross(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return vector (cross) product of two displacement vectors, as a vector in the coordinate system of th...
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
Definition: Rotation3D.h:65
TMarker * m
Definition: textangle.C:8
Scalar HesseDistance() const
Return the Hesse Distance (distance from the origin) of the plane or the d coefficient expressed in n...
Definition: Plane3D.h:174
Basic 3D Transformation class describing a rotation and then a translation The internal data are a 3D...
Definition: Transform3D.h:85
void GetCoordinates(Scalar &a, Scalar &b, Scalar &c) const
get internal data into 3 Scalar numbers
void Invert()
Invert the transformation in place.
static const double x1[5]
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
3D Vector based on the cartesian coordinates x,y,z in double precision
Definition: Vector3Dfwd.h:34
Namespace for new Math classes and functions.
Point operator()(const Point &p) const
Transformation operation for Position Vector in Cartesian coordinate For a Position Vector first a ro...
Definition: Transform3D.h:473
void SetComponents(IT begin, IT end)
Set the 12 matrix components given an iterator to the start of the desired data, and another to the e...
Definition: Transform3D.h:314
SVector< T, D > Unit(const SVector< T, D > &rhs)
Unit.
Definition: Functions.h:383
const Int_t n
Definition: legend1.C:16
Scalar Dot(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return the scalar (dot) product of two displacement vectors.