Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TEveTrans.cxx
Go to the documentation of this file.
1 // @(#)root/eve:$Id$
2 // Authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "TEveTrans.h"
13 #include "TEveUtil.h"
14 
15 #include "TBuffer.h"
16 #include "TClass.h"
17 #include "TMath.h"
18 
19 #include "Riostream.h"
20 
21 #include <cctype>
22 
23 #define F00 0
24 #define F01 4
25 #define F02 8
26 #define F03 12
27 
28 #define F10 1
29 #define F11 5
30 #define F12 9
31 #define F13 13
32 
33 #define F20 2
34 #define F21 6
35 #define F22 10
36 #define F23 14
37 
38 #define F30 3
39 #define F31 7
40 #define F32 11
41 #define F33 15
42 
43 /** \class TEveTrans
44 \ingroup TEve
45 TEveTrans is a 4x4 transformation matrix for homogeneous coordinates
46 stored internally in a column-major order to allow direct usage by
47 GL. The element type is Double32_t as statically the floats would
48 be precise enough but continuous operations on the matrix must
49 retain precision of column vectors.
50 
51 Cartan angles are stored in fA[1-3] (+z, -y, +x). They are
52 recalculated on demand.
53 
54 Direct element access (first two should be used with care):
55  - operator[i] direct access to elements, i:0->15
56  - CM(i,j) element 4*j + i; i,j:0->3 { CM ~ c-matrix }
57  - operator(i,j) element 4*(j-1) + i - 1 i,j:1->4
58 
59 Column-vector access:
60 USet Get/SetBaseVec(), Get/SetPos() and Arr[XYZT]() methods.
61 
62 For all methods taking the matrix indices:
63 1->X, 2->Y, 3->Z; 4->Position (if applicable). 0 reserved for time.
64 
65 Shorthands in method-names:
66 LF ~ LocalFrame; PF ~ ParentFrame; IP ~ InPlace
67 */
68 
69 ClassImp(TEveTrans);
70 
71 ////////////////////////////////////////////////////////////////////////////////
72 /// Default constructor.
73 
74 TEveTrans::TEveTrans() :
75  TObject(),
76  fA1(0), fA2(0), fA3(0), fAsOK(kFALSE),
77  fUseTrans (kTRUE),
78  fEditTrans(kFALSE),
79  fEditRotation(kTRUE),
80  fEditScale(kTRUE)
81 {
82  UnitTrans();
83 }
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 /// Constructor.
87 
88 TEveTrans::TEveTrans(const TEveTrans& t) :
89  TObject(),
90  fA1(t.fA1), fA2(t.fA2), fA3(t.fA3), fAsOK(t.fAsOK),
91  fUseTrans (t.fUseTrans),
92  fEditTrans(t.fEditTrans),
93  fEditRotation(kTRUE),
94  fEditScale(kTRUE)
95 {
96  SetTrans(t, kFALSE);
97 }
98 
99 ////////////////////////////////////////////////////////////////////////////////
100 /// Constructor.
101 
102 TEveTrans::TEveTrans(const Double_t arr[16]) :
103  TObject(),
104  fA1(0), fA2(0), fA3(0), fAsOK(kFALSE),
105  fUseTrans (kTRUE),
106  fEditTrans(kFALSE),
107  fEditRotation(kTRUE),
108  fEditScale(kTRUE)
109 {
110  SetFromArray(arr);
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 /// Constructor.
115 
116 TEveTrans::TEveTrans(const Float_t arr[16]) :
117  TObject(),
118  fA1(0), fA2(0), fA3(0), fAsOK(kFALSE),
119  fUseTrans (kTRUE),
120  fEditTrans(kFALSE),
121  fEditRotation(kTRUE),
122  fEditScale(kTRUE)
123 {
124  SetFromArray(arr);
125 }
126 
127 ////////////////////////////////////////////////////////////////////////////////
128 /// Reset matrix to unity.
129 
130 void TEveTrans::UnitTrans()
131 {
132  memset(fM, 0, 16*sizeof(Double_t));
133  fM[F00] = fM[F11] = fM[F22] = fM[F33] = 1;
134  fA1 = fA2 = fA3 = 0;
135  fAsOK = kTRUE;
136 }
137 
138 ////////////////////////////////////////////////////////////////////////////////
139 /// Reset matrix to zero, only the perspective scaling is set to w
140 /// (1 by default).
141 
142 void TEveTrans::ZeroTrans(Double_t w)
143 {
144  memset(fM, 0, 16*sizeof(Double_t));
145  fM[F33] = w;
146  fA1 = fA2 = fA3 = 0;
147  fAsOK = kFALSE;
148 }
149 
150 ////////////////////////////////////////////////////////////////////////////////
151 /// Reset rotation part of the matrix to unity.
152 
153 void TEveTrans::UnitRot()
154 {
155  memset(fM, 0, 12*sizeof(Double_t));
156  fM[F00] = fM[F11] = fM[F22] = 1;
157  fA1 = fA2 = fA3 = 0;
158  fAsOK = kTRUE;
159 }
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 /// Set matrix from another,
163 
164 void TEveTrans::SetTrans(const TEveTrans& t, Bool_t copyAngles)
165 {
166  memcpy(fM, t.fM, sizeof(fM));
167  if (copyAngles && t.fAsOK) {
168  fAsOK = kTRUE;
169  fA1 = t.fA1; fA2 = t.fA2; fA3 = t.fA3;
170  } else {
171  fAsOK = kFALSE;
172  }
173 }
174 
175 ////////////////////////////////////////////////////////////////////////////////
176 /// Set matrix from Double_t array.
177 
178 void TEveTrans::SetFromArray(const Double_t arr[16])
179 {
180  for(Int_t i=0; i<16; ++i) fM[i] = arr[i];
181  fAsOK = kFALSE;
182 }
183 
184 ////////////////////////////////////////////////////////////////////////////////
185 /// Set matrix from Float_t array.
186 
187 void TEveTrans::SetFromArray(const Float_t arr[16])
188 {
189  for(Int_t i=0; i<16; ++i) fM[i] = arr[i];
190  fAsOK = kFALSE;
191 }
192 
193 ////////////////////////////////////////////////////////////////////////////////
194 /// Setup the matrix as an elementary rotation.
195 /// Optimized versions of left/right multiplication with an elementary
196 /// rotation matrix are implemented in RotatePF/RotateLF.
197 /// Expects identity matrix.
198 
199 void TEveTrans::SetupRotation(Int_t i, Int_t j, Double_t f)
200 {
201  if(i == j) return;
202  TEveTrans& t = *this;
203  t(i,i) = t(j,j) = TMath::Cos(f);
204  Double_t s = TMath::Sin(f);
205  t(i,j) = -s; t(j,i) = s;
206  fAsOK = kFALSE;
207 }
208 
209 ////////////////////////////////////////////////////////////////////////////////
210 /// A function for creating a rotation matrix that rotates a vector called
211 /// "from" into another vector called "to".
212 /// Input : from[3], to[3] which both must be *normalized* non-zero vectors
213 /// Output: mtx[3][3] -- a 3x3 matrix in column-major form
214 ///
215 /// Authors: Tomas Möller, John Hughes
216 /// "Efficiently Building a Matrix to Rotate One Vector to Another"
217 /// Journal of Graphics Tools, 4(4):1-4, 1999
218 
219 void TEveTrans::SetupFromToVec(const TEveVector& from, const TEveVector& to)
220 {
221  static const float kFromToEpsilon = 0.000001f;
222 
223  ZeroTrans();
224 
225  Float_t e, f;
226  e = from.Dot(to);
227  f = (e < 0.0f) ? -e : e;
228 
229  if (f > 1.0f - kFromToEpsilon) /* "from" and "to"-vector almost parallel */
230  {
231  TEveVector u, v; /* temporary storage vectors */
232  TEveVector x; /* vector most nearly orthogonal to "from" */
233  Float_t c1, c2, c3; /* coefficients for later use */
234 
235  x.fX = (from.fX > 0.0f) ? from.fX : -from.fX;
236  x.fY = (from.fY > 0.0f) ? from.fY : -from.fY;
237  x.fZ = (from.fZ > 0.0f) ? from.fZ : -from.fZ;
238 
239  if (x.fX < x.fY)
240  {
241  if (x.fX < x.fZ) {
242  x.fX = 1.0f; x.fY = x.fZ = 0.0f;
243  } else {
244  x.fZ = 1.0f; x.fX = x.fY = 0.0f;
245  }
246  }
247  else
248  {
249  if (x.fY < x.fZ) {
250  x.fY = 1.0f; x.fX = x.fZ = 0.0f;
251  } else {
252  x.fZ = 1.0f; x.fX = x.fY = 0.0f;
253  }
254  }
255 
256  u.Sub(x, from);
257  v.Sub(x, to);
258 
259  c1 = 2.0f / u.Mag2();
260  c2 = 2.0f / v.Mag2();
261  c3 = c1 * c2 * u.Dot(v);
262 
263  for (int i = 0; i < 3; i++) {
264  for (int j = 0; j < 3; j++) {
265  CM(i, j) = - c1 * u[i] * u[j]
266  - c2 * v[i] * v[j]
267  + c3 * v[i] * u[j];
268  }
269  CM(i, i) += 1.0;
270  }
271  }
272  else /* the most common case, unless "from"="to", or "from"=-"to" */
273  {
274  TEveVector v = from.Cross(to);
275 
276  Float_t h, hvx, hvz, hvxy, hvxz, hvyz;
277  h = 1.0f/(1.0f + e);
278  hvx = h * v.fX;
279  hvz = h * v.fZ;
280  hvxy = hvx * v.fY;
281  hvxz = hvx * v.fZ;
282  hvyz = hvz * v.fY;
283 
284  CM(0, 0) = e + hvx * v.fX;
285  CM(0, 1) = hvxy - v.fZ;
286  CM(0, 2) = hvxz + v.fY;
287 
288  CM(1, 0) = hvxy + v.fZ;
289  CM(1, 1) = e + h * v.fY * v.fY;
290  CM(1, 2) = hvyz - v.fX;
291 
292  CM(2, 0) = hvxz - v.fY;
293  CM(2, 1) = hvyz + v.fX;
294  CM(2, 2) = e + hvz * v.fZ;
295  }
296 }
297 
298 ////////////////////////////////////////////////////////////////////////////////
299 /// Multiply from left: this = t * this.
300 
301 void TEveTrans::MultLeft(const TEveTrans& t)
302 {
303  Double_t buf[4];
304  Double_t* col = fM;
305  for(int c=0; c<4; ++c, col+=4) {
306  const Double_t* row = t.fM;
307  for(int r=0; r<4; ++r, ++row)
308  buf[r] = row[0]*col[0] + row[4]*col[1] + row[8]*col[2] + row[12]*col[3];
309  col[0] = buf[0]; col[1] = buf[1]; col[2] = buf[2]; col[3] = buf[3];
310  }
311  fAsOK = kFALSE;
312 }
313 
314 ////////////////////////////////////////////////////////////////////////////////
315 /// Multiply from right: this = this * t.
316 
317 void TEveTrans::MultRight(const TEveTrans& t)
318 {
319  Double_t buf[4];
320  Double_t* row = fM;
321  for(int r=0; r<4; ++r, ++row) {
322  const Double_t* col = t.fM;
323  for(int c=0; c<4; ++c, col+=4)
324  buf[c] = row[0]*col[0] + row[4]*col[1] + row[8]*col[2] + row[12]*col[3];
325  row[0] = buf[0]; row[4] = buf[1]; row[8] = buf[2]; row[12] = buf[3];
326  }
327  fAsOK = kFALSE;
328 }
329 
330 ////////////////////////////////////////////////////////////////////////////////
331 /// Copy, multiply from right and return product.
332 /// Avoid unless necessary.
333 
334 TEveTrans TEveTrans::operator*(const TEveTrans& t)
335 {
336  TEveTrans b(*this);
337  b.MultRight(t);
338  return b;
339 }
340 
341 ////////////////////////////////////////////////////////////////////////////////
342 /// Transpose 3x3 rotation sub-matrix.
343 
344 void TEveTrans::TransposeRotationPart()
345 {
346  Double_t x;
347  x = fM[F01]; fM[F01] = fM[F10]; fM[F10] = x;
348  x = fM[F02]; fM[F02] = fM[F20]; fM[F20] = x;
349  x = fM[F12]; fM[F12] = fM[F21]; fM[F21] = x;
350  fAsOK = kFALSE;
351 }
352 
353 ////////////////////////////////////////////////////////////////////////////////
354 /// Move in local-frame along axis with index ai.
355 
356 void TEveTrans::MoveLF(Int_t ai, Double_t amount)
357 {
358  const Double_t *col = fM + 4*--ai;
359  fM[F03] += amount*col[0]; fM[F13] += amount*col[1]; fM[F23] += amount*col[2];
360 }
361 
362 ////////////////////////////////////////////////////////////////////////////////
363 /// General move in local-frame.
364 
365 void TEveTrans::Move3LF(Double_t x, Double_t y, Double_t z)
366 {
367  fM[F03] += x*fM[0] + y*fM[4] + z*fM[8];
368  fM[F13] += x*fM[1] + y*fM[5] + z*fM[9];
369  fM[F23] += x*fM[2] + y*fM[6] + z*fM[10];
370 }
371 
372 ////////////////////////////////////////////////////////////////////////////////
373 /// Rotate in local frame. Does optimised version of MultRight.
374 
375 void TEveTrans::RotateLF(Int_t i1, Int_t i2, Double_t amount)
376 {
377  if(i1 == i2) return;
378  // Algorithm: TEveTrans a; a.SetupRotation(i1, i2, amount); MultRight(a);
379  // Optimized version:
380  const Double_t cos = TMath::Cos(amount), sin = TMath::Sin(amount);
381  Double_t b1, b2;
382  Double_t* row = fM;
383  --i1 <<= 2; --i2 <<= 2; // column major
384  for (int r=0; r<4; ++r, ++row) {
385  b1 = cos*row[i1] + sin*row[i2];
386  b2 = cos*row[i2] - sin*row[i1];
387  row[i1] = b1; row[i2] = b2;
388  }
389  fAsOK = kFALSE;
390 }
391 
392 ////////////////////////////////////////////////////////////////////////////////
393 /// Move in parent-frame along axis index ai.
394 
395 void TEveTrans::MovePF(Int_t ai, Double_t amount)
396 {
397  fM[F03 + --ai] += amount;
398 }
399 
400 ////////////////////////////////////////////////////////////////////////////////
401 /// General move in parent-frame.
402 
403 void TEveTrans::Move3PF(Double_t x, Double_t y, Double_t z)
404 {
405  fM[F03] += x;
406  fM[F13] += y;
407  fM[F23] += z;
408 }
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 /// Rotate in parent frame. Does optimised version of MultLeft.
412 
413 void TEveTrans::RotatePF(Int_t i1, Int_t i2, Double_t amount)
414 {
415  if(i1 == i2) return;
416  // Algorithm: TEveTrans a; a.SetupRotation(i1, i2, amount); MultLeft(a);
417 
418  // Optimized version:
419  const Double_t cos = TMath::Cos(amount), sin = TMath::Sin(amount);
420  Double_t b1, b2;
421  Double_t* col = fM;
422  --i1; --i2;
423  for(int c=0; c<4; ++c, col+=4) {
424  b1 = cos*col[i1] - sin*col[i2];
425  b2 = cos*col[i2] + sin*col[i1];
426  col[i1] = b1; col[i2] = b2;
427  }
428  fAsOK = kFALSE;
429 }
430 
431 ////////////////////////////////////////////////////////////////////////////////
432 /// Move in a's coord-system along axis-index ai.
433 
434 void TEveTrans::Move(const TEveTrans& a, Int_t ai, Double_t amount)
435 {
436  const Double_t* vec = a.fM + 4*--ai;
437  fM[F03] += amount*vec[0];
438  fM[F13] += amount*vec[1];
439  fM[F23] += amount*vec[2];
440 }
441 
442 ////////////////////////////////////////////////////////////////////////////////
443 /// General move in a's coord-system.
444 
445 void TEveTrans::Move3(const TEveTrans& a, Double_t x, Double_t y, Double_t z)
446 {
447  const Double_t* m = a.fM;
448  fM[F03] += x*m[F00] + y*m[F01] + z*m[F02];
449  fM[F13] += x*m[F10] + y*m[F11] + z*m[F12];
450  fM[F23] += x*m[F20] + y*m[F21] + z*m[F22];
451 }
452 
453 ////////////////////////////////////////////////////////////////////////////////
454 /// Rotate in a's coord-system, rotating base vector with index i1
455 /// into i2.
456 
457 void TEveTrans::Rotate(const TEveTrans& a, Int_t i1, Int_t i2, Double_t amount)
458 {
459  if(i1 == i2) return;
460  TEveTrans x(a);
461  x.Invert();
462  MultLeft(x);
463  RotatePF(i1, i2, amount);
464  MultLeft(a);
465  fAsOK = kFALSE;
466 }
467 
468 ////////////////////////////////////////////////////////////////////////////////
469 /// Set base-vector with index b.
470 
471 void TEveTrans::SetBaseVec(Int_t b, Double_t x, Double_t y, Double_t z)
472 {
473  Double_t* col = fM + 4*--b;
474  col[0] = x; col[1] = y; col[2] = z;
475  fAsOK = kFALSE;
476 }
477 
478 ////////////////////////////////////////////////////////////////////////////////
479 /// Set base-vector with index b.
480 
481 void TEveTrans::SetBaseVec(Int_t b, const TVector3& v)
482 {
483  Double_t* col = fM + 4*--b;
484  v.GetXYZ(col);
485  fAsOK = kFALSE;
486 }
487 
488 ////////////////////////////////////////////////////////////////////////////////
489 /// Get base-vector with index b.
490 
491 TVector3 TEveTrans::GetBaseVec(Int_t b) const
492 {
493  return TVector3(&fM[4*--b]);
494 }
495 
496 void TEveTrans::GetBaseVec(Int_t b, TVector3& v) const
497 {
498  // Get base-vector with index b.
499 
500  const Double_t* col = fM + 4*--b;
501  v.SetXYZ(col[0], col[1], col[2]);
502 }
503 
504 ////////////////////////////////////////////////////////////////////////////////
505 /// Set position (base-vec 4).
506 
507 void TEveTrans::SetPos(Double_t x, Double_t y, Double_t z)
508 {
509  fM[F03] = x; fM[F13] = y; fM[F23] = z;
510 }
511 
512 void TEveTrans::SetPos(Double_t* x)
513 {
514  // Set position (base-vec 4).
515  fM[F03] = x[0]; fM[F13] = x[1]; fM[F23] = x[2];
516 }
517 
518 void TEveTrans::SetPos(Float_t* x)
519 {
520  // Set position (base-vec 4).
521  fM[F03] = x[0]; fM[F13] = x[1]; fM[F23] = x[2];
522 }
523 
524 void TEveTrans::SetPos(const TEveTrans& t)
525 {
526  // Set position (base-vec 4).
527  const Double_t* m = t.fM;
528  fM[F03] = m[F03]; fM[F13] = m[F13]; fM[F23] = m[F23];
529 }
530 
531 ////////////////////////////////////////////////////////////////////////////////
532 /// Get position (base-vec 4).
533 
534 void TEveTrans::GetPos(Double_t& x, Double_t& y, Double_t& z) const
535 {
536  x = fM[F03]; y = fM[F13]; z = fM[F23];
537 }
538 
539 void TEveTrans::GetPos(Double_t* x) const
540 {
541  // Get position (base-vec 4).
542  x[0] = fM[F03]; x[1] = fM[F13]; x[2] = fM[F23];
543 }
544 
545 void TEveTrans::GetPos(Float_t* x) const
546 {
547  // Get position (base-vec 4).
548  x[0] = fM[F03]; x[1] = fM[F13]; x[2] = fM[F23];
549 }
550 
551 void TEveTrans::GetPos(TVector3& v) const
552 {
553  // Get position (base-vec 4).
554  v.SetXYZ(fM[F03], fM[F13], fM[F23]);
555 }
556 
557 TVector3 TEveTrans::GetPos() const
558 {
559  // Get position (base-vec 4).
560  return TVector3(fM[F03], fM[F13], fM[F23]);
561 }
562 
563 namespace
564 {
565 inline void clamp_angle(Float_t& a)
566 {
567  while(a < -TMath::TwoPi()) a += TMath::TwoPi();
568  while(a > TMath::TwoPi()) a -= TMath::TwoPi();
569 }
570 }
571 
572 void TEveTrans::SetRotByAngles(Float_t a1, Float_t a2, Float_t a3)
573 {
574  // Sets Rotation part as given by angles:
575  // a1 around z, -a2 around y, a3 around x.
576 
577  clamp_angle(a1); clamp_angle(a2); clamp_angle(a3);
578 
579  Double_t a, b, c, d, e, f;
580  a = TMath::Cos(a3); b = TMath::Sin(a3);
581  c = TMath::Cos(a2); d = TMath::Sin(a2); // should be -sin(a2) for positive direction
582  e = TMath::Cos(a1); f = TMath::Sin(a1);
583  Double_t ad = a*d, bd = b*d;
584 
585  fM[F00] = c*e; fM[F01] = -bd*e - a*f; fM[F02] = -ad*e + b*f;
586  fM[F10] = c*f; fM[F11] = -bd*f + a*e; fM[F12] = -ad*f - b*e;
587  fM[F20] = d; fM[F21] = b*c; fM[F22] = a*c;
588 
589  fA1 = a1; fA2 = a2; fA3 = a3;
590  fAsOK = kTRUE;
591 }
592 
593 ////////////////////////////////////////////////////////////////////////////////
594 /// Sets Rotation part as given by angles a1, a1, a3 and pattern pat.
595 /// Pattern consists of "XxYyZz" characters.
596 /// eg: x means rotate about x axis, X means rotate in negative direction
597 /// xYz -> R_x(a3) * R_y(-a2) * R_z(a1); (standard Gled representation)
598 /// Note that angles and pattern elements have inverted order!
599 ///
600 /// Implements Eulerian/Cardanian angles in a uniform way.
601 
602 void TEveTrans::SetRotByAnyAngles(Float_t a1, Float_t a2, Float_t a3,
603  const char* pat)
604 {
605  int n = strspn(pat, "XxYyZz"); if(n > 3) n = 3;
606  // Build Trans ... assign ...
607  Float_t a[] = { a3, a2, a1 };
608  UnitRot();
609  for(int i=0; i<n; i++) {
610  if(isupper(pat[i])) a[i] = -a[i];
611  switch(pat[i]) {
612  case 'x': case 'X': RotateLF(2, 3, a[i]); break;
613  case 'y': case 'Y': RotateLF(3, 1, a[i]); break;
614  case 'z': case 'Z': RotateLF(1, 2, a[i]); break;
615  }
616  }
617  fAsOK = kFALSE;
618 }
619 
620 ////////////////////////////////////////////////////////////////////////////////
621 /// Get Cardan rotation angles (pattern xYz above).
622 
623 void TEveTrans::GetRotAngles(Float_t* x) const
624 {
625  if(!fAsOK) {
626  Double_t sx, sy, sz;
627  GetScale(sx, sy, sz);
628  Double_t d = fM[F20]/sx;
629  if(d>1) d=1; else if(d<-1) d=-1; // Fix numerical errors
630  fA2 = TMath::ASin(d);
631  Double_t cos = TMath::Cos(fA2);
632  if(TMath::Abs(cos) > 8.7e-6) {
633  fA1 = TMath::ATan2(fM[F10], fM[F00]);
634  fA3 = TMath::ATan2(fM[F21]/sy, fM[F22]/sz);
635  } else {
636  fA1 = TMath::ATan2(fM[F10]/sx, fM[F11]/sy);
637  fA3 = 0;
638  }
639  fAsOK = kTRUE;
640  }
641  x[0] = fA1; x[1] = fA2; x[2] = fA3;
642 }
643 
644 ////////////////////////////////////////////////////////////////////////////////
645 /// Scale matrix. Translation part untouched.
646 
647 void TEveTrans::Scale(Double_t sx, Double_t sy, Double_t sz)
648 {
649  fM[F00] *= sx; fM[F10] *= sx; fM[F20] *= sx;
650  fM[F01] *= sy; fM[F11] *= sy; fM[F21] *= sy;
651  fM[F02] *= sz; fM[F12] *= sz; fM[F22] *= sz;
652 }
653 
654 ////////////////////////////////////////////////////////////////////////////////
655 /// Remove scaling, make all base vectors of unit length.
656 
657 Double_t TEveTrans::Unscale()
658 {
659  Double_t sx, sy, sz;
660  Unscale(sx, sy, sz);
661  return (sx + sy + sz)/3;
662 }
663 
664 ////////////////////////////////////////////////////////////////////////////////
665 /// Remove scaling, make all base vectors of unit length.
666 
667 void TEveTrans::Unscale(Double_t& sx, Double_t& sy, Double_t& sz)
668 {
669  GetScale(sx, sy, sz);
670  fM[F00] /= sx; fM[F10] /= sx; fM[F20] /= sx;
671  fM[F01] /= sy; fM[F11] /= sy; fM[F21] /= sy;
672  fM[F02] /= sz; fM[F12] /= sz; fM[F22] /= sz;
673 }
674 
675 ////////////////////////////////////////////////////////////////////////////////
676 /// Deduce scales from sizes of base vectors.
677 
678 void TEveTrans::GetScale(Double_t& sx, Double_t& sy, Double_t& sz) const
679 {
680  sx = TMath::Sqrt( fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20] );
681  sy = TMath::Sqrt( fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21] );
682  sz = TMath::Sqrt( fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22] );
683 }
684 
685 ////////////////////////////////////////////////////////////////////////////////
686 /// Set scaling.
687 
688 void TEveTrans::SetScale(Double_t sx, Double_t sy, Double_t sz)
689 {
690  sx /= TMath::Sqrt( fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20] );
691  sy /= TMath::Sqrt( fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21] );
692  sz /= TMath::Sqrt( fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22] );
693 
694  fM[F00] *= sx; fM[F10] *= sx; fM[F20] *= sx;
695  fM[F01] *= sy; fM[F11] *= sy; fM[F21] *= sy;
696  fM[F02] *= sz; fM[F12] *= sz; fM[F22] *= sz;
697 }
698 
699 ////////////////////////////////////////////////////////////////////////////////
700 /// Change x scaling.
701 
702 void TEveTrans::SetScaleX(Double_t sx)
703 {
704  sx /= TMath::Sqrt( fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20] );
705  fM[F00] *= sx; fM[F10] *= sx; fM[F20] *= sx;
706 }
707 
708 ////////////////////////////////////////////////////////////////////////////////
709 /// Change y scaling.
710 
711 void TEveTrans::SetScaleY(Double_t sy)
712 {
713  sy /= TMath::Sqrt( fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21] );
714  fM[F01] *= sy; fM[F11] *= sy; fM[F21] *= sy;
715 }
716 
717 ////////////////////////////////////////////////////////////////////////////////
718 /// Change z scaling.
719 
720 void TEveTrans::SetScaleZ(Double_t sz)
721 {
722  sz /= TMath::Sqrt( fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22] );
723  fM[F02] *= sz; fM[F12] *= sz; fM[F22] *= sz;
724 }
725 
726 ////////////////////////////////////////////////////////////////////////////////
727 /// Multiply vector in-place.
728 
729 void TEveTrans::MultiplyIP(TVector3& v, Double_t w) const
730 {
731  v.SetXYZ(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z() + fM[F03]*w,
732  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z() + fM[F13]*w,
733  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z() + fM[F23]*w);
734 }
735 
736 ////////////////////////////////////////////////////////////////////////////////
737 /// Multiply vector in-place.
738 
739 void TEveTrans::MultiplyIP(Double_t* v, Double_t w) const
740 {
741  Double_t r[3] = { v[0], v[1], v[2] };
742  v[0] = fM[F00]*r[0] + fM[F01]*r[1] + fM[F02]*r[2] + fM[F03]*w;
743  v[1] = fM[F10]*r[0] + fM[F11]*r[1] + fM[F12]*r[2] + fM[F13]*w;
744  v[2] = fM[F20]*r[0] + fM[F21]*r[1] + fM[F22]*r[2] + fM[F23]*w;
745 }
746 
747 ////////////////////////////////////////////////////////////////////////////////
748 /// Multiply vector in-place.
749 
750 void TEveTrans::MultiplyIP(Float_t* v, Double_t w) const
751 {
752  Double_t r[3] = { v[0], v[1], v[2] };
753  v[0] = fM[F00]*r[0] + fM[F01]*r[1] + fM[F02]*r[2] + fM[F03]*w;
754  v[1] = fM[F10]*r[0] + fM[F11]*r[1] + fM[F12]*r[2] + fM[F13]*w;
755  v[2] = fM[F20]*r[0] + fM[F21]*r[1] + fM[F22]*r[2] + fM[F23]*w;
756 }
757 
758 ////////////////////////////////////////////////////////////////////////////////
759 /// Multiply vector and return it.
760 
761 TVector3 TEveTrans::Multiply(const TVector3& v, Double_t w) const
762 {
763  return TVector3(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z() + fM[F03]*w,
764  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z() + fM[F13]*w,
765  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z() + fM[F23]*w);
766 }
767 
768 ////////////////////////////////////////////////////////////////////////////////
769 /// Multiply vector and fill output array vout.
770 
771 void TEveTrans::Multiply(const Double_t *vin, Double_t* vout, Double_t w) const
772 {
773  vout[0] = fM[F00]*vin[0] + fM[F01]*vin[1] + fM[F02]*vin[2] + fM[F03]*w;
774  vout[1] = fM[F10]*vin[0] + fM[F11]*vin[1] + fM[F12]*vin[1] + fM[F13]*w;
775  vout[2] = fM[F20]*vin[0] + fM[F21]*vin[1] + fM[F22]*vin[1] + fM[F23]*w;
776 }
777 
778 ////////////////////////////////////////////////////////////////////////////////
779 /// Rotate vector in-place. Translation is NOT applied.
780 
781 void TEveTrans::RotateIP(TVector3& v) const
782 {
783  v.SetXYZ(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z(),
784  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z(),
785  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z());
786 }
787 
788 ////////////////////////////////////////////////////////////////////////////////
789 /// Rotate vector in-place. Translation is NOT applied.
790 
791 void TEveTrans::RotateIP(Double_t* v) const
792 {
793  Double_t t[3] = { v[0], v[1], v[2] };
794 
795  v[0] = fM[F00]*t[0] + fM[F01]*t[1] + fM[F02]*t[2];
796  v[1] = fM[F10]*t[0] + fM[F11]*t[1] + fM[F12]*t[2];
797  v[2] = fM[F20]*t[0] + fM[F21]*t[1] + fM[F22]*t[2];
798 }
799 
800 ////////////////////////////////////////////////////////////////////////////////
801 /// Rotate vector in-place. Translation is NOT applied.
802 
803 void TEveTrans::RotateIP(Float_t* v) const
804 {
805  Double_t t[3] = { v[0], v[1], v[2] };
806 
807  v[0] = fM[F00]*t[0] + fM[F01]*t[1] + fM[F02]*t[2];
808  v[1] = fM[F10]*t[0] + fM[F11]*t[1] + fM[F12]*t[2];
809  v[2] = fM[F20]*t[0] + fM[F21]*t[1] + fM[F22]*t[2];
810 }
811 
812 ////////////////////////////////////////////////////////////////////////////////
813 /// Rotate vector and return the rotated vector. Translation is NOT applied.
814 
815 TVector3 TEveTrans::Rotate(const TVector3& v) const
816 {
817  return TVector3(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z(),
818  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z(),
819  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z());
820 }
821 
822 ////////////////////////////////////////////////////////////////////////////////
823 /// Norm 3-vector in column col.
824 
825 Double_t TEveTrans::Norm3Column(Int_t col)
826 {
827  Double_t* c = fM + 4*--col;
828  const Double_t l = TMath::Sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
829  c[0] /= l; c[1] /= l; c[2] /= l;
830  return l;
831 }
832 
833 ////////////////////////////////////////////////////////////////////////////////
834 /// Orto-norm 3-vector in column col with respect to column ref.
835 
836 Double_t TEveTrans::Orto3Column(Int_t col, Int_t ref)
837 {
838  Double_t* c = fM + 4*--col;
839  Double_t* rc = fM + 4*--ref;
840  const Double_t dp = c[0]*rc[0] + c[1]*rc[1] + c[2]*rc[2];
841  c[0] -= rc[0]*dp; c[1] -= rc[1]*dp; c[2] -= rc[2]*dp;
842  return dp;
843 }
844 
845 ////////////////////////////////////////////////////////////////////////////////
846 /// Orto-norm columns 1 to 3.
847 
848 void TEveTrans::OrtoNorm3()
849 {
850  Norm3Column(1);
851  Orto3Column(2,1); Norm3Column(2);
852  fM[F02] = fM[F10]*fM[F21] - fM[F11]*fM[F20];
853  fM[F12] = fM[F20]*fM[F01] - fM[F21]*fM[F00];
854  fM[F22] = fM[F00]*fM[F11] - fM[F01]*fM[F10];
855  // Cross-product faster than the following.
856  // Orto3Column(3,1); Orto3Column(3,2); Norm3Column(3);
857 }
858 
859 ////////////////////////////////////////////////////////////////////////////////
860 /// Invert matrix.
861 /// Copied from ROOT's TMatrixFCramerInv.
862 
863 Double_t TEveTrans::Invert()
864 {
865  static const TEveException eh("TEveTrans::Invert ");
866 
867  // Find all NECESSARY 2x2 dets: (18 of them)
868  const Double_t det2_12_01 = fM[F10]*fM[F21] - fM[F11]*fM[F20];
869  const Double_t det2_12_02 = fM[F10]*fM[F22] - fM[F12]*fM[F20];
870  const Double_t det2_12_03 = fM[F10]*fM[F23] - fM[F13]*fM[F20];
871  const Double_t det2_12_13 = fM[F11]*fM[F23] - fM[F13]*fM[F21];
872  const Double_t det2_12_23 = fM[F12]*fM[F23] - fM[F13]*fM[F22];
873  const Double_t det2_12_12 = fM[F11]*fM[F22] - fM[F12]*fM[F21];
874  const Double_t det2_13_01 = fM[F10]*fM[F31] - fM[F11]*fM[F30];
875  const Double_t det2_13_02 = fM[F10]*fM[F32] - fM[F12]*fM[F30];
876  const Double_t det2_13_03 = fM[F10]*fM[F33] - fM[F13]*fM[F30];
877  const Double_t det2_13_12 = fM[F11]*fM[F32] - fM[F12]*fM[F31];
878  const Double_t det2_13_13 = fM[F11]*fM[F33] - fM[F13]*fM[F31];
879  const Double_t det2_13_23 = fM[F12]*fM[F33] - fM[F13]*fM[F32];
880  const Double_t det2_23_01 = fM[F20]*fM[F31] - fM[F21]*fM[F30];
881  const Double_t det2_23_02 = fM[F20]*fM[F32] - fM[F22]*fM[F30];
882  const Double_t det2_23_03 = fM[F20]*fM[F33] - fM[F23]*fM[F30];
883  const Double_t det2_23_12 = fM[F21]*fM[F32] - fM[F22]*fM[F31];
884  const Double_t det2_23_13 = fM[F21]*fM[F33] - fM[F23]*fM[F31];
885  const Double_t det2_23_23 = fM[F22]*fM[F33] - fM[F23]*fM[F32];
886 
887  // Find all NECESSARY 3x3 dets: (16 of them)
888  const Double_t det3_012_012 = fM[F00]*det2_12_12 - fM[F01]*det2_12_02 + fM[F02]*det2_12_01;
889  const Double_t det3_012_013 = fM[F00]*det2_12_13 - fM[F01]*det2_12_03 + fM[F03]*det2_12_01;
890  const Double_t det3_012_023 = fM[F00]*det2_12_23 - fM[F02]*det2_12_03 + fM[F03]*det2_12_02;
891  const Double_t det3_012_123 = fM[F01]*det2_12_23 - fM[F02]*det2_12_13 + fM[F03]*det2_12_12;
892  const Double_t det3_013_012 = fM[F00]*det2_13_12 - fM[F01]*det2_13_02 + fM[F02]*det2_13_01;
893  const Double_t det3_013_013 = fM[F00]*det2_13_13 - fM[F01]*det2_13_03 + fM[F03]*det2_13_01;
894  const Double_t det3_013_023 = fM[F00]*det2_13_23 - fM[F02]*det2_13_03 + fM[F03]*det2_13_02;
895  const Double_t det3_013_123 = fM[F01]*det2_13_23 - fM[F02]*det2_13_13 + fM[F03]*det2_13_12;
896  const Double_t det3_023_012 = fM[F00]*det2_23_12 - fM[F01]*det2_23_02 + fM[F02]*det2_23_01;
897  const Double_t det3_023_013 = fM[F00]*det2_23_13 - fM[F01]*det2_23_03 + fM[F03]*det2_23_01;
898  const Double_t det3_023_023 = fM[F00]*det2_23_23 - fM[F02]*det2_23_03 + fM[F03]*det2_23_02;
899  const Double_t det3_023_123 = fM[F01]*det2_23_23 - fM[F02]*det2_23_13 + fM[F03]*det2_23_12;
900  const Double_t det3_123_012 = fM[F10]*det2_23_12 - fM[F11]*det2_23_02 + fM[F12]*det2_23_01;
901  const Double_t det3_123_013 = fM[F10]*det2_23_13 - fM[F11]*det2_23_03 + fM[F13]*det2_23_01;
902  const Double_t det3_123_023 = fM[F10]*det2_23_23 - fM[F12]*det2_23_03 + fM[F13]*det2_23_02;
903  const Double_t det3_123_123 = fM[F11]*det2_23_23 - fM[F12]*det2_23_13 + fM[F13]*det2_23_12;
904 
905  // Find the 4x4 det:
906  const Double_t det = fM[F00]*det3_123_123 - fM[F01]*det3_123_023 +
907  fM[F02]*det3_123_013 - fM[F03]*det3_123_012;
908 
909  if(det == 0) {
910  throw(eh + "matrix is singular.");
911  }
912 
913  const Double_t oneOverDet = 1.0/det;
914  const Double_t mn1OverDet = - oneOverDet;
915 
916  fM[F00] = det3_123_123 * oneOverDet;
917  fM[F01] = det3_023_123 * mn1OverDet;
918  fM[F02] = det3_013_123 * oneOverDet;
919  fM[F03] = det3_012_123 * mn1OverDet;
920 
921  fM[F10] = det3_123_023 * mn1OverDet;
922  fM[F11] = det3_023_023 * oneOverDet;
923  fM[F12] = det3_013_023 * mn1OverDet;
924  fM[F13] = det3_012_023 * oneOverDet;
925 
926  fM[F20] = det3_123_013 * oneOverDet;
927  fM[F21] = det3_023_013 * mn1OverDet;
928  fM[F22] = det3_013_013 * oneOverDet;
929  fM[F23] = det3_012_013 * mn1OverDet;
930 
931  fM[F30] = det3_123_012 * mn1OverDet;
932  fM[F31] = det3_023_012 * oneOverDet;
933  fM[F32] = det3_013_012 * mn1OverDet;
934  fM[F33] = det3_012_012 * oneOverDet;
935 
936  fAsOK = kFALSE;
937  return det;
938 }
939 
940 ////////////////////////////////////////////////////////////////////////////////
941 /// Stream an object of class TEveTrans.
942 
943 void TEveTrans::Streamer(TBuffer &R__b)
944 {
945  if (R__b.IsReading()) {
946  TEveTrans::Class()->ReadBuffer(R__b, this);
947  fAsOK = kFALSE;
948  } else {
949  TEveTrans::Class()->WriteBuffer(R__b, this);
950  }
951 }
952 
953 ////////////////////////////////////////////////////////////////////////////////
954 /// Print in reasonable format.
955 
956 void TEveTrans::Print(Option_t* /*option*/) const
957 {
958  const Double_t* row = fM;
959  for(Int_t i=0; i<4; ++i, ++row)
960  printf("%8.3f %8.3f %8.3f | %8.3f\n", row[0], row[4], row[8], row[12]);
961 }
962 
963 #include <iomanip>
964 
965 ////////////////////////////////////////////////////////////////////////////////
966 /// Print to std::ostream.
967 
968 std::ostream& operator<<(std::ostream& s, const TEveTrans& t)
969 {
970  s.setf(std::ios::fixed, std::ios::floatfield);
971  s.precision(3);
972  for(Int_t i=1; i<=4; i++)
973  for(Int_t j=1; j<=4; j++)
974  s << t(i,j) << ((j==4) ? "\n" : "\t");
975  return s;
976 }
977 
978 #include "TGeoMatrix.h"
979 #include "TBuffer3D.h"
980 
981 void TEveTrans::SetFrom(Double_t* carr)
982 {
983  // Initialize from array.
984 
985  fUseTrans = kTRUE;
986  memcpy(fM, carr, 16*sizeof(Double_t));
987  fAsOK = kFALSE;
988 }
989 
990 ////////////////////////////////////////////////////////////////////////////////
991 /// Initialize from TGeoMatrix.
992 
993 void TEveTrans::SetFrom(const TGeoMatrix& mat)
994 {
995  fUseTrans = kTRUE;
996  const Double_t *r = mat.GetRotationMatrix();
997  const Double_t *t = mat.GetTranslation();
998  Double_t *m = fM;
999  if (mat.IsScale())
1000  {
1001  const Double_t *s = mat.GetScale();
1002  m[0] = r[0]*s[0]; m[1] = r[3]*s[0]; m[2] = r[6]*s[0]; m[3] = 0;
1003  m[4] = r[1]*s[1]; m[5] = r[4]*s[1]; m[6] = r[7]*s[1]; m[7] = 0;
1004  m[8] = r[2]*s[2]; m[9] = r[5]*s[2]; m[10] = r[8]*s[2]; m[11] = 0;
1005  m[12] = t[0]; m[13] = t[1]; m[14] = t[2]; m[15] = 1;
1006  }
1007  else
1008  {
1009  m[0] = r[0]; m[1] = r[3]; m[2] = r[6]; m[3] = 0;
1010  m[4] = r[1]; m[5] = r[4]; m[6] = r[7]; m[7] = 0;
1011  m[8] = r[2]; m[9] = r[5]; m[10] = r[8]; m[11] = 0;
1012  m[12] = t[0]; m[13] = t[1]; m[14] = t[2]; m[15] = 1;
1013  }
1014  fAsOK = kFALSE;
1015 }
1016 
1017 ////////////////////////////////////////////////////////////////////////////////
1018 /// Set TGeoHMatrix mat.
1019 
1020 void TEveTrans::SetGeoHMatrix(TGeoHMatrix& mat)
1021 {
1022  Double_t *r = mat.GetRotationMatrix();
1023  Double_t *t = mat.GetTranslation();
1024  Double_t *s = mat.GetScale();
1025  if (fUseTrans)
1026  {
1027  mat.SetBit(TGeoMatrix::kGeoGenTrans);
1028  Double_t *m = fM;
1029  GetScale(s[0], s[1], s[2]);
1030  r[0] = m[0]/s[0]; r[3] = m[1]/s[0]; r[6] = m[2]/s[0]; m += 4;
1031  r[1] = m[0]/s[1]; r[4] = m[1]/s[1]; r[7] = m[2]/s[1]; m += 4;
1032  r[2] = m[0]/s[2]; r[5] = m[1]/s[2]; r[8] = m[2]/s[2]; m += 4;
1033  t[0] = m[0]; t[1] = m[1]; t[2] = m[2];
1034  }
1035  else
1036  {
1037  mat.ResetBit(TGeoMatrix::kGeoGenTrans);
1038  r[0] = 1; r[3] = 0; r[6] = 0;
1039  r[1] = 0; r[4] = 1; r[7] = 0;
1040  r[2] = 0; r[5] = 0; r[8] = 1;
1041  s[0] = s[1] = s[2] = 1;
1042  t[0] = t[1] = t[2] = 0;
1043  }
1044 }
1045 
1046 ////////////////////////////////////////////////////////////////////////////////
1047 /// Fill transformation part TBuffer3D core section.
1048 
1049 void TEveTrans::SetBuffer3D(TBuffer3D& buff)
1050 {
1051  buff.fLocalFrame = fUseTrans;
1052  if (fUseTrans) {
1053  // In phys-shape ctor the rotation part is transposed, due to
1054  // TGeo's convention for rotation matrix. So we have to transpose
1055  // it here, also.
1056  Double_t *m = buff.fLocalMaster;
1057  m[0] = fM[0]; m[1] = fM[4]; m[2] = fM[8]; m[3] = fM[3];
1058  m[4] = fM[1]; m[5] = fM[5]; m[6] = fM[9]; m[7] = fM[7];
1059  m[8] = fM[2]; m[9] = fM[6]; m[10] = fM[10]; m[11] = fM[11];
1060  m[12] = fM[12]; m[13] = fM[13]; m[14] = fM[14]; m[15] = fM[15];
1061  // Otherwise this would do:
1062  // memcpy(buff.fLocalMaster, fM, 16*sizeof(Double_t));
1063  }
1064 }
1065 
1066 ////////////////////////////////////////////////////////////////////////////////
1067 /// Test if the transformation is a scale.
1068 /// To be used by ROOT TGLObject descendants that potentially need to
1069 /// use GL_NORMALIZE.
1070 /// The low/high limits are expected to be squares of actual limits.
1071 ///
1072 /// Ideally this should be done by the TGLViewer [but is not].
1073 
1074 Bool_t TEveTrans::IsScale(Double_t low, Double_t high) const
1075 {
1076  if (!fUseTrans) return kFALSE;
1077  Double_t s;
1078  s = fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20];
1079  if (s < low || s > high) return kTRUE;
1080  s = fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21];
1081  if (s < low || s > high) return kTRUE;
1082  s = fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22];
1083  if (s < low || s > high) return kTRUE;
1084  return kFALSE;
1085 }