Index: Code/Common/itkSymmetricSecondRankTensor.h
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Common/itkSymmetricSecondRankTensor.h,v
retrieving revision 1.25
diff -u -r1.25 itkSymmetricSecondRankTensor.h
--- Code/Common/itkSymmetricSecondRankTensor.h	10 Mar 2008 22:48:13 -0000	1.25
+++ Code/Common/itkSymmetricSecondRankTensor.h	20 Feb 2009 18:36:17 -0000
@@ -164,11 +164,10 @@
   void ComputeEigenAnalysis( EigenValuesArrayType & eigenValues,
                              EigenVectorsMatrixType & eigenVectors ) const;
 
-  /** Pre-Multiply by a Matrix as ResultingTensor = Matrix * ThisTensor. */
-  Self PreMultiply( const MatrixType & m ) const;
-
-  /** Post-Multiply by a Matrix as ResultingTensor = ThisTensor * Matrix. */
-  Self PostMultiply( const MatrixType & m ) const;
+  /** Rotates the tensor by the matrix m and returns the result.  The
+  matrix must specify a pure rotation matrix. Assuming the current
+  tensor is D the result is m * D * m^T. */
+  Self Rotate( const MatrixType & m ) const;
 
 private:
 
Index: Code/Common/itkSymmetricSecondRankTensor.txx
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Common/itkSymmetricSecondRankTensor.txx,v
retrieving revision 1.14
diff -u -r1.14 itkSymmetricSecondRankTensor.txx
--- Code/Common/itkSymmetricSecondRankTensor.txx	10 Mar 2008 22:48:13 -0000	1.14
+++ Code/Common/itkSymmetricSecondRankTensor.txx	20 Feb 2009 18:36:17 -0000
@@ -350,58 +350,36 @@
 }
 
 /*
- * Pre-multiply the Tensor by a Matrix
+ * Rotate the Tensor by a rotation matrix
  */
 template<class T,unsigned int NDimension>
 SymmetricSecondRankTensor<T,NDimension>
 SymmetricSecondRankTensor<T,NDimension>
-::PreMultiply( const MatrixType & m ) const
+::Rotate( const MatrixType & m ) const
 {
-Self result;
-typedef typename NumericTraits<T>::AccumulateType  AccumulateType;
-for(unsigned int r=0; r<NDimension; r++)
-  {
-  for(unsigned int c=0; c<NDimension; c++)
+  typedef typename MatrixType::InternalMatrixType VnlMatrixType;
+  VnlMatrixType vnlthis;
+  for(unsigned int r=0; r<NDimension; r++)
     {
-    AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
-    for(unsigned int t=0; t<NDimension; t++)
+    for(unsigned int c=0; c<NDimension; c++)
       {
-      sum += m(r,t) * (*this)(t,c);
+        vnlthis(r,c) = (*this)(r,c);
       }
-    result(r,c) = static_cast<T>( sum );
     }
-  }
-return result;
-}
+  
+  vnlthis = m.GetVnlMatrix() * vnlthis * m.GetTranspose();
 
-/*
- * Post-multiply the Tensor by a Matrix
- */
-template<class T,unsigned int NDimension>
-SymmetricSecondRankTensor<T,NDimension>
-SymmetricSecondRankTensor<T,NDimension>
-::PostMultiply( const MatrixType & m ) const
-{
-Self result;
-typedef typename NumericTraits<T>::AccumulateType  AccumulateType;
-for(unsigned int r=0; r<NDimension; r++)
-  {
-  for(unsigned int c=0; c<NDimension; c++)
+  Self result;
+  for(unsigned int r=0; r<NDimension; r++)
     {
-    AccumulateType sum = NumericTraits<AccumulateType>::ZeroValue();
-    for(unsigned int t=0; t<NDimension; t++)
+    for(unsigned int c=0; c<NDimension; c++)
       {
-      sum += (*this)(r,t) * m(t,c);
+        result(r,c) = vnlthis(r,c);
       }
-    result(r,c) = static_cast<T>( sum );
     }
-  }
-return result;
+  return result;
 }
 
-
-
-
 /**
  * Print content to an ostream
  */
Index: Testing/Code/Common/itkSymmetricSecondRankTensorTest.cxx
===================================================================
RCS file: /cvsroot/Insight/Insight/Testing/Code/Common/itkSymmetricSecondRankTensorTest.cxx,v
retrieving revision 1.11
diff -u -r1.11 itkSymmetricSecondRankTensorTest.cxx
--- Testing/Code/Common/itkSymmetricSecondRankTensorTest.cxx	10 Mar 2008 22:48:13 -0000	1.11
+++ Testing/Code/Common/itkSymmetricSecondRankTensorTest.cxx	20 Feb 2009 18:36:17 -0000
@@ -24,6 +24,27 @@
 #include "itkImage.h"
 #include "itkImageRegionIterator.h"
 
+namespace {
+
+template<class T>
+int check_tensors_equal(const T& a, const T& b)
+{
+  for(unsigned int i = 0; i < 3; ++i)
+  {
+    for(unsigned int j = 0; j < 3; ++j)
+    {
+      if(fabs(a(i,j) - b(i,j)) > 1e-8)
+      {
+        return EXIT_FAILURE;
+      }
+    }
+  }
+  std::cout << "Check [PASSED]" << std::endl;
+  return EXIT_SUCCESS;
+}
+
+};
+
 int itkSymmetricSecondRankTensorTest(int, char* [] )
 { 
   // Test it all
@@ -429,10 +450,38 @@
     tensor3D(2,2) =  29.0;
 
     Double3DMatrixType matrix3D;
+    matrix3D(0,0) = 0.707106781186548;
+    matrix3D(0,1) = -0.707106781186548;
+    matrix3D(0,2) = 0;
+    matrix3D(1,0) = 0.707106781186548;
+    matrix3D(1,1) = 0.707106781186548;
+    matrix3D(1,2) = 0;
+    matrix3D(2,0) = 0;
+    matrix3D(2,1) = 0;
+    matrix3D(2,2) = 1.0;
+    
+    // Correct rotation of tensor3D by matrix3D
+    Double3DTensorType rottensor3D;
+    rottensor3D[0] =  21.0;
+    rottensor3D[1] =  -2.0; 
+    rottensor3D[2] =   4.949747468305833;
+    rottensor3D[3] =  21.0;
+    rottensor3D[4] =   4.949747468305833; 
+    rottensor3D[5] =  29.0;
+
+    std::cout << "Original: " << std::endl << tensor3D << std::endl;
+    std::cout << "Rotation: " << std::endl << matrix3D << std::endl;
+    std::cout << "Result:   " << std::endl << rottensor3D << std::endl;
 
-    Double3DTensorType result1 = tensor3D.PreMultiply( matrix3D );
-    Double3DTensorType result2 = tensor3D.PostMultiply( matrix3D );
-
+    Double3DTensorType result1 = tensor3D.Rotate( matrix3D );
+    if(check_tensors_equal(result1, rottensor3D))
+    {
+      std::cerr << "Tensor rotation [FAILED]" << std::endl;
+      std::cerr << "Correct: " << rottensor3D << std::endl;
+      std::cerr << "Found  : " << result1 << std::endl;
+      return EXIT_FAILURE;
+    }
+    
   } // end of Matrix * SymmetricSecondRankTensor test
 
 

Index: Code/Common/itkSymmetricSecondRankTensor.h
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Common/itkSymmetricSecondRankTensor.h,v
retrieving revision 1.26
diff -r1.26 itkSymmetricSecondRankTensor.h
172a184,202
>   /** Rotate by the provided matix. */
>   template<typename TMatrixValueType>
>   Self Rotate( const Matrix<TMatrixValueType, NDimension, NDimension> & m);
> 
>   template<typename TMatrixValueType>
>   Self Rotate( const vnl_matrix_fixed<TMatrixValueType, NDimension, NDimension> & m)
>   {
>     return this->Rotate( static_cast<Matrix<TMatrixValueType, NDimension, NDimension> >(m) );
>   }
> 
>   template<typename TMatrixValueType>
>   Self Rotate( const vnl_matrix<TMatrixValueType> & m)
>   {
>     return this->Rotate( static_cast<Matrix<TMatrixValueType> >(m) );
>   }
> 
Index: Code/Common/itkSymmetricSecondRankTensor.txx
===================================================================
RCS file: /cvsroot/Insight/Insight/Code/Common/itkSymmetricSecondRankTensor.txx,v
retrieving revision 1.15
diff -r1.15 itkSymmetricSecondRankTensor.txx
380a381,457
> /*
>  * Set the Tensor to a Rotated version of the current tensor.
>  * matrix * self * Transpose(matrix)
>  * */
> template<class T,unsigned int NDimension>
> template <typename TMatrixValueType>
> SymmetricSecondRankTensor<T,NDimension>
> SymmetricSecondRankTensor<T,NDimension>
> ::Rotate( const Matrix<TMatrixValueType, NDimension, NDimension> & m )
> {
>   Self result;
>   typedef Matrix<double, NDimension, NDimension> RotationMatrixType;
>   RotationMatrixType SCT; //self * Transpose(m)
>   for(unsigned int r=0; r<NDimension; r++)
>   {
>     for(unsigned int c=0; c<NDimension; c++)
>     {
>       double sum = 0.0;
>       for(unsigned int t=0; t<NDimension; t++)
>       {
>         sum += (*this)(r,t) * m(c,t);
>       }
>       SCT(r,c) = sum;
>     }
>   }
>   //self = m * sct;
>   for(unsigned int r=0; r<NDimension; r++)
>   {
>     for(unsigned int c=0; c<NDimension; c++)
>     {
>       double sum = 0.0;
>       for(unsigned int t=0; t<NDimension; t++)
>       {
>         sum += m(r,t) * SCT(t,c);
>       }
>       (result)(r,c) = static_cast<T>( sum );
>     }
>   }
>   return result;
> }
>