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:
== Test for the Diffusion Tensor pixel type ==
+
The DiffusionTensorPixel type derives from the itk::FixedArray<> instantiated for a fixed size of 6 elements. The rationale behind it is that we are dealing with symmetric tensors in this particular case, therefore for the sake of memory use we only store 6 components. The type however will also implement a matrix-like notation in order to provide a natural notation for tensor operations.
  
The following code is a typical test file for the functionalities implemented in the DiffusionTensorPixel type. It also exercises the use of this pixel type in an itk::Image instantiation.
+
A short version of the class essentials is below, a full compilable version is at the end of this page
  
<nowiki>
 
#if defined(_MSC_VER)
 
#pragma warning ( disable : 4786 )
 
#endif
 
 
#include <iostream>
 
 
#include "itkDiffusionTensorPixel.h"
 
#include "itkImage.h"
 
#include "itkImageRegionIterator.h"
 
 
   
 
   
  int itkDiffusionTensorPixelTest(int, char* [] )
+
  template < typename TComponent = float >
 +
class DiffusionTensorPixel: public FixedArray<TComponent,6>
 
  {
 
  {
   int i, j;
+
public:
 +
   typedef DiffusionTensorPixel  Self;
 +
  typedef FixedArray<TComponent, 6> SuperClass;
 +
  itkStaticConstMacro(Dimension, unsigned int, 6);
 +
  itkStaticConstMacro(Length, unsigned int, 6);
 +
  typedef FixedArray<TComponent, 6> BaseArray;
 
   
 
   
   // Test it all
+
   /** Array of eigen-values. */
 +
  typedef FixedArray<TComponent, 3> EigenValuesArray;
 +
  typedef Matrix<TComponent, 3, 3>  EigenVectorsMatrix;
 
   
 
   
   float val[6] = {1.8, 0.2, 0.5, 3.4, 2.0, 1.2};
+
   /** Matrix notation, in const and non-const forms. */
   itk::DiffusionTensorPixel<float> pixel(val);
+
   ValueType & operator()( unsigned int i, unsigned int j);
  unsigned char pixelInit0[6] = {255, 255, 255};
+
   const ValueType & operator()( unsigned int i, unsigned int j) const;
   unsigned char pixelInit1[6] = {255, 255, 244};
+
};
  itk::DiffusionTensorPixel<unsigned char> pixelArray[2];
 
  pixelArray[0] = pixelInit0;
 
  pixelArray[1] = pixelInit1;
 
 
   
 
   
  std::cout << "sizeof(pixel) = " << sizeof (pixel) << std::endl;
 
  if (sizeof(pixel) != 6 * sizeof(itk::DiffusionTensorPixel<float>::ComponentType))
 
    {
 
    std::cerr << "ERROR: sizeof(pixel) == " << sizeof(pixel) << " but is should be " << 6 * sizeof(itk::DiffusionTensorPixel<float>::ComponentType) << std::endl;
 
    return 1;
 
    }
 
  std::cout << "pixel.GetNumberOfComponents = " << pixel.GetNumberOfComponents() << std::endl;
 
  std::cout << "pixel.GetNthComponent()" << std::endl;
 
  for (i = 0; i < pixel.GetNumberOfComponents(); i++)
 
    {
 
    std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
 
    }
 
 
   
 
   
  pixel(0,0) = 11.0;
+
 
  pixel(0,1) = 21.0;
+
<br /> Here is the full version of the class. This code is compilable using ITK 2.0.
  pixel(0,2) = 15.0;
+
 
  pixel(1,0) = 11.0;
+
<nowiki>
  pixel(1,1) = 31.0;
+
#ifndef __itkDiffusionTensorPixel_h
  pixel(1,2) = 10.0;
+
#define __itkDiffusionTensorPixel_h
  pixel(2,0) = 11.0; // these three last element will overwrite its symmetric counterparts
 
  pixel(2,1) = 41.0;
 
  pixel(2,2) = 14.0;
 
 
   
 
   
  std::cout << "testing the pixel(i,j) APID" << std::endl;
+
// Undefine an eventual DiffusionTensorPixel macro
  for (i = 0; i < pixel.GetNumberOfComponents(); i++)
+
#ifdef DiffusionTensorPixel
    {
+
#undef DiffusionTensorPixel
    std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
+
#endif
    }
 
 
   
 
   
  std::cout << "pixel[0] = 111; pixel[1] = 222; pixel[2] = 333;" << std::endl;
+
#include <itkIndent.h>
  std::cout << "pixel[3] = 444; pixel[4] = 555; pixel[5] = 666;" << std::endl;
+
#include <itkFixedArray.h>
 +
#include <itkMatrix.h>
 +
#include "vnl/vnl_math.h"
 
   
 
   
  pixel[0] = 111; pixel[1] = 222; pixel[2] = 333;
+
namespace itk
  pixel[3] = 444; pixel[4] = 555; pixel[5] = 666;
+
{
 
   
 
   
  for (i = 0; i < pixel.GetNumberOfComponents(); i++)
+
/** \class DiffusionTensorPixel
    {
+
  * \brief Represent a diffusion tensor as used in DTI images.
    std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
+
  *
    }
+
  * This class implements a 3D symmetric tensor as it is used for representing
 +
  * diffusion of water molecules in Diffusion Tensor Images.
 +
  *
 +
  * Since DiffusionTensorPixel is a subclass of Array, you can access its components as:
 +
  *
 +
  * typedef itk::DiffusionTensorPixel< float >    TensorPixelType;
 +
  * TensorPixelType tensor;
 +
  *
 +
  *  tensor[0] = 1.233;
 +
  *  tensor[1] = 1.456;
 +
  *
 +
  * for convenience the indexed access is also available as
 +
  *
 +
  *  tensor(0,0) = 1.233;
 +
  *  tensor(2,0) = 1.233;
 +
  *
 +
  * The Tensor in principle represents a 3x3 matrix, but given that it is
 +
  * always symmetric the representation can be compacted into a 6-elements
 +
  * array that derives from the itk::FixedArray<T>
 +
  *
 +
  * \ingroup ImageObjects
 +
  */
 
   
 
   
  std::cout << "std::cout << pixel << std::endl;" << std::endl;
+
template < typename TComponent = float >
   std::cout << "\t" << pixel << std::endl;
+
class DiffusionTensorPixel: public FixedArray<TComponent,6>
 +
{
 +
public:
 +
  /** Standard class typedefs. */
 +
  typedef DiffusionTensorPixel  Self;
 +
   typedef FixedArray<TComponent, 6> SuperClass;
 
   
 
   
   for (j = 0; j < 2; j++)
+
   /** Dimension of the vector space. */
    {
+
  itkStaticConstMacro(Dimension, unsigned int, 6);
    std::cout << "pixelArray["<< j << "].GetNumberOfComponents = " << pixelArray[j].GetNumberOfComponents() << std::endl;
+
  itkStaticConstMacro(Length, unsigned int, 6);
    std::cout << "pixelArray[" << j << "].GetNthComponent()" << std::endl;
 
    for (i = 0; i < pixelArray[j].GetNumberOfComponents(); i++)
 
      {
 
      std::cout << "\tpixelArray[" << j << "].GetNthComponent(" << i << ") = " << static_cast<int>(pixelArray[j].GetNthComponent(i)) << std::endl;
 
      }
 
    }
 
 
   
 
   
   std::cout << "Testing arithmetic methods" << std::endl;
+
   /** Convenience typedefs. */
  itk::DiffusionTensorPixel< float > pa;
+
   typedef FixedArray<TComponent, 6> BaseArray;
   itk::DiffusionTensorPixel< float > pb;
 
 
   
 
   
   pa[0] = 1.25;
+
   /** Array of eigen-values. */
   pa[1] = 3.25;
+
   typedef FixedArray<TComponent, 3> EigenValuesArray;
  pa[2] = 5.25;
 
  pa[3] = 1.25;
 
  pa[4] = 3.25;
 
  pa[5] = 5.25;
 
 
   
 
   
   pb[0] = 1.55;
+
   /** Matrix of eigen-vectors. */
   pb[1] = 3.55;
+
   typedef Matrix<TComponent, 3, 3> EigenVectorsMatrix;
  pb[2] = 5.55;
 
  pb[3] = 1.55;
 
  pb[4] = 3.55;
 
  pb[5] = 5.55;
 
 
   
 
   
   itk::DiffusionTensorPixel< float > pc;
+
   /**  Define the component type. */
 +
  typedef TComponent ComponentType;
 +
  typedef typename SuperClass::ValueType ValueType;
 
   
 
   
   pc = pa + pb;
+
   /** Default constructor has nothing to do. */
   std::cout << "addition = " << pc << std::endl;
+
  DiffusionTensorPixel() {this->Fill(0);}
 +
   DiffusionTensorPixel (const ComponentType& r) { this->Fill(r); }
 
   
 
   
   pc = pa - pb;
+
   /** Pass-through constructor for the Array base class. */
   std::cout << "subtraction = " << pc << std::endl;
+
   DiffusionTensorPixel(const Self& r): BaseArray(r) {}
 +
  DiffusionTensorPixel(const ComponentType  r[6]): BaseArray(r) {}
 
   
 
   
   pc += pb;
+
   /** Pass-through assignment operator for the Array base class. */
   std::cout << "in-place addition = " << pc << std::endl;
+
  Self& operator= (const Self& r);
 +
   Self& operator= (const ComponentType r[6]);
 
   
 
   
   pc -= pb;
+
   /** Aritmetic operations between pixels. Return a new DiffusionTensorPixel. */
   std::cout << "in-place subtraction = " << pc << std::endl;
+
  Self operator+(const Self &vec) const;
 +
  Self operator-(const Self &vec) const;
 +
  const Self & operator+=(const Self &vec);
 +
   const Self & operator-=(const Self &vec);
 +
  Self operator*(const ComponentType &f) const;
 
   
 
   
  pc = pa * 3.2;
 
  std::cout << "product by scalar = " << pc << std::endl;
 
 
   
 
   
 +
  /** Return the number of components. */
 +
  static int GetNumberOfComponents(){ return 6;}
 
   
 
   
 +
  /** Return the value for the Nth component. */
 +
  ComponentType GetNthComponent(int c) const
 +
    { return this->operator[](c); }
 
   
 
   
 +
  /** Set the Nth component to v. */
 +
  void SetNthComponent(int c, const ComponentType& v)
 +
    {  this->operator[](c) = v; }
 
   
 
   
   /** Create an Image of tensors  */
+
   /** Matrix notation, in const and non-const forms. */
   typedef itk::DiffusionTensorPixel< float >  PixelType;
+
   ValueType & operator()( unsigned int i, unsigned int j);
   typedef itk::Image< PixelType, 3 >          ImageType;
+
   const ValueType & operator()( unsigned int i, unsigned int j) const;
 
   
 
   
  ImageType::Pointer dti = ImageType::New();
+
};
 
   
 
   
  ImageType::SizeType  size;
 
  ImageType::IndexType start;
 
  ImageType::RegionType region;
 
 
   
 
   
  region.SetIndex( start );
+
template< typename TComponent  >
  region.SetSize( size );
+
ITK_EXPORT std::ostream& operator<<(std::ostream& os,
 +
                                    const DiffusionTensorPixel<TComponent> & c);
 +
template< typename TComponent  >
 +
ITK_EXPORT std::istream& operator>>(std::istream& is,
 +
                                          DiffusionTensorPixel<TComponent> & c);
 
   
 
   
  dti->SetRegions( region );
+
} // end namespace itk
  dti->Allocate();
 
 
   
 
   
  ImageType::SpacingType spacing;
+
  #ifndef ITK_MANUAL_INSTANTIATION
  spacing[0] = 0.5;
+
  #include "itkDiffusionTensorPixel.txx"
  spacing[1] = 0.5;
+
  #endif
  spacing[2] = 1.5;
 
   
 
  ImageType::PointType origin;
 
  origin[0] = 25.5;
 
  origin[1] = 25.5;
 
  origin[2] = 27.5;
 
   
 
  dti->SetOrigin( origin );
 
  dti->SetSpacing( spacing );
 
 
  PixelType tensor;
 
 
  tensor[0] = 1.2;
 
  tensor[1] = 2.2;
 
  tensor[2] = 3.2;
 
  tensor[3] = 4.2;
 
  tensor[4] = 5.2;
 
  tensor[5] = 6.2;
 
 
  dti->FillBuffer( tensor );
 
 
  typedef itk::ImageRegionIterator< ImageType > IteratorType;
 
 
  IteratorType it( dti, region );
 
  it.GoToBegin();
 
 
  while( !it.IsAtEnd() )
 
    {
 
    it.Set( tensor );
 
    ++it;
 
    }
 
 
 
  return EXIT_SUCCESS;
 
  }
 
 
   
 
   
 +
#endif
 
  </nowiki>
 
  </nowiki>

Revision as of 13:31, 18 December 2006

Home < NAMIC Wiki:DTI:ITK

The DiffusionTensorPixel type derives from the itk::FixedArray<> instantiated for a fixed size of 6 elements. The rationale behind it is that we are dealing with symmetric tensors in this particular case, therefore for the sake of memory use we only store 6 components. The type however will also implement a matrix-like notation in order to provide a natural notation for tensor operations.

A short version of the class essentials is below, a full compilable version is at the end of this page


template < typename TComponent = float >
class DiffusionTensorPixel: public FixedArray<TComponent,6>
{
public:
  typedef DiffusionTensorPixel  Self;
  typedef FixedArray<TComponent, 6> SuperClass;
  itkStaticConstMacro(Dimension, unsigned int, 6);
  itkStaticConstMacro(Length, unsigned int, 6);
  typedef FixedArray<TComponent, 6> BaseArray;

  /** Array of eigen-values. */
  typedef FixedArray<TComponent, 3> EigenValuesArray;
  typedef Matrix<TComponent, 3, 3>  EigenVectorsMatrix;

  /** Matrix notation, in const and non-const forms. */
  ValueType & operator()( unsigned int i, unsigned int j);
  const ValueType & operator()( unsigned int i, unsigned int j) const;
};



Here is the full version of the class. This code is compilable using ITK 2.0.

 #ifndef __itkDiffusionTensorPixel_h
 #define __itkDiffusionTensorPixel_h
 
 // Undefine an eventual DiffusionTensorPixel macro
 #ifdef DiffusionTensorPixel
 #undef DiffusionTensorPixel
 #endif
 
 #include <itkIndent.h>
 #include <itkFixedArray.h>
 #include <itkMatrix.h>
 #include "vnl/vnl_math.h"
 
 namespace itk
 {
 
 /** \class DiffusionTensorPixel
  * \brief Represent a diffusion tensor as used in DTI images.
  *
  * This class implements a 3D symmetric tensor as it is used for representing
  * diffusion of water molecules in Diffusion Tensor Images.
  *
  * Since DiffusionTensorPixel is a subclass of Array, you can access its components as:
  *
  * typedef itk::DiffusionTensorPixel< float >    TensorPixelType;
  * TensorPixelType tensor;
  *
  *   tensor[0] = 1.233;
  *   tensor[1] = 1.456;
  *
  * for convenience the indexed access is also available as
  *
  *   tensor(0,0) = 1.233;
  *   tensor(2,0) = 1.233;
  *
  * The Tensor in principle represents a 3x3 matrix, but given that it is
  * always symmetric the representation can be compacted into a 6-elements
  * array that derives from the itk::FixedArray<T>
  *
  * \ingroup ImageObjects
  */
 
 template < typename TComponent = float >
 class DiffusionTensorPixel: public FixedArray<TComponent,6>
 {
 public:
   /** Standard class typedefs. */
   typedef DiffusionTensorPixel  Self;
   typedef FixedArray<TComponent, 6> SuperClass;
 
   /** Dimension of the vector space. */
   itkStaticConstMacro(Dimension, unsigned int, 6);
   itkStaticConstMacro(Length, unsigned int, 6);
 
   /** Convenience typedefs. */
   typedef FixedArray<TComponent, 6> BaseArray;
 
   /** Array of eigen-values. */
   typedef FixedArray<TComponent, 3> EigenValuesArray;
 
   /** Matrix of eigen-vectors. */
   typedef Matrix<TComponent, 3, 3> EigenVectorsMatrix;
 
   /**  Define the component type. */
   typedef TComponent ComponentType;
   typedef typename SuperClass::ValueType ValueType;
 
   /** Default constructor has nothing to do. */
   DiffusionTensorPixel() {this->Fill(0);}
   DiffusionTensorPixel (const ComponentType& r) { this->Fill(r); }
 
   /** Pass-through constructor for the Array base class. */
   DiffusionTensorPixel(const Self& r): BaseArray(r) {}
   DiffusionTensorPixel(const ComponentType  r[6]): BaseArray(r) {}
 
   /** Pass-through assignment operator for the Array base class. */
   Self& operator= (const Self& r);
   Self& operator= (const ComponentType r[6]);
 
   /** Aritmetic operations between pixels. Return a new DiffusionTensorPixel. */
   Self operator+(const Self &vec) const;
   Self operator-(const Self &vec) const;
   const Self & operator+=(const Self &vec);
   const Self & operator-=(const Self &vec);
   Self operator*(const ComponentType &f) const;
 
 
   /** Return the number of components. */
   static int GetNumberOfComponents(){ return 6;}
 
   /** Return the value for the Nth component. */
   ComponentType GetNthComponent(int c) const
     { return this->operator[](c); }
 
   /** Set the Nth component to v. */
   void SetNthComponent(int c, const ComponentType& v)
     {  this->operator[](c) = v; }
 
   /** Matrix notation, in const and non-const forms. */
   ValueType & operator()( unsigned int i, unsigned int j);
   const ValueType & operator()( unsigned int i, unsigned int j) const;
 
 };
 
 
 template< typename TComponent  >
 ITK_EXPORT std::ostream& operator<<(std::ostream& os,
                                     const DiffusionTensorPixel<TComponent> & c);
 template< typename TComponent  >
 ITK_EXPORT std::istream& operator>>(std::istream& is,
                                           DiffusionTensorPixel<TComponent> & c);
 
 } // end namespace itk
 
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "itkDiffusionTensorPixel.txx"
 #endif
 
 #endif