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:
−
== Implementation of the DiffusionTensorPixel type ==
+
== Test for the Diffusion Tensor pixel 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.
+
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.
−
 
 
−
<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.
 
−
 
 
−
<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.
 
−
 
 
−
<br /> Tensor operations are also open issues. In principle we should address
 
−
 
 
−
* Addition
 
−
* Subtraction
 
−
* Covariant product
 
−
* Tensor product with a vector
 
−
 
 
−
<br />
 
−
 
 
−
We probably should also provide operations for Tensor magnitude.
 
−
 
 
−
<br />
 
  
 
  <nowiki>
 
  <nowiki>
−
  #ifndef _itkDiffusionTensorPixel_txx
+
  #if defined(_MSC_VER)
−
  #define _itkDiffusionTensorPixel_txx
+
  #pragma warning ( disable : 4786 )
 +
#endif
 +
 +
#include <iostream>
 +
 
  #include "itkDiffusionTensorPixel.h"
 
  #include "itkDiffusionTensorPixel.h"
−
  #include "itkNumericTraits.h"
+
  #include "itkImage.h"
 +
#include "itkImageRegionIterator.h"
 
   
 
   
−
  namespace itk
+
  int itkDiffusionTensorPixelTest(int, char* [] )
 
  {
 
  {
 +
  int i, j;
 
   
 
   
−
/*
+
  // Test it all
−
  * Assignment Operator
 
−
  */
 
−
template<class T>
 
−
DiffusionTensorPixel<T>&
 
−
DiffusionTensorPixel<T>
 
−
::operator= (const Self& r)
 
−
{
 
−
  BaseArray::operator=(r);
 
−
  return *this;
 
−
}
 
 
   
 
   
 +
  float val[6] = {1.8, 0.2, 0.5, 3.4, 2.0, 1.2};
 +
  itk::DiffusionTensorPixel<float> pixel(val);
 +
  unsigned char pixelInit0[6] = {255, 255, 255};
 +
  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;
−
  * Assigment from a plain array
+
  if (sizeof(pixel) != 6 * sizeof(itk::DiffusionTensorPixel<float>::ComponentType))
−
  */
+
    {
−
template<class T>
+
    std::cerr << "ERROR: sizeof(pixel) == " << sizeof(pixel) << " but is should be " << 6 * sizeof(itk::DiffusionTensorPixel<float>::ComponentType) << std::endl;
−
DiffusionTensorPixel<T>&
+
    return 1;
−
DiffusionTensorPixel<T>
+
    }
−
::operator= (const ComponentType r[Dimension])
+
  std::cout << "pixel.GetNumberOfComponents = " << pixel.GetNumberOfComponents() << std::endl;
−
{
+
  std::cout << "pixel.GetNthComponent()" << std::endl;
−
  BaseArray::operator=(r);
+
  for (i = 0; i < pixel.GetNumberOfComponents(); i++)
−
  return *this;
+
    {
−
}
+
    std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
 +
    }
 
   
 
   
 +
  pixel(0,0) = 11.0;
 +
  pixel(0,1) = 21.0;
 +
  pixel(0,2) = 15.0;
 +
  pixel(1,0) = 11.0;
 +
  pixel(1,1) = 31.0;
 +
  pixel(1,2) = 10.0;
 +
  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;
−
/**
+
   for (i = 0; i < pixel.GetNumberOfComponents(); i++)
−
  * 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];
+
     std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
 
     }
 
     }
−
  return result;
 
−
}
 
 
   
 
   
 +
  std::cout << "pixel[0] = 111; pixel[1] = 222; pixel[2] = 333;" << std::endl;
 +
  std::cout << "pixel[3] = 444; pixel[4] = 555; pixel[5] = 666;" << std::endl;
 +
 +
  pixel[0] = 111; pixel[1] = 222; pixel[2] = 333;
 +
  pixel[3] = 444; pixel[4] = 555; pixel[5] = 666;
 
   
 
   
 +
  for (i = 0; i < pixel.GetNumberOfComponents(); i++)
 +
    {
 +
    std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
 +
    }
 
   
 
   
 +
  std::cout << "std::cout << pixel << std::endl;" << std::endl;
 +
  std::cout << "\t" << pixel << std::endl;
 
   
 
   
−
/**
+
   for (j = 0; j < 2; j++)
−
  * 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];
+
     std::cout << "pixelArray["<< j << "].GetNumberOfComponents = " << pixelArray[j].GetNumberOfComponents() << std::endl;
 +
    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;
 +
      }
 
     }
 
     }
−
  return result;
 
−
}
 
 
   
 
   
 +
  std::cout << "Testing arithmetic methods" << std::endl;
 +
  itk::DiffusionTensorPixel< float > pa;
 +
  itk::DiffusionTensorPixel< float > pb;
 +
 +
  pa[0] = 1.25;
 +
  pa[1] = 3.25;
 +
  pa[2] = 5.25;
 +
  pa[3] = 1.25;
 +
  pa[4] = 3.25;
 +
  pa[5] = 5.25;
 +
 +
  pb[0] = 1.55;
 +
  pb[1] = 3.55;
 +
  pb[2] = 5.55;
 +
  pb[3] = 1.55;
 +
  pb[4] = 3.55;
 +
  pb[5] = 5.55;
 +
 +
  itk::DiffusionTensorPixel< float > pc;
 +
 +
  pc = pa + pb;
 +
  std::cout << "addition = " << pc << std::endl;
 
   
 
   
 +
  pc = pa - pb;
 +
  std::cout << "subtraction = " << pc << std::endl;
 
   
 
   
−
/**
+
  pc += pb;
−
  * Returns a temporary copy of a vector
+
  std::cout << "in-place addition = " << pc << std::endl;
−
  */
 
−
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;
 
−
}
 
 
   
 
   
 +
  pc -= pb;
 +
  std::cout << "in-place subtraction = " << pc << std::endl;
 
   
 
   
 +
  pc = pa * 3.2;
 +
  std::cout << "product by scalar = " << pc << std::endl;
 
   
 
   
 
   
 
   
−
/**
 
−
  * 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;
 
−
}
 
 
   
 
   
 
   
 
   
 +
  /** Create an Image of tensors  */
 +
  typedef itk::DiffusionTensorPixel< float >  PixelType;
 +
  typedef itk::Image< PixelType, 3 >          ImageType;
 
   
 
   
 +
  ImageType::Pointer dti = ImageType::New();
 
   
 
   
 +
  ImageType::SizeType  size;
 +
  ImageType::IndexType start;
 +
  ImageType::RegionType region;
 
   
 
   
−
/**
+
  region.SetIndex( start );
−
  * Returns a temporary copy of a vector
+
   region.SetSize( size );
−
  */
 
−
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;
 
−
}
 
 
   
 
   
 +
  dti->SetRegions( region );
 +
  dti->Allocate();
 
   
 
   
−
/*
+
  ImageType::SpacingType spacing;
−
  * Matrix notation access to elements
+
   spacing[0] = 0.5;
−
  */
+
   spacing[1] = 0.5;
−
template<class T>
+
   spacing[2] = 1.5;
−
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];
 
−
}
 
 
   
 
   
 +
  ImageType::PointType origin;
 +
  origin[0] = 25.5;
 +
  origin[1] = 25.5;
 +
  origin[2] = 27.5;
 
   
 
   
 +
  dti->SetOrigin( origin );
 +
  dti->SetSpacing( spacing );
 
   
 
   
−
/*
+
   PixelType tensor;
−
  * 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];
 
−
}
 
 
   
 
   
 +
  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() )
−
  * 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]) << "  ";
+
     it.Set( tensor );
 +
    ++it;
 
     }
 
     }
−
  return os;
 
−
}
 
 
   
 
   
 
   
 
   
−
/**
+
   return EXIT_SUCCESS;
−
  * 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>
 
  </nowiki>

Revision as of 14:06, 18 December 2006

Home < NAMIC Wiki:DTI:ITK

Test for the Diffusion Tensor pixel type

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. 

 #if defined(_MSC_VER)
 #pragma warning ( disable : 4786 )
 #endif
 
 #include <iostream>
 
 #include "itkDiffusionTensorPixel.h"
 #include "itkImage.h"
 #include "itkImageRegionIterator.h"
 
 int itkDiffusionTensorPixelTest(int, char* [] )
 {
   int i, j;
 
   // Test it all
 
   float val[6] = {1.8, 0.2, 0.5, 3.4, 2.0, 1.2};
   itk::DiffusionTensorPixel<float> pixel(val);
   unsigned char pixelInit0[6] = {255, 255, 255};
   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;
   pixel(0,2) = 15.0;
   pixel(1,0) = 11.0;
   pixel(1,1) = 31.0;
   pixel(1,2) = 10.0;
   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;
   for (i = 0; i < pixel.GetNumberOfComponents(); i++)
     {
     std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
     }
 
   std::cout << "pixel[0] = 111; pixel[1] = 222; pixel[2] = 333;" << std::endl;
   std::cout << "pixel[3] = 444; pixel[4] = 555; pixel[5] = 666;" << std::endl;
 
   pixel[0] = 111; pixel[1] = 222; pixel[2] = 333;
   pixel[3] = 444; pixel[4] = 555; pixel[5] = 666;
 
   for (i = 0; i < pixel.GetNumberOfComponents(); i++)
     {
     std::cout << "\tpixel[" << i << "] = " << pixel.GetNthComponent(i) << std::endl;
     }
 
   std::cout << "std::cout << pixel << std::endl;" << std::endl;
   std::cout << "\t" << pixel << std::endl;
 
   for (j = 0; j < 2; j++)
     {
     std::cout << "pixelArray["<< j << "].GetNumberOfComponents = " << pixelArray[j].GetNumberOfComponents() << std::endl;
     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;
   itk::DiffusionTensorPixel< float > pa;
   itk::DiffusionTensorPixel< float > pb;
 
   pa[0] = 1.25;
   pa[1] = 3.25;
   pa[2] = 5.25;
   pa[3] = 1.25;
   pa[4] = 3.25;
   pa[5] = 5.25;
 
   pb[0] = 1.55;
   pb[1] = 3.55;
   pb[2] = 5.55;
   pb[3] = 1.55;
   pb[4] = 3.55;
   pb[5] = 5.55;
 
   itk::DiffusionTensorPixel< float > pc;
 
   pc = pa + pb;
   std::cout << "addition = " << pc << std::endl;
 
   pc = pa - pb;
   std::cout << "subtraction = " << pc << std::endl;
 
   pc += pb;
   std::cout << "in-place addition = " << pc << std::endl;
 
   pc -= pb;
   std::cout << "in-place subtraction = " << pc << std::endl;
 
   pc = pa * 3.2;
   std::cout << "product by scalar = " << pc << std::endl;
 
 
 
 
   /** Create an Image of tensors  */
   typedef itk::DiffusionTensorPixel< float >   PixelType;
   typedef itk::Image< PixelType, 3 >           ImageType;
 
   ImageType::Pointer dti = ImageType::New();
 
   ImageType::SizeType  size;
   ImageType::IndexType start;
   ImageType::RegionType region;
 
   region.SetIndex( start );
   region.SetSize( size );
 
   dti->SetRegions( region );
   dti->Allocate();
 
   ImageType::SpacingType spacing;
   spacing[0] = 0.5;
   spacing[1] = 0.5;
   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;
 }
Retrieved from "https://www.na-mic.org/w/index.php?title=NAMIC_Wiki:DTI:ITK&oldid=4530"