1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #include <cpExtensions/Algorithms/SecondRankDiffusionTensorToPolyData.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>
16 // -------------------------------------------------------------------------
18 typename cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
19 Self* cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
25 // -------------------------------------------------------------------------
27 void cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
28 SetInputData( const I* image )
30 this->m_ITKImage = image;
34 // -------------------------------------------------------------------------
36 cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
37 SecondRankDiffusionTensorToPolyData( )
38 : vtkPolyDataAlgorithm( )
40 this->SetNumberOfInputPorts( 0 );
43 // -------------------------------------------------------------------------
45 cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
46 ~SecondRankDiffusionTensorToPolyData( )
50 // -------------------------------------------------------------------------
52 int cpExtensions::Algorithms::SecondRankDiffusionTensorToPolyData< I >::
54 vtkInformation* request,
55 vtkInformationVector** inputVector,
56 vtkInformationVector* outputVector
59 if( this->m_ITKImage.IsNull( ) )
62 // get the info objects
63 vtkInformation* outInfo = outputVector->GetInformationObject( 0 );
64 vtkPolyData* polyData = vtkPolyData::SafeDownCast(
65 outInfo->Get(vtkDataObject::DATA_OBJECT())
68 itk::ImageRegionConstIteratorWithIndex< I > tensorIt( this->m_ITKImage, this->m_ITKImage->GetRequestedRegion() );
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 )
75 typename I::IndexType tensorIndex = tensorIt.GetIndex();
76 typename I::PixelType tensor = tensorIt.Get();
77 typename I::PointType ptensor;
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]);
84 * Definition included in cdifTypes.
85 * typedef itk::SymmetricSecondRankTensor< double > Diffusion2ndRankTensor;
87 //Diffusion2ndRankTensor Diffusion2ndRankTensor;
91 * SYMMETRIC SECOND RANK TENSOR
93 * The upper-right triangle of the matrix:
104 Diffusion2ndRankTensor = tensor[0];
105 Diffusion2ndRankTensor = tensor[1];
106 Diffusion2ndRankTensor = tensor[2];
107 Diffusion2ndRankTensor = tensor[3];
108 Diffusion2ndRankTensor = tensor[4];
109 Diffusion2ndRankTensor = tensor[5];
111 typename I::PixelType::EigenValuesArrayType evalues;
112 typename I::PixelType::EigenVectorsMatrixType evectors;
115 * Return an array containing EigenValues,
116 * and a matrix containing Eigenvectors.
118 tensor.ComputeEigenAnalysis(evalues, evectors);
121 std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen values: " << evalues[0] << ", " << evalues[1] << ", " << evalues[2] << std::endl;
123 std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen vectors: "
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;
129 evectors.GetVnlMatrix().scale_row(0, evalues[0]);
130 evectors.GetVnlMatrix().scale_row(1, evalues[1]);
131 evectors.GetVnlMatrix().scale_row(2, evalues[2]);
134 std::cout << "[C-DIFFUSION - NRRD Tensors Example] Eigen values x Eigen vectors : "
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;
140 vtkTensors->InsertTuple9(
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)
148 polyData->SetPoints( points );
149 vtkSmartPointer<vtkPointData> pointData = polyData->GetPointData();
150 pointData->SetTensors( vtkTensors );
154 * VISUALISATION CODE BELOW
157 // vtkSmartPointer<vtkCubeSource> cubeSource = vtkSmartPointer<vtkCubeSource>::New();
158 // cubeSource->Update();
160 // vtkSmartPointer<vtkTensorGlyph> tensorGlyph = vtkSmartPointer<vtkTensorGlyph>::New();
161 //#if VTK_MAJOR_VERSION <= 5
162 // tensorGlyph->SetInput(polyData);
164 // tensorGlyph->SetInputData(polyData);
166 // tensorGlyph->SetSourceConnection(cubeSource->GetOutputPort());
167 // tensorGlyph->ColorGlyphsOff();
168 // tensorGlyph->ThreeGlyphsOff();
169 // tensorGlyph->ExtractEigenvaluesOff();
170 // tensorGlyph->Update();
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 > >;