ANIMA  4.0
animaNLMeansPatientToGroupComparisonImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionIterator.h>
5 #include <itkImageRegionIteratorWithIndex.h>
6 #include <itkImageRegionConstIterator.h>
7 #include <itkImageRegionConstIteratorWithIndex.h>
8 #include <itkTimeProbe.h>
9 
11 #include <boost/math/distributions/fisher_f.hpp>
12 
13 namespace anima
14 {
15 
16 template <class PixelScalarType>
17 void
20 {
21  Superclass::BeforeThreadedGenerateData();
22 
23  // Checking consistency of the data and parameters
24 
25  unsigned int nbInputs = this->GetNumberOfIndexedInputs();
26  if (nbInputs <= 0)
27  itkExceptionMacro("Error: No inputs available... Exiting...");
28 }
29 
30 template <class PixelScalarType>
31 void
34 {
35  typedef itk::ImageRegionConstIteratorWithIndex< InputImageType > InIteratorType;
36  typedef itk::ImageRegionConstIterator < InputImageType > DatabaseIteratorType;
37  typedef itk::ImageRegionIterator < OutputImageType > OutRegionIteratorType;
38  typedef itk::ImageRegionConstIterator < MaskImageType > MaskRegionIteratorType;
39 
40  OutRegionIteratorType outIterator(this->GetOutput(0), outputRegionForThread);
41  OutRegionIteratorType outScoreIterator(this->GetOutput(1), outputRegionForThread);
42  OutRegionIteratorType outNumPatchesIterator(this->GetOutput(2), outputRegionForThread);
43 
44  MaskRegionIteratorType maskIterator (this->GetComputationMask(), outputRegionForThread);
45 
46  unsigned int numSamplesDatabase = m_DatabaseImages.size();
47  std::vector <DatabaseIteratorType> databaseIterators(numSamplesDatabase);
48  for (unsigned int k = 0;k < numSamplesDatabase;++k)
49  databaseIterators[k] = DatabaseIteratorType(m_DatabaseImages[k],outputRegionForThread);
50 
51  InputImageType *input = const_cast<InputImageType *> (this->GetInput());
52  InIteratorType patientIterator(input, outputRegionForThread);
53  VectorType patientSample;
54 
55  std::vector <VectorType> databaseSamples;
56  std::vector <double> databaseWeights;
57 
58  int maxAbsDisp = (int)floor((double)(m_SearchNeighborhood / m_SearchStepSize)) * m_SearchStepSize;
59 
61 
62  PatchSearcherType patchSearcher;
63  patchSearcher.SetPatchHalfSize(m_PatchHalfSize);
64  patchSearcher.SetSearchStepSize(m_SearchStepSize);
65  patchSearcher.SetMaxAbsDisp(maxAbsDisp);
66  patchSearcher.SetInputImage(input);
67  patchSearcher.SetBetaParameter(m_BetaParameter);
68  patchSearcher.SetWeightThreshold(m_WeightThreshold);
69  patchSearcher.SetMeanThreshold(m_MeanThreshold);
70  patchSearcher.SetVarianceThreshold(m_VarianceThreshold);
71  patchSearcher.SetDatabaseCovarianceDistanceAverage(m_DatabaseCovarianceDistanceAverage);
72  patchSearcher.SetDatabaseCovarianceDistanceStd(m_DatabaseCovarianceDistanceStd);
73  patchSearcher.SetDatabaseMeanDistanceAverage(m_DatabaseMeanDistanceAverage);
74  patchSearcher.SetDatabaseMeanDistanceStd(m_DatabaseMeanDistanceStd);
75  patchSearcher.SetDataMask(this->GetComputationMask());
76 
77  for (unsigned int k = 0;k < numSamplesDatabase;++k)
78  patchSearcher.AddComparisonImage(m_DatabaseImages[k]);
79 
80  while (!outIterator.IsAtEnd())
81  {
82  if (maskIterator.Get() == 0)
83  {
84  outIterator.Set(0.0);
85  outScoreIterator.Set(0.0);
86  outNumPatchesIterator.Set(0.0);
87  ++outIterator;
88  ++outScoreIterator;
89  ++maskIterator;
90  ++patientIterator;
91  ++outNumPatchesIterator;
92 
93  for (unsigned int k = 0;k < numSamplesDatabase;++k)
94  ++databaseIterators[k];
95 
96  continue;
97  }
98 
99  patchSearcher.UpdateAtPosition(patientIterator.GetIndex());
100 
101  databaseSamples = patchSearcher.GetDatabaseSamples();
102  databaseWeights = patchSearcher.GetDatabaseWeights();
103 
104  // Add center pixels
105  for (unsigned int k = 0;k < numSamplesDatabase;++k)
106  {
107  databaseSamples.push_back(databaseIterators[k].Get());
108  databaseWeights.push_back(1.0);
109  }
110 
111  double diffScore = 0;
112  patientSample = patientIterator.Get();
113  double pValue = 0;
114 
115  try
116  {
117  if (databaseSamples.size() > 1)
118  pValue = this->ComputeWeightedDistanceScore(patientSample,databaseWeights,databaseSamples,diffScore);
119  }
120  catch (...)
121  {
122  itkExceptionMacro("dying in final calculations");
123  }
124 
125  outIterator.Set(pValue);
126  outScoreIterator.Set(diffScore);
127  outNumPatchesIterator.Set(databaseSamples.size());
128 
129  this->IncrementNumberOfProcessedPoints();
130  ++outIterator;
131  ++outScoreIterator;
132  ++outNumPatchesIterator;
133  ++maskIterator;
134  ++patientIterator;
135 
136  for (unsigned int k = 0;k < numSamplesDatabase;++k)
137  ++databaseIterators[k];
138  }
139 }
140 
141 template <class PixelScalarType>
142 double
144 ::ComputeWeightedDistanceScore(itk::VariableLengthVector <double> &patientSample, std::vector < double > &databaseWeights,
145  std::vector < itk::VariableLengthVector <double> > &databaseSamples, double &diffScore)
146 {
147  diffScore = 0;
148 
149  unsigned int nsamples = databaseSamples.size();
150 
151  unsigned int ndim = patientSample.GetSize();
152 
153  if (nsamples <= ndim)
154  return 0;
155 
156  // Compute patient sample diff score
157  double sumWeights = 0;
158  itk::VariableLengthVector <double> databaseMean(ndim);
159  databaseMean.Fill(0);
160 
161  for (unsigned int i = 0;i < nsamples;++i)
162  sumWeights += databaseWeights[i];
163 
164  double sqrSumWeights = 0;
165  for (unsigned int i = 0;i < nsamples;++i)
166  {
167  databaseWeights[i] /= sumWeights;
168  sqrSumWeights += databaseWeights[i] * databaseWeights[i];
169  }
170 
171  for (unsigned int i = 0;i < nsamples;++i)
172  {
173  databaseMean += databaseWeights[i] * databaseSamples[i];
174  }
175 
176  vnl_matrix <double> covMatrixDatabase(ndim,ndim);
177  covMatrixDatabase.fill(0);
178 
179  for (unsigned int i = 0;i < nsamples;++i)
180  {
181  for (unsigned int j = 0;j < ndim;++j)
182  for (unsigned int k = j;k < ndim;++k)
183  covMatrixDatabase(j,k) += databaseWeights[i]*(databaseMean[j] - databaseSamples[i][j])*(databaseMean[k] - databaseSamples[i][k]);
184  }
185 
186  for (unsigned int j = 0;j < ndim;++j)
187  for (unsigned int k = j;k < ndim;++k)
188  {
189  covMatrixDatabase(j,k) *= 1.0 / (1.0 - sqrSumWeights);
190  if (j != k)
191  covMatrixDatabase(k,j) = covMatrixDatabase(j,k);
192  }
193 
194  vnl_matrix_inverse <double> invCovMatrixDatabase(covMatrixDatabase);
195  vnl_matrix <double> covInvDatabase = invCovMatrixDatabase.inverse();
196 
197  diffScore = 0;
198  for (unsigned int j = 0;j < ndim;++j)
199  for (unsigned int k = j;k < ndim;++k)
200  {
201  double factor = 2.0;
202  if (j == k)
203  factor = 1;
204 
205  diffScore += factor*(patientSample[j] - databaseMean[j])*(patientSample[k] - databaseMean[k])*covInvDatabase(j,k);
206  }
207 
208  double testScore = nsamples * (nsamples - ndim) * diffScore / ((nsamples *nsamples - 1) * ndim);
209  boost::math::fisher_f_distribution <double> fisherDist(ndim,nsamples - ndim);
210  double resVal = 1.0 - boost::math::cdf(fisherDist, testScore);
211 
212  diffScore = std::sqrt(diffScore);
213 
214  return resVal;
215 }
216 
217 } //end namespace anima
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE