ANIMA  4.0
animaDTIScalarMapsImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
2 
4 
5 #include <itkImageRegionConstIterator.h>
6 #include <itkNeighborhoodInnerProduct.h>
7 #include <itkImageRegionIterator.h>
8 #include <itkNeighborhoodAlgorithm.h>
9 #include <itkZeroFluxNeumannBoundaryCondition.h>
10 #include <itkOffset.h>
11 #include <itkProgressReporter.h>
12 #include <itkMatrix.h>
13 #include <itkSymmetricEigenAnalysis.h>
14 
15 #include <animaBaseTensorTools.h>
16 
17 namespace anima
18 {
19 
20 template < unsigned int ImageDimension>
22  Superclass()
23 {
24  this->SetNumberOfRequiredOutputs(4);
25  this->SetNthOutput(0, this->MakeOutput(0));
26  this->SetNthOutput(1, this->MakeOutput(1));
27  this->SetNthOutput(2, this->MakeOutput(2));
28  this->SetNthOutput(3, this->MakeOutput(3));
29 }
30 
34 template < unsigned int ImageDimension>
35 itk::DataObject::Pointer
36 DTIScalarMapsImageFilter< ImageDimension >::MakeOutput(itk::ProcessObject::DataObjectPointerArraySizeType idx)
37 {
38  return (OutputImageType::New()).GetPointer();
39 }
40 
41 template < unsigned int ImageDimension>
42 void
45 {
46  itk::ImageRegionConstIterator<TensorImageType> tensorIterator;
47  itk::ImageRegionIterator<ADCImageType> adcIterator;
48 
49  // Allocate output
50  typename InputImageType::ConstPointer tensorImage = this->GetInput();
51  typename ADCImageType::Pointer adcImage = this->GetOutput(0);
52 
53  TensorImageRegionType tensorRegionForThread;
54  tensorRegionForThread.SetIndex(outputRegionForThread.GetIndex());
55  tensorRegionForThread.SetSize(outputRegionForThread.GetSize());
56 
57  tensorIterator = itk::ImageRegionConstIterator<TensorImageType>(tensorImage, tensorRegionForThread);
58  adcIterator = itk::ImageRegionIterator<OutputImageType>(adcImage, outputRegionForThread);
59  tensorIterator.GoToBegin();
60  adcIterator.GoToBegin();
61 
62  itk::ImageRegionIterator<OutputImageType> faIterator;
63  typename FAImageType::Pointer faImage = this->GetOutput(1);
64  faIterator = itk::ImageRegionIterator<OutputImageType>(faImage, outputRegionForThread);
65  faIterator.GoToBegin();
66 
67  itk::ImageRegionIterator<OutputImageType> axialIterator;
68  typename FAImageType::Pointer axialImage = this->GetOutput(2);
69  axialIterator = itk::ImageRegionIterator<OutputImageType>(axialImage, outputRegionForThread);
70  axialIterator.GoToBegin();
71 
72  itk::ImageRegionIterator<OutputImageType> radialIterator;
73  typename FAImageType::Pointer radialImage = this->GetOutput(3);
74  radialIterator = itk::ImageRegionIterator<OutputImageType>(radialImage, outputRegionForThread);
75  radialIterator.GoToBegin();
76 
77  typedef vnl_matrix<double> TensorSymMatrixType;
78  typedef itk::Vector<double, 3> EigenVectorType;
79 
80  EigenVectorType eigenValue;
81  TensorSymMatrixType tensorSymMatrix(3,3);
82 
83  itk::SymmetricEigenAnalysis<TensorSymMatrixType, EigenVectorType> eigenAnalysis;
84  eigenAnalysis.SetDimension(3);
85 
86  while (!tensorIterator.IsAtEnd())
87  {
88  TensorVectorType tensor = tensorIterator.Get();
89  double a(tensor[0]), c(tensor[2]), f(tensor[5]), ADC(0);
90  ADC = (a+c+f) / 3.0;
91 
92  adcIterator.Set(ADC);
93 
94  anima::GetTensorFromVectorRepresentation(tensor, tensorSymMatrix,3,false);
95 
96  eigenAnalysis.ComputeEigenValues(tensorSymMatrix, eigenValue);
97 
98  double l1(eigenValue[2]), l2(eigenValue[1]), l3(eigenValue[0]), fa(1);
99  double num = std::sqrt ((l1 -l2) * (l1 -l2) + (l2 -l3) * (l2 -l3) + (l3 - l1) * (l3 - l1));
100  double den = std::sqrt (l1*l1 + l2*l2 + l3*l3);
101 
102  if (den == 0)
103  fa = 0;
104  else
105  fa = std::sqrt(0.5) * (num / den);
106 
107  faIterator.Set(fa);
108 
109  axialIterator.Set(l1);
110  radialIterator.Set((l2+l3) / 2.0);
111 
112  this->IncrementNumberOfProcessedPoints();
113  ++tensorIterator;
114  ++adcIterator;
115  ++faIterator;
116  ++axialIterator;
117  ++radialIterator;
118  }
119 }
120 
121 } // end of namespace anima
Implements a class to handle thread number in a dynamic way for multithreaded methods needing thread ...
OutputImageType::RegionType OutputImageRegionType
itk::DataObject::Pointer MakeOutput(itk::ProcessObject::DataObjectPointerArraySizeType idx) ITK_OVERRIDE
Applies an variance filter to an image.
TensorImageType::RegionType TensorImageRegionType
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
void GetTensorFromVectorRepresentation(const itk::VariableLengthVector< T1 > &vector, vnl_matrix< T2 > &tensor, unsigned int tensDim=0, bool scale=false)