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> 20 template <
unsigned int ImageDimension>
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));
34 template <
unsigned int ImageDimension>
35 itk::DataObject::Pointer
38 return (OutputImageType::New()).GetPointer();
41 template <
unsigned int ImageDimension>
46 itk::ImageRegionConstIterator<TensorImageType> tensorIterator;
47 itk::ImageRegionIterator<ADCImageType> adcIterator;
50 typename InputImageType::ConstPointer tensorImage = this->GetInput();
51 typename ADCImageType::Pointer adcImage = this->GetOutput(0);
54 tensorRegionForThread.SetIndex(outputRegionForThread.GetIndex());
55 tensorRegionForThread.SetSize(outputRegionForThread.GetSize());
57 tensorIterator = itk::ImageRegionConstIterator<TensorImageType>(tensorImage, tensorRegionForThread);
58 adcIterator = itk::ImageRegionIterator<OutputImageType>(adcImage, outputRegionForThread);
59 tensorIterator.GoToBegin();
60 adcIterator.GoToBegin();
62 itk::ImageRegionIterator<OutputImageType> faIterator;
63 typename FAImageType::Pointer faImage = this->GetOutput(1);
64 faIterator = itk::ImageRegionIterator<OutputImageType>(faImage, outputRegionForThread);
65 faIterator.GoToBegin();
67 itk::ImageRegionIterator<OutputImageType> axialIterator;
68 typename FAImageType::Pointer axialImage = this->GetOutput(2);
69 axialIterator = itk::ImageRegionIterator<OutputImageType>(axialImage, outputRegionForThread);
70 axialIterator.GoToBegin();
72 itk::ImageRegionIterator<OutputImageType> radialIterator;
73 typename FAImageType::Pointer radialImage = this->GetOutput(3);
74 radialIterator = itk::ImageRegionIterator<OutputImageType>(radialImage, outputRegionForThread);
75 radialIterator.GoToBegin();
77 typedef vnl_matrix<double> TensorSymMatrixType;
78 typedef itk::Vector<double, 3> EigenVectorType;
80 EigenVectorType eigenValue;
81 TensorSymMatrixType tensorSymMatrix(3,3);
83 itk::SymmetricEigenAnalysis<TensorSymMatrixType, EigenVectorType> eigenAnalysis;
84 eigenAnalysis.SetDimension(3);
86 while (!tensorIterator.IsAtEnd())
89 double a(tensor[0]), c(tensor[2]), f(tensor[5]), ADC(0);
96 eigenAnalysis.ComputeEigenValues(tensorSymMatrix, eigenValue);
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);
105 fa = std::sqrt(0.5) * (num / den);
109 axialIterator.Set(l1);
110 radialIterator.Set((l2+l3) / 2.0);
112 this->IncrementNumberOfProcessedPoints();
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)
TensorImageType::PixelType TensorVectorType