ANIMA  4.0
animaVectorImagePatchStatistics.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionConstIteratorWithIndex.h>
5 #include <itkSymmetricEigenAnalysis.h>
6 
7 namespace anima
8 {
9 
10 template <class T1, class T2, unsigned int Dimension>
11 unsigned int computePatchMeanAndCovariance(const itk::VectorImage <T1, Dimension> *inputImage, const itk::ImageRegion<Dimension> &patchRegion,
12  itk::VariableLengthVector <T2> &patchMean, vnl_matrix <T2> &patchCov)
13 {
14  unsigned int ndim = inputImage->GetNumberOfComponentsPerPixel();
15  if (patchMean.GetSize() != ndim)
16  patchMean.SetSize(ndim);
17  patchMean.Fill(0);
18 
19  patchCov.set_size(ndim,ndim);
20  patchCov.fill(0);
21 
22  typedef itk::VectorImage <T1, Dimension> VectorImageType;
23  typedef itk::ImageRegionConstIteratorWithIndex< VectorImageType > InIteratorType;
24 
25  InIteratorType imageIt(inputImage,patchRegion);
26 
27  unsigned int numPixels = 0;
28  while (!imageIt.IsAtEnd())
29  {
30  patchMean += imageIt.Get();
31  ++numPixels;
32  ++imageIt;
33  }
34 
35  patchMean /= numPixels;
36  imageIt.GoToBegin();
37 
38  itk::VariableLengthVector <T2> patchDiff(ndim);
39  patchDiff.Fill(0);
40 
41  while (!imageIt.IsAtEnd())
42  {
43  for (unsigned int i = 0;i < ndim;++i)
44  patchDiff[i] = patchMean[i] - imageIt.Get()[i];
45 
46  for (unsigned int i = 0;i < ndim;++i)
47  for (unsigned int j = i;j < ndim;++j)
48  {
49  double sqDiff = patchDiff[i] * patchDiff[j];
50  patchCov(i,j) += sqDiff;
51  }
52 
53  ++imageIt;
54  }
55 
56  patchCov /= (numPixels - 1.0);
57  for (unsigned int i = 0;i < ndim;++i)
58  for (unsigned int j = i+1;j < ndim;++j)
59  patchCov(j,i) = patchCov(i,j);
60 
61  return numPixels;
62 }
63 
64 template <class T1, class T2, unsigned int Dimension>
65 void computeAverageLocalCovariance(vnl_matrix <T2> &resVariance, itk::VectorImage <T1, Dimension> *inputImage,
66  itk::Image<unsigned char, Dimension> *maskImage, const itk::ImageRegion<Dimension> &averagingRegion,
67  int localNeighborhood)
68 {
69  //Generalization of what Pierrick Coupe was using in his NL-Means on scalar images
70 
71  typedef itk::VectorImage <T1, Dimension> VectorImageType;
72  typedef itk::ImageRegionConstIteratorWithIndex< VectorImageType > InIteratorType;
73  typedef itk::ImageRegionConstIterator < itk::Image<unsigned char, Dimension> > MaskRegionIteratorType;
74 
75  typedef itk::ImageRegion<Dimension> ImageRegionType;
76  typedef typename VectorImageType::IndexType ImageIndexType;
77  typedef itk::VariableLengthVector <T2> VectorType;
78 
79  ImageRegionType largestRegion = inputImage->GetLargestPossibleRegion();
80 
81  MaskRegionIteratorType maskIterator (maskImage, averagingRegion);
82  InIteratorType dataIterator (inputImage, averagingRegion);
83 
84  unsigned int ndim = inputImage->GetNumberOfComponentsPerPixel();
85  resVariance.set_size(ndim,ndim);
86  resVariance.fill(0);
87 
88  VectorType averageLocalSignal(ndim), diffSignal(ndim);
89  VectorType baseSignal(ndim);
90  ImageIndexType baseIndex;
91 
92  unsigned int numEstimations = 0;
93 
94  while (!maskIterator.IsAtEnd())
95  {
96  if (maskIterator.Get() == 0)
97  {
98  ++maskIterator;
99  ++dataIterator;
100  continue;
101  }
102 
103  baseSignal = dataIterator.Get();
104  baseIndex = dataIterator.GetIndex();
105  averageLocalSignal = - baseSignal;
106 
107  ImageRegionType diffRegion;
108  for (unsigned int i = 0;i < Dimension;++i)
109  {
110  diffRegion.SetIndex(i,std::max(0,(int)baseIndex[i] - localNeighborhood));
111  diffRegion.SetSize(i,std::min((unsigned int)(largestRegion.GetSize()[i] - 1),(unsigned int)(baseIndex[i] + localNeighborhood)) - diffRegion.GetIndex(i) + 1);
112  }
113 
114  unsigned int numLocalPixels = diffRegion.GetSize()[0];
115  for (unsigned int i = 1;i < Dimension;++i)
116  numLocalPixels *= diffRegion.GetSize()[i];
117 
118  InIteratorType localIterator(inputImage,diffRegion);
119 
120  while (!localIterator.IsAtEnd())
121  {
122  averageLocalSignal += localIterator.Get();
123  ++localIterator;
124  }
125 
126  averageLocalSignal /= numLocalPixels;
127  diffSignal = sqrt(numLocalPixels / (numLocalPixels + 1.0)) * (baseSignal - averageLocalSignal);
128 
129  for (unsigned int i = 0;i < ndim;++i)
130  for (unsigned int j = i;j < ndim;++j)
131  {
132  double tmpVal = diffSignal[i] * diffSignal[j];
133  resVariance(i,j) += tmpVal;
134  if (i != j)
135  resVariance(j,i) += tmpVal;
136  }
137 
138  ++numEstimations;
139  ++maskIterator;
140  ++dataIterator;
141  }
142 
143  // Now divide by number of estimations and compute average variance
144  resVariance /= numEstimations;
145 }
146 
147 template <class T> double VectorCovarianceTest(vnl_matrix <T> &logRefPatchCov, vnl_matrix <T> &movingPatchCov)
148 {
149  // Here comes a test on the distance between variances
150  unsigned int ndim = logRefPatchCov.rows();
151 
152  vnl_matrix <double> eVec(ndim,ndim);
153  vnl_diag_matrix <double> eVals(ndim);
154  itk::SymmetricEigenAnalysis <vnl_matrix <T>, vnl_diag_matrix <double>, vnl_matrix <double> > EigenAnalysis(ndim);
155 
156  EigenAnalysis.ComputeEigenValuesAndVectors(movingPatchCov, eVals, eVec);
157 
158  for (unsigned int i = 0;i < ndim;++i)
159  eVals[i] = log(eVals[i]);
160 
161  vnl_matrix <double> logMoving = eVec.transpose() * eVals * eVec;
162 
163  double varsDist = 0;
164  for (unsigned int i = 0;i < ndim;++i)
165  for (unsigned int j = i;j < ndim;++j)
166  {
167  if (i == j)
168  varsDist += (logRefPatchCov(i,j) - logMoving(i,j)) * (logRefPatchCov(i,j) - logMoving(i,j));
169  else
170  varsDist += 2.0 * (logRefPatchCov(i,j) - logMoving(i,j)) * (logRefPatchCov(i,j) - logMoving(i,j));
171  }
172 
173  varsDist = std::sqrt(varsDist);
174 
175  return varsDist;
176 }
177 
178 template <class T>
179 double VectorMeansTest(itk::VariableLengthVector <T> &refPatchMean, itk::VariableLengthVector <T> &movingPatchMean,
180  const unsigned int &refPatchNumElts, const unsigned int &movingPatchNumElts,
181  vnl_matrix <T> &refPatchCov, vnl_matrix <T> &movingPatchCov)
182 {
183  // Here comes a test on the distance between variances
184  unsigned int ndim = refPatchMean.GetNumberOfElements();
185 
186  itk::VariableLengthVector <T> meansDiff;
187 
188  vnl_matrix <double> S_pooled(ndim,ndim), S_pooled_inv(ndim,ndim);
189 
190  S_pooled = refPatchCov * refPatchNumElts + movingPatchCov * movingPatchNumElts;
191  S_pooled /= (refPatchNumElts + movingPatchNumElts - 2.0);
192 
193  vnl_matrix_inverse <double> spoolInverter(S_pooled);
194  S_pooled_inv = spoolInverter.inverse();
195 
196  meansDiff = refPatchMean - movingPatchMean;
197  double distMeans = 0;
198 
199  for (unsigned int j = 0;j < ndim;++j)
200  for (unsigned int k = j;k < ndim;++k)
201  {
202  if (j != k)
203  distMeans += 2.0 * S_pooled_inv(j,k) * meansDiff[j] * meansDiff[k];
204  else
205  distMeans += S_pooled_inv(j,k) * meansDiff[j] * meansDiff[k];
206  }
207 
208  distMeans = std::sqrt(distMeans);
209 
210  return distMeans;
211 }
212 
213 } // 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)
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 computeAverageLocalCovariance(vnl_matrix< T2 > &resVariance, itk::VectorImage< T1, Dimension > *inputImage, itk::Image< unsigned char, Dimension > *maskImage, const itk::ImageRegion< Dimension > &averagingRegion, int localNeighborhood)
Noise estimation for a patch of a vector image.