ANIMA  4.0
animaGaussianREMEstimator.hxx
Go to the documentation of this file.
2 
3 namespace anima
4 {
5 
6 template <typename TInputImage, typename TMaskImage>
8 {
9  std::vector<vnl_matrix<GaussianREMEstimator::NumericType> > inverseCovariance;
10  inverseCovariance.reserve(this->m_GaussianModel.size());
11  double *determinantCovariance = new double[this->m_GaussianModel.size()];
12 
13  if(!this->m_APosterioriProbability.empty())
14  this->m_APosterioriProbability.clear();
15  if(!this->m_JointHistogram.empty())
16  this->m_JointHistogram.clear();
17 
18  GaussianFunctionType::CovarianceMatrixType *covari = new GaussianFunctionType::CovarianceMatrixType[this->m_GaussianModel.size()];
19  for(unsigned int i = 0; i < this->m_GaussianModel.size(); i++)
20  {
21  covari[i] = (this->m_GaussianModel[i])->GetCovariance();
22  }
23 
24  //1. We calculate covariance inverse and determinant
25  for(unsigned int i = 0; i < this->m_GaussianModel.size(); i++)
26  {
27  determinantCovariance[i] = vnl_determinant(covari[i].GetVnlMatrix());
28  if(std::fabs(determinantCovariance[i]) < 1e-12)
29  {
30  return false;
31  }
32  inverseCovariance.push_back(covari[i].GetInverse());
33  }
34 
35  Histogram::iterator histoIt;
36  GaussianFunctionType::CovarianceMatrixType intensities(this->m_OriginalJointHistogram.begin()->first.size(),1);
37  std::vector<double> probas(this->m_GaussianModel.size());
38 
39  ResidualMap residualMap;
40  double numberOfPixels=0;
41 
42  GaussianFunctionType::CovarianceMatrixType *xi = new GaussianFunctionType::CovarianceMatrixType[probas.size()];
43  for(unsigned int i = 0; i < probas.size(); i++)
44  {
45  GaussianFunctionType::MeanVectorType mu = (this->m_GaussianModel[i])->GetMean();
46  (xi[i]).SetSize(mu.Size(),1);
47  for(unsigned int j = 0; j < mu.Size(); j++)
48  {
49  (xi[i])(j,0) = mu[j];
50  }
51  }
52 
53  GaussianFunctionType::CovarianceMatrixType x,xT;
54  for(histoIt = this->m_OriginalJointHistogram.begin(); histoIt != this->m_OriginalJointHistogram.end(); ++histoIt)
55  {
56  //We set the intensities of the histogram in a vector
57  for(unsigned int i = 0; i < histoIt->first.size(); i++)
58  {
59  intensities(i,0) = static_cast<double>(histoIt->first[i]);
60  }
61 
62  double concentrationValue = 0.0;
63 
64  //Calculate the exponential term and the minimum of them
65  for(unsigned int i = 0; i < probas.size(); i++)
66  {
67  x = intensities - xi[i];
68  xT = (x.GetTranspose());
69  vnl_matrix<GaussianREMEstimator::NumericType> result = xT.GetVnlMatrix() * inverseCovariance[i] * x.GetVnlMatrix();
70  probas[i] = result(0,0)/2;
71  concentrationValue += this->m_Alphas[i] * std::exp(-probas[i]) / std::sqrt( determinantCovariance[i] );
72  }
73 
74  //Add value to the GenericContainer
75  // We are storing the value inside of the exponential( it will be use in expectation)
76  //this->m_APosterioriProbability.insert(GenericContainer::value_type(histoIt->first,probas));
77  //Fill residualmap with probabilities of the mixed gaussian (probability = constant * concentrationvalue)
78  residualMap.insert(ResidualMap::value_type(std::log(concentrationValue),histoIt->first));
79  numberOfPixels+=histoIt->second;
80  }
81 
82  ResidualMap::iterator it = residualMap.begin();
83 
84  //number of rejected pixels
85  double numberOfRejections = this->m_RejectionRatio * numberOfPixels;
86  double rejected = 0;
87  this->m_JointHistogram = this->m_OriginalJointHistogram;
88 
89  for(it=residualMap.begin(); it != residualMap.end(); ++it)
90  {
91  if(rejected >= numberOfRejections)
92  break;
93  histoIt = this->m_JointHistogram.find(it->second);
94  if(histoIt == this->m_JointHistogram.end())
95  {
96  return false;
97  }
98  double actual = histoIt->second;
99  if(actual+rejected >= numberOfRejections)
100  {
101  //We pass the limit...we get only some points of this Intensities
102  histoIt->second = actual+rejected-numberOfRejections;
103  break;
104  }
105  else
106  {
107  //We don't pass the limit... we eliminate this Intensities
108  this->m_JointHistogram.erase(histoIt);
109  rejected += actual;
110  }
111  }
112 
113  delete[] determinantCovariance;
114  determinantCovariance = NULL;
115 
116  delete[] xi;
117  xi = NULL;
118 
119  return true;
120 }
121 
122 template <typename TInputImage, typename TMaskImage>
124 ::PrintSolution(std::vector<double> alphas, std::vector<GaussianFunctionType::Pointer> model)
125 {
126  unsigned int nbTissus = alphas.size();
127  unsigned int nbModalities = (model[0])->GetMean().Size();
128  for(unsigned int i = 0; i < nbTissus; i++)
129  {
130  std::cout << "* Class: " << i << std::endl;
131  std::cout << " Alpha: " << alphas[i] << std::endl;
132  GaussianFunctionType::MeanVectorType mu = (model[i])->GetMean();
133  std::cout << " Mean: " << std::endl;
134  for(unsigned int j = 0; j < nbModalities; j++)
135  {
136  std::cout << " " << mu[j] << std::endl;
137  }
138 
139  GaussianFunctionType::CovarianceMatrixType covar = (model[i])->GetCovariance();
140  std::cout << " covar: " << std::endl;
141  for(unsigned int k = 0; k < nbModalities; k++)
142  {
143  std::cout << " " << covar(k,0) << " " << covar(k,1) << " " << covar(k,2) << std::endl;
144  }
145  }
146  return 0;
147 }
148 
149 template <typename TInputImage, typename TMaskImage>
151 {
152  this->createJointHistogram();
153 
154  this->m_OriginalJointHistogram = this->m_JointHistogramInitial;
155 
156  unsigned int iter = 0; //number of current iterations
157  double distance = 0.0;
158 
159  this->m_Likelihood = 0.0;
160 
161  if( this->m_StremMode )
162  this->m_JointHistogram = this->m_OriginalJointHistogram;
163  else
164  {
165  if( !this->concentration() )
166  {
167  this->m_Likelihood = 0.0;
168  return;
169  }
170  }
171 
172  do
173  {
174  int emIter = 0; // iterations between concentration steps
175  do
176  {
177  this->m_Likelihood = this->expectation();
178 
179  if( this->m_Likelihood >= 0.0 )
180  {
181  this->m_Likelihood = 0.0;
182  return;
183  }
184 
185  std::vector<GaussianFunctionType::Pointer> newModel;
186  std::vector<double> newAlphas;
187  if( !this->maximization(newModel,newAlphas) )
188  {
189  this->m_Likelihood = 0.0;
190  return;
191  }
192 
193  distance = this->computeDistance(newModel);
194 
195  this->m_GaussianModel = newModel;
196  this->m_Alphas = newAlphas;
197 
198  if( this->m_Verbose )
199  {
200  std::cout << "current iteration / max iteration between concentration step: "<< emIter << " / " << this->m_MaxIterationsConc << std::endl;
201  std::cout << "current distance / minimum distance: "<< distance << " / " << this->m_ModelMinDistance << std::endl;
202  std::cout << "Current model: " << std::endl;
203  this->PrintSolution(this->m_Alphas, this->m_GaussianModel);
204  std::cout << std::endl;
205 
206  double ratio = this->m_ModelMinDistance / distance;
207  if( ratio > 1 )
208  {
209  ratio = 1;
210  }
211  this->UpdateProgress(ratio);
212  std::cout << std::endl;
213  }
214 
215  emIter++;
216 
217  }while((distance > this->m_ModelMinDistance) && emIter < this->m_MaxIterationsConc);
218 
219  if( !this->concentration() )
220  {
221  this->m_Likelihood = 0.0;
222  return;
223  }
224 
225  iter++;
226  }while((distance > this->m_ModelMinDistance) && iter < this->m_MaxIterations);
227 
228  this->m_Likelihood = this->expectation();
229 }
230 
231 }
double xi(const T &k)
int PrintSolution(std::vector< double > alphas, std::vector< GaussianFunctionType::Pointer > model)
virtual void Update() ITK_OVERRIDE
std::map< double, std::vector< unsigned short > > ResidualMap