labels.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */

#ifndef PRIMATOLOGIST_PROBA_LABELS_H
#define PRIMATOLOGIST_PROBA_LABELS_H

#include <primatologist/proba/proba.h>
#include <cartodata/volume/volume.h>
#include <map>
#include <vector>

namespace aims {
namespace proba {

  namespace internal {
    const int TO_REMOVE = -3000;
  }

  template <typename ProbaType>
  void prepare( carto::VolumeRef<ProbaType> & proba,
                carto::VolumeRef<ProbaType> & nbprior,
                std::vector<int>            & labels,
                std::map<int,int>           & indices,
                const std::vector<int>      & ignore,
                int                           nclass = 0 )
  {
    int new_size;
    if( proba.get() )
      new_size = proba.getSizeT();
    else
      new_size = nclass;

    // initialize labels
    int prev_size = labels.size();
    labels.resize(new_size);
    for( int t = prev_size; t < new_size; ++t )
      labels[t] = t;

    // initialize indices
    std::multimap<int,int> pre_indices;
    for( int t = 0; t < labels.size(); ++t )
      pre_indices.insert( std::pair<int,int>(labels[t], t) );

    // manage ignored labels
    std::vector<int>::const_iterator ig;
    for( ig = ignore.begin(); ig != ignore.end(); ++ig )
    {
      std::multimap<int,int>::const_iterator i;
      std::pair<std::multimap<int,int>::const_iterator, std::multimap<int,int>::const_iterator> ii;
      ii = pre_indices.equal_range(*ig);
      for( i = ii.first; i != ii.second; ++i )
      {
        labels[i->second] = internal::TO_REMOVE;
      }
    }

    // compute new size
    new_size = labels.size();
    std::set<int> seen;
    for( int t = 0; t < labels.size(); ++t )
    {
      if( labels[t] == internal::TO_REMOVE )
        --new_size;
      else
      {
        if( seen.count(labels[t]) )
          --new_size;
        else
          seen.insert(labels[t]);
      }
    }

    if( new_size == labels.size() )
      return;

    // compute new proba
    if( proba.get() )
    {
      carto::VolumeRef<ProbaType> new_proba( proba.getSizeX(), proba.getSizeY(),
                                      proba.getSizeZ(), new_size );
      new_proba->copyHeaderFrom( proba.header() );
      new_proba->fill(0.);

      std::vector<int> old_labels = labels;
      labels.resize(0);
      pre_indices.clear();
      for( int t = 0; t < old_labels.size(); ++t )
      {
        if( old_labels[t] != internal::TO_REMOVE )
        {
          if( ! indices.count(old_labels[t]) )
          {
            indices[old_labels[t]] = labels.size();
            labels.push_back( old_labels[t] );
          }
          pre_indices.insert( std::pair<int,int>(indices[old_labels[t]], t) );
          for( int z = 0; z < proba.getSizeZ(); ++z )
          for( int y = 0; y < proba.getSizeY(); ++y )
          for( int x = 0; x < proba.getSizeX(); ++x )
            new_proba(x, y, z, indices[old_labels[t]]) += proba(x, y, z, t);
        }
      }

      proba = new_proba;
    }

    if( nbprior.get() )
    {
      std::multimap<int,int>::const_iterator i, j;
      std::pair<std::multimap<int,int>::const_iterator, std::multimap<int,int>::const_iterator> ii, jj;
      carto::VolumeRef<ProbaType> new_nbprior(new_size, new_size);
      new_nbprior->fill(0.);
      for( int t = 0; t < new_size; ++t ) {
        ii = pre_indices.equal_range(t);
        for( i = ii.first;  i != ii.second; ++i ) {
          for( int u = 0; u < new_size; ++u ) {
            jj = pre_indices.equal_range(u);
            for( j = jj.first; j != jj.second; ++j )
            {
              new_nbprior(t, u) += nbprior(i->second, j->second);
            }
          }
        }
      }

      double sum;
      for( int t = 0; t < new_size; ++t ) {
        sum = 0.;
        for( int u = 0; u < new_size; ++u )
          sum += new_nbprior(t, u);
        for( int u = 0; u < new_size; ++u )
          new_nbprior(t, u) /= sum;
      }

      nbprior = new_nbprior;
    }
  }


  template <typename ProbaType>
  void addProbaCatch( carto::VolumeRef<ProbaType> & proba,
                      carto::VolumeRef<ProbaType> & nbprior,
                      double probacatch = 1. )
  {
    if( probacatch >= 1. )
      probacatch = 0.01;

    if( proba.get() ) {

      int nclass_old = proba.getSizeT();
      int nclass_new = nclass_old + 1;

      proba->reallocate( proba.getSizeX(), proba.getSizeY(), proba.getSizeZ(),
                         proba.getSizeT() + 1, true );
      proba::normalize( proba );

      for( int z = 0; z < proba.getSizeZ(); ++z )
      for( int y = 0; y < proba.getSizeY(); ++y )
      for( int x = 0; x < proba.getSizeX(); ++x )
      {
        proba(x, y, z, proba.getSizeT() - 1) = probacatch;
      }

    }

    if( nbprior.get() ) {

      int nclass_old = proba.getSizeT();
      int nclass_new = nclass_old + 1;

      // pre normalize
      for( int j = 0; j < nbprior.getSizeY(); ++j )
      {
        double sum = 0.;
        for( int i = 0; i < nbprior.getSizeX(); ++i )
          sum += nbprior(i, j);
        for( int i = 0; i < nbprior.getSizeX(); ++i ) {
          if( sum > 0. )
            nbprior(i, j) /= sum; // / (1. - probacatch);
        }
      }

      // save data
      carto::VolumeRef<ProbaType> old_nbprior( new carto::Volume<ProbaType>( *nbprior ) );
      // reallocate
      nbprior->reallocate( nclass_new, nclass_new );
      for( int j = 0; j < nclass_old; ++j )
      for( int i = 0; i < nclass_old; ++i )
        nbprior(i, j) = old_nbprior(i, j);

      for( int i = 0; i < nclass_old; ++i ) {
        nbprior( nclass_old, i ) = 0;
        nbprior( i, nclass_old ) = .25 / nclass_old;
      }
      nbprior( nclass_old, nclass_old ) = 0.75;

    }
  }

} // namespace proba
} // namespace aims

#endif // PRIMATOLOGIST_PROBA_LABELS_H