ANIMA  4.0
animaLocalPatchMeanDistanceImageFilter.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 <unsigned int> numPixels(numSamplesDatabase, 0);
46  std::vector <VectorType> meanVectors(numSamplesDatabase, patchMean);
47  std::vector <CovarianceType> varianceVector(numSamplesDatabase, baseVar);
48  OutputImageRegionType tmpBlockRegion;
49 
50  unsigned int numDistances = numSamplesDatabase * (numSamplesDatabase + 1) / 2 - numSamplesDatabase;
51 
52  InputImageIndexType curIndex;
53  OutputImageRegionType largestRegionOut = this->GetOutput(0)->GetLargestPossibleRegion();
54 
55  while (!maskIterator.IsAtEnd())
56  {
57  if (maskIterator.Get() == 0)
58  {
59  outMeanIterator.Set(0.0);
60  outStdIterator.Set(0.0);
61 
62  ++outMeanIterator;
63  ++outStdIterator;
64  ++maskIterator;
65  continue;
66  }
67 
68  curIndex = maskIterator.GetIndex();
69 
70  for (unsigned int i = 0;i < 3;++i)
71  {
72  tmpBlockRegion.SetIndex(i,std::max(0,(int)curIndex[i] - (int)m_PatchHalfSize));
73  tmpBlockRegion.SetSize(i,std::min((unsigned int)(largestRegionOut.GetSize()[i] - 1),(unsigned int)(curIndex[i] + m_PatchHalfSize)) - tmpBlockRegion.GetIndex(i) + 1);
74  }
75 
76  for (unsigned int i = 0;i < numSamplesDatabase;++i)
77  numPixels[i] = anima::computePatchMeanAndCovariance(this->GetInput(i),tmpBlockRegion,meanVectors[i],varianceVector[i]);
78 
79  double meanDist = 0;
80  double varDist = 0;
81  for (unsigned int i = 0;i < numSamplesDatabase;++i)
82  for (unsigned int j = i+1;j < numSamplesDatabase;++j)
83  {
84  double tmpDist = anima::VectorMeansTest(meanVectors[i], meanVectors[j], numPixels[i], numPixels[j],
85  varianceVector[i], varianceVector[j]);
86  meanDist += tmpDist;
87  varDist += tmpDist * tmpDist;
88  }
89 
90  varDist /= numDistances;
91  meanDist /= numDistances;
92  varDist -= meanDist * meanDist;
93  varDist *= numDistances / (numDistances - 1.0);
94 
95  outMeanIterator.Set(meanDist);
96  outStdIterator.Set(std::sqrt(varDist));
97 
98  this->IncrementNumberOfProcessedPoints();
99  ++outMeanIterator;
100  ++outStdIterator;
101  ++maskIterator;
102  }
103 }
104 
105 } // end 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 VectorMeansTest(itk::VariableLengthVector< T > &refPatchMean, itk::VariableLengthVector< T > &movingPatchMean, const unsigned int &refPatchNumElts, const unsigned int &movingPatchNumElts, vnl_matrix< T > &refPatchCov, vnl_matrix< T > &movingPatchCov)
Test if vector means are different (returns distance)
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE