/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/utility/image_analysis/image_analysis_compute_intensity_gradient_histogram.h

Go to the documentation of this file.

#ifndef OPENTISSUE_UTILITY_IMAGE_ANALYSIS_IMAGE_ANALYSIS_COMPUTE_INTENSITY_GRADIENT_HISTOGRAM_H
#define OPENTISSUE_UTILITY_IMAGE_ANALYSIS_IMAGE_ANALYSIS_COMPUTE_INTENSITY_GRADIENT_HISTOGRAM_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/containers/grid/grid.h>
#include <OpenTissue/core/containers/grid/util/grid_transform.h>
#include <OpenTissue/core/math/math_basic_types.h>

#include <boost/bind.hpp>

#include <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>

namespace OpenTissue
{

  namespace image_analysis
  {

    namespace detail
    {
      template<typename grid_iterator>
      class GridGradientMagnitudeFunctor
      {
      public:

        typename grid_iterator::math_types::real_type 
          operator()( grid_iterator const & iter )
        {
          return OpenTissue::grid::gradient_magnitude(iter);
        }

      };
    }



    template< typename image_type, typename histogram_type >
    void compute_intensity_gradient_histogram(
      image_type const & U
      , unsigned short const & cnt_intensity_buckets
      , unsigned short const & cnt_grad_mag_buckets
      , histogram_type & histogram
      , OpenTissue::grid::Grid<size_t, OpenTissue::math::default_math_types>  & index_grid
      )
    {
      typedef typename image_type::const_index_iterator        intensity_iterator;
      typedef typename image_type::value_type                  intensity_type;

      typedef double                                              real_type;

      typedef OpenTissue::grid::Grid<real_type,math::default_math_types> grid_type;

      typedef typename grid_type::index_iterator                   grid_iterator;
      typedef OpenTissue::math::Vector3<real_type>                vector3_type;

      //--- Clear histogram
      {
        for(size_t g=0;g<cnt_grad_mag_buckets;++g)
          for(size_t i=0;i<cnt_intensity_buckets;++i)
            histogram[i][g] = 0;
      }

      //--- Allocate temporary grid to hold intermediate results
      grid_type tmp;
      tmp.create(U.min_coord(),U.max_coord(),U.I(),U.J(),U.K());

      //--- Allocate index grid
      index_grid.create(U.min_coord(),U.max_coord(),U.I(),U.J(),U.K());

      //--- Generate a gradient magnitude grid
      //std::transform(
      //    U.begin()
      //  , U.end()
      //  , tmp.begin()
      //  , boost::bind( OpenTissue::grid::gradient_magnitude<intensity_iterator> , _1 )
      //  );


      // 2005-09-18 KE: We need a full type instanciation of template
      // function in order to get a function pointer for the bind-function. Yrgh!!!
      OpenTissue::grid::transform(
        U.begin()
        , U.end()
        , tmp.begin()
        //        , boost::bind<real_type>( OpenTissue::grid::gradient_magnitude<intensity_iterator> , _1 )  // works only on .NET 2005
        //        , boost::bind( OpenTissue::grid::gradient_magnitude<intensity_iterator> , _1 )   // works on gcc if function is renamed
        , detail::GridGradientMagnitudeFunctor<intensity_iterator>()  // Ugly workaround, it works:-(
        );

      //--- Compute histogram buckets widths
      double  min_val       = 1.0*min(tmp);
      double  max_val       = 1.0*max(tmp);
      double  min_intensity = 1.0*min(U);
      double  max_intensity = 1.0*max(U);
      double  i_width = (max_intensity-min_intensity)/cnt_intensity_buckets;
      double  g_width  = (max_val-min_val)/cnt_grad_mag_buckets;

      //--- Iterate over voxels and map them into the histogram buckets they belong
      //--- Note we could have computed the gradient maginitude on the fly and
      //--- thus saved allocated space for the temporary grid.
      intensity_iterator Ibegin = U.begin();
      intensity_iterator Iend   = U.end();
      intensity_iterator i;
      grid_iterator       Gbegin = tmp.begin();
      grid_iterator       Gend   = tmp.end();
      grid_iterator       g;
      for(i=Ibegin,g=Gbegin;i!=Iend;++i,++g)
      {
        size_t i_bucket = static_cast<size_t>( std::floor((*i - min_intensity)/i_width) );
        size_t g_bucket = static_cast<size_t>( std::floor((*g - min_val)/g_width) );
        if(i_bucket==cnt_intensity_buckets) --i_bucket; //--- Better safe than sorry:-)
        if(g_bucket==cnt_grad_mag_buckets) --g_bucket;
        assert(i_bucket<cnt_intensity_buckets);
        assert(g_bucket<cnt_grad_mag_buckets);

        index_grid(i.get_index()) =  static_cast<size_t>((i_bucket & 0x00FF) << 16)  | (g_bucket & 0x00FF);
        ++(histogram[i_bucket][g_bucket]);
      }
    }

  }// end namespace image_analysis

} // namespace OpenTissue

//  OPENTISSUE_UTILITY_IMAGE_ANALYSIS_IMAGE_ANALYSIS_COMPUTE_INTENSITY_GRADIENT_HISTOGRAM_H
#endif