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

From NAMIC Wiki
Jump to: navigation, search
m (Update from Wiki)
m (Update from Wiki)
 
(6 intermediate revisions by the same user not shown)
Line 1: Line 1:
== Test for the Diffusion Tensor pixel type ==
+
This page host a discussion on the implementation details of a Tensor class.
  
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.
+
== Tensors and Traits ==
  
<nowiki>
+
Torsten Rohlfing for SRI International Neuroscience Program proposed the following for ITK tensor processing:
#if defined(_MSC_VER)
+
 
#pragma warning ( disable : 4786 )
+
* Separate the operations which depend on the tensor size from those that do not.
#endif
+
** The former are in tensor traits, the latter are in tensor class.
+
* Filter classes operate on different tensor pixel type images (by templating).
#include <iostream>
+
** 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]]
#include "itkDiffusionTensorPixel.h"
+
** Contains functions that do not benefit from knowledge of the tensor size, e.g., GetVnlMatrix()
#include "itkImage.h"
+
* [[NAMIC_Wiki:DTI:ITK-SymmetricTensorTraits::Header|itkSymmetricTensorTraits.h]]
#include "itkImageRegionIterator.h"
+
** Operations that depend on the tensor size, e.g., eigenvalue computation.
+
* itkTensorToFractionalAnisotropyImageFilter.{h,txx}
int itkDiffusionTensorPixelTest(int, char* [] )
+
** Example filter that operates on tensor data.
{
+
** Operates on any tensor class that implements T::ComputeEigenvalues()
  int i, j;
+
* DiffusionTensorToFractionalAnisotropy.cxx
+
** Example application. Reads a tensor image from a raw data file, computes the FA image, and writes the result.
  // Test it all
+
 
+
=== A First version of the SymmetricSecondRankTensor class committed to ITK ===
  float val[6] = {1.8, 0.2, 0.5, 3.4, 2.0, 1.2};
+
 
  itk::DiffusionTensorPixel<float> pixel(val);
+
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.
  unsigned char pixelInit0[6] = {255, 255, 255};
+
 
  unsigned char pixelInit1[6] = {255, 255, 244};
+
The class has been committed into ITK (May 2 2005) and it is expected to evolve in the CVS repository.
  itk::DiffusionTensorPixel<unsigned char> pixelArray[2];
+
 
  pixelArray[0] = pixelInit0;
+
=== A first version of SymmetricEigenAnalysis class committed to ITK ===
  pixelArray[1] = pixelInit1;
+
 
+
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.
  std::cout << "sizeof(pixel) = " << sizeof (pixel) << std::endl;
+
 
  if (sizeof(pixel) != 6 * sizeof(itk::DiffusionTensorPixel<float>::ComponentType))
+
=== A Hessian filter uses the Symmetric second rank tensor class ===
    {
+
 
    std::cerr << "ERROR: sizeof(pixel) == " << sizeof(pixel) << " but is should be " << 6 * sizeof(itk::DiffusionTensorPixel<float>::ComponentType) << std::endl;
+
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.
    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;
 
}
 
 
</nowiki>
 

Latest revision as of 14:06, 18 December 2006

Home < NAMIC Wiki:DTI:ITK

This page host a discussion on the implementation details of a Tensor class.

Tensors and Traits

Torsten Rohlfing for SRI International Neuroscience Program proposed the following for ITK tensor processing:

  • Separate the operations which depend on the tensor size from those that do not.
    • 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.
  • itkSymmetricTensor.h
    • Contains functions that do not benefit from knowledge of the tensor size, e.g., GetVnlMatrix()
  • 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

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.

The class has been committed into ITK (May 2 2005) and it is expected to evolve in the CVS repository.

A first version of SymmetricEigenAnalysis class committed to ITK

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.

A Hessian filter uses the Symmetric second rank tensor class

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.