primatologist-gpl  5.1.2
proba_d.h
Go to the documentation of this file.
1 /* Copyright (C) 2000-2013 CEA
2  *
3  * This software and supporting documentation were developed by
4  * bioPICSEL
5  * CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
6  * 18, route du Panorama
7  * 92265 Fontenay-aux-Roses
8  * France
9  */
10 
11 #ifndef PRIMATOLOGIST_PROBA_PROBA_D_H
12 #define PRIMATOLOGIST_PROBA_PROBA_D_H
13 
14 #include <primatologist/proba/proba.h>
15 #include <primatologist/utility/volume.h>
16 #include <aims/math/mathelem.h>
17 #include <aims/vector/vector.h>
18 #include <cartodata/volume/volume.h>
19 #include <algorithm>
20 #include <cmath>
21 #include <exception>
22 #include <vector>
23 
24 //----------------------------------------------------------------------------
25 // probabilities: various methods
26 //----------------------------------------------------------------------------
27 
28 namespace aims {
29 namespace proba {
30 
31  //--- attenuate ------------------------------------------------------------
32 
33  double attenuate( double proba, float alpha, int n )
34  {
35  if( alpha == 1. )
36  return proba;
37  else if( alpha == 0. )
38  return 1. / n;
39  else
40  return proba * (double)alpha + ( 1. - (double)alpha ) / n;
41  }
42 
43  template <typename P, typename M>
44  carto::VolumeRef<P> &
45  attenuate( carto::VolumeRef<P> & proba,
46  float alpha,
47  const carto::VolumeRef<M> & mask )
48  {
49  if( alpha != 1. )
50  {
51  int n = proba.getSizeT();
52  for( int z = 0; z < proba.getSizeZ(); ++z )
53  for( int y = 0; y < proba.getSizeY(); ++y )
54  for( int x = 0; x < proba.getSizeX(); ++x )
55  if( !mask.get() || mask( x, y, z ) != 0 )
56  for( int t = 0; t < proba.getSizeT(); ++t )
57  proba(x, y, z, t) = attenuate( proba(x, y, z, t), alpha, n );
58  }
59  return proba;
60  }
61 
62  template <typename P>
63  carto::VolumeRef<P> &
64  attenuate( carto::VolumeRef<P> & proba,float alpha )
65  {
66  return attenuate( proba, alpha, vol::empty<int>() );
67  }
68 
69  template <typename P, typename M>
70  carto::VolumeRef<P>
71  newAttenuate( const carto::VolumeRef<P> & proba,
72  float alpha,
73  const carto::VolumeRef<M> & mask )
74  {
75  carto::VolumeRef<P> copy = proba.copy();
76  attenuate( copy, alpha, mask );
77  return copy;
78  }
79 
80  template <typename P>
81  carto::VolumeRef<P>
82  newAttenuate( const carto::VolumeRef<P> & proba, float alpha )
83  {
84  carto::VolumeRef<P> copy = proba.copy();
85  attenuate( copy, alpha );
86  return copy;
87  }
88 
89  template <typename P, typename M>
90  carto::VolumeRef<P> &
91  attenuateZ( carto::VolumeRef<P> & proba,
92  float alpha,
93  const carto::VolumeRef<M> & mask )
94  {
95  if( alpha != 1. )
96  {
97  int n = proba.getSizeT();
98  double sum;
99  for( int z = 0; z < proba.getSizeZ(); ++z )
100  for( int y = 0; y < proba.getSizeY(); ++y )
101  for( int x = 0; x < proba.getSizeX(); ++x )
102  if( !mask.get() || mask( x, y, z ) != 0 )
103  {
104  sum = 0.;
105  for( int t = 0; t < proba.getSizeT(); ++t ) {
106  proba( x, y, z, t ) = std::pow( (double)proba( x, y, z, t ), (double)alpha );
107  sum += proba( x, y, z, t );
108  }
109  if( sum > 0 )
110  for( int t = 0; t < proba.getSizeT(); ++t )
111  proba( x, y, z, t ) /= sum;
112  else
113  for( int t = 0; t < proba.getSizeT(); ++t )
114  proba( x, y, z, t ) = 1./n;
115  }
116  }
117  return proba;
118  }
119 
120  template <typename P>
121  carto::VolumeRef<P> &
122  attenuateZ( carto::VolumeRef<P> & proba, float alpha )
123  {
124  return attenuateZ( proba, alpha, vol::empty<int>() );
125  }
126 
127  template <typename P, typename M>
128  carto::VolumeRef<P>
129  newAttenuateZ( const carto::VolumeRef<P> & proba,
130  float alpha,
131  const carto::VolumeRef<M> & mask )
132  {
133  carto::VolumeRef<P> copy = proba.copy();
134  attenuateZ( copy, alpha, mask );
135  return copy;
136  }
137 
138  template <typename P>
139  carto::VolumeRef<P>
140  newAttenuateZ( const carto::VolumeRef<P> & proba, float alpha )
141  {
142  carto::VolumeRef<P> copy = proba.copy();
143  attenuateZ( copy, alpha );
144  return copy;
145  }
146 
147  //--- normalize ------------------------------------------------------------
148 
149  template <typename P, typename M>
150  carto::VolumeRef<P> &
151  normalize( carto::VolumeRef<P> & proba,
152  const carto::VolumeRef<M> & mask )
153  {
154  double sum;
155  int n = proba.getSizeT();
156  for( int z = 0; z < proba.getSizeZ(); ++z )
157  for( int y = 0; y < proba.getSizeY(); ++y )
158  for( int x = 0; x < proba.getSizeX(); ++x )
159  if( !mask.get() || mask( x, y, z ) != 0 )
160  {
161  sum = 0.;
162  for( int t = 0; t < proba.getSizeT(); ++t )
163  sum += proba( x, y, z, t );
164  if( sum > 0 )
165  for( int t = 0; t < proba.getSizeT(); ++t )
166  proba( x, y, z, t ) /= sum;
167  else
168  for( int t = 0; t < proba.getSizeT(); ++t )
169  proba( x, y, z, t ) = 1./n;
170  }
171  return proba;
172  }
173 
174  template <typename P>
175  carto::VolumeRef<P> &
176  normalize( carto::VolumeRef<P> & proba )
177  {
178  return normalize( proba, vol::empty<int>() );
179  }
180 
181  template <typename P, typename M>
182  carto::VolumeRef<P>
183  newNormalize( const carto::VolumeRef<P> & proba,
184  const carto::VolumeRef<M> & mask )
185  {
186  carto::VolumeRef<P> copy = proba.copy();
187  normalize( copy, mask );
188  return copy;
189  }
190 
191  template <typename P>
192  carto::VolumeRef<P>
193  newNormalize( const carto::VolumeRef<P> & proba )
194  {
195  carto::VolumeRef<P> copy = proba.copy();
196  normalize( copy );
197  return copy;
198  }
199 
200  //--- pdf ------------------------------------------------------------------
201 
202  // TODO: multivariate case
203  template <typename P, typename I, typename M, typename Distrib>
204  carto::VolumeRef<P> &
205  pdf( const std::vector<Distrib> & param,
206  const carto::VolumeRef<I> & values,
207  carto::VolumeRef<P> & output,
208  const carto::VolumeRef<M> & mask )
209  {
210  output.reallocate( values.getSizeX(), values.getSizeY(),
211  values.getSizeZ(), param.size() );
212  for( int z = 0; z < output.getSizeZ(); ++z )
213  for( int y = 0; y < output.getSizeY(); ++y )
214  for( int x = 0; x < output.getSizeX(); ++x )
215  if( !mask.get() || mask(x, y, z) != 0 )
216  {
217  for( int t = 0; t < output.getSizeT(); ++t )
218  output(x, y, z, t) = (P)param[t].pdf( (double)values(x, y, z) );
219  }
220  return output;
221  }
222 
223  template <typename P, typename I, typename Distrib>
224  carto::VolumeRef<P> &
225  pdf( const std::vector<Distrib> & param,
226  const carto::VolumeRef<I> & values,
227  carto::VolumeRef<P> & output )
228  {
229  return pdf( param, values, output, vol::empty<int>() );
230  }
231 
232  template <typename P, typename I, typename M, typename Distrib>
233  carto::VolumeRef<P>
234  newPdf( const std::vector<Distrib> & param,
235  const carto::VolumeRef<I> & values,
236  const carto::VolumeRef<M> & mask )
237  {
238  carto::VolumeRef<P> output;
239  output.copyHeaderFrom( values.header() );
240  pdf( param, values, output, mask );
241  return output;
242  }
243 
244  template <typename P, typename I, typename Distrib>
245  carto::VolumeRef<P>
246  newPdf( const std::vector<Distrib> & param,
247  const carto::VolumeRef<I> & values )
248  {
249  return newPdf<P>( param, values, vol::empty<int>() );
250  }
251 
252  //--- bayes : stationary ---------------------------------------------------
253 
254  template <typename OP, typename P1, typename P2, typename M>
255  carto::VolumeRef<OP> &
256  bayes( const carto::VolumeRef<P1> & conditional,
257  const std::vector<P2> & prior,
258  carto::VolumeRef<OP> & output,
259  const carto::VolumeRef<M> & mask )
260  {
261  if( conditional.getSizeT() != prior.size() )
262  throw std::invalid_argument( "aims::proba::bayes: conditional and prior "
263  "should have similar number of classes" );
264 
265  output.reallocate( conditional.getSizeX(), conditional.getSizeY(),
266  conditional.getSizeZ(), conditional.getSizeT() );
267  double sum;
268  for( int z = 0; z < output.getSizeZ(); ++z )
269  for( int y = 0; y < output.getSizeY(); ++y )
270  for( int x = 0; x < output.getSizeX(); ++x )
271  if( !mask.get() || mask(x, y, z) != 0 )
272  {
273  sum = 0;
274  for( int t = 0; t < output.getSizeT(); ++t )
275  output(x, y, z, t) = conditional( x, y, z, t ) * prior[t];
276  }
277  return normalize( output, mask );
278  }
279 
280  template <typename OP, typename P1, typename P2>
281  carto::VolumeRef<OP> &
282  bayes( const carto::VolumeRef<P1> & conditional,
283  const std::vector<P2> & prior,
284  carto::VolumeRef<OP> & output)
285  {
286  return bayes( conditional, prior, output, vol::empty<int>() );
287  }
288 
289  template <typename OP, typename P1, typename P2, typename M>
290  carto::VolumeRef<OP>
291  newBayes( const carto::VolumeRef<P1> & conditional,
292  const std::vector<P2> & prior,
293  const carto::VolumeRef<M> & mask )
294  {
295  carto::VolumeRef<OP> output;
296  output.copyHeaderFrom( conditional.header() );
297  return bayes( conditional, prior, output, mask );
298  }
299 
300  template <typename OP, typename P1, typename P2>
301  carto::VolumeRef<OP>
302  newBayes( const carto::VolumeRef<P1> & conditional,
303  const std::vector<P2> & prior )
304  {
305  return newBayes<OP>( conditional, prior, vol::empty<int>() );
306  }
307 
308  //--- bayes : non-stationary -----------------------------------------------
309 
310  template <typename OP, typename P1, typename P2, typename M>
311  carto::VolumeRef<OP> &
312  bayes( const carto::VolumeRef<P1> & conditional,
313  const carto::VolumeRef<P2> & prior,
314  carto::VolumeRef<OP> & output,
315  const carto::VolumeRef<M> & mask )
316  {
317  if( conditional.getSizeX() != prior.getSizeX() ||
318  conditional.getSizeY() != prior.getSizeY() ||
319  conditional.getSizeZ() != prior.getSizeZ() ||
320  conditional.getSizeT() != prior.getSizeT() )
321  throw std::invalid_argument( "aims::proba::bayes: conditional and prior "
322  "volumes should have similar dimensions" );
323 
324  output.reallocate( prior.getSizeX(), prior.getSizeY(),
325  prior.getSizeZ(), prior.getSizeT() );
326  double sum;
327  for( int z = 0; z < output.getSizeZ(); ++z )
328  for( int y = 0; y < output.getSizeY(); ++y )
329  for( int x = 0; x < output.getSizeX(); ++x )
330  if( !mask.get() || mask(x, y, z) != 0 )
331  {
332  sum = 0;
333  for( int t = 0; t < output.getSizeT(); ++t )
334  output(x, y, z, t) = conditional( x, y, z, t ) * prior( x, y, z, t );
335  }
336  return normalize( output, mask );
337  }
338 
339  template <typename OP, typename P1, typename P2>
340  carto::VolumeRef<OP> &
341  bayes( const carto::VolumeRef<P1> & conditional,
342  const carto::VolumeRef<P2> & prior,
343  carto::VolumeRef<OP> & output)
344  {
345  return bayes( conditional, prior, output, vol::empty<int>() );
346  }
347 
348  template <typename OP, typename P1, typename P2, typename M>
349  carto::VolumeRef<OP>
350  newBayes( const carto::VolumeRef<P1> & conditional,
351  const carto::VolumeRef<P2> & prior,
352  const carto::VolumeRef<M> & mask )
353  {
354  carto::VolumeRef<OP> output;
355  output.copyHeaderFrom( prior.header() );
356  return bayes( conditional, prior, output, mask );
357  }
358 
359  template <typename OP, typename P1, typename P2>
360  carto::VolumeRef<OP>
361  newBayes( const carto::VolumeRef<P1> & conditional,
362  const carto::VolumeRef<P2> & prior )
363  {
364  return newBayes<OP>( conditional, prior, vol::empty<int>() );
365  }
366 
367  //--- robust bayes : helper ------------------------------------------------
368 
369  namespace internal {
370  template <typename P>
371  struct CompPairLabelProba
372  {
373  bool operator() ( const std::pair<int,P> & a, const std::pair<int,P> & b ) const
374  {
375  return a.second > b.second;
376  }
377  };
378  }
379 
380  //--- robust bayes : stationary --------------------------------------------
381 
382  template <typename OP, typename P1, typename P2, typename M>
383  carto::VolumeRef<OP> &
384  robustBayes( int n_best,
385  const carto::VolumeRef<P1> & conditional,
386  const std::vector<P2> & prior,
387  carto::VolumeRef<OP> & output,
388  const carto::VolumeRef<M> & mask )
389  {
390  if( conditional.getSizeT() != prior.size() )
391  throw std::invalid_argument( "aims::proba::bayes: conditional and prior "
392  "should have similar number of classes" );
393 
394  output.reallocate( conditional.getSizeX(), conditional.getSizeY(),
395  conditional.getSizeZ(), conditional.getSizeT() );
396  std::vector<std::pair<int,P2> > sorted( prior.size() );
397  for( int t = 0; t < prior.size(); ++t ) {
398  sorted[t].first = t;
399  sorted[t].second = prior[t];
400  }
401  std::sort( sorted.begin(), sorted.end(), internal::CompPairLabelProba<P2>() );
402  double sum;
403  for( int z = 0; z < output.getSizeZ(); ++z )
404  for( int y = 0; y < output.getSizeY(); ++y )
405  for( int x = 0; x < output.getSizeX(); ++x )
406  if( !mask.get() || mask(x, y, z) != 0 )
407  {
408  sum = 0;
409  for( int t = 0; t < output.getSizeT() && t < n_best; ++t )
410  output( x, y, z, sorted[t].first ) = conditional( x, y, z, sorted[t].first ) * sorted[t].second;
411  for( int t = n_best; t < output.getSizeT(); ++t )
412  output( x, y, z, sorted[t].first ) = 0.;
413  }
414  return normalize( output, mask );
415  }
416 
417  template <typename OP, typename P1, typename P2>
418  carto::VolumeRef<OP> &
419  robustBayes( int n_best,
420  const carto::VolumeRef<P1> & conditional,
421  const std::vector<P2> & prior,
422  carto::VolumeRef<OP> & output)
423  {
424  return robustBayes( n_best, conditional, prior, output, vol::empty<int>() );
425  }
426 
427  template <typename OP, typename P1, typename P2, typename M>
428  carto::VolumeRef<OP>
429  newRobustBayes( int n_best,
430  const carto::VolumeRef<P1> & conditional,
431  const std::vector<P2> & prior,
432  const carto::VolumeRef<M> & mask )
433  {
434  carto::VolumeRef<OP> output;
435  output.copyHeaderFrom( conditional.header() );
436  return robustBayes( n_best, conditional, prior, output, mask );
437  }
438 
439  template <typename OP, typename P1, typename P2>
440  carto::VolumeRef<OP>
441  newRobustBayes( int n_best,
442  const carto::VolumeRef<P1> & conditional,
443  const std::vector<P2> & prior )
444  {
445  return newRobustBayes<OP>( n_best, conditional, prior, vol::empty<int>() );
446  }
447 
448  //--- robust bayes : non-stationary ----------------------------------------
449 
450  template <typename OP, typename P1, typename P2, typename M>
451  carto::VolumeRef<OP> &
452  robustBayes( int n_best,
453  const carto::VolumeRef<P1> & conditional,
454  const carto::VolumeRef<P2> & prior,
455  carto::VolumeRef<OP> & output,
456  const carto::VolumeRef<M> & mask )
457  {
458  if( conditional.getSizeX() != prior.getSizeX() ||
459  conditional.getSizeY() != prior.getSizeY() ||
460  conditional.getSizeZ() != prior.getSizeZ() ||
461  conditional.getSizeT() != prior.getSizeT() )
462  throw std::invalid_argument( "aims::proba::bayes: conditional and prior "
463  "volumes should have similar dimensions" );
464 
465  output.reallocate( prior.getSizeX(), prior.getSizeY(),
466  prior.getSizeZ(), prior.getSizeT() );
467  double sum;
468  std::vector<std::pair<int,P2> > sorted( prior.getSizeT() );
469  for( int z = 0; z < output.getSizeZ(); ++z )
470  for( int y = 0; y < output.getSizeY(); ++y )
471  for( int x = 0; x < output.getSizeX(); ++x )
472  if( !mask.get() || mask(x, y, z) != 0 )
473  {
474  for( int t = 0; t < output.getSizeT(); ++t ) {
475  sorted[t].first = t;
476  sorted[t].second = prior( x, y, z, t );
477  }
478  std::sort( sorted.begin(), sorted.end(), internal::CompPairLabelProba<P2>() );
479  sum = 0;
480  for( int t = 0; t < output.getSizeT() && t < n_best; ++t )
481  output( x, y, z, sorted[t].first ) = conditional( x, y, z, sorted[t].first ) * sorted[t].second;
482  for( int t = n_best; t < output.getSizeT(); ++t )
483  output( x, y, z, sorted[t].first ) = 0.;
484  }
485  return normalize( output, mask );
486  }
487 
488  template <typename OP, typename P1, typename P2>
489  carto::VolumeRef<OP> &
490  robustBayes( int n_best,
491  const carto::VolumeRef<P1> & conditional,
492  const carto::VolumeRef<P2> & prior,
493  carto::VolumeRef<OP> & output)
494  {
495  return robustBayes( n_best, conditional, prior, output, vol::empty<int>() );
496  }
497 
498  template <typename OP, typename P1, typename P2, typename M>
499  carto::VolumeRef<OP>
500  newRobustBayes( int n_best,
501  const carto::VolumeRef<P1> & conditional,
502  const carto::VolumeRef<P2> & prior,
503  const carto::VolumeRef<M> & mask )
504  {
505  carto::VolumeRef<OP> output;
506  output.copyHeaderFrom( prior.header() );
507  return robustBayes( n_best, conditional, prior, output, mask );
508  }
509 
510  template <typename OP, typename P1, typename P2>
511  carto::VolumeRef<OP>
512  newRobustBayes( int n_best,
513  const carto::VolumeRef<P1> & conditional,
514  const carto::VolumeRef<P2> & prior )
515  {
516  return newRobustBayes<OP>( n_best, conditional, prior, vol::empty<int>() );
517  }
518 
519  //--- maximum likelihood ---------------------------------------------------
520 
521  template <typename C, typename P, typename M>
522  carto::VolumeRef<C> &
523  maximumLikelihood( const carto::VolumeRef<P> & proba,
524  carto::VolumeRef<C> & output,
525  const carto::VolumeRef<M> & mask )
526  {
527  output.reallocate( proba.getSizeX(), proba.getSizeY(),
528  proba.getSizeZ(), 1 );
529  double max, current;
530  for( int z = 0; z < proba.getSizeZ(); ++z )
531  for( int y = 0; y < proba.getSizeY(); ++y )
532  for( int x = 0; x < proba.getSizeX(); ++x )
533  if( !mask.get() || mask(x, y, z) != 0 )
534  {
535  max = -1;
536  for( int t = 0; t < proba.getSizeT(); ++t )
537  {
538  current = proba( x, y, z, t );
539  if( current > max )
540  {
541  max = current;
542  output( x, y, z ) = t;
543  }
544  }
545  }
546  return output;
547  }
548 
549  template <typename C, typename P>
550  carto::VolumeRef<C> &
551  maximumLikelihood( const carto::VolumeRef<P> & proba,
552  carto::VolumeRef<C> & output )
553  {
554  return maximumLikelihood( proba, output, vol::empty<int>() );
555  }
556 
557  template <typename C, typename P, typename M>
558  carto::VolumeRef<C>
559  newMaximumLikelihood( const carto::VolumeRef<P> & proba,
560  const carto::VolumeRef<M> & mask )
561  {
562  carto::VolumeRef<C> output;
563  output.copyHeaderFrom( proba.header() );
564  return maximumLikelihood( proba, output, mask );
565  }
566 
567  template <typename C, typename P>
568  carto::VolumeRef<C>
569  newMaximumLikelihood( const carto::VolumeRef<P> & proba )
570  {
571  return newMaximumLikelihood<C>( proba, vol::empty<int>() );
572  }
573 
574  //--- maximum a posteriori -------------------------------------------------
575 
576  template <typename C, typename P1, typename P2, typename M>
577  carto::VolumeRef<C> &
578  maximumAPosteriori( const carto::VolumeRef<P1> & conditional,
579  const carto::VolumeRef<P2> & prior,
580  carto::VolumeRef<C> & output,
581  const carto::VolumeRef<M> & mask )
582  {
583  if( conditional.getSizeX() != prior.getSizeX() ||
584  conditional.getSizeY() != prior.getSizeY() ||
585  conditional.getSizeZ() != prior.getSizeZ() ||
586  conditional.getSizeT() != prior.getSizeT() )
587  throw std::invalid_argument( "aims::proba::map: conditional and prior "
588  "volumes should have similar dimensions" );
589 
590  output.reallocate( conditional.getSizeX(), conditional.getSizeY(),
591  conditional.getSizeZ(), 1 );
592  double max, current;
593  for( int z = 0; z < conditional.getSizeZ(); ++z )
594  for( int y = 0; y < conditional.getSizeY(); ++y )
595  for( int x = 0; x < conditional.getSizeX(); ++x )
596  if( !mask.get() || mask(x, y, z) != 0 )
597  {
598  max = -1;
599  for( int t = 0; t < conditional.getSizeT(); ++t )
600  {
601  current = conditional( x, y, z, t ) * prior( x, y, z, t );
602  if( current > max )
603  {
604  max = current;
605  output( x, y, z ) = t;
606  }
607  }
608  }
609  return output;
610  }
611 
612  template <typename C, typename P1, typename P2>
613  carto::VolumeRef<C> &
614  maximumAPosteriori( const carto::VolumeRef<P1> & conditional,
615  const carto::VolumeRef<P2> & prior,
616  carto::VolumeRef<C> & output )
617  {
618  return maximumAPosteriori( conditional, prior, output, vol::empty<int>() );
619  }
620 
621  template <typename C, typename P1, typename P2, typename M>
622  carto::VolumeRef<C>
623  newMaximumAPosteriori( const carto::VolumeRef<P1> & conditional,
624  const carto::VolumeRef<P2> & prior,
625  const carto::VolumeRef<M> & mask )
626  {
627  carto::VolumeRef<C> output;
628  output.copyHeaderFrom( conditional.header() );
629  return maximumAPosteriori( conditional, prior, output, mask );
630  }
631 
632  template <typename C, typename P1, typename P2>
633  carto::VolumeRef<C>
634  newMaximumAPosteriori( const carto::VolumeRef<P1> & conditional,
635  const carto::VolumeRef<P2> & prior )
636  {
637  return newMaximumAPosteriori<C>( conditional, prior, vol::empty<int>() );
638  }
639 
640  //--- robust maximum a posteriori ------------------------------------------
641 
642  template <typename C, typename P1, typename P2, typename M>
643  carto::VolumeRef<C> &
644  robustMaximumAPosteriori( int n_best,
645  const carto::VolumeRef<P1> & conditional,
646  const carto::VolumeRef<P2> & prior,
647  carto::VolumeRef<C> & output,
648  const carto::VolumeRef<M> & mask )
649  {
650  if( conditional.getSizeX() != prior.getSizeX() ||
651  conditional.getSizeY() != prior.getSizeY() ||
652  conditional.getSizeZ() != prior.getSizeZ() ||
653  conditional.getSizeT() != prior.getSizeT() )
654  throw std::invalid_argument( "aims::proba::map: conditional and prior "
655  "volumes should have similar dimensions" );
656 
657  output.reallocate( conditional.getSizeX(), conditional.getSizeY(),
658  conditional.getSizeZ(), 1 );
659  double max, current;
660  std::vector<std::pair<int,P2> > sorted( prior.getSizeT() );
661  for( int z = 0; z < conditional.getSizeZ(); ++z )
662  for( int y = 0; y < conditional.getSizeY(); ++y )
663  for( int x = 0; x < conditional.getSizeX(); ++x )
664  if( !mask.get() || mask(x, y, z) != 0 )
665  {
666  for( int t = 0; t < output.getSizeT(); ++t ) {
667  sorted[t].first = t;
668  sorted[t].second = prior( x, y, z, t );
669  }
670  std::sort( sorted.begin(), sorted.end(), internal::CompPairLabelProba<P2>() );
671  max = -1;
672  for( int t = 0; t < output.getSizeT() && t < n_best; ++t )
673  {
674  current = conditional( x, y, z, sorted[t].first ) * sorted[t].second;
675  if( current > max )
676  {
677  max = current;
678  output( x, y, z ) = t;
679  }
680  }
681  }
682  return output;
683  }
684 
685  template <typename C, typename P1, typename P2>
686  carto::VolumeRef<C> &
687  robustMaximumAPosteriori( int n_best,
688  const carto::VolumeRef<P1> & conditional,
689  const carto::VolumeRef<P2> & prior,
690  carto::VolumeRef<C> & output )
691  {
692  return robustMaximumAPosteriori( n_best, conditional, prior, output, vol::empty<int>() );
693  }
694 
695  template <typename C, typename P1, typename P2, typename M>
696  carto::VolumeRef<C>
697  newRobustMaximumAPosteriori( int n_best,
698  const carto::VolumeRef<P1> & conditional,
699  const carto::VolumeRef<P2> & prior,
700  const carto::VolumeRef<M> & mask )
701  {
702  carto::VolumeRef<C> output;
703  output.copyHeaderFrom( conditional.header() );
704  return robustMaximumAPosteriori( n_best, conditional, prior, output, mask );
705  }
706 
707  template <typename C, typename P1, typename P2>
708  carto::VolumeRef<C>
709  newRobustMaximumAPosteriori( int n_best,
710  const carto::VolumeRef<P1> & conditional,
711  const carto::VolumeRef<P2> & prior )
712  {
713  return newRobustMaximumAPosteriori<C>( n_best, conditional, prior, vol::empty<int>() );
714  }
715 
716  //--- log-likelihood -------------------------------------------------------
717 
718  template <typename P, typename M>
719  double
720  logLikelihoodSum( const carto::VolumeRef<P> & conditional,
721  const carto::VolumeRef<M> & mask )
722  {
723  double ll = 0.,
724  sum;
725  for( int z = 0; z < conditional.getSizeZ(); ++z )
726  for( int y = 0; y < conditional.getSizeY(); ++y )
727  for( int x = 0; x < conditional.getSizeX(); ++x )
728  if( !mask.get() || mask(x, y, z) != 0 )
729  {
730  sum = 0.;
731  for( int t = 0; t < conditional.getSizeT(); ++t )
732  if( std::isfinite( conditional( x, y, z, t ) ) )
733  sum += conditional( x, y, z, t );
734  if( sum > 0 )
735  ll += std::log( sum );
736  }
737  return ll;
738  }
739 
740  template <typename P>
741  double
742  logLikelihoodSum( const carto::VolumeRef<P> & proba)
743  {
744  return logLikelihoodSum( proba, vol::empty<int>() );
745  }
746 
747 } // namespace proba
748 } // namespace aims
749 
750 #endif // PRIMATOLOGIST_PROBA_PROBA_D_H
max
float max(float x, float y)
aims::proba::newRobustBayes
carto::VolumeRef< OP > newRobustBayes(int n_best, const carto::VolumeRef< P1 > &conditional, const std::vector< P2 > &prior, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:429
aims::proba::normalize
carto::VolumeRef< P > & normalize(carto::VolumeRef< P > &proba, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:151
aims::proba::attenuate
double attenuate(double proba, float alpha, int n)
Definition: proba_d.h:33
aims::proba::newPdf
carto::VolumeRef< P > newPdf(const std::vector< Distrib > &param, const carto::VolumeRef< I > &values, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:234
aims::proba::logLikelihoodSum
double logLikelihoodSum(const carto::VolumeRef< P > &conditional, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:720
aims::proba::newNormalize
carto::VolumeRef< P > newNormalize(const carto::VolumeRef< P > &proba, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:183
aims::proba::newAttenuateZ
carto::VolumeRef< P > newAttenuateZ(const carto::VolumeRef< P > &proba, float alpha, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:129
aims::proba::newRobustMaximumAPosteriori
carto::VolumeRef< C > newRobustMaximumAPosteriori(int n_best, const carto::VolumeRef< P1 > &conditional, const carto::VolumeRef< P2 > &prior, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:697
aims::proba::newAttenuate
carto::VolumeRef< P > newAttenuate(const carto::VolumeRef< P > &proba, float alpha, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:71
aims::proba::maximumAPosteriori
carto::VolumeRef< C > & maximumAPosteriori(const carto::VolumeRef< P1 > &conditional, const carto::VolumeRef< P2 > &prior, carto::VolumeRef< C > &output, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:578
aims::proba::pdf
carto::VolumeRef< P > & pdf(const std::vector< Distrib > &param, const carto::VolumeRef< I > &values, carto::VolumeRef< P > &output, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:205
aims::proba::newMaximumAPosteriori
carto::VolumeRef< C > newMaximumAPosteriori(const carto::VolumeRef< P1 > &conditional, const carto::VolumeRef< P2 > &prior, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:623
aims::proba::newMaximumLikelihood
carto::VolumeRef< C > newMaximumLikelihood(const carto::VolumeRef< P > &proba, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:559
aims::proba::robustMaximumAPosteriori
carto::VolumeRef< C > & robustMaximumAPosteriori(int n_best, const carto::VolumeRef< P1 > &conditional, const carto::VolumeRef< P2 > &prior, carto::VolumeRef< C > &output, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:644
aims::proba::robustBayes
carto::VolumeRef< OP > & robustBayes(int n_best, const carto::VolumeRef< P1 > &conditional, const std::vector< P2 > &prior, carto::VolumeRef< OP > &output, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:384
aims::proba::maximumLikelihood
carto::VolumeRef< C > & maximumLikelihood(const carto::VolumeRef< P > &proba, carto::VolumeRef< C > &output, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:523
aims::proba::newBayes
carto::VolumeRef< OP > newBayes(const carto::VolumeRef< P1 > &conditional, const std::vector< P2 > &prior, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:291
aims::proba::attenuateZ
carto::VolumeRef< P > & attenuateZ(carto::VolumeRef< P > &proba, float alpha, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:91
aims::proba::bayes
carto::VolumeRef< OP > & bayes(const carto::VolumeRef< P1 > &conditional, const std::vector< P2 > &prior, carto::VolumeRef< OP > &output, const carto::VolumeRef< M > &mask)
Definition: proba_d.h:256
aims::vol::sum
carto::DataTypeTraits< T >::LongType sum(const carto::VolumeRef< T > &in, const carto::VolumeRef< M > &mask)
Definition: volume_d.h:163
aims::proba::internal::CompPairLabelProba::operator()
bool operator()(const std::pair< int, P > &a, const std::pair< int, P > &b) const
Definition: proba_d.h:373