From 2e275ad1c6bc24012f8bc0ff593a797a13154b99 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (laptop)" <j.wuttke@fz-juelich.de>
Date: Mon, 25 Mar 2013 21:41:35 +0100
Subject: [PATCH] broken intermediate state

---
 Core/Geometry/inc/Transform3D.h   | 129 +++++------------------
 Core/Geometry/src/Transform3D.cpp | 166 ++++++++++++++----------------
 2 files changed, 103 insertions(+), 192 deletions(-)

diff --git a/Core/Geometry/inc/Transform3D.h b/Core/Geometry/inc/Transform3D.h
index afcd75f5daf..6e406f30e32 100644
--- a/Core/Geometry/inc/Transform3D.h
+++ b/Core/Geometry/inc/Transform3D.h
@@ -43,8 +43,7 @@ namespace Geometry {
 
 //! Transformations of 3D geometrical objects (points, vectors, normals).
 
-//! Uses a 4x3 double-precision transform matrix to rotate, translate,
-//! reflect and scale.
+//! Uses a 3x3 double-precision transform matrix to rotate vectors.
 //!
 //! Identity transformation:
 //!   Transform3D::Identity   - global identity transformation;
@@ -69,30 +68,6 @@ namespace Geometry {
 //!   RotateY3D(ang)           - rotation around Y-axis;
 //!   RotateZ3D(ang)           - rotation around Z-axis;
 //!
-//! Translations:
-//!   Translate3D(v)           - translation given by CLHEP::Hep3Vector "v";
-//!   Translate3D(dx,dy,dz)    - translation on vector (dx,dy,dz);
-//!   TraslateX3D(dx)          - translation along X-axis;
-//!   TraslateY3D(dy)          - translation along Y-axis;
-//!   TraslateZ3D(dz)          - translation along Z-axis;
-//!
-//! Reflections:
-//!   Reflect3D(a,b,c,d)       - reflection in the plane a*x+b*y+c*z+d=0;
-//!   Reflect3D(normal,p)      - reflection in the plane going through "p"
-//!                              and whose normal is equal to "normal";
-//!   ReflectX3D(a)            - reflect X in the plane x=a (default a=0);
-//!   ReflectY3D(a)            - reflect Y in the plane y=a (default a=0);
-//!   ReflectZ3D(a)            - reflect Z in the plane z=a (default a=0);
-//!
-//! Scalings:
-//!   Scale3D(sx,sy,sz)        - general scaling with factors "sx","sy","sz"
-//!                                 along X, Y and Z;
-//!   Scale3D(s)               - scaling with constant factor "s" along all 
-//!                                 directions;
-//!   ScaleX3D(sx)             - scale X;
-//!   ScaleY3D(sy)             - scale Y;
-//!   ScaleZ3D(sz)             - scale Z;
-//!
 //! Inverse transformation:
 //!   m.inverse() or           - returns inverse transformation;
 //!
@@ -100,44 +75,10 @@ namespace Geometry {
 //!   m3 = m2 * m1             - it is relatively slow in comparison with
 //!                              transformation of a vector. Use parenthesis
 //!                              to avoid this operation (see example below);
-//! Transformation of point:
-//!   p2 = m * p1
-//!
-//! Transformation of vector:
-//!   v2 = m * v1
-//!
-//! Transformation of normal:
-//!   n2 = m * n1
-//!
-//! The following table explains how different transformations affect
-//! point, vector and normal. "+" means affect, "-" means do not affect,
-//! "*" meas affect but in different way than "+" 
-//!
-//!                     Point  Vector  Normal
-//!      -------------+-------+-------+-------
-//!       Rotation    !   +   !   +   !   +
-//!       Translation !   +   !   -   !   -
-//!       Reflection  !   +   !   +   !   *
-//!       Scaling     !   +   !   +   !   *
-//!      -------------+-------+-------+-------
-//!
-//! Example of the usage:
-//!
-//!   Transform3D   m1, m2, m3;
-//!   BasicVector3D v2, v1(0,0,0);
-//!
-//!   m1 = Rotate3D(angle, Vector3D(1,1,1));
-//!   m2 = Translate3D(dx,dy,dz);
-//!   m3 = m1.inverse();
-//!
-//!   v2 = m3*(m2*(m1*v1));
 //!
 //! Several specialized classes are derived from it:
 //!
-//! TranslateX3D, TranslateY3D, TranslateZ3D, Translate3D,<br>
-//! RotateX3D,    RotateY3D,    RotateZ3D,    Rotate3D,   <br>
-//! ScaleX3D,     ScaleY3D,     ScaleZ3D,     Scale3D,    <br>
-//! ReflectX3D,   ReflectY3D,   ReflectZ3D,   Reflect3D.
+//! RotateX3D,    RotateY3D,    RotateZ3D,    Rotate3D.
 //!
 //! The idea behind these classes is to provide some additional constructors
 //! for Transform3D, they normally should not be used as separate classes.
@@ -153,25 +94,25 @@ namespace Geometry {
 class Transform3D {
   protected:
     // 4x3  Transformation Matrix
-    double xx_, xy_, xz_, dx_,     
-           yx_, yy_, yz_, dy_,
-           zx_, zy_, zz_, dz_;
+    double xx_, xy_, xz_,     
+           yx_, yy_, yz_,
+           zx_, zy_, zz_;
 
     //! Protected constructor.
-    Transform3D(double XX, double XY, double XZ, double DX,
-                double YX, double YY, double YZ, double DY,
-                double ZX, double ZY, double ZZ, double DZ)
-        : xx_(XX), xy_(XY), xz_(XZ), dx_(DX),
-          yx_(YX), yy_(YY), yz_(YZ), dy_(DY),
-          zx_(ZX), zy_(ZY), zz_(ZZ), dz_(DZ) {}
+    Transform3D(double XX, double XY, double XZ,
+                double YX, double YY, double YZ,
+                double ZX, double ZY, double ZZ)
+        : xx_(XX), xy_(XY), xz_(XZ),
+          yx_(YX), yy_(YY), yz_(YZ),
+          zx_(ZX), zy_(ZY), zz_(ZZ) {}
 
     //! Sets transformation matrix.
-    void setTransform(double XX, double XY, double XZ, double DX,
-                      double YX, double YY, double YZ, double DY,
-                      double ZX, double ZY, double ZZ, double DZ) {
-        xx_ = XX; xy_ = XY; xz_ = XZ; dx_ = DX;
-        yx_ = YX; yy_ = YY; yz_ = YZ; dy_ = DY;
-        zx_ = ZX; zy_ = ZY; zz_ = ZZ; dz_ = DZ;
+    void setTransform(double XX, double XY, double XZ,
+                      double YX, double YY, double YZ,
+                      double ZX, double ZY, double ZZ) {
+        xx_ = XX; xy_ = XY; xz_ = XZ;
+        yx_ = YX; yy_ = YY; yz_ = YZ;
+        zx_ = ZX; zy_ = ZY; zz_ = ZZ;
     }
 
   public:
@@ -190,28 +131,20 @@ class Transform3D {
 
     //! Default constructor - sets the Identity transformation.
     Transform3D()
-        : xx_(1), xy_(0), xz_(0), dx_(0),
-          yx_(0), yy_(1), yz_(0), dy_(0),
-          zx_(0), zy_(0), zz_(1), dz_(0) {}
+        : xx_(1), xy_(0), xz_(0),
+          yx_(0), yy_(1), yz_(0),
+          zx_(0), zy_(0), zz_(1) {}
   
-    //! Constructor: transformation of basis (assumed - no reflection).
-    Transform3D(const Point3D<double>& fr0,
-                const Point3D<double>& fr1,
-                const Point3D<double>& fr2,
-                const Point3D<double>& to0,
-                const Point3D<double>& to1,
-                const Point3D<double>& to2);
-
     //! Copy constructor.
     Transform3D(const Transform3D& m)
-        : xx_(m.xx_), xy_(m.xy_), xz_(m.xz_), dx_(m.dx_),
-          yx_(m.yx_), yy_(m.yy_), yz_(m.yz_), dy_(m.dy_),
-          zx_(m.zx_), zy_(m.zy_), zz_(m.zz_), dz_(m.dz_) {}
+        : xx_(m.xx_), xy_(m.xy_), xz_(m.xz_),
+          yx_(m.yx_), yy_(m.yy_), yz_(m.yz_),
+          zx_(m.zx_), zy_(m.zy_), zz_(m.zz_) {}
 
     //! Destructor. 
     //! Virtual for now as some persistency mechanism needs that,
     //! in future releases this might go away again.
-    virtual ~Transform3D() { /* nop */ }
+    virtual ~Transform3D() {}
 
     //! Returns object of the helper class for C-style subscripting r[i][j]
     inline const Transform3D_row operator [] (int) const; 
@@ -237,24 +170,18 @@ class Transform3D {
     double zy() const { return zy_; }
     //! Gets zz-element of the transformation matrix.
     double zz() const { return zz_; }
-    //! Gets dx-element of the transformation matrix.
-    double dx() const { return dx_; }
-    //! Gets dy-element of the transformation matrix.
-    double dy() const { return dy_; }
-    //! Gets dz-element of the transformation matrix.
-    double dz() const { return dz_; }
     
     //! Assignment.
     Transform3D& operator=(const Transform3D&m) {
-        xx_= m.xx_; xy_= m.xy_; xz_= m.xz_; dx_= m.dx_;
-        yx_= m.yx_; yy_= m.yy_; yz_= m.yz_; dy_= m.dy_;
-        zx_= m.zx_; zy_= m.zy_; zz_= m.zz_; dz_= m.dz_;
+        xx_= m.xx_; xy_= m.xy_; xz_= m.xz_;
+        yx_= m.yx_; yy_= m.yy_; yz_= m.yz_;
+        zx_= m.zx_; zy_= m.zy_; zz_= m.zz_;
         return *this;
     }
 
     //! Sets the Identity transformation.
     void setIdentity() { 
-        xy_= xz_= dx_= yx_= yz_= dy_= zx_= zy_= dz_= 0; xx_= yy_= zz_= 1;
+        xy_= xz_= yx_= yz_= zx_= zy_= 0; xx_= yy_= zz_= 1;
     }
     
     //! Returns the inverse transformation.
diff --git a/Core/Geometry/src/Transform3D.cpp b/Core/Geometry/src/Transform3D.cpp
index cd1d825facc..b0d6c02b809 100644
--- a/Core/Geometry/src/Transform3D.cpp
+++ b/Core/Geometry/src/Transform3D.cpp
@@ -37,22 +37,18 @@ double Transform3D::operator () (int i, int j) const
         if (j == 0) { return xx_; }
         if (j == 1) { return xy_; }
         if (j == 2) { return xz_; }
-        if (j == 3) { return dx_; }
     } else if (i == 1) {
         if (j == 0) { return yx_; }
         if (j == 1) { return yy_; }
         if (j == 2) { return yz_; }
-        if (j == 3) { return dy_; }
     } else if (i == 2) {
         if (j == 0) { return zx_; }
         if (j == 1) { return zy_; }
         if (j == 2) { return zz_; }
-        if (j == 3) { return dz_; }
     } else if (i == 3) {
         if (j == 0) { return 0.0; }
         if (j == 1) { return 0.0; }
         if (j == 2) { return 0.0; }
-        if (j == 3) { return 1.0; }
     }
     std::cerr << "Transform3D subscripting: bad indeces "
               << "(" << i << "," << j << ")" << std::endl;
@@ -64,90 +60,15 @@ double Transform3D::operator () (int i, int j) const
 Transform3D Transform3D::operator*(const Transform3D & b) const
 {
     return Transform3D
-            (xx_*b.xx_+xy_*b.yx_+xz_*b.zx_, xx_*b.xy_+xy_*b.yy_+xz_*b.zy_,
-             xx_*b.xz_+xy_*b.yz_+xz_*b.zz_, xx_*b.dx_+xy_*b.dy_+xz_*b.dz_+dx_,
-             yx_*b.xx_+yy_*b.yx_+yz_*b.zx_, yx_*b.xy_+yy_*b.yy_+yz_*b.zy_,
-             yx_*b.xz_+yy_*b.yz_+yz_*b.zz_, yx_*b.dx_+yy_*b.dy_+yz_*b.dz_+dy_,
-             zx_*b.xx_+zy_*b.yx_+zz_*b.zx_, zx_*b.xy_+zy_*b.yy_+zz_*b.zy_,
-             zx_*b.xz_+zy_*b.yz_+zz_*b.zz_, zx_*b.dx_+zy_*b.dy_+zz_*b.dz_+dz_);
-}
-
-//! Constructor: transformation of basis (assumed - no reflection).
-
-//! Transformation from coordinate system defined by its origin "fr0"
-//! and two axes "fr0->fr1", "fr0->fr2" to another coordinate system
-//! "to0", "to0->to1", "to0->to2"
-//!
-//! @author E. Chernyaev 1996-2003
-
-Transform3D::Transform3D(const Point3D<double> & fr0,
-                         const Point3D<double> & fr1,
-                         const Point3D<double> & fr2,
-                         const Point3D<double> & to0,
-                         const Point3D<double> & to1,
-                         const Point3D<double> & to2)
-    : xx_(1), xy_(0), xz_(0), dx_(0),
-      yx_(0), yy_(1), yz_(0), dy_(0),
-      zx_(0), zy_(0), zz_(1), dz_(0)
-{
-    Vector3D<double> x1,y1,z1, x2,y2,z2;
-    x1 = (fr1 - fr0).unit();
-    y1 = (fr2 - fr0).unit();
-    x2 = (to1 - to0).unit();
-    y2 = (to2 - to0).unit();
-    
-    // -- Check angles --
-    
-    double cos1, cos2;
-    cos1 = x1.dot(y1);
-    cos2 = x2.dot(y2);
-    
-    if (std::abs(1.0-cos1) <= 0.000001 || std::abs(1.0-cos2) <= 0.000001) {
-        std::cerr << "Transform3D: zero angle between axes" << std::endl;
-        setIdentity();
-    } else {
-        if (std::abs(cos1-cos2) > 0.000001) {
-            std::cerr << "Transform3D: angles between axes are not equal"
-                      << std::endl;
-        }
-
-        // -- Find rotation matrix --
-
-        z1 = (x1.cross(y1)).unit();
-        y1  = z1.cross(x1);
-
-        z2 = (x2.cross(y2)).unit();
-        y2  = z2.cross(x2);
-
-        double detxx =  (y1.y()*z1.z() - z1.y()*y1.z());
-        double detxy = -(y1.x()*z1.z() - z1.x()*y1.z());
-        double detxz =  (y1.x()*z1.y() - z1.x()*y1.y());
-        double detyx = -(x1.y()*z1.z() - z1.y()*x1.z());
-        double detyy =  (x1.x()*z1.z() - z1.x()*x1.z());
-        double detyz = -(x1.x()*z1.y() - z1.x()*x1.y());
-        double detzx =  (x1.y()*y1.z() - y1.y()*x1.z());
-        double detzy = -(x1.x()*y1.z() - y1.x()*x1.z());
-        double detzz =  (x1.x()*y1.y() - y1.x()*x1.y());
-
-        double txx = x2.x()*detxx + y2.x()*detyx + z2.x()*detzx;
-        double txy = x2.x()*detxy + y2.x()*detyy + z2.x()*detzy;
-        double txz = x2.x()*detxz + y2.x()*detyz + z2.x()*detzz;
-        double tyx = x2.y()*detxx + y2.y()*detyx + z2.y()*detzx;
-        double tyy = x2.y()*detxy + y2.y()*detyy + z2.y()*detzy;
-        double tyz = x2.y()*detxz + y2.y()*detyz + z2.y()*detzz;
-        double tzx = x2.z()*detxx + y2.z()*detyx + z2.z()*detzx;
-        double tzy = x2.z()*detxy + y2.z()*detyy + z2.z()*detzy;
-        double tzz = x2.z()*detxz + y2.z()*detyz + z2.z()*detzz;
-
-        // -- set transformation --
-
-        double dx1 = fr0.x(), dy1 = fr0.y(), dz1 = fr0.z();
-        double dx2 = to0.x(), dy2 = to0.y(), dz2 = to0.z();
-
-        setTransform(txx, txy, txz, dx2-txx*dx1-txy*dy1-txz*dz1,
-                     tyx, tyy, tyz, dy2-tyx*dx1-tyy*dy1-tyz*dz1,
-                     tzx, tzy, tzz, dz2-tzx*dx1-tzy*dy1-tzz*dz1);
-    }
+        (xx_*b.xx_+xy_*b.yx_+xz_*b.zx_,
+         xx_*b.xy_+xy_*b.yy_+xz_*b.zy_,
+         xx_*b.xz_+xy_*b.yz_+xz_*b.zz_, 
+         yx_*b.xx_+yy_*b.yx_+yz_*b.zx_,
+         yx_*b.xy_+yy_*b.yy_+yz_*b.zy_,
+         yx_*b.xz_+yy_*b.yz_+yz_*b.zz_, 
+         zx_*b.xx_+zy_*b.yx_+zz_*b.zx_,
+         zx_*b.xy_+zy_*b.yy_+zz_*b.zy_,
+         zx_*b.xz_+zy_*b.yz_+zz_*b.zz_ );
 }
 
 //! Returns the inverse transformation.
@@ -172,9 +93,72 @@ Transform3D Transform3D::inverse() const
     double detzy = (xx_*yz_ - xz_*yx_)*det;
     double detzz = (xx_*yy_ - xy_*yx_)*det;
     return Transform3D
-            (detxx, -detyx,  detzx, -detxx*dx_+detyx*dy_-detzx*dz_,
-             -detxy,  detyy, -detzy,  detxy*dx_-detyy*dy_+detzy*dz_,
-             detxz, -detyz,  detzz, -detxz*dx_+detyz*dy_-detzz*dz_);
+            (detxx, -detyx,  detzx,
+             -detxy,  detyy, -detzy,
+             detxz, -detyz,  detzz );
+}
+
+//! Difference between corresponding matrix elements less than tolerance?
+
+bool Transform3D::isNear(const Transform3D & t, double tolerance) const
+{
+    return ( (std::abs(xx_ - t.xx_) <= tolerance) &&
+             (std::abs(xy_ - t.xy_) <= tolerance) &&
+             (std::abs(xz_ - t.xz_) <= tolerance) &&
+             (std::abs(yx_ - t.yx_) <= tolerance) &&
+             (std::abs(yy_ - t.yy_) <= tolerance) &&
+             (std::abs(yz_ - t.yz_) <= tolerance) &&
+             (std::abs(zx_ - t.zx_) <= tolerance) &&
+             (std::abs(zy_ - t.zy_) <= tolerance) &&
+             (std::abs(zz_ - t.zz_) <= tolerance) );
+}
+
+//! Test for equality.
+
+bool Transform3D::operator==(const Transform3D & t) const
+{
+    return (this == &t) ? true :
+                          (xx_==t.xx_ && xy_==t.xy_ && xz_==t.xz_ &&
+                           yx_==t.yx_ && yy_==t.yy_ && yz_==t.yz_ &&
+                           zx_==t.zx_ && zy_==t.zy_ && zz_==t.zz_ );
+}
+
+
+//! Construct rotation by angle a around axis p1->p2.
+//!
+//! @author E. Chernyaev 1996
+
+Rotate3D::Rotate3D(double a,
+                   const Point3D<double> & p1,
+                   const Point3D<double> & p2)
+        : Transform3D()
+{
+    if (a == 0) return;
+
+    double cx = p2.x()-p1.x(), cy = p2.y()-p1.y(), cz = p2.z()-p1.z();
+    double ll = std::sqrt(cx*cx + cy*cy + cz*cz);
+    if (ll == 0) {
+        std::cerr << "Rotate3D: zero axis" << std::endl;
+    } else {
+        double cosa = std::cos(a), sina = std::sin(a);
+        cx /= ll; cy /= ll; cz /= ll;
+
+        double txx = cosa + (1-cosa)*cx*cx;
+        double txy =        (1-cosa)*cx*cy - sina*cz;
+        double txz =        (1-cosa)*cx*cz + sina*cy;
+
+        double tyx =        (1-cosa)*cy*cx + sina*cz;
+        double tyy = cosa + (1-cosa)*cy*cy;
+        double tyz =        (1-cosa)*cy*cz - sina*cx;
+
+        double tzx =        (1-cosa)*cz*cx - sina*cy;
+        double tzy =        (1-cosa)*cz*cy + sina*cx;
+        double tzz = cosa + (1-cosa)*cz*cz;
+
+        setTransform(txx, txy, txz,
+                     tyx, tyy, tyz,
+                     tzx, tzy, tzz);
+    }
 }
 
 }  // namespace Geometry
-- 
GitLab