ANIMA  4.0
animaLocalPatchCovarianceDistanceImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionIteratorWithIndex.h>
5 #include <itkSymmetricEigenAnalysis.h>
6 
8 
9 namespace anima
10 {
11 
12 template <class PixelScalarType>
13 void
16 {
17  Superclass::BeforeThreadedGenerateData();
18 
19  // Checking consistency of the data and parameters
20 
21  unsigned int nbInputs = this->GetNumberOfIndexedInputs();
22  if (nbInputs <= 1)
23  itkExceptionMacro("Error: Not enough inputs available... Exiting...");
24 }
25 
26 template <class PixelScalarType>
27 void
30 {
31  typedef itk::ImageRegionIteratorWithIndex< OutputImageType > OutRegionIteratorType;
32  typedef itk::ImageRegionConstIteratorWithIndex < MaskImageType > MaskRegionIteratorType;
33 
34  OutRegionIteratorType outMeanIterator(this->GetOutput(0), outputRegionForThread);
35  OutRegionIteratorType outStdIterator(this->GetOutput(1), outputRegionForThread);
36 
37  MaskRegionIteratorType maskIterator (this->GetComputationMask(), outputRegionForThread);
38 
39  unsigned int numSamplesDatabase = this->GetNumberOfIndexedInputs();
40 
41  unsigned int ndim = this->GetInput(0)->GetNumberOfComponentsPerPixel();
42  CovarianceType baseVar(ndim,ndim,0);
43  VectorType patchMean(ndim);
44 
45  std::vector <CovarianceType> varianceVector(numSamplesDatabase, baseVar);
46  std::vector <CovarianceType> logVarianceVector(numSamplesDatabase, baseVar);
47  OutputImageRegionType tmpBlockRegion;
48 
49  unsigned int numDistances = numSamplesDatabase * (numSamplesDatabase + 1) / 2 - numSamplesDatabase;
50 
51  CovarianceType eVec(ndim,ndim);
52  vnl_diag_matrix <double> eVals(ndim);
53  itk::SymmetricEigenAnalysis <vnl_matrix <double>, vnl_diag_matrix <double>, vnl_matrix <double> > EigenAnalysis(ndim);
54 
55  InputImageIndexType curIndex;
56  OutputImageRegionType largestRegionOut = this->GetOutput(0)->GetLargestPossibleRegion();
57 
58  while (!maskIterator.IsAtEnd())
59  {
60  if (maskIterator.Get() == 0)
61  {
62  outMeanIterator.Set(0.0);
63  outStdIterator.Set(0.0);
64 
65  ++outMeanIterator;
66  ++outStdIterator;
67  ++maskIterator;
68  continue;
69  }
70 
71  curIndex = maskIterator.GetIndex();
72 
73  for (unsigned int i = 0;i < 3;++i)
74  {
75  tmpBlockRegion.SetIndex(i,std::max(0,(int)curIndex[i] - (int)m_PatchHalfSize));
76  tmpBlockRegion.SetSize(i,std::min((unsigned int)(largestRegionOut.GetSize()[i] - 1),(unsigned int)(curIndex[i] + m_PatchHalfSize)) - tmpBlockRegion.GetIndex(i) + 1);
77  }
78 
79  for (unsigned int i = 0;i < numSamplesDatabase;++i)
80  {
81  anima::computePatchMeanAndCovariance(this->GetInput(i),tmpBlockRegion,patchMean,varianceVector[i]);
82  EigenAnalysis.ComputeEigenValuesAndVectors(varianceVector[i], eVals, eVec);
83 
84  for (unsigned int j = 0;j < ndim;++j)
85  eVals[j] = log(eVals[j]);
86 
87  logVarianceVector[i] = eVec.transpose() * eVals * eVec;
88  }
89 
90  double meanDist = 0;
91  double varDist = 0;
92  for (unsigned int i = 0;i < numSamplesDatabase;++i)
93  for (unsigned int j = i+1;j < numSamplesDatabase;++j)
94  {
95  double tmpDist = anima::VectorCovarianceTest(logVarianceVector[i], varianceVector[j]);
96  meanDist += tmpDist;
97  varDist += tmpDist * tmpDist;
98  }
99 
100  varDist /= numDistances;
101  meanDist /= numDistances;
102  varDist -= meanDist * meanDist;
103  varDist *= numDistances / (numDistances - 1.0);
104 
105  outMeanIterator.Set(meanDist);
106  outStdIterator.Set(sqrt(varDist));
107 
108  this->IncrementNumberOfProcessedPoints();
109  ++outMeanIterator;
110  ++outStdIterator;
111  ++maskIterator;
112  }
113 }
114 
115 } // end of namespace anima
unsigned int computePatchMeanAndCovariance(const itk::VectorImage< T1, Dimension > *inputImage, const itk::ImageRegion< Dimension > &patchRegion, itk::VariableLengthVector< T2 > &patchMean, vnl_matrix< T2 > &patchCov)
Computes the average and covariance matrix from a patch of a vector image.
double VectorCovarianceTest(vnl_matrix< T > &logRefPatchCov, vnl_matrix< T > &movingPatchCov)
Test if covariance matrices are different (returns distance)
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE