Difference between revisions of "NAMIC Wiki:DTI:ITK"

From NAMIC Wiki
Jump to: navigation, search
m (Update from Wiki)
m (Update from Wiki)
Line 1: Line 1:
−
This page host a discussion on the implementation details of a Tensor class.
+
== Implementation of the DiffusionTensorPixel type ==
  
−
== Tensors and Traits ==
+
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.
  
−
Torsten Rohlfing for SRI International Neuroscience Program proposed the following for ITK tensor processing:
+
<br /> 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.
  
−
* Separate the operations which depend on the tensor size from those that do not.
+
<br /> 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.
−
** The former are in tensor traits, the latter are in tensor class.
 
−
* Filter classes operate on different tensor pixel type images (by templating).
 
−
** The tensor classes can be more or less specialized. In particular, the filters will be able to support more general tensor classes that may be implemented later, as long as they provide the necessary interface functions for any given filter.
 
−
* [[NAMIC_Wiki:DTI:ITK-SymmetricTensorPixelType:Header|itkSymmetricTensor.h]]
 
−
** Contains functions that do not benefit from knowledge of the tensor size, e.g., GetVnlMatrix()
 
−
* [[NAMIC_Wiki:DTI:ITK-SymmetricTensorTraits::Header|itkSymmetricTensorTraits.h]]
 
−
** Operations that depend on the tensor size, e.g., eigenvalue computation.
 
−
* itkTensorToFractionalAnisotropyImageFilter.{h,txx}
 
−
** Example filter that operates on tensor data.
 
−
** Operates on any tensor class that implements T::ComputeEigenvalues()
 
−
* DiffusionTensorToFractionalAnisotropy.cxx
 
−
** Example application. Reads a tensor image from a raw data file, computes the FA image, and writes the result.
 
  
−
=== A First version of the SymmetricSecondRankTensor class committed to ITK ===
+
<br /> Tensor operations are also open issues. In principle we should address
  
−
Given that Tensor can have any rank, from zero (scalar), one (vector), two (matrices) to N. It seemed appropriate to specify that this particular class was representing a symmetric tensor of second rank.
+
* Addition
 +
* Subtraction
 +
* Covariant product
 +
* Tensor product with a vector
  
−
The class has been committed into ITK (May 2 2005) and it is expected to evolve in the CVS repository.
+
<br />
  
−
=== A first version of SymmetricEigenAnalysis class committed to ITK ===
+
We probably should also provide operations for Tensor magnitude.
  
−
Serves as a thread safe alternative to vnl_symmetric_eigensystem, which calls netlib C routines is not thread safe. Please use this class in any multi-threaded filters.
+
<br />
  
−
=== A Hessian filter uses the Symmetric second rank tensor class ===
+
<nowiki>
−
 
+
#ifndef _itkDiffusionTensorPixel_txx
−
A filter for computing the Hessian of an image was also committed to ITK. This filter illustrates the first use of the SymmetricSecondRankTensor class as pixel type of an image.
+
#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
 +
</nowiki>

Revision as of 14:06, 18 December 2006

Home < 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
Retrieved from "https://www.na-mic.org/w/index.php?title=NAMIC_Wiki:DTI:ITK&oldid=4529"