A.I.M.S algorithms


ppca_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 
36 #ifndef AIMS_MATH_PPCA_D_H
37 #define AIMS_MATH_PPCA_D_H
38 
40 #include <aims/data/data.h>
41 #include <aims/math/svd.h>
42 #include <aims/math/ppca.h>
43 #include <aims/math/mathelem.h>
44 #include <aims/math/gausslu.h>
45 #include <aims/io/writer.h>
47 
48 
49 namespace aims
50 {
51 
52  template <typename T>
53  void
54  ProbabilisticPcaElement::doIt( const std::list< Point3d >& selectedPoints,
55  const AimsData<T>& data, double distanceRef )
56  {
57  if( selectedPoints.size() <= (unsigned) data.dimT() )
58  {
59  _valid = false ;
60  _computed = true ;
61  return ;
62  }
63  _valid = true ;
64  AimsFastAllocationData<T> individuals( std::max(int(selectedPoints.size()), 1),
65  data.dimT() );
66  int i = 0 ;
67  std::list< Point3d >::const_iterator iter( selectedPoints.begin() ),
68  last( selectedPoints.end() ) ;
69  Point3df mean ;
70 
71  while( iter != last )
72  {
73  for(int t = 0 ; t < data.dimT() ; ++t )
74  individuals( i, t ) = T( data( (*iter)[0], (*iter)[1], (*iter)[2],
75  t ) ) ;
76  ++i ; ++iter ;
77  //std::cout << "Indiv " << i << " = " << *iter << std::endl ;
78  mean += Point3df( (*iter)[0], (*iter)[1], (*iter)[2] ) ;
79  }
80  //std::cout << "Center = " << mean / float( individuals.dimX() )
81  // << std::endl ;
82  doIt( individuals, distanceRef ) ;
83  }
84 
85  template <typename T>
86  void
88  double distanceRef )
89  {
90  _distanceRef = distanceRef ;
91 // AimsData<double> CiT(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
92 // AimsData<double> invCiT(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
93 // AimsData<double> invCi2T(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
94 // AimsData<double> Wi2T(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
95 // AimsData<double> matVarCovT(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
96 // AimsData<double> CiD(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
97 // AimsData<double> invCiD(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
98 // AimsData<double> invCi2D(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
99 // AimsData<double> Wi2D(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
100 // AimsData<double> matVarCovD(matriceIndiv.dimY(), matriceIndiv.dimY()) ;
101 
102  _x = AimsData<double>( matriceIndiv.dimY() ) ;
103  _xT = AimsData<double>( 1, matriceIndiv.dimY() ) ;
104  /* int ignoreVPs = 0 ; */
105  std::cout << "toto" << std::endl << std::endl ;
106 
107  if( matriceIndiv.dimX() <= matriceIndiv.dimY() ){
108  _valid = false ;
109  _computed = true ;
110  return ;
111  }
112  _valid = true ;
113 
114  //----------------- Just for testing to separate kidney core and cortex
115  int ignoreVPs = 0 ;
116 
117  _computed = true ;
118  int nbFrame = matriceIndiv.dimY() ;
119 
120  if( matriceIndiv.dimY() < nbFrame )
121  std::cerr << "Beware : not enough data to evaluate svd" << std::endl ;
122 
123  AimsFastAllocationData<double> centeredIndivMatrix( matriceIndiv.dimX(), matriceIndiv.
124  dimY() ) ;
125  _mean = AimsFastAllocationData<double>(matriceIndiv.dimY()) ;
126  _mean = 0. ;
127  T val ;
128  for( int ind = 0 ; ind < matriceIndiv.dimX() ; ++ind )
129  for( int t = 0 ; t < matriceIndiv.dimY() ; ++t )
130  {
131  val = matriceIndiv( ind, t ) ;
132  _mean[t] += val ;
133  }
134 
135  //TEMP
136  double meanMean = 0. ;
137  //END TEMP
138 
139  for( int t = 0 ; t < nbFrame ; ++t )
140  {
141  _mean[t] /= centeredIndivMatrix.dimX() ;
142  for( int ind = 0 ; ind < centeredIndivMatrix.dimX() ; ++ind ){
143  centeredIndivMatrix( ind, t ) = matriceIndiv(ind, t) - _mean[t] ;
144  meanMean += centeredIndivMatrix( ind, t ) ;
145  }
146  }
147  std::cout << "Mean centered mean = " << meanMean << " (should be 0) " << std::endl ;
148 
149  /* std::cout << "Mean : " ; */
150  /* for( int i = 0 ; i < _mean.dimX() ; ++ i ){ */
151  /* std::cout << _mean(i) << " " ; */
152  /* } */
153  /* std::cout << std::endl << _meanMean << std::endl ; */
154 
155  Writer< AimsData< double> > wri02( "centeredIndiv.ima" ) ;
156  wri02.write(centeredIndivMatrix) ;
157 
158  // Matrice des correlations
159  AimsFastAllocationData<double> matVarCov(nbFrame, nbFrame) ;
160  int x, y ;
161  ForEach2d( matVarCov, x, y )
162  {
163  for(int k=0; k < centeredIndivMatrix.dimX() ;++k)
164  matVarCov(x, y) += centeredIndivMatrix(k, x)
165  * centeredIndivMatrix(k, y);
166  matVarCov(x, y) /= centeredIndivMatrix.dimX() - 1 ;
167  }
168 
169  Writer< AimsData< double> > wriMVC( "matVarCov.ima" ) ;
170  wriMVC.write(matVarCov) ;
171 
172  // Decomposition SVD
173  AimsSVD< double > svd;
175  _EValues = svd.doit( matVarCov );
176 // Writer< AimsData< double> > wri04( "eValAvant.ima" ) ;
177 // wri04.write(_EValues) ;
178 
179  Writer< AimsData< double> > wriMVC2( "matVarCov2.ima" ) ;
180  wriMVC2.write((matVarCov.clone().transpose()).cross( _EValues.cross(matVarCov) ) ) ;
181 
182  Writer< AimsData< double> > wriMVC3( "matVarCov3.ima" ) ;
183  wriMVC3.write((matVarCov).cross( _EValues.cross(matVarCov.clone().transpose()) ) ) ;
184 
185  svd.sort(matVarCov, _EValues) ;
186 // Writer< AimsData< double> > wri05( "eValApres.ima" ) ;
187 // wri05.write(_EValues) ;
188 
190  for( int i = 0 ; i < nbFrame ; ++i )
191  for( int j = 0 ; j < _nbOfSignificantEV ; ++j )
192  _EVect(i, j) = matVarCov(i, j) ;
193 
194  _Sigma2 = 0. ;
195  for( int i = _nbOfSignificantEV + ignoreVPs ; i < nbFrame ; ++i )
196  _Sigma2 += _EValues(i, i) ;
197 
198  _Sigma2 /= double(nbFrame - ( _nbOfSignificantEV + ignoreVPs ) ) ;
199 
200  for( int c = 0 ; c < _EValues.dimX() ; ++c )
201  std::cout << "\tLamda " << c+1 << "EValue = " << _EValues(c, c) ;
202  std::cout <<"\tSigma2 = " << _Sigma2 << std::endl ;
203 
204 // std::cout << "Lamda 1 = " << _EValues(0, 0) << "Lamda 2 = "
205 // << _EValues(1, 1)
206 // <<"\tSigma2 = " << _Sigma2 << std::endl ;
207 
209  for( int i = 0; i < _nbOfSignificantEV ; ++i )
210  _Wi(i, i) = sqrt( _EValues(i) - _Sigma2 ) ;
211 
212  std::cout << "_Wi.dimZ = " << _Wi.dimZ()
213  << "\t_Wi.dimT = " << _Wi.dimT() << std::endl ;
214 
215  std::cout << "_EVect.dimZ = " << _EVect.dimZ()
216  << "\t_EVect.dimT = " << _EVect.dimT() << std::endl ;
217 
218  _Wi = _EVect.cross( _Wi ) ;
219 // Writer< AimsData< double> > wri0( "Wi.ima" ) ;
220 // wri0.write(_Wi) ;
221 
222  AimsFastAllocationData<double> Ci = _Wi.cross( _Wi.clone().transpose() ) ;
223 // Writer< AimsData< double > > wri1( "Wi2.ima" ) ;
224 // wri1.write(Ci) ;
225 
226  ForEach2d( matVarCov, x, y )
227  if( x < y ){
228  double val = 0.5 * ( Ci(x, y) + Ci(y, x) ) ;
229  Ci(x, y) = Ci(y, x) = val ;
230  }
231 
232  for( int i = 0; i < nbFrame ; ++i )
233  Ci(i, i) += _Sigma2 ;
234 
235 // Writer< AimsData< double> > wri2( "Ci.ima" ) ;
236 // wri2.write(Ci) ;
237 
238  //-------------------------------------------------------------------------
239  double s = 0.00000001 ;
240  AimsFastAllocationData< double > u = Ci.clone(); // because it is replaced on SVD output
242 
243  AimsSVD< double > svd2;
244  AimsFastAllocationData< double > w = svd2.doit( u, &v );
245 
246  // Set the singular values lower than threshold to zero
247  double ts = s * w.maximum();
248  int i, n = w.dimX();
249  _detCi = 1. ;
250  short unkept = 0 ;
251  for ( i=0; i<n; i++ )
252  {
253  _detCi *= w( i, i ) / _Sigma2 ;
254  if ( w( i, i ) > ts ) w( i, i ) = 1.0f / w( i, i );
255  else {
256  w( i, i ) = 0.0f;
257  ++unkept ;
258  }
259  }
260 
261  std::cout << "Unkept values = " << unkept << std::endl ;
262 
264  _invCi = v.cross( w.cross( u.transpose() ) ) ;
265 
266 // Writer< AimsData< double> > wriI( "invCi.ima" ) ;
267 // wriI.write(_invCi) ;
268 
269  //TEMP
270  AimsData<float> ciFloat( Ci.dimX(), Ci.dimY() ) ;
271  ForEach2d(ciFloat, x, y)
272  ciFloat(x, y) = float( Ci(x, y) ) ;
273 
274 // Writer< AimsData< float> > wriI2( "invCi2.ima" ) ;
275 // wriI2.write( AimsInversionLU(ciFloat) ) ;
276 
277 // Writer< AimsData< float> > wriI3( "invCi3.ima" ) ;
278 // wriI3.write( AimsInversionLU( AimsInversionLU(ciFloat) ) ) ;
279  // END TEMP
280 
281  std::cout << "_invCiInside.dimZ = " << _invCi.dimZ()
282  << "\t_invCiInside.dimT = " << _invCi.dimT() << std::endl ;
283 
284 // _normFactor = 1. / ( pow( 2*M_PI, nbFrame/2. ) * pow( _detCi, 0.5 )
285 // * pow( _Sigma2, nbFrame/2. ) ) ;
286  _normFactor = 1. / ( pow(_detCi, 0.5 ) * pow( _Sigma2, nbFrame/2. ) ) ;
287 
288  std::cout << "-0.5 * log( _detCi ) = " << -0.5 * log( _detCi ) << "\t log(_PIj) = " << log(_PIj) << std::endl ;
289  //x << "- (nbFrame/2) * log(_Sigma2) = " << - (nbFrame/2) * log(_Sigma2) << std::endl ;
290 
291  _lnAddFactor = log( 1. / sqrt(_detCi) / pow( _Sigma2, nbFrame/2. ) ) ;
292  //_lnAddFactor = -0.5 * log( _detCi ) /*+ log(_PIj)*/ - (nbFrame/2) * log(_Sigma2) ;
293 
294 
295  //_lnAddFactor = log( _detCi / ( _PIj * _PIj ) ) + nbFrame * log( _Sigma2 ) ;
296 
297  std::cout << "Det Ci = " << _detCi << "\tNorm Factor = " << _normFactor
298  << "\tLn Add Factor = " << _lnAddFactor << "\tCorresp factor = " << exp(_lnAddFactor)<< std::endl ;
299 
300 // AimsData<float> invCi2( AimsInversionLU(ciFloat) ) ;
301 // ForEach2d(matVarCov, x, y){
302 // CiT(y, x) = Ci(x, y) ;
303 // CiD(x, y) = Ci(x, y) - Ci(y, x) ;
304 // invCiT(y, x) = _invCi(x, y) ;
305 // invCiD(x, y) = _invCi(x, y) - _invCi(y, x) ;
306 // invCi2T(y, x) = invCi2(x, y) ;
307 // invCi2D(x, y) = invCi2(x, y) - invCi2(y, x) ;
308 // // Wi2T(y, x) = Wi2(x, y) ;
309 // // Wi2D(x, y) = Wi2(x, y) - Wi2(y, x) ;
310 // }
311 // Writer< AimsData<double> > w1("matVarCovT.ima") ;
312 // w1.write(matVarCovT) ;
313 // Writer< AimsData<double> > w2("matVarCovD.ima") ;
314 // w2.write(matVarCovD) ;
315 // Writer< AimsData<double> > w3("CiT.ima") ;
316 // w3.write(CiT) ;
317 // Writer< AimsData<double> > w4("CiD.ima") ;
318 // w4.write(CiD) ;
319 // Writer< AimsData<double> > w5("invCiT.ima") ;
320 // w5.write(invCiT) ;
321 // Writer< AimsData<double> > w6("invCiD.ima") ;
322 // w6.write(invCiD) ;
323 // Writer< AimsData<double> > w7("invCi2T.ima") ;
324 // w7.write(invCi2T) ;
325 // Writer< AimsData<double> > w8("invCi2D.ima") ;
326 // w8.write(invCi2D) ;
327 
328  ForEach2d( matVarCov, x, y )
329  if( x < y ){
330  double val = 0.5 * ( _invCi(x, y) + _invCi(y, x) ) ;
331  _invCi(x, y) = _invCi(y, x) = val ;
332  }
333 // ForEach2d( matVarCov, x, y )
334 // invCiD(x, y) = _invCi(x, y) - _invCi(y, x) ;
335 
336 // Writer< AimsData<double> > w9("invCiSymD.ima") ;
337 // w9.write(invCiD) ;
338 
339  }
340 
341  template <class T>
342  double
344  double & meanNorm ) const
345  {
346  int ignoreVPs = 0 ;
347 
348  int nbFrame = matriceIndiv.dimY() ;
349 
350  if( matriceIndiv.dimY() < nbFrame )
351  std::cerr << "Beware : not enough data to evaluate svd" << std::endl ;
352 
353  //-------------------TMP-------------------------------
354  //AimsData<T> indMatCopy(matriceIndiv.clone()) ;
355  //Writer< AimsData<T> > indWri("indiv.ima") ;
356  //indWri.write(indMatCopy) ;
357 
358  //-------------------TMP-------------------------------
359 
360  AimsFastAllocationData<double> centeredIndivMatrix( matriceIndiv.dimX(),
361  matriceIndiv.dimY() ) ;
362  AimsFastAllocationData<double> mean(matriceIndiv.dimY()) ;
363 
364  for( int ind = 0 ; ind < matriceIndiv.dimX() ; ++ind )
365  for( int t = 0 ; t < matriceIndiv.dimY() ; ++t )
366  mean[t] += matriceIndiv( ind, t ) ;
367 
368 
369  meanNorm = 0. ;
370  for( int t = 0 ; t < nbFrame ; ++t ){
371  mean[t] /= centeredIndivMatrix.dimX() ;
372  meanNorm += mean[t] * mean[t] ;
373  for( int ind = 0 ; ind < centeredIndivMatrix.dimX() ; ++ind ){
374  centeredIndivMatrix( ind, t ) = matriceIndiv(ind, t) - mean[t] ;
375  }
376  }
377  meanNorm /= mean.dimX() ;
378 
379 
380  // Matrice des correlations
381  AimsFastAllocationData<double> matVarCov(nbFrame, nbFrame) ;
382  int x, y ;
383  ForEach2d( matVarCov, x, y )
384  {
385  for(int k=0; k < centeredIndivMatrix.dimX() ;++k)
386  matVarCov(x, y) += centeredIndivMatrix(k, x)
387  * centeredIndivMatrix(k, y);
388  matVarCov(x, y) /= centeredIndivMatrix.dimX() - 1 ;
389  }
390 
391  //Writer< AimsFastAllocationData<T> > indWri("indiv.ima") ;
392  //indWri.write(matriceIndiv) ;
393 
394  //Writer< AimsData<double> > matVarWri("matVarCov.ima") ;
395  //matVarWri.write(matVarCov) ;
396 
397  // Decomposition SVD
398  AimsSVD< double > svd;
400  AimsFastAllocationData<double> EValues = svd.doit( matVarCov );
401 
402  svd.sort(matVarCov, EValues) ;
403 
404 
405  double Sigma2 = 0. ;
406  for( int i = _nbOfSignificantEV + ignoreVPs ; i < nbFrame ; ++i )
407  Sigma2 += EValues(i, i) ;
408 
409  Sigma2 /= double(nbFrame - ( _nbOfSignificantEV + ignoreVPs ) ) ;
410 
411  return Sigma2 ;
412  }
413 
414 
415  template <class T>
417  const std::vector<
418  std::list< Point3d > >&
419  classes,
420  int nbOfSignificantEigenValues,
421  const std::vector<double>& PIj ) :
422  _classes(classes), _data(data), _PIj(PIj),
423  _nbOfSignificantEigenValues(nbOfSignificantEigenValues)
424  {
425  //std::cout << "ProbabilisticPca<T>::ProbabilisticPca classes size "
426  // << _classes.size() << std::endl ;
427  if( PIj.size() != 0 )
428  if( PIj.size() == classes.size() )
429  _PIj = PIj ;
430  else
431  throw std::runtime_error( "Inconsistant data : number of classes and "
432  "number of classes prior probabilities are "
433  "different !") ;
434  else
435  {
436  _PIj.reserve( _classes.size() ) ;
437  for( unsigned int c = 0 ; c < _classes.size() ; ++c )
438  _PIj.push_back( 1./_classes.size() ) ;
439  }
440 
441  _discrElements.reserve( _classes.size() ) ;
442  _distanceRef = 0. ;
443 
444  // ------------------ TEST ---------------------
445  AimsFastAllocationData<double> eigenValues( _data.dimT(), _classes.size() ) ;
446  // ------------------ TEST ---------------------
447 
448 
449  for( unsigned int c = 0 ; c < _classes.size() ; ++c ){
450  if( int(_classes[c].size()) > _data.dimT() ) {
451  ProbabilisticPcaElement el( _nbOfSignificantEigenValues, _PIj[c] ) ;
452  if( c == 0 ){
453  el.doIt( _classes[c], _data ) ;
454  _distanceRef = log( el.normFactor() ) ;
455  std::cout << "Noise Ref = " << _distanceRef << std::endl ;
456  }
457  el.doIt( _classes[c], _data, _distanceRef ) ;
458 
459  // ------------------ TEST ---------------------
460  for( int t = 0 ; t < _data.dimT() ; ++t )
461  eigenValues( t, c ) = el.eigenValues()(t) ;
462  // ------------------ TEST ---------------------
463 
464  std::cout << "Class " << c << " done !" << std::endl ;
465  _discrElements.push_back( el ) ;
466  } else {
467  std::cerr << "Not enough data for element " << c << std::endl ;
468  }
469  }
470 
471  // ------------------ TEST ---------------------
472  Writer< AimsData<double> > eigenValuesWri( "classesEigenValues.ima") ;
473  eigenValuesWri.write(eigenValues) ;
474  std::cout << "Probabilistic Pca initization completed" << std::endl ;
475  //std::cout << "ProbabilisticPca<T>::ProbabilisticPca end classes size "
476  // << _classes.size() << std::endl ;
477 
478  }
479 
480  template <class T>
481  double
482  ProbabilisticPca<T>::weight( double norm2, int c, float tolerance )
483  { if( tolerance * norm2 > _discrElements[c].meanNorm() )
484  return 1. ;
485  else
486  return 2 * pow( 0.5, _discrElements[c].meanNorm()
487  / (tolerance * tolerance * norm2 ) ) ;
488  }
489 
490  template <class T>
491  double
492  ProbabilisticPca<T>::wienerFilter( double sigma2, double norm2,
493  double factor )
494  {
495  double res = 1 - factor * factor * sigma2 / norm2 ;
496  if( res < 0. )
497  return 0. ;
498  else
499  return res ;
500  }
501 
502  template <class T>
503  std::vector<double>
505  std::vector<double>& maxProbaByClass )
506  {
507  //std::cout << "ProbabilisticPca<T>::andersonScores classes size "
508  // << _classes.size() << std::endl ;
509 
510  std::vector<double> res( _discrElements.size() ) ;
511  //std::cout << "Nb of classes "<< _discrElements.size()<< std::endl ;
512  double scoreSum = 0. ;
513  for( unsigned int c = 0 ; c < _discrElements.size() ; ++c )
514  {
515  //std::cout << "Class " << c << std::endl ;
516 
517  res[c] = _discrElements[c].posteriorProbability( x, 1 ) * _PIj[c] ;
518  if( res[c] > maxProbaByClass[c] )
519  maxProbaByClass[c] = res[c] ;
520  scoreSum += res[c] ;
521  }
522 
523  for( unsigned int c = 0 ; c < _discrElements.size() ; ++c )
524  {
525  res[c] /= scoreSum ;
526  //std::cout << "\t" << res[c] ;
527  }
528  //std::cout << std::endl ;
529 
530  return res ;
531  }
532 
533  template <class T>
534  std::vector<double>
536  double,
537  std::vector<double>&
538  maxProbaByClass )
539  {
540  std::vector<double> res( _discrElements.size() ) ;
541  /* std::vector<double> normFactor ( _classes.size() ) ; */
542  /* for( int i = 0 ; i < _classes.size() ; ++i ) */
543  /* normFactor[i] = pow( _discrElements[c].noiseVariance()
544  / _discrElements[0].noiseVariance(), x.dimX() / 2. ) ; */
545 
546  for( unsigned int c = 0 ; c < _discrElements.size() ; ++c )
547  {
548  res[c] = _discrElements[c].posteriorProbability( x, 1 ) ;
549  if( res[c] > maxProbaByClass[c] )
550  maxProbaByClass[c] = res[c] ;
551  /* * *//* wienerFilter( _discrElements[c].noiseVariance(),
552  norm2, 5. ) ; */
553  /* _discrElements[c].normalizedPosteriorProbability( x, 1, norm2 )
554  *//* * weight(norm2, c, 3) */ ;
555  }
556  return res ;
557  }
558 
559  template <class T>
560  int
562  {
563  double lnP ;
564  double lnPMin = _discrElements[0].lnPosteriorProbability( x ) ;
565  int indMin = 0 ;
566  for( unsigned int c = 0 ; c < _discrElements.size() ; ++c )
567  {
568  lnP = _discrElements[c].lnPosteriorProbability( x ) ;
569  //std::cout << "Classe " << c << " : lnP = "<< lnP << std::endl ;
570  if( lnP < lnPMin ){
571  lnPMin = lnP ;
572  indMin = c ;
573  }
574  }
575 
576  //std::cout << "lnPMin = "<< lnPMin << std::endl ;
577  return indMin ;
578  }
579 
580 
581  template <class T>
582  bool
584  const AimsData<byte>& mask,
585  AimsData<short>& segmented )
586  {
587  int x, y, z ;
588  if( dynamicImage.dimT() != _data.dimT() )
589  return false ;
590  if( segmented.dimX() != _data.dimX() || segmented.dimY() != _data.dimY()
591  || segmented.dimZ() != _data.dimZ() )
592  {
593  segmented = AimsFastAllocationData<short> ( _data.dimX(), _data.dimY(),
594  _data.dimZ() ) ;
595  segmented.setSizeXYZT( _data.sizeX(), _data.sizeY(), _data.sizeZ(),
596  1.0 ) ;
597  }
598 
599  AimsFastAllocationData<double> indiv( dynamicImage.dimT() ) ;
600  ForEach3d( dynamicImage, x, y, z )
601  {
602  if( mask(x, y, z) )
603  {
604  //std::cout << "x : " << x << ", y : " << y << ", z : " << z
605  // << std::endl ;
606  for( int t = 0 ; t < _data.dimT() ; ++t )
607  indiv(t) = dynamicImage(x, y, z, t ) ;
608  segmented(x, y, z) = affectedTo( indiv ) ;
609  //std::cout << " affected To " << segmented(x, y, z) << std::endl ;
610  }
611  }
612 
613  std::cout << "Probabilistic Pca classification completed" << std::endl ;
614 
615  return true ;
616  }
617 
618  template <class T>
619  bool
621  const AimsData<byte>& mask,
622  AimsData<float>& fuzzySegmented,
623  double /* probaThreshold */,
624  double /* andersonScoreThreshold*/,
625  const AimsData<double>&
626  indivPriorProbabilities )
627  {
628  bool noPriorProba =
629  ( indivPriorProbabilities.dimX( ) == 1
630  && indivPriorProbabilities.dimY( ) == 1
631  && indivPriorProbabilities.dimZ( ) == 1 ) ;
632  int count = 0 ;
633  double uniformPriorProba ;
634  int x, y, z ;
635  if( noPriorProba )
636  {
637  ForEach3d( mask, x, y, z )
638  if( mask( x, y, z ) )
639  ++count ;
640  uniformPriorProba = 1. / count ;
641  }
642  if( dynamicImage.dimT() != _data.dimT() )
643  return false ;
644  if( fuzzySegmented.dimX() != _data.dimX()
645  || fuzzySegmented.dimY() != _data.dimY() ||
646  fuzzySegmented.dimZ() != _data.dimZ()
647  || fuzzySegmented.dimT() != int(_classes.size() ) )
648  {
649  fuzzySegmented = AimsFastAllocationData<float> ( _data.dimX(), _data.dimY(),
650  _data.dimZ(), _classes.size() ) ;
651  fuzzySegmented.setSizeXYZT( _data.sizeX(), _data.sizeY(),
652  _data.sizeZ(), 1.0 ) ;
653  }
654 
655  std::vector<double> probabilityRepartition ;
656  std::vector<double> maxProbaByClass( _discrElements.size(), 0. ) ;
657  std::vector<double> max( _discrElements.size(), 0. ) ;
658  AimsFastAllocationData<double> indiv( dynamicImage.dimT() ) ;
659  AimsFastAllocationData<double> maxOnClass( dynamicImage.dimX(), dynamicImage.dimY(),
660  dynamicImage.dimZ() ) ;
661 
662  ForEach3d( dynamicImage, x, y, z )
663  {
664  if( mask(x, y, z) )
665  {
666  double norm = 0. ;
667  /* std::cout << "x : " << x << ", y : " << y << ", z : " << z
668  << std::endl ; */
669  for( int t = 0 ; t < _data.dimT() ; ++t )
670  {
671  indiv(t) = dynamicImage(x, y, z, t ) ;
672  norm += indiv(t) * indiv(t);
673  }
674  /* probabilityRepartition = posteriorProbabilities( indiv, norm,
675  maxProbaByClass ) ; */
676  probabilityRepartition = andersonScores(indiv, norm,
677  maxProbaByClass) ;
678  for( int c = 0 ; c < fuzzySegmented.dimT() ; ++c )
679  {
680  fuzzySegmented(x, y, z, c) = probabilityRepartition[c] ;
681  if( probabilityRepartition[c] > max[c] )
682  max[c] = probabilityRepartition[c] ;
683  if( probabilityRepartition[c] > maxOnClass(x, y, z) )
684  maxOnClass(x, y, z) = probabilityRepartition[c] ;
685  }
686  }
687  }
688 
689  std::cout << "Max on c = " ;
690  for( int c = 0 ; c < fuzzySegmented.dimT() ; ++c )
691  std::cout << max[c] << " " << std::endl ;
692 
693  Writer< AimsData< float> > fuzzySegmentedW( "fuzzy.ima" ) ;
694  fuzzySegmentedW.write( fuzzySegmented ) ;
695 
696  // float f = 0.000001 ;
697 
698  /* ForEach3d(fuzzySegmented, x, y, z ) */
699  /* if( mask(x, y, z) ) */
700  /* for( int c = 0 ; c < fuzzySegmented.dimT() ; ++c ){ */
701  /* if( fuzzySegmented(x, y, z, c) < f * max[c] ) */
702  /* fuzzySegmented(x, y, z, c) = 0. ; */
703 
704  /* fuzzySegmented(x, y, z, c) /= max[c] ; */
705  /* //std::cout << fuzzySegmented(x, y, z, c) << " " << std::endl ; */
706  /* } */
707  //std::cout << std::endl ;
708 
709 
710  Writer< AimsData< float> > pbByClassW( "probaByClass.ima" ) ;
711  pbByClassW.write( fuzzySegmented ) ;
712 
713  std::cout << "Probabilistic Pca : fuzzy classification completed"
714  << std::endl ;
715 
716  return true ;
717  }
718 
719 
720  template <class T>
721  double
723  {
724  std::vector<double> res( _discrElements.size() ) ;
725  //cout << "Nb of classes "<< _discrElements.size()<< endl ;
726  double scoreSum = 0. ;
727  for( unsigned int c = 0 ; c < _discrElements.size() ; ++c )
728  scoreSum = _discrElements[c].posteriorProbability( x, 1 ) * _PIj[c] ;
729 
730  return scoreSum ;
731  }
732 
733 }
734 
735 #endif
736 
737 
738 
739 
AimsFastAllocationData< double > _EValues
Definition: ppca.h:116
double pX(const AimsData< double > &x)
Definition: ppca_d.h:722
const AimsData< double > & eigenValues()
Definition: ppca.h:87
AimsData< T > & transpose()
float max(float x, float y)
Definition: thickness.h:97
BucketMap< Void > * mask(const BucketMap< Void > &src, const BucketMap< Void > &m, bool intersect=true)
void sort(AimsData< T > &, AimsData< T > &, AimsData< T > *v=NULL)
sort the U and V matrices and the W vector in decreasing order
const AimsData< double > & mean() const
Definition: ppca.h:77
bool classification(const AimsData< T > &dynamicImage, const AimsData< byte > &mask, AimsData< short > &segmented)
Definition: ppca_d.h:583
std::vector< double > andersonScores(const AimsData< double > &x, double px, std::vector< double > &maxProbaByClass)
Definition: ppca_d.h:504
AimsData< T > doit(AimsData< T > &, AimsData< T > *v=NULL)
Singular Value Decomposition.
int dimX() const
#define ForEach3d(thing, x, y, z)
int dimT() const
#define ForEach2d(thing, x, y)
std::vector< double > posteriorProbabilities(const AimsData< double > &x, double px, std::vector< double > &maxProbaByClass)
Definition: ppca_d.h:535
void setReturnType(SVDReturnType rt)
Definition: svd.h:73
Point3df cross(const Point3df &A, const Point3df &B)
Definition: projection.h:234
T maximum() const
virtual bool write(const T &obj, bool ascii=false, const std::string *format=0)
void doIt(const AimsData< T > &individuals, double distanceRef=0.)
Definition: ppca_d.h:87
AimsFastAllocationData< double > _mean
Definition: ppca.h:108
ProbabilisticPca(const AimsData< T > &data, const std::vector< std::list< Point3d > > &classes, int nbOfSignificantEigenValues, const std::vector< double > &PIj=std::vector< double >())
Definition: ppca_d.h:416
AimsData< T > clone() const
AimsData< T > cross(const AimsData< T > &other)
int affectedTo(const AimsData< double > &x)
Definition: ppca_d.h:561
AimsFastAllocationData< double > _x
Definition: ppca.h:111
AimsFastAllocationData< double > _xT
Definition: ppca.h:112
Definition: svd.h:55
bool fuzzyClassification(const AimsData< T > &dynamicImage, const AimsData< byte > &mask, AimsData< float > &fuzzySegmented, double thresholdOnMaxPercentage=0., double andersonScoreThreshold=0.2, const AimsData< double > &indivPriorProbabilities=aims::AimsFastAllocationData< double >())
Definition: ppca_d.h:620
int dimZ() const
AimsFastAllocationData< double > _Wi
Definition: ppca.h:114
void setSizeXYZT(float sizex=1.0f, float sizey=1.0f, float sizez=1.0f, float sizet=1.0f)
int dimY() const
AIMSDATA_API float norm(const Tensor &thing)
AimsFastAllocationData< double > _invCi
Definition: ppca.h:115
double noiseVariance() const
Definition: ppca.h:78
AimsFastAllocationData< double > _EVect
Definition: ppca.h:117
double normFactor() const
Definition: ppca.h:90