A.I.M.S algorithms


ppca.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 #ifndef AIMS_MATH_PPCA_H
36 #define AIMS_MATH_PPCA_H
37 
38 #include <aims/def/general.h>
40 
41 namespace aims
42 {
43 
45  {
46  public:
47  ProbabilisticPcaElement( int nbOfSignificativeEV, double PIj = 1. ) ;
49 
50  template <class T>
51  void doIt( const AimsData<T>& individuals, double distanceRef = 0. ) ;
52 
53  template <class T>
54  void doIt( const std::list< Point3d>& selectedPoints,
55  const AimsData<T>& data, double distanceRef = 0. ) ;
56 
57  // set a priori class probability
58  void setPIj( double PIj )
59  {
60  _PIj = PIj ;
61  if( _computed )
62  _lnAddFactor = log( _detCi / ( _PIj * _PIj ) ) ;
63  }
64 
65  inline double distance( const AimsData<double>& individual ) ;
66 
67  inline double posteriorProbability( const AimsData<double>& individual,
68  double pX ) ;
69 
70 // double normalizedPosteriorProbability( const AimsData<double>& individual,
71 // double pX, double normX ) ;
72 
73  /* for comparison purposes only, because x prior probability is not taken
74  into account */
75  inline double lnPosteriorProbability( const AimsData<double>& individual ) ;
76 
77  const AimsData<double>& mean() const { return _mean ; }
78  double noiseVariance() const
79  { return _Sigma2 ; }
80 
81  template <class T>
82  double noiseVariance( const AimsData<T>& individuals,
83  double& normMean ) const ;
84 
85  double meanNorm() { return _meanNorm ; }
86 
87  const AimsData<double>& eigenValues() { return _EValues ; }
89 
90  double normFactor() const { return _normFactor ; }
91 
92  protected:
93  static double * _exp ;
94  inline double exponential( double x )
95  {
96  if( x < -750 )
97  return 0. ;
98  int approx = int( x * 10.) ;
99  double res = (x * 10. - approx) ;
100  return ( 1. - res ) * _exp[approx] + res * _exp[approx+1] ;
101  }
102  double _distanceRef ;
103  bool _computed ;
104  bool _valid ;
106  double _PIj ;
107 
109  double _meanMean ;
110  double _meanNorm ;
113 
118  double _Sigma2 ;
119 
120  double _detCi ;
121  double _normFactor ;
122  double _lnAddFactor ;
123  };
124 
125 
126  template <class T>
128  {
129  public:
130  ProbabilisticPca( const AimsData<T>& data,
131  const std::vector< std::list <Point3d> >& classes,
132  int nbOfSignificantEigenValues,
133  const std::vector<double>& PIj
134  = std::vector<double>() ) ;
136 
137  std::vector<double>
138  posteriorProbabilities( const AimsData<double>& x, double px,
139  std::vector<double>& maxProbaByClass ) ;
140  std::vector<double>
141  andersonScores( const AimsData<double>& x, double px,
142  std::vector<double>& maxProbaByClass ) ;
143  int affectedTo( const AimsData<double>& x ) ;
144 
145  bool classification( const AimsData<T>& dynamicImage,
146  const AimsData<byte>& mask,
147  AimsData<short>& segmented ) ;
148  bool fuzzyClassification( const AimsData<T>& dynamicImage,
149  const AimsData<byte>& mask,
150  AimsData<float>& fuzzySegmented,
151  double thresholdOnMaxPercentage = 0.,
152  double andersonScoreThreshold = 0.2,
153  const AimsData<double>& indivPriorProbabilities
155 
156 
157  float posteriorProbability( const AimsData<double>& x, float pX, unsigned int classNb )
158  {
159  if( classNb >= _discrElements.size() )
160  throw std::runtime_error("Class number exceeds number of classes") ;
161  else
162  return _discrElements[classNb].posteriorProbability(x, pX) ;
163  }
164 
165  double pX(const AimsData<double>& x ) ;
166 
167  short nbOfClasses() const { return _discrElements.size() ; }
168 
169  private:
170  double weight( double norm2, int classe, float tolerance ) ;
171  double wienerFilter( double sigma2, double norm2, double factor ) ;
172  const std::vector< std::list< Point3d > >& _classes ;
173  const AimsData<T>& _data ;
174  double _distanceRef ;
175  std::vector<double> _PIj ;
176  int _nbOfSignificantEigenValues ;
177  std::vector<ProbabilisticPcaElement> _discrElements ;
178  } ;
179 }
180 
181 double
183 {
184  if( ! _computed )
185  {
186  std::cerr << "ProbabilisticPcaElement::distance : " << std::endl ;
187  throw std::runtime_error( "You must doIt first, parameter not yet "
188  "computed !" ) ;
189  }
190 
191  if( !_valid ){
192  std::cerr << "Invalid ppca ! "<< std::endl ;
193  return 10000000. ;
194  }
195 
196 // for( int t = 0 ; t < indiv.dimX() ; ++t ){
197 // _xT(1, t) = indiv(t) - _mean(t) ;
198 // _x(t) = indiv(t) - _mean(t) ;
199 // }
200 
201  double distance = 0. ;
202  for(int i = 0 ; i < x.dimX() ; ++i ){
203  double sum = 0. ;
204  for( int j = 0 ; j < x.dimX() ; ++j )
205  sum += _invCi(i, j) * ( x(j) - _mean(j) ) ;
206  distance += ( x(i) - _mean(i) ) * sum ;
207  }
208 // double distance = _xT.cross( _invCi.cross( _x ) )(0, 0) ;
209 
210  if( distance < 0. )
211  std::cout << "Distance = " << distance << std::endl ;
212 
213  return distance - _distanceRef ;
214 }
215 
216 
217 double
219  double pX )
220 {
221  if( ! _computed )
222  {
223  std::cerr << "ProbabilisticPcaElement::posteriorProbability : "
224  << std::endl ;
225  throw std::runtime_error( "You must doIt first, parameter not yet "
226  "computed !") ;
227  }
228 
229  if( !_valid )
230  return 0. ;
231 
232 
233  double distance = 0. ;
234  for(int i = 0 ; i < x.dimX() ; ++i ){
235  double sum = 0. ;
236  for( int j = 0 ; j < x.dimX() ; ++j )
237  sum += _invCi(i, j) * ( x(j) - _mean(j) ) ;
238  distance += ( x(i) - _mean(i) ) * sum ;
239  }
240 
241  if( distance < 0. )
242  std::cout << "Distance = " << distance << std::endl ;
243 
244 // double distance = _xT.cross( _invCi.cross( _x ) )(0, 0) ;
245 // double distance = xCentered.clone( ).transpose().cross
246 // ( _invCi.cross( xCentered ) )(0, 0) ;
247  int index = int( 0.5 * (distance - _distanceRef) * 100.
248  + 0.5 ) ;
249  if( index > 100000 )
250  return 0. ;
251  else if( index < 0. )
252  return _normFactor / _exp[ -index ] * _PIj / pX ;
253 
254  return _normFactor * _exp[ index ] * _PIj / pX ;
255 }
256 
257 double
259 {
260  if( ! _computed ){
261  std::cerr << "ProbabilisticPcaElement::lnPosteriorProbability : "
262  << std::endl ;
263  throw std::runtime_error( "You must doIt first, parameter not yet "
264  "computed !") ;
265  }
266 
267  if( !_valid )
269 
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 += _invCi(i, j) * ( x(j) - _mean(j) ) ;
276  distance += ( x(i) - _mean(i) ) * sum ;
277  }
278 // double distance = _xT.cross( _invCi.cross( _x ) )(0, 0) ;
279 
280  if( distance < 0. )
281  std::cout << "Distance = " << distance << std::endl ;
282 
283 // double distance = xCentered.clone( ).transpose().cross
284 // ( _invCi.cross( xCentered ) )(0, 0) ;
285  //std::cout << "distance = " << distance << "\tNormFactor = "
286  // << _lnAddFactor << std::endl ;
287  return ( _distanceRef - distance /*/ (x.dimX() - _nbOfSignificantEV )*/ + _lnAddFactor ) ;
288 }
289 
290 #endif
291 
292 
AimsFastAllocationData< double > _EValues
Definition: ppca.h:116
float posteriorProbability(const AimsData< double > &x, float pX, unsigned int classNb)
Definition: ppca.h:157
double pX(const AimsData< double > &x)
Definition: ppca_d.h:722
const AimsData< double > & eigenValues()
Definition: ppca.h:87
double exponential(double x)
Definition: ppca.h:94
BucketMap< Void > * mask(const BucketMap< Void > &src, const BucketMap< Void > &m, bool intersect=true)
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
int dimX() const
static double * _exp
Definition: ppca.h:93
std::vector< double > posteriorProbabilities(const AimsData< double > &x, double px, std::vector< double > &maxProbaByClass)
Definition: ppca_d.h:535
double lnPosteriorProbability(const AimsData< double > &individual)
Definition: ppca.h:258
short nbOfClasses() const
Definition: ppca.h:167
ProbabilisticPcaElement(int nbOfSignificativeEV, double PIj=1.)
const AimsData< double > & selectedEigenVectors()
Definition: ppca.h:88
double distance(const AimsData< double > &individual)
Definition: ppca.h:182
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
void setPIj(double PIj)
Definition: ppca.h:58
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
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
AimsFastAllocationData< double > _Wi
Definition: ppca.h:114
AimsFastAllocationData< double > _invCi
Definition: ppca.h:115
double noiseVariance() const
Definition: ppca.h:78
AimsFastAllocationData< double > _EVect
Definition: ppca.h:117
double posteriorProbability(const AimsData< double > &individual, double pX)
Definition: ppca.h:218
virtual ~ProbabilisticPcaElement()
Definition: ppca.h:48
double normFactor() const
Definition: ppca.h:90
static _Tp max()