ShapeStatistics.hpp

// This file is part of the imaging2 class library.
//
// University of Innsbruck, Infmath Imaging, 2009.
// http://infmath.uibk.ac.at
//
// All rights reserved.


#ifndef SHAPE_SHAPESTATISTICS_H
#define SHAPE_SHAPESTATISTICS_H

#include <shape/ShapeInterface.hpp>
#include <lapack/linear_algebra.hpp>
#include <core/distribution_utilities.hpp>

#include <map>
#include <boost/shared_ptr.hpp>


namespace imaging
{
  namespace shape_statistics_impl
  {
    void constructor(const std::vector<const ShapeInterface*> shape_ptrs, ShapeInterface & mean_shape, ublas::matrix<float_t> & covariance);
    void compute_statistics(const std::map<std::size_t, float_t> & manual_modes, const ublas::matrix<float_t> & covariance, ublas::matrix<float_t> & root_of_covariance);
  }

  template <class shape_t>
  class ShapeStatistics
  {
    const static std::size_t SHAPE_DIMENSION = shape_t::SHAPE_DIMENSION;
    
    ublas::matrix<float_t> _covariance;
    ublas::matrix<float_t> _root_of_covariance;
    
    shape_t _mean_shape;
    
    std::map<std::size_t, float_t> _manual_modes;
     
  public:
  
    ShapeStatistics(const shape_t & mean_shape) : 
      _mean_shape(mean_shape),
      _covariance(ublas::identity_matrix<float_t>(mean_shape.dimension())),
      _root_of_covariance(ublas::identity_matrix<float_t>(mean_shape.dimension())) {}
    
    ShapeStatistics(const std::vector<shape_t> & shapes)
    {
      std::vector<const ShapeInterface *> shape_ptrs(shapes.size());
      for(size_t i = 0; i < shapes.size(); ++i)
        shape_ptrs[i] = &(shapes[i]);
      
      shape_statistics_impl::constructor(shape_ptrs, _mean_shape, _covariance);
    }
    
    ShapeStatistics(const std::vector< boost::shared_ptr<shape_t> > & shapes)
    {
      std::vector<const ShapeInterface *> shape_ptrs(shapes.size());
      for(size_t i = 0; i < shapes.size(); ++i)
        shape_ptrs[i] = shapes[i].get();
      
      shape_statistics_impl::constructor(shape_ptrs, _mean_shape, _covariance);
    }
    
    void set_mode_manually(std::size_t mode, float_t deviation)
    {
      _manual_modes.insert(std::pair<std::size_t, float_t>(mode, deviation));
    }
    
    void compute_statistics()
    {
      shape_statistics_impl::compute_statistics(_manual_modes, _covariance, _root_of_covariance);
    }
    
    void shape_vector(const ublas::vector<float_t> & coefficients, ublas::vector<float_t> & vector) const
    {
      vector.resize(dimension());
      vector = prod(_root_of_covariance, coefficients);
    }
    
    void shape_vector(const ublas::vector<float_t> & coefficients, ublas::vector<float_t> & vector, float_t & squared_distance) const
    {
      vector.resize(dimension());
      vector = prod(_root_of_covariance, coefficients);
      
      squared_distance = 0.0;
      for(std::size_t i = 0; i < coefficients.size(); ++i)
        if(_manual_modes.count(i) == 0)
          squared_distance += square(coefficients(i));
    }
  
    void shape_sample(const ublas::vector<float_t> & coefficients, shape_t & out_shape) const
    {
      float_t squared_distance;
      shape_sample(coefficients, out_shape, squared_distance);
    }

    void shape_sample(const ublas::vector<float_t> & coefficients, shape_t & out_shape, float_t & squared_distance) const
    { 
      ublas::vector<float_t> sample_vector(_mean_shape.dimension());
      shape_vector(coefficients, sample_vector, squared_distance);
  
      _mean_shape.exponential(sample_vector, out_shape);
    }  
    
    void gradient(const ublas::vector<float_t> & coefficients, ublas::vector<float_t> & gradient) const
    { 
      gradient.resize(_mean_shape.dimension());
      
      gradient = 2 * coefficients;
      
      for(std::map<std::size_t, float_t>::const_iterator iter = _manual_modes.begin(); iter != _manual_modes.end(); ++iter)
        gradient(iter->first) = 0.0;
    }
    
    void random_shape_sample(float_t sample_radius, shape_t & out_shape, float_t & squared_distance) const
    {
      ublas::vector<float_t> coefficients(_covariance.size1());
      
      for(std::size_t i = 0; i < coefficients.size(); ++i)
      {
        if(_manual_modes.count(i) > 0)
          coefficients(i) = symmetric_uniform_distribution();
        else
          coefficients(i) = sample_radius * normal_distribution();
      }
          
      shape_sample(coefficients, out_shape, squared_distance);
    }
    
    const shape_t & mean_shape() const { return _mean_shape; }
    
    shape_t & mode_shape(std::size_t mode, float_t offset, shape_t & out_shape) const
    {
      float_t distance_from_mean;
      
      ublas::vector<float_t> coefficients(_covariance.size1());
      
      coefficients = ublas::scalar_vector<float_t>(coefficients.size(), 0.0);
        
      coefficients(mode) = offset;
      
      shape_sample(coefficients, out_shape, distance_from_mean);
      
      return out_shape;
    }
    
    size_t dimension() const
    {
      return _mean_shape.dimension();
    }
    
  };
  
}


#endif

---