A.I.M.S algorithms


discriminantanalysis_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 
36 #include <aims/math/svd.h>
37 #include <aims/math/mathelem.h>
38 #include <aims/math/random.h>
39 #include <aims/math/gausslu.h>
41 #include <aims/io/writer.h>
42 
43 namespace aims
44 {
45 
47  _significantEV(significantEV), _computed(false), _dataScaleFactor(1.),
48  _probaScaleFactor(1.), _PIj(PIj), _detVarCov(1.),
49  _normFactor(1.), _lnAddFactor(0.)
50 {
51 
52 }
53 
54 template <typename T>
55 void
56 DiscriminantAnalysisElement::doIt( const std::list< Point3d >& selectedPoints,
57  const AimsData<T>& data )
58 {
59  aims::AimsFastAllocationData<T> individuals( std::max(int(selectedPoints.size()), 1), data.dimT() ) ;
60  int i = 0 ;
61  std::list< Point3d >::const_iterator iter( selectedPoints.begin() ), last( selectedPoints.end() ) ;
62  Point3df mean ;
63 
64  _indivPosition = std::vector<Point3d>( selectedPoints.size() ) ;
65  while( iter != last ){
66  _indivPosition[ i ] = *iter ;
67  for(int t = 0 ; t < data.dimT() ; ++t )
68  individuals( i, t ) = T(data( (*iter)[0], (*iter)[1], (*iter)[2], t ) ) ;
69  ++i ; ++iter ;
70  //cout << "Indiv " << i << " = " << *iter << endl ;
71  mean += Point3df( (*iter)[0], (*iter)[1], (*iter)[2] ) ;
72  }
73  //cout << "Center = " << mean / float( individuals.dimX() ) << endl ;
74  doIt( individuals ) ;
75 }
76 
77 template <typename T>
78 void
80 {
81  _computed = true ;
82  int nbFrame = individuals.dimY() ;
83 
84  aims::AimsFastAllocationData<double> centeredIndivMatrix( individuals.dimX(), individuals.dimY() ) ;
86  _mean = 0 ;
87  T val ;
88  for( int ind = 0 ; ind < individuals.dimX() ; ++ind )
89  for( int t = 0 ; t < individuals.dimY() ; ++t ){
90  val = individuals( ind, t ) ;
91  _mean[t] += val ;
92  }
93 
94  //TEMP
95  double meanMean = 0. ;
96  //END TEMP
97 
98  for( int t = 0 ; t < nbFrame ; ++t ){
99  _mean[t] = _mean[t] / centeredIndivMatrix.dimX() ;
100  for( int ind = 0 ; ind < centeredIndivMatrix.dimX() ; ++ind ){
101  centeredIndivMatrix( ind, t ) = individuals( ind, t ) - _mean[t] ;
102  meanMean += centeredIndivMatrix( ind, t ) ;
103  }
104  }
105 
106  // Matrice des correlations
107  aims::AimsFastAllocationData<double> matVarCov(centeredIndivMatrix.dimY(), centeredIndivMatrix.dimY() ) ;
108  int x1, y1;
109  ForEach2d( matVarCov, x1, y1 )
110  {
111  for(int k=0; k < centeredIndivMatrix.dimX() ;++k)
112  matVarCov(x1, y1) += centeredIndivMatrix(k, x1) * centeredIndivMatrix(k, y1);
113  matVarCov(x1, y1) /= centeredIndivMatrix.dimX() - 1 ;
114  }
115 
116  //Writer< AimsData<double> > wriMatVarCov("matVarCov.ima") ;
117  //wriMatVarCov.write( matVarCov ) ;
118 
119  aims::AimsFastAllocationData< double > u = matVarCov.clone(); // because it is replaced on SVD output
120  aims::AimsFastAllocationData< double > v( matVarCov.dimY(), matVarCov.dimY() );
121 
122  AimsSVD< double > svd2;
123  aims::AimsFastAllocationData< double > w = svd2.doit( u, &v );
124  //svd2.sort(u, w, &v) ;
125 
126  AimsData<double> invW( w.dimX(), w.dimY() ) ;
127 
128  // JUST FOR SOCIETY FOR MOLECULAR IMAGING
129 // int indMax = 0 ;
130  double sig2 = 0. ;
131 // for( int i = 1 ; i < w.dimX() ; ++i )
132 // if( w(i, i) > w(indMax, indMax) )
133 // indMax = i ;
134  for( int i = _significantEV ; i < w.dimX() ; ++i ){
135  sig2 += w(i, i) ;
136  }
137  sig2 /= w.dimX() - _significantEV ;
138 
139  std::cout << "signif EV = " << _significantEV << " sig2 = " << sig2 << std::endl ;
140 
141  double lnSig2 = log(sig2) ;
142  double sqrtSig2 = sqrt(sig2) ;
143  int i, n = w.dimX();
144  double lnDetVarCov = 0 ;
145  _detVarCov = 1. ;
146  for ( i = 0; i < _significantEV ; i++ )
147  {
148  lnDetVarCov += 0.5 * log( w( i, i ) ) ;
149  _detVarCov *= sqrt(w( i, i )) ;
150  invW( i, i ) = 1.0f / w( i, i );
151  }
152  for ( i = _significantEV ; i<n; i++ )
153  {
154  lnDetVarCov += 0.5 * lnSig2 ;
155  _detVarCov *= sqrtSig2 ;
156  w( i, i ) = sig2 ;
157  invW( i, i ) = 1.0f / sig2 ;
158  }
159 
160  // JUST FOR SOCIETY FOR MOLECULAR IMAGING END
161 
162  // Set the singular values lower than threshold to zero
163 // double ts = s * w.maximum();
164 // int i, n = w.dimX();
165 // _detVarCov = 1. ;
166 // for ( i=0; i<n; i++ )
167 // {
168 // _detVarCov *= w( i, i ) ;
169 // if ( w( i, i ) > ts ) w( i, i ) = 1.0f / w( i, i );
170 // else w( i, i ) = 0.0f;
171 // }
172 
173 
174  //cout << "Det Var Cov = " << _detVarCov << endl ;
175  if( _detVarCov != 0. ){
176  _invVarCov = (u.clone().transpose()).cross( invW.cross( u ) ) ;
177  _normFactor = 1. / _detVarCov ;
178  _lnAddFactor = -lnDetVarCov + log( _PIj * _PIj ) ;
179 
180  std::cout << "normFactor = " << _normFactor << "\t_lnAddFactor = " << _lnAddFactor << std::endl ;
181  ForEach2d( matVarCov, x1, y1 )
182  if( x1 < y1 ){
183  double val = 0.5 * ( _invVarCov(x1, y1) + _invVarCov(y1, x1) ) ;
184  _invVarCov(x1, y1) = _invVarCov(y1, x1) = val ;
185  }
186  } else{
187  std::cerr << "Bad var cov matrix !" << std::endl ;
188  throw std::runtime_error( "Bad var cov matrix !") ;
189  }
190 }
191 
192 double
194  double ) const
195 {
196  if( ! _computed ){
197  std::cerr << "You must doIt first, parameter not yet computed !"
198  << std::endl ;
199 
200  throw std::runtime_error( "You must doIt first, parameter not yet computed !") ;
201  }
202 // aims::AimsFastAllocationData<double> xCentered(x.dimX()) ;
203 // for( int t = 0 ; t < x.dimX() ; ++t )
204 // xCentered(t) = ( x(t) - _mean(t) ) * _dataScaleFactor ;
205 
206 
207  double distance = 0. ;
208  for(int i = 0 ; i < x.dimX() ; ++i ){
209  double sum = 0. ;
210  for( int j = 0 ; j < x.dimX() ; ++j )
211  sum += _invVarCov(i, j) * ( x(j) - _mean(j) ) ;
212  //double term = ( x(i) - _mean(i) ) * sum ;
213  distance += ( x(i) - _mean(i) ) * sum ;
214  }
215 
216  return _normFactor * exp( -distance ) ;
217 }
218 
219 double
221 {
222  if( ! _computed ){
223  std::cerr <<"You must doIt first, parameter not yet computed !"
224  << std::endl ;
225  throw std::runtime_error( "You must doIt first, parameter not yet computed !") ;
226  }
227 // cout << "DiscriminantAnalysisElement::distance " << endl ;
228 // cout << "inv var cov size = " << _invVarCov.dimX() << " , " << _invVarCov.dimY() << endl ;
229  double distance = 0. ;
230  for(int i = 0 ; i < x.dimX() ; ++i ){
231 // cout << "i = " << i << endl ;
232  double sum = 0. ;
233  for( int j = 0 ; j < x.dimX() ; ++j ){
234 // cout << "\tj = " << j << endl ;
235  sum += _invVarCov(i, j) * ( x(j) - _mean(j) ) ;
236  }
237  //double term = ( x(i) - _mean(i) ) * sum ;
238  distance += ( x(i) - _mean(i) ) * sum ;
239  }
240 // cout << "DiscriminantAnalysisElement::distance end " << endl ;
241 
242  return distance ;
243 }
244 
245 double
247 {
248  if( ! _computed ){
249  std::cerr << "You must doIt first, parameter not yet computed !"
250  << std::endl ;
251  throw std::runtime_error( "You must doIt first, parameter not yet computed !") ;
252  }
253  double distance = 0. ;
254  for(int i = 0 ; i < x.dimX() ; ++i ){
255  double sum = 0. ;
256  for( int j = 0 ; j < x.dimX() ; ++j )
257  sum += _invVarCov(i, j) * ( x(j) - _mean(j) ) ;
258  //double term = ( x(i) - _mean(i) ) * sum ;
259  distance += ( x(i) - _mean(i) ) * sum ;
260  }
261 
262 // double distanceTest = 0. ;
263 // AimsData<double> Ux( x.dimX() ) ;
264 // for(int i = 0 ; i < _U.dimX() ; ++i )
265 // for( int j = 0 ; j < _U.dimY() ; ++j )
266 // Ux(i) += _U(j, i) * (x(j) - _mean(i) ) ;
267 
268 // for( int i = 0 ; i < Ux.dimX() ; ++i )
269 // distanceTest += Ux(i) * Ux(i) * _invEigenValues[i] ;
270 
271 // double distance = 0. ;
272 // for(int i = 0 ; i < x.dimX() ; ++i ){
273 // double sum = 0. ;
274 // for( int j = 0 ; j < x.dimX() ; ++j )
275 // sum += _invVarCov(i, j) * ( x(j) - _mean(j) ) ;
276 // distance += ( x(i) - _mean(i) ) * sum ;
277 // cout << i << "\txi = " << x(i) << "\tmean i = " << _mean(i) <<
278 // " diff = " << x(i) - _mean(i) << " dist = " << distance << endl ;
279 // }
280 
281  return -distance + _lnAddFactor ;
282 }
283 
284 
285 const AimsData<double>&
287 {
288  if( ! _computed ){
289  std::cerr <<"You must doIt first, parameter not yet computed !"
290  << std::endl ;
291  throw std::runtime_error( "You must doIt first, parameter not yet computed !") ;
292  }
293  return _mean ;
294 }
295 
296 template <class T>
298  const AimsData<T>& data,
299  const std::vector< std::list< Point3d > >& classes,
300  int significantEV, const std::vector<double>& PIj )
301  : _significantEV(significantEV <= 0 ? data.dimT() : significantEV),
302  _classes(classes), _data(data), _PIj(PIj)
303 {
304  if( PIj.size() != 0 )
305  if( PIj.size() == classes.size() )
306  _PIj = PIj ;
307  else{
308  std::cerr << "Inconsistant data : number of classes and number of classes prior probabilities are different !" << std::endl ;
309 
310  throw std::runtime_error("Inconsistant data : number of classes and number of classes prior probabilities are different !") ;
311  }
312  else{
313  _PIj.reserve( _classes.size() ) ;
314  for( unsigned int c = 0 ; c < _classes.size() ; ++c )
315  _PIj.push_back( 1./_classes.size() ) ;
316  }
317  //--------------------------------TEMP---------------------------------------------
318 /* AimsData<double> indiv(12000, 3) ; */
319 /* Point3dd mean( 56., 32., 100.) ; */
320 /* Point3dd sigma( 20., 2., 10.) ; */
321 
322 /* for( int ind = 0 ; ind < indiv.dimX() ; ++ind ) */
323 /* for( int t = 0 ; t < indiv.dimY() ; ++t ) */
324 /* indiv(ind, t) = GaussianRandom ( mean[t], sigma[t] ) ; */
325 
326  //--------------------------------TEMP---------------------------------------------
327  _discrElements.reserve( _classes.size() ) ;
328  for( unsigned int c = 0 ; c < _classes.size() ; ++c ){
329  DiscriminantAnalysisElement el( significantEV, _PIj[c] ) ;
330  el.doIt( _classes[c], _data ) ;
331 /* el.doIt( indiv ) ; */
332 
333  std::cout << "Class " << c << " done !" << std::endl ;
334 
335  _discrElements.push_back( el ) ;
336  }
337  std::cout << "Discriminant analysis initization completed" << std::endl ;
338 }
339 
340 template <class T>
341 std::vector<double>
343 {
344  std::vector<double> res( _classes.size() ) ;
345  for( unsigned int c = 0 ; c < _classes.size() ; ++c )
346  res[c] = _discrElements[c].posteriorProbability( x, pX ) ;
347  return res ;
348 }
349 
350 template <class T>
351 std::vector<double>
353 {
354  std::vector<double> res( _classes.size() ) ;
355  double sum = 0. ;
356  for( unsigned int c = 0 ; c < _classes.size() ; ++c ){
357  res[c] = _discrElements[c].posteriorProbability( x, 1. ) ;
358  sum += res[c] ;
359  }
360 
361  for( unsigned int c = 0 ; c < _classes.size() ; ++c )
362  res[c] /= sum ;
363 
364  return res ;
365 }
366 
367 template <class T>
368 int
370 {
371  double lnP ;
372  double lnPMax = _discrElements[0].lnPosteriorProbability( x ) ;
373  int indMax = 0 ;
374  for( unsigned int c = 0 ; c < _classes.size() ; ++c ){
375  lnP = _discrElements[c].lnPosteriorProbability( x ) ;
376  if( lnP > lnPMax ){
377  lnPMax = lnP ;
378  indMax = c ;
379  }
380  }
381 
382  return indMax ;
383 }
384 
385 
386 template <class T>
387 bool
389  const AimsData<byte>& mask,
390  AimsData<short>& segmented )
391 {
392  int x, y, z ;
393  if( dynamicImage.dimT() != _data.dimT() )
394  return false ;
395  if( segmented.dimX() != _data.dimX() && segmented.dimY() != _data.dimY()
396  && segmented.dimZ() != _data.dimZ() )
397  {
398  segmented = aims::AimsFastAllocationData<short> ( _data.dimX(), _data.dimY(), _data.dimZ() );
399  segmented.setSizeXYZT( _data.sizeX(), _data.sizeY(),
400  _data.sizeZ(), 1.0 ) ;
401  }
402 
403  aims::AimsFastAllocationData<double> indiv( dynamicImage.dimT() ) ;
404  ForEach3d( dynamicImage, x, y, z ){
405  if( mask(x, y, z) ) {
406  std::cout << "x : " << x << ", y : " << y << ", z : " << z << std::endl ;
407  for( int t = 0 ; t < _data.dimT() ; ++t )
408  indiv(t) = dynamicImage(x, y, z, t ) ;
409  segmented(x, y, z) = affectedTo( indiv ) ;
410  }
411  }
412 
413  std::cout << "Discriminant analysis classification completed" << std::endl ;
414 
415  return true ;
416 }
417 
418 template <class T>
419 bool
421  const AimsData<byte>& mask,
422  AimsData<float>& fuzzySegmented,
423  const AimsData<double> &
424  indivPriorProbabilities )
425 {
426  bool noPriorProba
427  = ( indivPriorProbabilities.dimX( ) == 1
428  && indivPriorProbabilities.dimY( ) == 1
429  && indivPriorProbabilities.dimZ( ) == 1 );
430  int count = 0 ;
431  double uniformPriorProba ;
432 
433  int x, y, z ;
434  if( noPriorProba )
435  {
436  ForEach3d( mask, x, y, z )
437  if( mask( x, y, z ) )
438  ++count ;
439  uniformPriorProba = 1. / count ;
440  }
441  if( dynamicImage.dimT() != _data.dimT() )
442  return false ;
443  if( fuzzySegmented.dimX() != _data.dimX()
444  && fuzzySegmented.dimY() != _data.dimY()
445  && fuzzySegmented.dimZ() != _data.dimZ() &&
446  fuzzySegmented.dimT() != int(_classes.size() ) )
447  {
448  fuzzySegmented = aims::AimsFastAllocationData<float> ( _data.dimX(), _data.dimY(),
449  _data.dimZ(), _classes.size() ) ;
450  fuzzySegmented.setSizeXYZT( _data.sizeX(), _data.sizeY(), _data.sizeZ(),
451  1.0 ) ;
452  }
453 
454  std::vector<double> probabilityRepartition ;
455  aims::AimsFastAllocationData<double> indiv( dynamicImage.dimT() ) ;
456  ForEach3d( dynamicImage, x, y, z ){
457  if( mask(x, y, z) ) {
458  //cout << "x : " << x << ", y : " << y << ", z : " << z << endl ;
459  for( int t = 0 ; t < _data.dimT() ; ++t )
460  indiv(t) = dynamicImage(x, y, z, t ) ;
461  /* probabilityRepartition = andersonScores( indiv ) ; */
462  probabilityRepartition = posteriorProbabilities( indiv, 1. ) ;
463 
464  for( int c = 0 ; c < fuzzySegmented.dimT() ; ++c )
465  fuzzySegmented(x, y, z, c) = probabilityRepartition[c] ;
466  }
467  }
468 
469  std::cout << "Discriminant analysis : fuzzy classification completed"
470  << std::endl ;
471 
472  return true ;
473 }
474 
475 }
476 
477 
const AimsData< double > & mean() const
float max(float x, float y)
Definition: thickness.h:97
BucketMap< Void > * mask(const BucketMap< Void > &src, const BucketMap< Void > &m, bool intersect=true)
std::vector< double > andersonScores(const AimsData< double > &x)
void doIt(const AimsData< T > &individuals)
int dimX() const
#define ForEach3d(thing, x, y, z)
int dimT() const
#define ForEach2d(thing, x, y)
Point3df cross(const Point3df &A, const Point3df &B)
Definition: projection.h:234
double posteriorProbability(const AimsData< double > &individual, double pX) const
bool fuzzyClassification(const AimsData< T > &dynamicImage, const AimsData< byte > &mask, AimsData< float > &fuzzySegmented, const AimsData< double > &indivPriorProbabilities=AimsData< double >())
std::vector< double > posteriorProbabilities(const AimsData< double > &x, double px)
AimsData< T > clone() const
int affectedTo(const AimsData< double > &x)
double lnPosteriorProbability(const AimsData< double > &individual) const
DiscriminantAnalysisElement(int significantEV=-1, double PIj=1.)
DiscriminantAnalysis(const AimsData< T > &data, const std::vector< std::list< Point3d > > &classes, int significantEV=-1, const std::vector< double > &PIj=std::vector< double >())
Definition: svd.h:55
int dimZ() const
bool classification(const AimsData< T > &dynamicImage, const AimsData< byte > &mask, AimsData< short > &segmented)
AimsFastAllocationData< double > _invVarCov
void setSizeXYZT(float sizex=1.0f, float sizey=1.0f, float sizez=1.0f, float sizet=1.0f)
int dimY() const
double distance(const AimsData< double > &x) const
AimsFastAllocationData< double > _mean