active_contour.h

Go to the documentation of this file.

//
// Active contour 
//
// Bala Amavasai (bala@amavasai.org)
// Sat Apr  28 10:21:18 BST 2001 
//
// Original C code by bpa & fbc

#ifndef ACTIVE_CONTOUR_H
#define ACTIVE_CONTOUR_H

/* #pragma interface  */

#include <iostream>
#include <fstream>
#include <list>
#include <cstring>
#include <cmath>
#include <values.h>
#include "cartesian.h"
#include "image.h"
#include "image_fileoutput.h"
#include "image_funcs.h"

namespace mimas {
/* correlation tracking method extends image */
template<typename T> 
class activeContour : public image <T>
{
  protected:
    int numPolyPoints; /* number of points on polygon */

  private:
    // void findRectTemplate(Cartesian<int> &topLeft, Cartesian<int> &bottomRight);
    int pointInPolygon(int x, int y);
    image<double> imageEnergy;
    double eucDist(double x1, double y1, double x2, double y2);
    void computeImageEnergy(const const_image_ref<T> &inputImage);
    
  public:
    void pushPolyCoord(int xcoord, int ycoord);
    void showPolyCoords(void);
    void clearPolyCoords(void);
    void sortPolyCoords(void);
    void showPolygonOnImage( image_ref<T> &someImage);
    
    void setPolyCircle(int xcenter, int ycenter, int circleRadius, int numPoints);
    int adaptContour( const const_image_ref<T> &inputImage );
      std::list< Cartesian<int> > adaptContourList( const const_image_ref<T> &trackingImage );
      std::list< Cartesian<int> > polyCoords;
    
    // variables that can be changed by the user
    int maxTotalPointsMoved;
    int maxNumIterations;
    double alpha, beta, gamma;
#ifndef NDEBUG
    bool debug;
#endif
      bool smooth;
    int showPolygonOnImage_graylevel, showPolygonOnImage_thickness;

    // constructor
    activeContour()
    {
      numPolyPoints=0;
      maxTotalPointsMoved=4000;
      maxNumIterations=1000;
      alpha=1.0;
      beta=1.0;
      gamma=5.0;
#ifndef NDEBUG
      debug=false;
#endif
      smooth=false; // smoothing is off by default
      showPolygonOnImage_graylevel=127;
      showPolygonOnImage_thickness=1;
    }
    
    // destructor
    ~activeContour()
    {
      //polyCoords.clear();
    }
};


template<typename T> 
double activeContour<T>::eucDist(double x1, double y1, double x2, double y2)
{
  return sqrt(((x1-x2)*(x1-x2))+((y1-y2)*(y1-y2)));
}


template<typename T> 
int activeContour<T>::adaptContour( const const_image_ref<T> &inputImage )
{
  int j, k, jmin, kmin;
  int pointcount, ptsmoved, totptsmoved=0, itercount=0;
  double Esnake, Emin, E, Eimage, aveDist;
  double Econt[9], Econtmax, Econtmin;
  double Ecurv[9], Ecurvmax, Ecurvmin;
  
  std::list< Cartesian<int> >::iterator curpoint, prevpoint, nextpoint;

  computeImageEnergy(inputImage);

  do
  { 
    /* compute average distance between points */
    aveDist=0.0;
    for(curpoint=polyCoords.begin(), prevpoint=--polyCoords.end() ; curpoint != polyCoords.end(); prevpoint=curpoint,++curpoint)
    {
      aveDist+=eucDist((*curpoint).x,(*curpoint).y,(*prevpoint).x,(*prevpoint).y);
    } 
    aveDist=aveDist/(double)numPolyPoints;

    /* cycle through the snake point list */
    ptsmoved=0;
    Esnake=0.0;

//    for(nextpoint=++polyCoords.begin(), curpoint=polyCoords.begin(), prevpoint=--polyCoords.end() ; 
//      curpoint != polyCoords.end(); prevpoint=curpoint, curpoint=nextpoint, nextpoint++)
    
    prevpoint=--polyCoords.end();
    nextpoint=++polyCoords.begin();
    for (curpoint=polyCoords.begin(); curpoint != polyCoords.end(); prevpoint++, curpoint++, nextpoint++)
    {
        

      if (prevpoint==polyCoords.end()) prevpoint=polyCoords.begin();
      if (nextpoint==polyCoords.end()) nextpoint=polyCoords.begin();
      
      /* calculate Econt and Ecurv */
      for (j=-1, pointcount=0, Econtmax=Ecurvmax=-FLT_MAX, Econtmin=Ecurvmin=FLT_MAX; j<=1; j++)
        for (k=-1; k<=1; k++, pointcount++)
        { 
          /* continuity energy */   
          Econt[pointcount]=fabs(aveDist-
                eucDist( (*prevpoint).x, (*prevpoint).y,
                  (*curpoint).x + k, (*curpoint).y + j ) );

          if (Econt[pointcount]>Econtmax) Econtmax=Econt[pointcount];
          if (Econt[pointcount]<Econtmin) Econtmin=Econt[pointcount];
          
          /* curvature energy */  
          Ecurv[pointcount]=
            ((*prevpoint).x - 2.0 * ((*curpoint).x + k) + (*nextpoint).x) *
            ((*prevpoint).x - 2.0 * ((*curpoint).x + k) + (*nextpoint).x) +
            ((*prevpoint).y - 2.0 * ((*curpoint).y + j) + (*nextpoint).y) *
            ((*prevpoint).y - 2.0 * ((*curpoint).y + j) + (*nextpoint).y);

          if (Ecurv[pointcount]>Ecurvmax) Ecurvmax=Ecurv[pointcount];
          if (Ecurv[pointcount]<Ecurvmin) Ecurvmin=Ecurv[pointcount];
        }
        
        
      /* calculate Emin */
      for (j=-1, pointcount=0, jmin=kmin=0, Emin=FLT_MAX; j<=1; j++)
        for (k=-1; k<=1; k++, pointcount++)
        {
          if (Econtmax != Econtmin) Econt[pointcount]=(Econt[pointcount]-Econtmin)/(Econtmax-Econtmin);
          if (Ecurvmax != Ecurvmin) Ecurv[pointcount]=(Ecurv[pointcount]-Ecurvmin)/(Ecurvmax-Ecurvmin);
          Eimage = imageEnergy.getPixel((*curpoint).x + k, (*curpoint).y + j);
            
          E = alpha * Econt[pointcount] + beta * Ecurv[pointcount] + gamma * Eimage;
              
          if (E < Emin)
          {  Emin = E; jmin = j; kmin = k; }
        }
        
      /* update total snake energy */
      Esnake += Emin;

      /* update snake list point position */
      if ((jmin!=0) || (kmin!=0))
      {
        (*curpoint).y = (*curpoint).y + jmin;
        (*curpoint).x = (*curpoint).x + kmin;
        ptsmoved++; totptsmoved++;
      }     
    }
#ifndef NDEBUG
    if (debug) std::cerr << itercount << "  " << totptsmoved   << std::endl;
#endif


    // If debugging then dump one frame per iteration
#ifndef NDEBUG
    if (debug)
    {
        image<int> debugimage( inputImage );
    
      showPolygonOnImage(debugimage);
                  char framefilename[80];
        sprintf(framefilename,"frame%.3d.pgm", itercount);
      std::ofstream file( framefilename, std::ios::binary );
      file << debugimage;
      debugimage.clear();
    }
#endif

  } while ((++itercount < maxNumIterations) && (totptsmoved < maxTotalPointsMoved) && (ptsmoved > 0));
    
  // ensure that snake is within image boundary
  std::list< Cartesian<int> > ::iterator m;
  for(m=polyCoords.begin(); m != polyCoords.end(); m++)
  {
    if ((*m).x<0) (*m).x=0;
    if ((*m).x>=inputImage.getWidth()) (*m).x=inputImage.getWidth()-1;
    if ((*m).y<0) (*m).y=0;
    if ((*m).y>=inputImage.getHeight()) (*m).y=inputImage.getHeight()-1;
  }
      


  return (itercount);
}


template<typename T> 
void activeContour<T>::setPolyCircle(int xcenter, int ycenter, int circleRadius, int numPoints)
{
  double angle;
  int count;
    

  for (angle=0.0, count=0; angle<360.0; angle+=(360.0/numPoints))
  {
    /* random circle
    snakelist[count][YCOORD]=(int)(ycenter-circleRadius*sin(angle*M_PI/180.0))+(5.0*drand48());
    snakelist[count][XCOORD]=(int)(xcenter+circleRadius*cos(angle*M_PI/180.0))+(5.0*drand48());
    */
        
    pushPolyCoord((int)(xcenter+circleRadius*cos(angle*M_PI/180.0)),
      (int)(ycenter-circleRadius*sin(angle*M_PI/180.0)));
    
    if (++count==numPoints) break;
  }
}


template<typename T>
void activeContour<T>::computeImageEnergy(const const_image_ref<T> &inputImage)
{
   if (smooth)
     imageEnergy = edgeSobelSqr( softenHeavy( inputImage ) );
   else
     imageEnergy = edgeSobelSqr( inputImage );

   double
     minval = min_val( imageEnergy ),
     maxval = max_val( imageEnergy );

  /* this is a hack. see Williams and Shah paper */
  if ((maxval-minval)<5.0) maxval=minval+5.0;
  
  /* normalise image pixels between 0 and 1 */
   imageEnergy = ( imageEnergy - minval ) * ( 1.0 / ( maxval - minval ) );
}




template<typename T>
void activeContour<T>::pushPolyCoord(int xcoord, int ycoord)
{
  //polyCoords.push_front(xcoord,ycoord);
  
  Cartesian<int> tempcoord;
  tempcoord.x=xcoord;
  tempcoord.y=ycoord;

  polyCoords.push_back(tempcoord);

  numPolyPoints++;
}


template<typename T>
void activeContour<T>::showPolyCoords(void)
{
  std::list< Cartesian<int> > ::iterator i;

  for(i=polyCoords.begin(); i != polyCoords.end(); ++i) 
    std::cout << *i << std::endl; 
}

template<typename T>
void activeContour<T>::clearPolyCoords(void)
{
  polyCoords.clear();
}


template<typename T>
void activeContour<T>::sortPolyCoords(void)
{
  polyCoords.sort();
}

/*
template<typename T>
void activeContour<T>::findRectTemplate(Cartesian<int> &topLeft, Cartesian<int> &bottomRight)
{
  std::list< Cartesian<int> > ::iterator i;

  int xmax=-1, ymax=-1, xmin=INT_MAX, ymin=INT_MAX;

  for(i=polyCoords.begin(); i != polyCoords.end(); ++i)
  {
    xmax=( (xmax<(*i).x) ? (*i).x:xmax );
    ymax=( (ymax<(*i).y) ? (*i).y:ymax );
    xmin=( (xmin>(*i).x) ? (*i).x:xmin );
    ymin=( (ymin>(*i).y) ? (*i).y:ymin );   
  }
  templateRows=xmax-xmin;
  templateCols=ymax-ymin;

  topLeft.x=xmin; topLeft.y=ymin;
  bottomRight.x=xmax; bottomRight.y=ymax;

  offsetForFirstCoord.x=(*polyCoords.begin()).x-xmin;
  offsetForFirstCoord.y=(*polyCoords.begin()).y-ymin;
} */


template<typename T>
int activeContour<T>::pointInPolygon(int x, int y)
{
  std::list< Cartesian<int> > ::iterator i,j;
  int c = 0;

  // ugly logic in the following line but I can't see how to achieve it otherwise!
  for(i=polyCoords.begin(), j=--polyCoords.end() ; i != polyCoords.end(); j=i,++i)
  {
    if (((((*i).y<=y) && (y<(*j).y)) ||
      (((*j).y<=y) && (y<(*i).y))) &&
      (x < ((*j).x - (*i).x) * (y - (*i).y) / ((*j).y - (*i).y) + (*i).x))
      c = !c;
  }

  return c;
}


template<typename T>
void activeContour<T>::showPolygonOnImage( image_ref<T> &someImage)
{
  std::list< Cartesian<int> > ::iterator i,j;
  for(i=polyCoords.begin(), j=--polyCoords.end() ; i != polyCoords.end(); j=i,++i)
  if (showPolygonOnImage_thickness == 1)
    drawLine(someImage, (*i).x, (*i).y, (*j).x, (*j).y, showPolygonOnImage_graylevel);
  else
    drawThickLine(someImage, (*i).x, (*i).y, (*j).x, (*j).y, showPolygonOnImage_graylevel, 5);
}


// Run adaptContour but return list of points
template<typename T>
std::list< Cartesian<int> > activeContour<T>::adaptContourList( const const_image_ref<T> &trackingImage )
{
  std::list< Cartesian<int> > ::iterator i;
  std::list< Cartesian<int> >  outCoordList;
  Cartesian<int> tempcoord;

  // call the snake algorithm
  adaptContour(trackingImage);

  for(i=polyCoords.begin(); i != polyCoords.end(); i++)
  {
    tempcoord.x=(*i).x;
    tempcoord.y=(*i).y;
    outCoordList.push_back(tempcoord);      
  }

/*
  for(k=polyCoords.begin(); k != polyCoords.end(); ++k) 
  {
    Cartesian<int> tempcoord;

    tempcoord.x=(*k).x-(*polyCoords.begin()).x+outCoord.x+offsetForFirstCoord.x;
    tempcoord.y=(*k).y-(*polyCoords.begin()).y+outCoord.y+offsetForFirstCoord.y;

    outCoordList.push_back(tempcoord);      
  }
*/
  
  return outCoordList;
}   

}

#endif

---