primatologist-gpl 6.0.4
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
28namespace aims {
29namespace 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
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::empty
carto::VolumeRef< T > empty()
Definition volume_d.h:22
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