/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/math/interpolation/interpolation_variational_interpolator.h

Go to the documentation of this file.

#ifndef OPENTISSUE_CORE_MATH_INTERPOLATION_INTERPOLATION_VARIATIONAL_INTERPOLATOR_H
#define OPENTISSUE_CORE_MATH_INTERPOLATION_INTERPOLATION_VARIATIONAL_INTERPOLATOR_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/core/math/big/big_lu.h>

#include <boost/lambda/lambda.hpp>

#include <cmath>
#include <algorithm> // for copy
#include <iterator>  // for back_inserter

#ifndef NDEBUG
#include <fstream>
#endif

namespace OpenTissue
{
  namespace interpolation
  {

    namespace detail
    {

      template<typename vector_type>
      inline typename vector_type::value_type  radial_basis_function(vector_type const & x) 
      {
        using std::log;
        typedef typename vector_type::value_type  real_type;

        real_type dot_prod = x*x;
        real_type norm     = sqrt(dot_prod);
        if(norm>0)
          return dot_prod * log( norm );
        return norm;
      }


      template<typename small_vector_type>
      class ImplicitFunction
      {
      public:

        typedef typename small_vector_type::value_type  real_type;

        typedef typename ublas::matrix<real_type>       big_matrix_type;
        typedef typename ublas::vector<real_type>       big_vector_type;
        typedef typename big_matrix_type::size_type     size_type;

        typedef small_vector_type              return_type;

      protected:

        big_vector_type                  m_x;      
        big_vector_type                  m_b;      
        big_matrix_type                  m_A;      

        size_type                        m_N;      
        size_type                        m_dim;    
        size_type                        m_M;      
        std::vector<small_vector_type>   m_c;      

      public:

        template<typename iterator>
          void init( iterator point_begin,iterator point_end, iterator normal_begin,iterator normal_end)
        {
          real_type   k = 0.01;  // TODO: KE 2006-0606 Should depend on geometry size! But I am lazy

          using std::distance;

          size_type  cnt_points  = distance(point_begin,point_end);
          size_type  cnt_normals = distance(normal_begin,normal_end);

          assert(cnt_points==cnt_normals || !"implicit_function::init(): The number of normals and points did not match");

          m_dim = point_begin->size();

          m_N  = cnt_points+cnt_normals;
          m_M  = m_N + m_dim + 1;

          std::copy(point_begin, point_end, std::back_inserter( m_c ) );

          // TODO: Rewrite into a one-liner?
          // hd 2006-06-09 - Using transform and boost::lambda, this is possible.
          //                 However, I couldn't convince boost to deduce the return type, 
          //                 so I had to help it with ret<small_vector_type> that looks
          //                 a bit clumsy.
          //                 Further, I don't know if we are allowed to use using namespace
          //                 like this without polluting the global namespace? Dropping the 
          //                 useing stuff, however, would make the one-liner even uglier...
          {
            using namespace boost::lambda;
            std::transform( point_begin, point_end, normal_begin, std::back_inserter( m_c ), 
              ret<small_vector_type>(_1 + ret<small_vector_type>(_2*k)) );
          }
          //iterator p = point_begin;
          //iterator n = normal_begin;
          //for(;n!=normal_end;++p,++n)
          //  m_c.push_back( (*p) + (*n)*k );

          m_x.resize( m_M );
          m_b.resize( m_M );
          m_A.resize( m_M, m_M );

          m_A.clear();
          m_b.clear();
          m_x.clear();

          //--- take care of b-vector
          for(unsigned int i=0;i<m_M;++i)
          {
            if(i<cnt_points)
              m_b(i) = 0;
            else if (i<m_N)
              m_b(i) = 1;
            else
              m_b(i) = 0;
          }
          //--- take care of phi-part
          for(unsigned int i=0;i<m_N;++i)
            for(unsigned int j=i;j<m_N;++j)
            {
              real_type phi_ij = radial_basis_function( m_c[i] - m_c[j] );
              m_A(i,j) = m_A(j,i)  = phi_ij;
            }
            //--- take care of polynomial part
            for(unsigned int i=0;i<m_N;++i)
            {
              m_A(i, m_N) = m_A(m_N , i) = 1;
              for(unsigned int j=0;j<m_dim;++j)
                m_A(i,j+m_N+1) = m_A(j+m_N+1,i)  = m_c[i](j);
            }
            math::big::lu(m_A,m_x,m_b);
        }

      public:

        real_type operator()( small_vector_type const & x ) const
        {
          real_type value = 0;

          for(unsigned int j=0;j<m_N;++j)
            value += m_x(j)*radial_basis_function(   x - m_c[j] );

          value += m_x(m_N);
          for(unsigned int j=0;j<m_dim;++j)
            value += m_x(m_N + j +1)*(x(j)); 

          return value;
        }

      };

    }// namespace detail


    template<typename iterator>
    inline   detail::ImplicitFunction<typename iterator::value_type> variational_interpolator(
      iterator point_begin
      , iterator point_end
      , iterator normal_begin
      , iterator normal_end
      )
    {
      typedef typename iterator::value_type                 small_vector_type;//--- Hmmm, use iterator_traits for this!
      typedef detail::ImplicitFunction<small_vector_type>   implicit_function_type;
      implicit_function_type  f;
      f.init(point_begin,point_end,normal_begin,normal_end);
      return f;
    }

  } // namespace interpolation

} // namespace OpenTissue

//OPENTISSUE_CORE_MATH_INTERPOLATION_INTERPOLATION_VARIATIONAL_INTERPOLATOR_H
#endif