contourLine.h

Go to the documentation of this file.

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//

#ifndef BALL_DATATYPE_CONTOURLINE_H
#define BALL_DATATYPE_CONTOURLINE_H

#ifndef BALL_COMMON_H
# include <BALL/common.h>
#endif

#ifndef BALL_DATATYPE_REGULARDATA2D_H
# include <BALL/DATATYPE/regularData2D.h>
#endif

#include <vector>

namespace BALL
{
  // First I define some macros needed for the marching cube-algorithm. 
  // The names come from the number associated with the different corners of the square.
  #define INTERPOL12 { \
            vec = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x, act_cell_y)));\
            d1  = from[act_cell_x + act_cell_y*(number_of_cells_x+1)];\
            d2  = from[act_cell_x + 1 + act_cell_y*(number_of_cells_x+1)];\
            vec2 = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x + 1, act_cell_y + 1)));\
            slope = (d2 - d1) / (vec2.x - vec.x);\
            vec.x += (threshold - d1)/slope;\
            data_.push_back(vec);\
  } 

  #define INTERPOL18 { \
            vec = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x, act_cell_y)));\
            d1  = from[act_cell_x + act_cell_y*(number_of_cells_x+1)];\
            d2  = from[act_cell_x + (act_cell_y+1)*(number_of_cells_x+1)];\
            vec2 = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x, act_cell_y+1)));\
            slope = (d2 - d1) / (vec2.y - vec.y);\
            vec.y += (threshold - d1)/slope;\
            data_.push_back(vec);\
  }

  #define INTERPOL24 {  \
            vec = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x+1, act_cell_y)));\
            d1  = from[act_cell_x+1 + act_cell_y*(number_of_cells_x+1)];\
            d2  = from[act_cell_x+1 + (act_cell_y+1)*(number_of_cells_x+1)];\
            vec2 = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x+1, act_cell_y+1)));\
            slope = (d2 - d1) / (vec2.y - vec.y);\
            vec.y += (threshold - d1)/slope;\
            data_.push_back(vec);\
  }

  // is it vec.x += or vec.y += ...?
  #define INTERPOL48 {  \
            vec = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x+1, act_cell_y+1)));\
            d1  = from[act_cell_x+1 + (act_cell_y+2)*(number_of_cells_x+1)];\
            d2  = from[act_cell_x   + (act_cell_y+1)*(number_of_cells_x+1)];\
            vec2 = from.getCoordinates(from.getClosestIndex(Vector2(act_cell_x, act_cell_y+1)));\
            slope = (d2 - d1) / (vec2.x - vec.x);\
            vec.x += (threshold - d1)/slope;\
            data_.push_back(vec);\
  }


  template <typename T>  
  class TContourLine
  {
    public:

      
      typedef Vector2 PointType;
      
      typedef std::vector<PointType> VectorType;

      
      TContourLine(T height = 0);

      TContourLine(const TContourLine& copyTContourLine);

      virtual ~TContourLine();

      void createContourLine(TRegularData2D<T>& from);

      void interpol12();
      void interpol18();
      void interpol24();
      void interpol48();

      
      const TContourLine& operator = (const TContourLine& assigTContourLine);

      virtual void clear();


      bool operator == (const TContourLine& compTContourLine) const;



      bool getNextPoint(PointType &p);

      void resetCounter();


      //  private:
      T height_;
      VectorType data_;
      typename VectorType::iterator it_;
      Position index_;
    };

    typedef TContourLine<float> ContourLine;

    template <typename T>
      TContourLine<T>::TContourLine(T height)
      : height_(height), 
        index_(0)
    {
    }

    template <typename T>
    TContourLine<T>::~TContourLine()
    {
    }

    template <typename T>
    TContourLine<T>::TContourLine(const TContourLine<T>& from)
      : height_(from.height_),
        data_(from.data_),
        it_(from.it_),
        index_(from.index_)
    {
    }

    template <typename T>
    void TContourLine<T>::clear()
    {
      data_.clear();
      it_=data_.begin();
      index_ = 0;
    }

    template <typename T>
    const TContourLine<T>& TContourLine<T>::operator = (const TContourLine<T>& data)
    {
      data_ = data.data_;
      height_ = data.height_;
      it_ = data.it_;
      index_ = data.index_;

      return *this;
    }

    template <typename T>
    bool TContourLine<T>:: operator == (const TContourLine<T>& data) const
    {
      return ((height_    == data.height_)
               && (data_  == data.data_)
               && (it_    == data.it_)
               && (index_ == data.index_));
    }

    template <typename T>
    void TContourLine<T>::createContourLine(TRegularData2D<T>& from)
    {
      // This function uses a "marching cubes"-style algorithm to determine the contour-lines.
      //Size number_of_cells;
      Size number_of_cells_x;
      Size number_of_cells_y;
      Position act_cell_x;
      Position act_cell_y;
      PointType vec, vec2;
      double d1, d2, slope;
      double threshold = height_;

      number_of_cells_x = (Size) from.getSize().x - 1;
      number_of_cells_y = (Size) from.getSize().y - 1;
      
      for (act_cell_y = 0; act_cell_y < number_of_cells_y; act_cell_y++)
      {
        for (act_cell_x = 0; act_cell_x < number_of_cells_x; act_cell_x++)
        {
          // First we have to find out the topology of the actual square.
          int topology = 0;
          
          if (from[act_cell_x + act_cell_y * (number_of_cells_x+1)] > threshold)
          {
            topology |= 1;
          }
          if (from[act_cell_x + 1 + act_cell_y * (number_of_cells_x+1)] > threshold)
          {
            topology |= 2;
          }
          if (from[act_cell_x + 1 + (act_cell_y + 1)*(number_of_cells_x + 1)] > threshold)
          {
            topology |= 4;
          }
          if (from[act_cell_x + (act_cell_y + 1) * (number_of_cells_x + 1)] > threshold)
          {
            topology |= 8;
          }
          // now we can use this information to compute the contour-line.
          switch (topology)
          {
            // no cut of contour-line here
            case 0  :
            case 15 : break;
              
            // Line from upper left to lower right
            case 1  : 
            case 14 : INTERPOL18
                      INTERPOL12
                      break;
            
            case 4  :
            case 11 : INTERPOL48
                      INTERPOL24
                      break;
            
            // Line from upper right to lower left
            case 2  :
            case 13 : INTERPOL12
                      INTERPOL24
                      break;

            case 8  :
            case 7  : INTERPOL18
                      INTERPOL48
                      break;
            
            // Line through the middle (upwards)
            case 9  :
            case 6  : INTERPOL12
                      INTERPOL48
                      break;

            // Line through the middle (left to right)
            case 3  :
            case 12 : INTERPOL18
                      INTERPOL24
                      break;

            // Two lines from upper right to lower left
            case 10 : INTERPOL18
                      INTERPOL12
                      INTERPOL48
                      INTERPOL24
                      break;
          
            // Two lines from upper left to lower right
            case 5  : INTERPOL12
                      INTERPOL24
                      INTERPOL18
                      INTERPOL48
                      break;
          };
        }
      }
      index_ = 0;
      it_ = data_.begin();
    }

    template <typename T>
    bool TContourLine<T>::getNextPoint(typename TContourLine<T>::PointType &p)
    {
      if (index_ < data_.size()) 
      {
        p = *it_;
        index_++;
        it_++;
        return true;
      } 
      else 
      {
        return false;
      }
    }

    template <typename T>
    void TContourLine<T>::resetCounter()
    {
      it_ = data_.begin();
      index_ = 0;
    }
}
#endif