From 8c79e07343f4a9927ef4b6d670dbfe1059f99d2d Mon Sep 17 00:00:00 2001
From: eap <eap@opencascade.com>
Date: Wed, 17 Dec 2008 15:57:35 +0000
Subject: [PATCH] MEDMEM Industrialization 2008 +    void reverseApply(double*
 destPt, const double* srcPt) const;

---
 src/INTERP_KERNEL/TetraAffineTransform.cxx | 45 ++++++++++++++++++++++
 src/INTERP_KERNEL/TetraAffineTransform.hxx |  4 ++
 2 files changed, 49 insertions(+)

diff --git a/src/INTERP_KERNEL/TetraAffineTransform.cxx b/src/INTERP_KERNEL/TetraAffineTransform.cxx
index e018176a3..6a57bf49a 100644
--- a/src/INTERP_KERNEL/TetraAffineTransform.cxx
+++ b/src/INTERP_KERNEL/TetraAffineTransform.cxx
@@ -35,6 +35,9 @@ namespace INTERP_KERNEL
             _linearTransform[3*j + i] = (pts[i+1])[j] - (pts[0])[j];
           }
       }
+
+    // remember _linearTransform for the reverse transformation
+    memcpy( _backLinearTransform, _linearTransform, 9*sizeof(double));
     
     calculateDeterminant();
     
@@ -127,6 +130,48 @@ namespace INTERP_KERNEL
       }
   }
 
+  /**
+   * Calculates the reverse transform of point srcPt and stores the result in destPt.
+   * If destPt == srcPt, then srcPt is overwritten safely.
+   *
+   * @param destPt  double[3] in which to store the transformed point
+   * @param srcPt   double[3] containing coordinates of points to be transformed
+   */
+  void TetraAffineTransform::reverseApply(double* destPt, const double* srcPt) const
+  {
+    double* dest = destPt;
+    
+    // are we self-allocating ?
+    const bool selfAllocation = (destPt == srcPt);
+    
+    if(selfAllocation)
+      {
+        // alloc temporary memory
+        dest = new double[3];
+       
+        LOG(6, "Info : Self-affectation in TetraAffineTransform::reverseApply");
+      }
+    
+    for(int i = 0 ; i < 3 ; ++i)
+      {
+        // matrix - vector multiplication
+        dest[i] = _backLinearTransform[3*i] * srcPt[0] + _backLinearTransform[3*i + 1] * srcPt[1] + _backLinearTransform[3*i + 2] * srcPt[2];
+
+        // translation
+        dest[i] -= _translation[i];
+      }
+
+    if(selfAllocation)
+      {
+        // copy result back to destPt
+        for(int i = 0 ; i < 3 ; ++i)
+          {
+            destPt[i] = dest[i];
+          }
+        delete[] dest;
+      }
+  }
+
   /**
    * Returns the determinant of the linear part A of the transform.
    *
diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx
index 835bdbb28..138e53875 100644
--- a/src/INTERP_KERNEL/TetraAffineTransform.hxx
+++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx
@@ -21,6 +21,8 @@ namespace INTERP_KERNEL
 
     void apply(double* destPt, const double* srcPt) const;
 
+    void reverseApply(double* destPt, const double* srcPt) const;
+
     double determinant() const;
 
     void dump() const;
@@ -49,6 +51,8 @@ namespace INTERP_KERNEL
     /// The determinant of the matrix A is calculated at the construction of the object and cached.
     double _determinant;
 
+    /// 3x3 matrix AT is transposed A matrix used for y -> AT(y - b) transformation
+    double _backLinearTransform[9];
   };
 
 
-- 
2.39.2