ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 
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 
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
PositionVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > Point
Definition: Transform3D.h:91
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
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)
Class describing a generic position vector (point) in 3 dimensions.
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 GetCoordinates(Scalar &a, Scalar &b, Scalar &c) const
get internal data into 3 Scalar numbers
Scalar Dot(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return the scalar (dot) product of two displacement vectors.
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
static const double x2[5]
Vector Normal() const
Return normal vector to the plane as Cartesian DisplacementVector.
Definition: Plane3D.h:166
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
int d
Definition: tornado.py:11
TThread * t[5]
Definition: threadsh1.C:13
Transform3D()
Default constructor (identy rotation) + zero translation.
Definition: Transform3D.h:105
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
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
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
ROOT::R::TRInterface & r
Definition: Object.C:4
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
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
Definition: Rotation3D.h:65
TMarker * m
Definition: textangle.C:8
Basic 3D Transformation class describing a rotation and then a translation The internal data are a 3D...
Definition: Transform3D.h:85
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
void Invert()
Invert the transformation in place.
static const double x1[5]
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
3D Vector based on the cartesian coordinates x,y,z in double precision
Definition: Vector3Dfwd.h:34
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...
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