+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#include <cpExtensions/Algorithms/SecondRankDiffusionTensorToPolyData.h>
+#include <itkImage.h>
+#include <itkSymmetricSecondRankTensor.h>
+#include <itkImageRegionConstIteratorWithIndex.h>
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
+#include <vtkPoints.h>
+#include <vtkPointData.h>
+#include <vtkSmartPointer.h>
+#include <vtkDoubleArray.h>
+
+// -------------------------------------------------------------------------
+template< class I >
+typename cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
+Self* cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
+New( )
+{
+ return( new Self );
+}
+
+// -------------------------------------------------------------------------
+template< class I >
+void cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
+SetInputData( const I* image )
+{
+ this->m_ITKImage = image;
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class I >
+cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
+SecondRankDiffusionTensorToPolyData( )
+ : vtkPolyDataAlgorithm( )
+{
+ this->SetNumberOfInputPorts( 0 );
+}
+
+// -------------------------------------------------------------------------
+template< class I >
+cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
+~SecondRankDiffusionTensorToPolyData( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class I >
+int cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
+RequestData(
+ vtkInformation* request,
+ vtkInformationVector** inputVector,
+ vtkInformationVector* outputVector
+ )
+{
+ if( this->m_ITKImage.IsNull( ) )
+ return( 0 );
+
+ // get the info objects
+ vtkInformation* outInfo = outputVector->GetInformationObject( 0 );
+ vtkPolyData* polyData = vtkPolyData::SafeDownCast(
+ outInfo->Get(vtkDataObject::DATA_OBJECT())
+ );
+
+ itk::ImageRegionConstIteratorWithIndex< I > tensorIt( this->m_ITKImage, this->m_ITKImage->GetRequestedRegion() );
+
+ vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+ vtkSmartPointer<vtkDoubleArray> vtkTensors = vtkSmartPointer<vtkDoubleArray>::New();
+ vtkTensors->SetNumberOfComponents(9);
+ for ( tensorIt.GoToBegin(); !tensorIt.IsAtEnd(); ++tensorIt )
+ {
+ typename I::IndexType tensorIndex = tensorIt.GetIndex();
+ typename I::PixelType tensor = tensorIt.Get();
+ typename I::PointType ptensor;
+
+ // TODO: std::cout << "\n Index [" << tensorIndex[0] << ", " << tensorIndex[1] << ", " << tensorIndex[2] << "]" << std::endl;
+ this->m_ITKImage->TransformIndexToPhysicalPoint( tensorIndex, ptensor );
+ vtkIdType vtkPointID = points->InsertNextPoint(ptensor[0], ptensor[1], ptensor[2]);
+
+ /**
+ * Definition included in cdifTypes.
+ * typedef itk::SymmetricSecondRankTensor< double > Diffusion2ndRankTensor;
+ */
+ //Diffusion2ndRankTensor Diffusion2ndRankTensor;
+
+
+ /**
+ * SYMMETRIC SECOND RANK TENSOR
+ *
+ * The upper-right triangle of the matrix:
+ *
+ * | 0 1 2 |
+ * | X 3 4 |
+ * | X X 5 |
+ *
+ * | xx xy xz |
+ * | X yy yz |
+ * | X X zz |
+ */
+ /*
+ Diffusion2ndRankTensor = tensor[0];
+ Diffusion2ndRankTensor = tensor[1];
+ Diffusion2ndRankTensor = tensor[2];
+ Diffusion2ndRankTensor = tensor[3];
+ Diffusion2ndRankTensor = tensor[4];
+ Diffusion2ndRankTensor = tensor[5];
+ */
+ typename I::PixelType::EigenValuesArrayType evalues;
+ typename I::PixelType::EigenVectorsMatrixType evectors;
+
+ /**
+ * Return an array containing EigenValues,
+ * and a matrix containing Eigenvectors.
+ */
+ tensor.ComputeEigenAnalysis(evalues, evectors);
+
+ /*
+ std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen values: " << evalues[0] << ", " << evalues[1] << ", " << evalues[2] << std::endl;
+
+ std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen vectors: "
+
+ << "(" << evectors(0,0) << ", " << evectors(0,1) << ", " << evectors(0,2) << "), "
+ << "(" << evectors(1,0) << ", " << evectors(1,1) << ", " << evectors(1,2) << "), "
+ << "(" << evectors(2,0) << ", " << evectors(2,1) << ", " << evectors(2,2) << ")" << std::endl;
+ */
+ evectors.GetVnlMatrix().scale_row(0, evalues[0]);
+ evectors.GetVnlMatrix().scale_row(1, evalues[1]);
+ evectors.GetVnlMatrix().scale_row(2, evalues[2]);
+
+ /*
+ std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen values x Eigen vectors : "
+
+ << "(" << evectors(0,0) << ", " << evectors(0,1) << ", " << evectors(0,2) << "), "
+ << "(" << evectors(1,0) << ", " << evectors(1,1) << ", " << evectors(1,2) << "), "
+ << "(" << evectors(2,0) << ", " << evectors(2,1) << ", " << evectors(2,2) << ")" << std::endl;
+ */
+ vtkTensors->InsertTuple9(
+ vtkPointID,
+ evectors(0,0), evectors(0,1), evectors(0,2),
+ evectors(1,0), evectors(1,1), evectors(1,2),
+ evectors(2,0), evectors(2,1), evectors(2,2)
+ );
+ } // rof
+
+ polyData->SetPoints( points );
+ vtkSmartPointer<vtkPointData> pointData = polyData->GetPointData();
+ pointData->SetTensors( vtkTensors );
+ return( 1 );
+
+ /**
+ * VISUALISATION CODE BELOW
+ */
+
+// vtkSmartPointer<vtkCubeSource> cubeSource = vtkSmartPointer<vtkCubeSource>::New();
+// cubeSource->Update();
+//
+// vtkSmartPointer<vtkTensorGlyph> tensorGlyph = vtkSmartPointer<vtkTensorGlyph>::New();
+//#if VTK_MAJOR_VERSION <= 5
+// tensorGlyph->SetInput(polyData);
+//#else
+// tensorGlyph->SetInputData(polyData);
+//#endif
+// tensorGlyph->SetSourceConnection(cubeSource->GetOutputPort());
+// tensorGlyph->ColorGlyphsOff();
+// tensorGlyph->ThreeGlyphsOff();
+// tensorGlyph->ExtractEigenvaluesOff();
+// tensorGlyph->Update();
+// .
+// .
+// .
+
+}
+
+// -------------------------------------------------------------------------
+template class cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< itk::Image< itk::SymmetricSecondRankTensor< float, 3 >, 3 > >;
+template class cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< itk::Image< itk::SymmetricSecondRankTensor< double, 3 >, 3 > >;
+
+// eof - $RCSfile$