ANIMA  4.0
animaClassificationStrategy.hxx
Go to the documentation of this file.
2 
3 
4 namespace anima
5 {
6 template <typename TInputImage, typename TMaskImage>
8 {
9  if( this->m_RandomInitializer.IsNull() )
10  {
11  std::cerr <<" -- Error in ClassificationStrategy: RandomInitializer is null"<<std::endl;
12  return;
13  }
14 
15  m_ListGaussianModels.clear();
16  m_ListAlphas.clear();
17 
18  m_ListGaussianModels.resize(this->m_NumberOfEstimators.size());
19  m_ListAlphas.resize(this->m_NumberOfEstimators.size());
20 
21  // First iteration with RandomInitialization
22  unsigned int errors = 0;
23 
24  for(unsigned int iter = 0; iter < this->m_NumberOfEstimators[0]; iter++)
25  {
26  this->m_RandomInitializer->Update();
27 
28  std::vector<GaussianFunctionType::Pointer> model = this->m_RandomInitializer->GetInitialization();
29  std::vector<double> alphas = this->m_RandomInitializer->GetAlphas();
30 
31  this->m_Estimator->SetInitialGaussianModel(model);
32  this->m_Estimator->SetInitialAlphas(alphas);
33  this->m_Estimator->Update();
34  double likelihood = m_Estimator->GetLikelihood();
35 
36  if( (!m_EM_Mode && likelihood > 0.0) || (m_EM_Mode && likelihood < 0.0) )
37  {
38  std::vector<GaussianFunctionType::Pointer> result = this->m_Estimator->GetGaussianModel();
39  std::vector<double> resultAlpha = this->m_Estimator->GetAlphas();
40 
41  if(m_EM_Mode)
42  {
43  likelihood=-likelihood; //we change the sign to have the best solution in the first position
44  }
45 
46  std::map< double, std::vector<GaussianFunctionType::Pointer> >::iterator mapIt = this->m_ListGaussianModels[0].begin();
47  bool foundSame = false;
48 
49  for(mapIt = this->m_ListGaussianModels[0].begin(); mapIt != this->m_ListGaussianModels[0].end(); ++mapIt)
50  {
51  if(sameModel(mapIt->second, result))
52  {
53  foundSame = true;
54  break;
55  }
56 
57  }
58 
59  if(!foundSame)
60  {
61  this->m_ListGaussianModels[0].insert(std::map<double, std::vector<GaussianFunctionType::Pointer> >::value_type(likelihood,result));
62  this->m_ListAlphas[0].insert(std::map<double,std::vector<double> >::value_type(likelihood,resultAlpha));
63  }
64 
65  }
66  else
67  {
68  iter--;
69  errors++;
70  if(errors > 10 * this->m_NumberOfEstimators[0])
71  break;
72  }
73 
74 
75  double ratio = static_cast<double>(iter+1) / static_cast<double>(this->m_NumberOfEstimators[0]);
76  this->UpdateProgress(ratio);
77  }
78 
79 
80  // if several estimators
81  for(unsigned int step = 1; step < this->m_NumberOfEstimators.size(); step++)
82  {
83  if(this->m_ListGaussianModels[step-1].size() == 0)
84  {
85  std::cout<<"-- ERROR in ClassificationStrategy: No solution found in step " << step-1 <<std::endl;
86  return;
87  }
88 
89  std::map<double,std::vector<GaussianFunctionType::Pointer> >::iterator mapIt = this->m_ListGaussianModels[step-1].begin();
90  for(unsigned int iter = 0 ; mapIt != this->m_ListGaussianModels[step-1].end() && iter< this->m_NumberOfEstimators[step]; iter++, ++mapIt)
91  {
92  this->m_Estimator->SetInitialGaussianModel(mapIt->second);
93 
94  m_Estimator->Update();
95  double likelihood = m_Estimator->GetLikelihood();
96 
97  if((m_EM_Mode && likelihood < 0.0) || (!m_EM_Mode && likelihood > 0.0))
98  {
99  std::vector<GaussianFunctionType::Pointer> result = m_Estimator->GetGaussianModel();
100  std::vector<double> resultAlpha = m_Estimator->GetAlphas();
101 
102  mapIt = this->m_ListGaussianModels[step].begin();
103  bool foundSame = false;
104  for(mapIt = this->m_ListGaussianModels[step].begin(); mapIt != this->m_ListGaussianModels[step].end(); ++mapIt)
105  {
106  if(sameModel(mapIt->second, result))
107  {
108  foundSame=true;
109  break;
110  }
111  }
112 
113  if(m_EM_Mode)
114  likelihood=-likelihood; //we change sign to have best m_Solutions first
115 
116  if(!foundSame)
117  this->m_ListGaussianModels[step].insert(std::map<double,std::vector<GaussianFunctionType::Pointer> >::value_type(likelihood,result));
118  this->m_ListAlphas[step].insert(std::map<double,std::vector<double> >::value_type(likelihood,resultAlpha));
119  }
120  else
121  {
122  //error
123  iter--;
124  errors++;
125  if(errors > 10 * this->m_NumberOfEstimators[step])
126  break;
127  }
128  }
129  }
130 }
131 
132 template <typename TInputImage, typename TMaskImage>
133 bool ClassificationStrategy<TInputImage,TMaskImage>::GetSolutionMap(std::map< double, std::vector<GaussianFunctionType::Pointer> > &solution,
134  std::map< double, std::vector<double> > &solutionAlpha, int step)
135 {
136  if(this->m_ListGaussianModels.size() == 0 || static_cast<int>(this->m_ListGaussianModels.size()) <= step)
137  return false;
138 
139  if(step < 0 )
140  {
141  step=this->m_ListGaussianModels.size()-1;
142  }
143 
144  if(this->m_ListGaussianModels[step].size() == 0 )
145  return false;
146 
147  solution = this->m_ListGaussianModels[step];
148  solutionAlpha = this->m_ListAlphas[step];
149  return true;
150 }
151 
152 
153 template <typename TInputImage, typename TMaskImage>
154 void ClassificationStrategy<TInputImage,TMaskImage>::SetStrategy( std::vector< unsigned int >& ems,std::vector<unsigned int> &iters )
155 {
156  this->m_NumberOfEstimators=ems;
157  this->m_NumberOfIterations=iters;
158 
159  if(this->m_NumberOfEstimators.size()!= this->m_NumberOfIterations.size())
160  {
161  while( this->m_NumberOfEstimators.size() > this->m_NumberOfIterations.size() )
162  {
163  this->m_NumberOfEstimators.pop_back();
164  }
165  while( this->m_NumberOfEstimators.size() < this->m_NumberOfIterations.size() )
166  {
167  this->m_NumberOfIterations.pop_back();
168  }
169  }
170 
171 }
172 
173 template <typename TInputImage, typename TMaskImage>
174 bool ClassificationStrategy<TInputImage,TMaskImage>::sameModel( std::vector<GaussianFunctionType::Pointer> &mod1, std::vector<GaussianFunctionType::Pointer> &mod2)
175 {
176  unsigned int NbClasses = mod1.size();
177  unsigned int NbDimension = (mod1[0])->GetMean().Size();
178 
179  std::vector<int> comparisonVector(NbClasses,-1);
180  unsigned int params = NbDimension + NbDimension*NbDimension;
181 
182  double *p1 = new double[params];
183  double *p2 = new double[params];
184 
185  for(unsigned int i = 0; i < NbClasses; i++)
186  {
187  for(unsigned int j = 0; j < NbClasses; j++)
188  {
189  if((mod1[i])->GetMean().Size() != (mod2[j])->GetMean().Size())
190  {
191  comparisonVector[i] = -1;
192  continue;
193  }
194 
195  unsigned int t = 0;
196  for(unsigned int l = 0; l < NbDimension; l++)
197  {
198  p1[t] = static_cast<double>((mod1[i])->GetMean()[l]);
199  p2[t] = static_cast<double>((mod2[j])->GetMean()[l]);
200  t++;
201  }
202 
203  for(unsigned int l = 0; l < NbDimension; l++)
204  {
205  for(unsigned int k = l; k < NbDimension; k++)
206  {
207  p1[t] = static_cast<double>((mod1[i])->GetCovariance()[l][k]);
208  p2[t] = static_cast<double>((mod2[j])->GetCovariance()[l][k]);
209  t++;
210  }
211  }
212 
213  double distance = 0.0;
214  for(unsigned int k = 0; k < params; k++)
215  {
216  distance+=p1[k]-p2[k];
217  }
218 
219  distance/=static_cast<double>(params);
220  if(distance < 1e-2)
221  {
222  comparisonVector[i]=j;
223  break;
224  }
225  }
226  }
227 
228  delete[] p1;
229  delete[] p2;
230 
231  for (unsigned int i=0; i < NbClasses; i++)
232  {
233  if(comparisonVector[i]==-1)
234  return false;
235  for(unsigned int j=i+1;j < NbClasses; j++)
236  if( comparisonVector[i]==comparisonVector[j])
237  return false;
238  }
239 
240  return true;
241 }
242 
243 
244 }
bool sameModel(std::vector< GaussianFunctionType::Pointer > &mod1, std::vector< GaussianFunctionType::Pointer > &mod2)
compare two models
void SetStrategy(std::vector< unsigned int > &ems, std::vector< unsigned int > &iters)
bool GetSolutionMap(std::map< double, std::vector< GaussianFunctionType::Pointer > > &solution, std::map< double, std::vector< double > > &solutionAlpha, int step=-1)
return the solutions map for each step
void Update() ITK_OVERRIDE
executes the strategy