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

Go to the documentation of this file.

#ifndef OPENTISSUE_UTILITY_IMAGE_ANALYSIS_IMAGE_ANALYSIS_COMPUTE_DEFAULT_LABELS_H
#define OPENTISSUE_UTILITY_IMAGE_ANALYSIS_IMAGE_ANALYSIS_COMPUTE_DEFAULT_LABELS_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 <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>

#include <OpenTissue/core/containers/grid/grid.h>
#include <OpenTissue/utility/image_analysis/image_analysis_compute_intensity_gradient_histogram.h>
#include <OpenTissue/core/math/math_basic_types.h>

namespace OpenTissue
{
  namespace image_analysis
  {
    template<typename image_type>
    void compute_default_labels(image_type & U)
    {
      using std::max;
      using std::log;

      OpenTissue::grid::Grid<size_t,math::default_math_types> index_grid;
      size_t histogram[512][128];

      compute_intensity_gradient_histogram(U,512,128,histogram,index_grid);

      //--- Take the logarithm of histogram
      double log_histogram[512][128];
      for(size_t j=0;j<128;++j)
        for(size_t i=0;i<512;++i)
          log_histogram[i][j] = log(1.0*histogram[i][j]+1.0);

      //--- Find maximum value in log histogram
      double max_value = 0;
      for(size_t j=0;j<128;++j)
        for(size_t i=0;i<512;++i)
          max_value = max(max_value,log_histogram[i][j]);

      //--- Write out log histogram as a matlab matrix, great for inspection.
      //--- Try examine the matrix with matlab's imaging toolbox
      std::ofstream file("histogram.m");
      file << " H = [ ";
      for(size_t j=0;j<128;++j)
      {
        for(size_t i=0;i<512;++i)
          file << " " << log_histogram[i][j];
        if(j<127)
          file << ";" << std::endl;
      }
      file << "];" << std::endl;
      file << "imagesc(H);" << std::endl;
      file << "colormap(gray);" << std::endl;
      file.flush();
      file.close();

      //--- Now use histogram to search for features and regions in the image and label them.
      size_t labels[512][128]; //--- Think of this as a color or transfer function.
      for(size_t j=0;j<128;++j)
        for(size_t i=0;i<512;++i)
          labels[i][j] = static_cast<size_t>( log_histogram[i][j] );

      //--- Use index map to assign labels to each voxel of U
      for(size_t x=0; x<U.I(); ++x)
        for(size_t y=0; y<U.J(); ++y)
          for(size_t z=0; z<U.K(); ++z)
          {
            size_t idx = index_grid(x,y,z);
            size_t i_idx = (idx >> 16) & 0x00FF;
            size_t g_idx = idx & 0x00FF;
            U(x,y,z) = labels[i_idx][g_idx];
          }

          std::cout << "default_labelling done..." << std::endl;
    }


  }// end namespace image_analysis

} // namespace OpenTissue

//  OPENTISSUE_UTILITY_IMAGE_ANALYSIS_IMAGE_ANALYSIS_COMPUTE_DEFAULT_LABELS_H
#endif