lib/cpExtensions/Algorithms/SecondRankDiffusionTensorToPolyData.cxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #include <cpExtensions/Algorithms/SecondRankDiffusionTensorToPolyData.h>
   6 #include <itkImage.h>
   7 #include <itkSymmetricSecondRankTensor.h>
   8 #include <itkImageRegionConstIteratorWithIndex.h>
   9 #include <vtkInformation.h>
  10 #include <vtkInformationVector.h>
  11 #include <vtkPoints.h>
  12 #include <vtkPointData.h>
  13 #include <vtkSmartPointer.h>
  14 #include <vtkDoubleArray.h>
  15
  16 // -------------------------------------------------------------------------
  17 template< class I >
  18 typename cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
  19 Self* cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
  20 New( )
  21 {
  22   return( new Self );
  23 }
  24
  25 // -------------------------------------------------------------------------
  26 template< class I >
  27 void cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
  28 SetInputData( const I* image )
  29 {
  30   this->m_ITKImage = image;
  31   this->Modified( );
  32 }
  33
  34 // -------------------------------------------------------------------------
  35 template< class I >
  36 cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
  37 SecondRankDiffusionTensorToPolyData( )
  38   : vtkPolyDataAlgorithm( )
  39 {
  40   this->SetNumberOfInputPorts( 0 );
  41 }
  42
  43 // -------------------------------------------------------------------------
  44 template< class I >
  45 cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
  46 ~SecondRankDiffusionTensorToPolyData( )
  47 {
  48 }
  49
  50 // -------------------------------------------------------------------------
  51 template< class I >
  52 int cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
  53 RequestData(
  54   vtkInformation* request,
  55   vtkInformationVector** inputVector,
  56   vtkInformationVector* outputVector
  57   )
  58 {
  59   if( this->m_ITKImage.IsNull( ) )
  60     return( 0 );
  61
  62   // get the info objects
  63   vtkInformation* outInfo = outputVector->GetInformationObject( 0 );
  64   vtkPolyData* polyData = vtkPolyData::SafeDownCast(
  65     outInfo->Get(vtkDataObject::DATA_OBJECT())
  66     );
  67
  68   itk::ImageRegionConstIteratorWithIndex< I > tensorIt( this->m_ITKImage, this->m_ITKImage->GetRequestedRegion() );
  69
  70   vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
  71   vtkSmartPointer<vtkDoubleArray> vtkTensors = vtkSmartPointer<vtkDoubleArray>::New();
  72   vtkTensors->SetNumberOfComponents(9);
  73   for ( tensorIt.GoToBegin(); !tensorIt.IsAtEnd(); ++tensorIt )
  74   {
  75     typename I::IndexType tensorIndex = tensorIt.GetIndex();
  76     typename I::PixelType tensor = tensorIt.Get();
  77     typename I::PointType ptensor;
  78
  79     // TODO: std::cout << "\n Index [" << tensorIndex[0] << ", " << tensorIndex[1] << ", " << tensorIndex[2] << "]" << std::endl;
  80     this->m_ITKImage->TransformIndexToPhysicalPoint( tensorIndex, ptensor );
  81     vtkIdType vtkPointID = points->InsertNextPoint(ptensor[0], ptensor[1], ptensor[2]);
  82
  83     /**
  84      * Definition included in cdifTypes.
  85      * typedef itk::SymmetricSecondRankTensor< double > Diffusion2ndRankTensor;
  86      */
  87     //Diffusion2ndRankTensor Diffusion2ndRankTensor;
  88
  89
  90     /**
  91      * SYMMETRIC SECOND RANK TENSOR
  92      *
  93      * The upper-right triangle of the matrix:
  94      *
  95      *       | 0  1  2  |
  96      *       | X  3  4  |
  97      *       | X  X  5  |
  98      *
  99      *       | xx  xy  xz  |
 100      *       | X   yy  yz  |
 101      *       | X   X   zz  |
 102      */
 103     /*
 104       Diffusion2ndRankTensor = tensor[0];
 105       Diffusion2ndRankTensor = tensor[1];
 106       Diffusion2ndRankTensor = tensor[2];
 107       Diffusion2ndRankTensor = tensor[3];
 108       Diffusion2ndRankTensor = tensor[4];
 109       Diffusion2ndRankTensor = tensor[5];
 110     */
 111     typename I::PixelType::EigenValuesArrayType    evalues;
 112     typename I::PixelType::EigenVectorsMatrixType  evectors;
 113
 114     /**
 115      * Return an array containing EigenValues,
 116      * and a matrix containing Eigenvectors.
 117      */
 118     tensor.ComputeEigenAnalysis(evalues, evectors);
 119
 120     /*
 121       std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen values: " << evalues[0] << ", " << evalues[1] << ", " << evalues[2] << std::endl;
 122
 123       std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen vectors: "
 124
 125       << "(" << evectors(0,0) << ", " << evectors(0,1) << ", " << evectors(0,2) << "), "
 126       << "(" << evectors(1,0) << ", " << evectors(1,1) << ", " << evectors(1,2) << "), "
 127       << "(" << evectors(2,0) << ", " << evectors(2,1) << ", " << evectors(2,2) << ")" << std::endl;
 128     */
 129     evectors.GetVnlMatrix().scale_row(0, evalues[0]);
 130     evectors.GetVnlMatrix().scale_row(1, evalues[1]);
 131     evectors.GetVnlMatrix().scale_row(2, evalues[2]);
 132
 133     /*
 134       std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen values x Eigen vectors : "
 135
 136       << "(" << evectors(0,0) << ", " << evectors(0,1) << ", " << evectors(0,2) << "), "
 137       << "(" << evectors(1,0) << ", " << evectors(1,1) << ", " << evectors(1,2) << "), "
 138       << "(" << evectors(2,0) << ", " << evectors(2,1) << ", " << evectors(2,2) << ")" << std::endl;
 139     */
 140     vtkTensors->InsertTuple9(
 141       vtkPointID,
 142       evectors(0,0), evectors(0,1), evectors(0,2),
 143       evectors(1,0), evectors(1,1), evectors(1,2),
 144       evectors(2,0), evectors(2,1), evectors(2,2)
 145       );
 146   } // rof
 147
 148   polyData->SetPoints( points );
 149   vtkSmartPointer<vtkPointData> pointData = polyData->GetPointData();
 150   pointData->SetTensors( vtkTensors );
 151   return( 1 );
 152
 153   /**
 154    * VISUALISATION CODE BELOW
 155    */
 156
 157 //    vtkSmartPointer<vtkCubeSource> cubeSource = vtkSmartPointer<vtkCubeSource>::New();
 158 //    cubeSource->Update();
 159 //
 160 //    vtkSmartPointer<vtkTensorGlyph> tensorGlyph = vtkSmartPointer<vtkTensorGlyph>::New();
 161 //#if VTK_MAJOR_VERSION <= 5
 162 //    tensorGlyph->SetInput(polyData);
 163 //#else
 164 //    tensorGlyph->SetInputData(polyData);
 165 //#endif
 166 //    tensorGlyph->SetSourceConnection(cubeSource->GetOutputPort());
 167 //    tensorGlyph->ColorGlyphsOff();
 168 //    tensorGlyph->ThreeGlyphsOff();
 169 //    tensorGlyph->ExtractEigenvaluesOff();
 170 //    tensorGlyph->Update();
 171 //    .
 172 //    .
 173 //    .
 174
 175 }
 176
 177 // -------------------------------------------------------------------------
 178 template class cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< itk::Image< itk::SymmetricSecondRankTensor< float, 3 >, 3 > >;
 179 template class cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< itk::Image< itk::SymmetricSecondRankTensor< double, 3 >, 3 > >;
 180
 181 // eof - $RCSfile$