NAMIC Wiki:DTI:ITK

Implementation of the DiffusionTensorPixel type

This page contains the code implementing the functionalities of the Diffusion Tensor pixel type. Given that the class is templated over the component type, the code is to be stored in a file with extension .txx. This is an initial draft of the suggested implementation.

Note that in particular it is missing the methods for returning EigenVectors and EigenValues. It is likely that we will be able to use a 3D specific implementation of eigen analysis in order to improve performance.

Another topic for discussion is the proper error management of out-of-bound calls. For example, when a user request the element tensor(5,3). The current suggested implementation is simply returning the element tensor(0,0), however it is arguable whether we should throw an exception and make the user aware of the problem.

Tensor operations are also open issues. In principle we should address

• Addition
• Subtraction
• Covariant product
• Tensor product with a vector

We probably should also provide operations for Tensor magnitude.

 #ifndef _itkDiffusionTensorPixel_txx
 #define _itkDiffusionTensorPixel_txx
 #include "itkDiffusionTensorPixel.h"
 #include "itkNumericTraits.h"
 
 namespace itk
 {
 
 /*
  * Assignment Operator
  */
 template<class T>
 DiffusionTensorPixel<T>&
 DiffusionTensorPixel<T>
 ::operator= (const Self& r)
 {
   BaseArray::operator=(r);
   return *this;
 }
 
 
 /*
  * Assigment from a plain array
  */
 template<class T>
 DiffusionTensorPixel<T>&
 DiffusionTensorPixel<T>
 ::operator= (const ComponentType r[Dimension])
 {
   BaseArray::operator=(r);
   return *this;
 }
 
 
 
 /**
  * Returns a temporary copy of a vector
  */
 template<class T>
 DiffusionTensorPixel<T>
 DiffusionTensorPixel<T>
 ::operator+(const Self & r) const
 {
   Self result;
   for( unsigned int i=0; i<Dimension; i++)
     {
     result[i] = (*this)[i] + r[i];
     }
   return result;
 }
 
 
 
 
 /**
  * Returns a temporary copy of a vector
  */
 template<class T>
 DiffusionTensorPixel<T>
 DiffusionTensorPixel<T>
 ::operator-(const Self & r) const
 {
   Self result;
   for( unsigned int i=0; i<Dimension; i++)
     {
     result[i] = (*this)[i] - r[i];
     }
   return result;
 }
 
 
 
 /**
  * Returns a temporary copy of a vector
  */
 template<class T>
 const DiffusionTensorPixel<T> &
 DiffusionTensorPixel<T>
 ::operator+=(const Self & r)
 {
   for( unsigned int i=0; i<Dimension; i++)
     {
     (*this)[i] += r[i];
     }
   return *this;
 }
 
 
 
 
 /**
  * Returns a temporary copy of a vector
  */
 template<class T>
 const DiffusionTensorPixel<T> &
 DiffusionTensorPixel<T>
 ::operator-=(const Self & r)
 {
   for( unsigned int i=0; i<Dimension; i++)
     {
     (*this)[i] -= r[i];
     }
   return *this;
 }
 
 
 
 
 
 /**
  * Returns a temporary copy of a vector
  */
 template<class T>
 DiffusionTensorPixel<T>
 DiffusionTensorPixel<T>
 ::operator*(const ComponentType & r) const
 {
   Self result;
   for( unsigned int i=0; i<Dimension; i++)
     {
     result[i] = (*this)[i] * r;
     }
   return result;
 }
 
 
 /*
  * Matrix notation access to elements
  */
 template<class T>
 const typename DiffusionTensorPixel<T>::ValueType &
 DiffusionTensorPixel<T>
 ::operator()(unsigned int i, unsigned int j) const
 {
   unsigned int k = i * 3 + j;
   if( k >= 5 )
     {
     return  (*this)[0];
     }
   return (*this)[k];
 }
 
 
 
 /*
  * Matrix notation access to elements
  */
 template<class T>
 typename DiffusionTensorPixel<T>::ValueType &
 DiffusionTensorPixel<T>
 ::operator()(unsigned int i, unsigned int j)
 {
   unsigned int k = i * 3 + j;
   if( k >= 5 )
     {
     return  (*this)[0];
     }
   return (*this)[k];
 }
 
 
 
 
 
 /**
  * Print content to an ostream
  */
 template<class TComponent>
 std::ostream &
 operator<<(std::ostream& os,const DiffusionTensorPixel<TComponent> & c )
 {
   for(unsigned int i=0; i<c.GetNumberOfComponents(); i++)
     {
     os <<  static_cast<typename NumericTraits<TComponent>::PrintType>(c[i]) << "  ";
     }
   return os;
 }
 
 
 /**
  * Read content from an istream
  */
 template<class TComponent>
 std::istream &
 operator>>(std::istream& is, DiffusionTensorPixel<TComponent> & c )
 {
   TComponent red;
   TComponent green;
   TComponent blue;
   for(unsigned int i=0; i<is.GetNumberOfComponents(); i++)
     {
     is >> si[i];
     }
   return is;
 }
 
 } // end namespace itk
 
 #endif