normal_d.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_NORMAL_D_H
#define PRIMATOLOGIST_PROBA_NORMAL_D_H
 
 
//--- primatologist ----------------------------------------------------------
#include <primatologist/proba/normal.h>
//--- aims -------------------------------------------------------------------
#include <aims/math/mathelem.h>
//--- cartodata --------------------------------------------------------------
#include <cartodata/volume/volume.h>
//--- cartobase --------------------------------------------------------------
#include <cartobase/exception/assert.h>
//--- std --------------------------------------------------------------------
#include <cmath>
#include <vector>
//----------------------------------------------------------------------------
 
 
namespace aims {
 
  //==========================================================================
  //
  //                      NORMAL ESTIMATION
  //
  //==========================================================================
 
  template <typename DataType, typename WeightType, typename MaskType>
  NormalEstimation<DataType,WeightType,MaskType>::NormalEstimation():
    _distrib(1),
    _data((carto::Volume<DataType>*)0),
    _weight((carto::Volume<WeightType>*)0),
    _mask((carto::Volume<MaskType>*)0),
    _bessel(false),
    _robust(false),
    _multi(false),
    _vs_data(4,1.),
    _vs_mask(4,1.)
  {}
 
  template <typename DataType, typename WeightType, typename MaskType>
  NormalEstimation<DataType,WeightType,MaskType>::~NormalEstimation()
  {}
 
 
  //--- set volumes ----------------------------------------------------------
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::setData(
    const carto::VolumeRef<DataType> & data )
  {
    _data = data;
    if( _data.get() ) {
      _vs_data.assign(4, 1.);
      _data.header().getProperty("voxel_size", _vs_data);
    }
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::setWeight(
    const carto::VolumeRef<WeightType> & weight )
  {
    _weight = weight;
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::setMask(
    const carto::VolumeRef<MaskType> & mask )
  {
    _mask = mask;
    if( _mask.get() ) {
      _vs_mask.assign(4, 1.);
      _mask.header().getProperty("voxel_size", _vs_mask);
    }
  }
 
  //--- options --------------------------------------------------------------
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::setBessel( bool bessel )
  {
    _bessel = bessel;
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::setRobust( bool robust )
  {
    _robust = robust;
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::setMultiWeights( bool multi )
  {
    _multi = multi;
  }
 
  //--- results --------------------------------------------------------------
 
  template <typename DataType, typename WeightType, typename MaskType>
  double NormalEstimation<DataType,WeightType,MaskType>::getMu() const
  {
    return _distrib[0].getMu();
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  double NormalEstimation<DataType,WeightType,MaskType>::getSigma() const
  {
    return _distrib[0].getSigma();
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  NormalDistribution
  NormalEstimation<DataType,WeightType,MaskType>::getDistribution() const
  {
    return _distrib[0];
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  std::vector<NormalDistribution>
  NormalEstimation<DataType,WeightType,MaskType>::getDistributionVector() const
  {
    return _distrib;
  }
 
  //--- helper --------------------------------------------------------------
 
  template <typename DataType, typename WeightType, typename MaskType>
  inline bool NormalEstimation<DataType,WeightType,MaskType>::inMask(
    long x, long y, long z, long t )
  {
    if( !_mask.get() )
      return true;
 
    if(x)
      x = (long)((double)x * _vs_data[0] / _vs_mask[0]);
    if(y)
      y = (long)((double)y * _vs_data[1] / _vs_mask[1]);
    if(z)
      z = (long)((double)z * _vs_data[2] / _vs_mask[2]);
    if(t)
      t = (long)((double)t * _vs_data[3] / _vs_mask[3]);
 
    return 0 <= x && x < _mask.getSizeX() &&
           0 <= y && y < _mask.getSizeY() &&
           0 <= z && z < _mask.getSizeZ() &&
           0 <= t && ( ( t < _mask.getSizeT()  && _mask(x, y, z, t) ) ||
                       ( _mask.getSizeT() == 1  && _mask(x, y, z) ) );
  }
 
  //--- fit ------------------------------------------------------------------
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::fit()
  {
    if( _multi )
    {
      if( _weight.get() )
        _distrib.resize( _weight.getSizeT() );
      else
        _distrib.resize(1);
      for( int t = 0; t < ( _weight.get() ? _weight.getSizeT() : 1 ); ++t )
        fitOne( _weight.view( Point4di(0, 0, 0, t),
                              Point4di( _weight.getSizeX(),
                                        _weight.getSizeY(),
                                        _weight.getSizeZ(), 1 ) ),
                _distrib[t] );
    }
    else
    {
      _distrib.resize(1);
      fitOne( _weight, _distrib[0] );
    }
  }
 
  template <typename DataType, typename WeightType, typename MaskType>
  void NormalEstimation<DataType,WeightType,MaskType>::fitOne(
    const carto::VolumeRef<WeightType> & volweight,
    NormalDistribution & distrib )
  {
    if( _weight.get() ) {
      ASSERT( _data.getSizeX() <= volweight.getSizeX() );
      ASSERT( _data.getSizeY() <= volweight.getSizeY() );
      ASSERT( _data.getSizeZ() <= volweight.getSizeZ() );
      ASSERT( _data.getSizeT() <= volweight.getSizeT() );
    }
 
    double data, weight, sum_weight = 0.;
    long count = 0;
    double mu = 0.;
    double sigma = 0.;
 
    for( int t = 0; t < _data.getSizeT(); ++t )
    for( int z = 0; z < _data.getSizeZ(); ++z )
    for( int y = 0; y < _data.getSizeY(); ++y )
    for( int x = 0; x < _data.getSizeX(); ++x )
    if( inMask(x, y, z, t) )
    {
      data   = (double)_data(x, y, z, t);
      weight = ( _weight ? (double)volweight(x, y, z, t) : 1. );
      sum_weight += weight;
      mu         += weight * data;
      sigma      += weight * sqr( data );
      ++count;
    }
 
    mu /= sum_weight;
    sigma /= sum_weight;
    sigma -= sqr( mu );
    sigma  = std::sqrt( sigma );
 
    if( _bessel ) {
      double mean_weight = sum_weight / count;
      sigma *= sum_weight / ( sum_weight - mean_weight );
    }
 
    if( _robust )
    {
      double prev_mu = mu, prev_sigma = sigma;
      mu = 0.;
      sigma = 0.;
      sum_weight = 0.;
 
      for( int t = 0; t < _data.getSizeT(); ++t )
      for( int z = 0; z < _data.getSizeZ(); ++z )
      for( int y = 0; y < _data.getSizeY(); ++y )
      for( int x = 0; x < _data.getSizeX(); ++x )
      if( inMask(x, y, z, t) )
      {
        data = (double)_data(x, y, z, t);
        if( std::abs( prev_mu - data ) < prev_sigma * 3. )
        {
          weight = ( _weight ? (double)volweight(x, y, z, t) : 1. );
          sum_weight += weight;
          mu         += weight * data;
          sigma      += weight * aims::sqr( data );
        }
      }
 
      mu /= sum_weight;
      sigma /= sum_weight;
      sigma -= aims::sqr( mu );
      sigma  = std::sqrt( sigma );
 
      if( _bessel ) {
        double mean_weight = sum_weight /
                            ( _data.getSizeT() * _data.getSizeZ() *
                              _data.getSizeY() * _data.getSizeX() );
        sigma *= sum_weight / ( sum_weight - mean_weight );
      }
    }
 
    distrib.setMu( mu );
    distrib.setSigma( sigma );
  }
 
  //==========================================================================
  //
  //                      NORMAL ESTIMATION
  //
  //==========================================================================
 
#define instantiate_normalestim_3( DataType, WeightType, MaskType )          \
  template class NormalEstimation<DataType,WeightType,MaskType>;
 
#define instantiate_normalestim_2( DataType, WeightType )                    \
  instantiate_normalestim_3( DataType, WeightType, int8_t )                  \
  instantiate_normalestim_3( DataType, WeightType, uint8_t )                 \
  instantiate_normalestim_3( DataType, WeightType, int16_t )                 \
  instantiate_normalestim_3( DataType, WeightType, uint16_t )                \
  instantiate_normalestim_3( DataType, WeightType, int32_t )                 \
  instantiate_normalestim_3( DataType, WeightType, uint32_t )                \
  instantiate_normalestim_3( DataType, WeightType, int64_t )                 \
  instantiate_normalestim_3( DataType, WeightType, uint64_t )                \
  instantiate_normalestim_3( DataType, WeightType, float )                   \
  instantiate_normalestim_3( DataType, WeightType, double )
 
#define instantiate_normalestim_1( DataType )                                \
  instantiate_normalestim_2( DataType, int8_t )                              \
  instantiate_normalestim_2( DataType, uint8_t )                             \
  instantiate_normalestim_2( DataType, int16_t )                             \
  instantiate_normalestim_2( DataType, uint16_t )                            \
  instantiate_normalestim_2( DataType, int32_t )                             \
  instantiate_normalestim_2( DataType, uint32_t )                            \
  instantiate_normalestim_2( DataType, int64_t )                             \
  instantiate_normalestim_2( DataType, uint64_t )                            \
  instantiate_normalestim_2( DataType, float )                               \
  instantiate_normalestim_2( DataType, double )
 
} // namespace aims
 
#endif // PRIMATOLOGIST_PROBA_NORMAL_D_H