ANIMA  4.0
animaGaussianEMEstimator.hxx
Go to the documentation of this file.
2 
3 namespace anima
4 {
5 
6 template <typename TInputImage, typename TMaskImage>
8 {
9  this->SetNthInput(0, const_cast<TMaskImage*>(mask));
10 }
11 
12 template <typename TInputImage, typename TMaskImage>
14 {
15  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
16  m_IndexImage1=m_NbInputs;
17  m_NbInputs++;
18 }
19 
20 template <typename TInputImage, typename TMaskImage>
22 {
23  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
24  m_IndexImage2=m_NbInputs;
25  m_NbInputs++;
26 }
27 
28 template <typename TInputImage, typename TMaskImage>
30 {
31  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
32  m_IndexImage3=m_NbInputs;
33  m_NbInputs++;
34 }
35 
36 template <typename TInputImage, typename TMaskImage>
38 {
39  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
40  m_IndexImage4=m_NbInputs;
41  m_NbInputs++;
42 }
43 
44 template <typename TInputImage, typename TMaskImage>
46 {
47  this->SetNthInput(m_NbInputs, const_cast<TInputImage*>(image));
48  m_IndexImage5=m_NbInputs;
49  m_NbInputs++;
50 }
51 
52 template <typename TInputImage, typename TMaskImage>
53 typename TMaskImage::ConstPointer GaussianEMEstimator<TInputImage,TMaskImage>::GetMask()
54 {
55  return static_cast< const TMaskImage * >
56  ( this->itk::ProcessObject::GetInput(0) );
57 }
58 
59 template <typename TInputImage, typename TMaskImage>
61 {
62  return static_cast< const TInputImage * >
63  ( this->itk::ProcessObject::GetInput(m_IndexImage1) );
64 }
65 
66 template <typename TInputImage, typename TMaskImage>
68 {
69  return static_cast< const TInputImage * >
70  ( this->itk::ProcessObject::GetInput(m_IndexImage2) );
71 }
72 
73 template <typename TInputImage, typename TMaskImage>
75 {
76  return static_cast< const TInputImage * >
77  ( this->itk::ProcessObject::GetInput(m_IndexImage3) );
78 }
79 
80 template <typename TInputImage, typename TMaskImage>
82 {
83  return static_cast< const TInputImage * >
84  ( this->itk::ProcessObject::GetInput(m_IndexImage4) );
85 }
86 
87 template <typename TInputImage, typename TMaskImage>
89 {
90  return static_cast< const TInputImage * >
91  ( this->itk::ProcessObject::GetInput(m_IndexImage5) );
92 }
93 
94 template <typename TInputImage, typename TMaskImage>
96 {
97  m_ImagesVector.clear();
98  m_JointHistogramInitial.clear();
99 
100  if(m_IndexImage1 < m_nbMaxImages){m_ImagesVector.push_back(this->GetInputImage1());}
101  if(m_IndexImage2 < m_nbMaxImages){m_ImagesVector.push_back(this->GetInputImage2());}
102  if(m_IndexImage3 < m_nbMaxImages){m_ImagesVector.push_back(this->GetInputImage3());}
103  if(m_IndexImage4 < m_nbMaxImages){m_ImagesVector.push_back(this->GetInputImage4());}
104  if(m_IndexImage5 < m_nbMaxImages){m_ImagesVector.push_back(this->GetInputImage5());}
105 
106  unsigned int histoDimension = m_ImagesVector.size();
107  std::vector<InputConstIteratorType> ImagesVectorIt;
108  for ( unsigned int i = 0; i < m_ImagesVector.size(); i++ )
109  {
110  InputConstIteratorType It(m_ImagesVector[i],m_ImagesVector[i]->GetLargestPossibleRegion() );
111  ImagesVectorIt.push_back(It);
112  }
113  MaskConstIteratorType MaskIt (this->GetMask(), this->GetMask()->GetLargestPossibleRegion() );
114  Histogram::iterator it;
115  while (!MaskIt.IsAtEnd())
116  {
117  if(MaskIt.Get()!=0)
118  {
119  Intensities value(histoDimension);
120  for(unsigned int m = 0; m < histoDimension; m++ )
121  {
122  if(static_cast<MeasureType>(ImagesVectorIt[m].Get()) < 0 )
123  value[m] = 0;
124  if(static_cast<MeasureType>(ImagesVectorIt[m].Get()) > static_cast<double>(std::numeric_limits<MeasureType>::max()))
125  value[m]= std::numeric_limits<MeasureType>::max();
126  else
127  value[m] = static_cast<MeasureType>(ImagesVectorIt[m].Get());
128  }
129  it = m_JointHistogramInitial.find(value);
130  if(it == m_JointHistogramInitial.end())
131  {
132  m_JointHistogramInitial.insert(Histogram::value_type(value,1));
133  }
134  else
135  {
136  (it->second)++;
137  }
138  }
139  for ( unsigned int i = 0; i < histoDimension; i++ )
140  {
141  ++ImagesVectorIt[i];
142  }
143  ++MaskIt;
144  }
145 }
146 
147 template <typename TInputImage, typename TMaskImage>
149 {
150  unsigned int nbClasses = m_GaussianModel.size();
151  this->m_APosterioriProbability.clear();
152 
153  //1. We calculate the inverse of the covariance and the determinant;
154  GaussianFunctionType::CovarianceMatrixType *inverseCovariance = new GaussianFunctionType::CovarianceMatrixType[nbClasses];
155  double *determinantCovariance = new double[nbClasses];
156  for(unsigned int i = 0 ; i < nbClasses; i++)
157  {
158  GaussianFunctionType::CovarianceMatrixType covar = (m_GaussianModel[i])->GetCovariance();
159 
160  determinantCovariance[i] = vnl_determinant(covar.GetVnlMatrix());
161  if(std::abs(determinantCovariance[i]) < 1e-12)
162  {
163  if(inverseCovariance != NULL)
164  {
165  delete[] inverseCovariance;
166  inverseCovariance=NULL;
167  }
168  if(determinantCovariance != NULL)
169  {
170  delete[] determinantCovariance;
171  determinantCovariance= NULL;
172  }
173  return 1.0;
174  }
175  inverseCovariance[i] = covar.GetInverse();
176  }
177 
178  //2. We calculate the a posterirori probability
179  Histogram::iterator histoIt;
180  GaussianFunctionType::CovarianceMatrixType intensities(this->m_JointHistogram.begin()->first.size(),1);
181  std::vector<double> probas(nbClasses);
182 
183  GaussianFunctionType::CovarianceMatrixType *xi = new GaussianFunctionType::CovarianceMatrixType[nbClasses];
184  unsigned int matrixSize = 0;
185  for(unsigned int i = 0; i < nbClasses; i++)
186  {
187  GaussianFunctionType::MeanVectorType mu = (m_GaussianModel[i])->GetMean();
188  matrixSize = mu.Size();
189  for(unsigned int j = 0; j < mu.Size(); j++)
190  {
191  (xi[i]).SetSize(mu.Size(),1);
192  (xi[i])(j,0) = mu[j];
193  }
194  }
195 
196  GaussianFunctionType::CovarianceMatrixType x;
197 
198  for(histoIt = this->m_JointHistogram.begin(); histoIt != this->m_JointHistogram.end(); ++histoIt)
199  {
200  //We set the intensities of the histogram in a vector
201  for(unsigned int i = 0; i < histoIt->first.size(); i++)
202  {
203  intensities(i,0) = static_cast<double>(histoIt->first[i]);
204  }
205 
206  // Calculate probability a posteriori for each pixel
207  // To eliminate problems with too small numbers we are going to substract in the exponetial
208  // the minimum found to at least have one "significant" value (equivalent to multiply the whole for a constant)
209  // Afterwards the a posteriory probability is normalize so this constant is eliminated
210  double minExpoTerm = 1e10;
211  double sumProba = 0.0;
212 
213  //Calculate the exponential term and the minimum of them
214  for(unsigned int i = 0; i < nbClasses; i++)
215  {
216  x = intensities - xi[i];
217 
218  double result = 0;
219  for (unsigned int j = 0; j < matrixSize; ++j)
220  {
221  result += inverseCovariance[i](j,j) * x(j,0) * x(j,0);
222  for (unsigned int k = j+1; k < matrixSize; ++k)
223  result += 2 * inverseCovariance[i](j,k) * x(j,0) * x(k,0);
224  }
225  probas[i] = result;
226 
227  if( minExpoTerm > probas[i])
228  {
229  minExpoTerm = probas[i];
230  }
231  }
232 
233  //Calculate numerator and denominator of a posteriori probability
234  for(unsigned int i = 0; i < probas.size(); i++)
235  {
236  // Constant * conditional probability of pixel knwoing gaussian i, last value multiply by the gaussian proportion
237  probas[i] = m_Alphas[i] * ((std::exp(0.5* (minExpoTerm - probas[i])) / std::sqrt(std::fabs(determinantCovariance[i]))));
238  // adding all classes to have the denominator
239  sumProba += probas[i];
240  }
241 
242  // This division eliminate the constant ment before
243  for(unsigned int i = 0; i < probas.size(); i++)
244  {
245  probas[i] /= sumProba;
246  }
247 
248  //Add probability to the GenericContainer
249  this->m_APosterioriProbability.insert(GenericContainer::value_type(histoIt->first, probas));
250  }
251 
252  double likelihoodValue = this->likelihood(inverseCovariance, determinantCovariance);
253 
254  // memory free
255  if(inverseCovariance != NULL)
256  {
257  delete[] inverseCovariance;
258  inverseCovariance = NULL;
259  }
260  if(determinantCovariance != NULL)
261  {
262  delete[] determinantCovariance;
263  determinantCovariance = NULL;
264  }
265 
266  delete[] xi;
267  xi = NULL;
268 
269  return likelihoodValue;
270 }
271 
272 template <typename TInputImage, typename TMaskImage>
273 bool GaussianEMEstimator<TInputImage,TMaskImage>::maximization(std::vector<GaussianFunctionType::Pointer> &newModel, std::vector<double> &newAlphas)
274 {
275  GenericContainer::iterator it;
276  Histogram::iterator histoIt;
277 
278  unsigned int numberOfClasses = m_GaussianModel.size();
279  unsigned int dimensions = this->m_JointHistogram.begin()->first.size();
280  double numberOfPixels = 0;
281 
282  // initializations for estimations
283  double *mixedProportions = new double[numberOfClasses];
284  double **means = new double*[numberOfClasses];
285  for (unsigned int i = 0; i < numberOfClasses; ++i)
286  {
287  means[i] = new double [dimensions];
288  }
289 
290  std::vector<GaussianFunctionType::CovarianceMatrixType> covariances(numberOfClasses, GaussianFunctionType::CovarianceMatrixType(dimensions,dimensions));
291 
292  for(unsigned int i = 0; i < numberOfClasses;i++)
293  {
294  mixedProportions[i] = 0.0;
295  for(unsigned int j = 0; j < dimensions; j++)
296  {
297  means[i][j] = 0.0;
298  for(unsigned int k = 0; k < dimensions; k++)
299  covariances[i](j,k)=0.0;
300  }
301  }
302 
303  //Mixing proportions and gaussian means
304  for(it = this->m_APosterioriProbability.begin(),histoIt = this->m_JointHistogram.begin(); it != this->m_APosterioriProbability.end(); ++it, ++histoIt) //for all histogram
305  {
306  numberOfPixels += static_cast<double>(histoIt->second);//counting all pixels
307  for(unsigned int i = 0; i < numberOfClasses; i++)
308  {
309  mixedProportions[i] += it->second[i] * histoIt->second; //mixedProportions [A posteriori probability] * [occurrences]
310 
311  for(unsigned int j = 0; j < dimensions; j++)
312  means[i][j] += it->second[i] * histoIt->second * histoIt->first[j]; // means: [A posteriori probability] * [occurrences]*[intensity]
313  }
314  }
315 
316  for(unsigned int i = 0; i < numberOfClasses; i++)
317  {
318  // normalization of means by sum( [A posteriori probability] * [occurrences])
319  for(unsigned int j = 0; j < dimensions; j++)
320  means[i][j] /= mixedProportions[i];
321  }
322 
323  // Covariance matrix for gaussians
324  for(it = this->m_APosterioriProbability.begin(), histoIt = this->m_JointHistogram.begin() ; it != this->m_APosterioriProbability.end(); ++it, ++histoIt) //for all histogram
325  {
326  for(unsigned int i = 0; i < numberOfClasses; i++)
327  for(unsigned int j = 0; j < dimensions; j++)
328  for(unsigned int k = j; k < dimensions; k++)
329  covariances[i](j,k) += it->second[i] * histoIt->second * ( histoIt->first[j] - means[i][j] ) * ( histoIt->first[k] - means[i][k] ); //[post proba] [occurrences] ([intensity]-[mean])^2
330  }
331 
332  double mix_pro;
333  for(unsigned int i = 0; i < numberOfClasses; i++)
334  {
335  for(unsigned int j = 0; j < dimensions; j++)
336  {
337  mix_pro = mixedProportions[i];
338  covariances[i](j,j) /= mixedProportions[i];
339  for(unsigned int k = j+1; k < dimensions; k++)
340  {
341  covariances[i](j,k) /= mixedProportions[i];
342  covariances[i](k,j) = covariances[i](j,k);
343  }
344  }
345  mixedProportions[i] /= static_cast<double>(numberOfPixels); // normalization of proportions by [numberOfPixels]
346  }
347 
348 
349  //storing values in an appropiate class
350  newModel.clear();
351  int *sort = new int[numberOfClasses]; //sorting in increasing order the means[0]
352  for (unsigned int i = 0; i < numberOfClasses;i++)
353  {
354  sort[i] =-1;
355  double minValue = 1e32;
356  for(unsigned int j = 0; j < numberOfClasses;j++)
357  {
358  bool used = false;
359  //we check j wasn't used yet
360  for(unsigned int k = 0; k < i; k++)
361  {
362  if(j == static_cast<unsigned int>(sort[k]))
363  {
364  used=true;
365  break;
366  }
367  }
368  // if not used we get the min
369  if(!used && means[j][0] < minValue)
370  {
371  minValue = means[j][0];
372  sort[i] = j;
373  }
374  }
375  }
376 
377  // storing classes in first image order
378  for (unsigned int i = 0; i < numberOfClasses;i++)
379  {
380  if(sort[i]==-1)
381  {
382  return false;
383  }
384 
385  GaussianFunctionType::MeanVectorType mu(dimensions);
386  for(unsigned int j = 0; j < dimensions; j++)
387  {
388  mu[j] = means[sort[i]][j];
389  }
390 
391  GaussianFunctionType::Pointer tmp = GaussianFunctionType::New();
392  tmp->SetMean(mu);
393  tmp->SetCovariance(covariances[sort[i]]);
394  newAlphas.push_back(mixedProportions[sort[i]]);
395  newModel.push_back(tmp);
396  }
397 
398  delete[] sort;
399  for (unsigned int i = 0; i < numberOfClasses; ++i)
400  {
401  delete[] (means[i]);
402  }
403  delete[] mixedProportions;
404 
405  return true;
406 }
407 
408 template <typename TInputImage, typename TMaskImage>
410 {
411  this->createJointHistogram();
412  this->m_JointHistogram = this->m_JointHistogramInitial;
413  unsigned int iter = 0; //number of current iterations
414  double distance = 0.0;
415  m_Likelihood = 0.0;
416 
417  do
418  {
419  m_Likelihood = this->expectation();
420  if( m_Likelihood >= 0.0 )
421  {
422  m_Likelihood = 0.0;
423  return;
424  }
425 
426  std::vector<GaussianFunctionType::Pointer> newModel;
427  std::vector<double> newAlphas;
428 
429  if( !this->maximization(newModel,newAlphas) )
430  {
431  m_Likelihood = 0.0;
432  return;
433  }
434  distance = this->computeDistance(newModel);
435  m_GaussianModel = newModel;
436  m_Alphas = newAlphas;
437  iter++;
438 
439  }while((distance > this->m_ModelMinDistance) && iter < this->m_MaxIterations);
440  m_Likelihood = this->expectation();
441 }
442 
443 template <typename TInputImage, typename TMaskImage>
444 double GaussianEMEstimator<TInputImage,TMaskImage>::likelihood(GaussianFunctionType::CovarianceMatrixType *invCovariance, double *detCovariance)
445 {
446  double likelihoodValue = 0.0;
447  unsigned int nbClasses = m_GaussianModel.size();
448  GaussianFunctionType::CovarianceMatrixType *inverseCovariance = NULL;
449  double *determinantCovariance = NULL;
450  if(invCovariance != NULL && detCovariance !=NULL)
451  {
452  inverseCovariance = invCovariance;
453  determinantCovariance = detCovariance;
454  }
455  else
456  {
457  inverseCovariance = new GaussianFunctionType::CovarianceMatrixType[nbClasses];
458  determinantCovariance = new double[nbClasses];
459 
460  //1. We calculate covariance inverse and determinant
461  for(unsigned int i = 0 ; i < nbClasses; i++)
462  {
463  GaussianFunctionType::CovarianceMatrixType covar = (m_GaussianModel[i])->GetCovariance();
464  determinantCovariance[i] = vnl_determinant(covar.GetVnlMatrix());
465  if(std::fabs(determinantCovariance[i]) < 1e-9)
466  {
467  return likelihoodValue;
468  }
469  inverseCovariance[i] = covar.GetInverse();
470  }
471  }
472 
473  //2. We calculate the a posteriori probability
474  Histogram::iterator histoIt;
475  GenericContainer::iterator it = this->m_APosterioriProbability.begin();
476 
477  GaussianFunctionType::CovarianceMatrixType intensities(this->m_JointHistogram.begin()->first.size(),1);
478 
479  for(histoIt = this->m_JointHistogram.begin(); histoIt != this->m_JointHistogram.end(); ++histoIt, ++it)
480  {
481  //We set the intensities of the histogram in a vector
482  for(unsigned int i = 0; i < histoIt->first.size(); i++)
483  {
484  intensities(i,0) = static_cast<double> (histoIt->first[i]);
485  }
486 
487  unsigned int maxIndex = 0;
488  double maxPostProba = 0.0;
489  //we look for the max post proba to resolve de ecuation
490  for(unsigned int i = 0; i < nbClasses; i++)
491  {
492  if(it->second[i] > maxPostProba)
493  {
494  maxPostProba = it->second[i];
495  maxIndex = i;
496  }
497  }
498 
499  GaussianFunctionType::CovarianceMatrixType x,xT,result;
500  GaussianFunctionType::MeanVectorType mu = (m_GaussianModel[maxIndex])->GetMean();
501  for(unsigned int i = 0; i < mu.Size(); i++)
502  {
503  x.SetSize(mu.Size(),1);
504  x(i,0) = mu[i];
505  }
506  x = intensities - x;
507  xT = (x.GetTranspose());
508  result = xT * inverseCovariance[maxIndex] * x;
509  double proba = result(0,0);
510 
511  likelihoodValue += histoIt->second * ( -proba/2.0 - std::log(std::sqrt(pow(2*M_PI,static_cast<int>(this->m_JointHistogram.begin()->first.size())) *
512  std::fabs(determinantCovariance[maxIndex])))
513  + std::log(m_Alphas[maxIndex]/maxPostProba));
514  }
515 
516  //Cleaning pointers if created in this function
517  if((invCovariance == NULL)&&(inverseCovariance!=NULL))
518  {
519  delete[] inverseCovariance;
520  inverseCovariance=NULL;
521  }
522  if((detCovariance ==NULL)&&(determinantCovariance!=NULL))
523  {
524  delete[] determinantCovariance;
525  determinantCovariance=NULL;
526  }
527 
528  return likelihoodValue;
529 }
530 
531 
532 template <typename TInputImage, typename TMaskImage>
533 double GaussianEMEstimator<TInputImage,TMaskImage>::computeDistance(std::vector<GaussianFunctionType::Pointer> &newModel)
534 {
535  unsigned int nbClasses = m_GaussianModel.size();
536 
537  if(newModel.size() != nbClasses)
538  return 1e9; //arbitrary high value
539 
540  double criterium = 0.0;
541 
542  unsigned int numberOfGaussians = (newModel[0])->GetMean().Size();
543 
544  for(unsigned int i = 0; i < m_GaussianModel.size(); i++)
545  {
546  GaussianFunctionType::MeanVectorType newMu = (m_GaussianModel[i])->GetMean();
547  GaussianFunctionType::MeanVectorType oldMu = (newModel[i])->GetMean();
548 
549  if(newMu.Size()!=oldMu.Size())
550  return 1e9;
551 
552  for(unsigned int j = 0; j < numberOfGaussians; j++)
553  criterium += std::fabs( oldMu[j] - newMu[j] );
554  }
555 
556  return criterium/(static_cast<double>(nbClasses * numberOfGaussians));
557 }
558 
559 }
std::vector< MeasureType > Intensities
void SetInputImage1(const TInputImage *image)
TInputImage::ConstPointer GetInputImage1()
void SetMask(const TMaskImage *mask)
itk::ImageRegionConstIterator< MaskImageType > MaskConstIteratorType
virtual void Update() ITK_OVERRIDE
TInputImage::ConstPointer GetInputImage2()
double xi(const T &k)
void SetInputImage5(const TInputImage *image)
double likelihood(GaussianFunctionType::CovarianceMatrixType *invCovariance=NULL, double *detCovariance=NULL)
virtual bool maximization(std::vector< GaussianFunctionType::Pointer > &newModel, std::vector< double > &newAlphas)
void SetInputImage4(const TInputImage *image)
itk::ImageRegionConstIterator< InputImageType > InputConstIteratorType
TInputImage::ConstPointer GetInputImage5()
TInputImage::ConstPointer GetInputImage3()
TInputImage::ConstPointer GetInputImage4()
double computeDistance(std::vector< GaussianFunctionType::Pointer > &newModel)
void SetInputImage2(const TInputImage *image)
TMaskImage::ConstPointer GetMask()
void SetInputImage3(const TInputImage *image)