canny.h

Go to the documentation of this file.

// WARNING! DO NOT USE THIS CLASS
// This class has been replaced with the much better canny_tool class
// canny will be deprecated in the next release

#ifndef CANNY_H
#define CANNY_H

#include <iostream> 
#include <fstream>
#include <cmath>
#include "cartesian.h"
#include "image.h"
#include "image_draw.h"

namespace mimas {
  
// $Header: /cvs/mimas2/include/canny.h,v 1.1.1.1 2005/08/09 15:37:45 engmb Exp $

/* Scale floating point magnitudes and angles to 8 bits */
// #define ORI_SCALE 40.0
#define MAG_SCALE 20.0

template<typename T>
class canny : public object
            //perhaps public virtual drawable one day
{
protected:
    float norm (float x, float y);
    float gauss(float x, float sigma);
    float dGauss (float x, float sigma);
    float meanGauss (float x, float sigma);
    void separable_convolution( const image< int > &im,
                                const std::vector< float > &gau,
                                image< float > &smx,
                                image< float > &smy );
    void dxy_separable_convolution( const std::vector< float > &dgau,
                                    const image< float > &smx,
                                    const image< float > &smy,
                                    image< float > &dx,
                                    image< float > &dy );
    void nonmax_suppress( const image< float > &dx,
                          const image< float > &dy,
                          image< float > &mag );
    void estimate_thresh( const image< float > &mag, int fringe );
    float distancePointToPoint(int ax, int ay, int bx, int by);
    float distancePointToLine(int ax, int ay, int bx, int by, int px, int py);
    void hysteresis( const image< float > &mag, image< int > &im,
                     int fringe );
    int trace( const image< float > &mag, image< int > &im, int i, int j);

    /*
    image<int> ori, mag, im;
    image<int> traced;
    image<float> smx, smy, dx, dy; */


public:
    float ratio; /* Fraction of pixels that should be above the HIGH threshold */
    float sd; /* Gaussian standard deviation */
    int highthres, lowthres; /* high and low thresholds */
    bool autothreshold; /* autothresholding */

    void setCannyPixHighThresRatio(double tempval);
    double getCannyPixHighThresRatio(void);
    void setCannyStdDev(double tempval);
    double getCannyStdDev(void);
    void setCannyHighThres(int tempval);
    int getCannyHighThres(void);
    void setCannyLowThres(int tempval);
    int getCannyLowThres(void);
    void setCannyAutoThres(bool tempval);
    bool getCannyAutoThres(void);
    void getEdges(image<T> &inputImage, image<T> &outimage);
    void getEdges(image<T> &image);


    canny()
    {
        ratio=0.1;
        sd=1.0;
        highthres=127;
        lowthres=63;
        autothreshold=true;
    }
};


template<typename T>
void canny<T>::setCannyPixHighThresRatio(double tempval)
{
    ratio=tempval;
}

template<typename T>
double canny<T>::getCannyPixHighThresRatio(void)
{
    return ratio;
}

template<typename T>
void canny<T>::setCannyStdDev(double tempval)
{
    sd=tempval;
}

template<typename T>
double canny<T>::getCannyStdDev(void)
{
    return sd;
}

template<typename T>
void canny<T>::setCannyHighThres(int tempval)
{
    highthres=tempval;
}

template<typename T>
int canny<T>::getCannyHighThres(void)
{
    return highthres;
}

template<typename T>
void canny<T>::setCannyLowThres(int tempval)
{
    lowthres=tempval;
}

template<typename T>
int canny<T>::getCannyLowThres(void)
{
    return lowthres;
}


template<typename T>
void canny<T>::setCannyAutoThres(bool tempval)
{
    autothreshold=tempval;
}

template<typename T>
bool canny<T>::getCannyAutoThres(void)
{
    return autothreshold;
}

template<typename T>
float canny<T>::norm (float x, float y)
{
    return (float) sqrt ( (double)(x*x + y*y) );
}

template<typename T>
void canny<T>::getEdges(image<T> &inputImage, image<T> &outimage)
{
    std::vector< float > gau, dgau;

    image< int > im( inputImage );

    /* Create a Gaussian and a derivative of Gaussian filter mask */
    while ( true ) {
        int i = gau.size();
        gau.push_back( meanGauss ((float)i, sd) );
        if (gau[gau.size()-1] < 0.005)
          break;
        dgau.push_back( dGauss ((float)i, sd) );
    }

    //cout << "Smoothing with a Gaussian (width = %d) ..." << std::endl;
    /* Convolution of source image with a Gaussian in X and Y directions  */
    image< float > smx, smy;
    separable_convolution( im, gau, smx, smy );

    image< float > dx, dy;
    /* Now convolve smoothed data with a derivative */
    //cout << "Convolution with the derivative of a Gaussian..." << std::endl;
    dxy_separable_convolution( dgau, smx, smy, dx, dy );

    image< float > mag; mag.init( im.getWidth(), im.getHeight() );

    /* Non-maximum suppression - edge pixels should be a local max */
    nonmax_suppress( dx, dy, mag );
    hysteresis( mag, im, gau.size() / 2 );

    outimage.init(im.getWidth(),im.getHeight());
    for ( int i=1; i<im.getHeight()-1; i++)
        for ( int j=1; j<im.getWidth()-1; j++)
          outimage.pixel(j,i)=im.pixel(j,i);
    
    // set boundaries to zero due to an artefact of this implement of Canny
    drawLine(outimage,0,0,outimage.getWidth()-1,0,T(0));
    drawLine(outimage, outimage.getWidth()-1,0,outimage.getWidth()-1,outimage.getHeight()-1,T(0));
    drawLine(outimage, 0,outimage.getHeight()-1, outimage.getWidth()-1,
     outimage.getHeight()-1,T(0));
    drawLine(outimage,0,0,0,outimage.getHeight()-1,T(0));
    

    //dx.clear();
    //dy.clear();
}

template<typename T>
void canny<T>::getEdges(image<T> &image)
{
    getEdges(image,image);
}
  
 /*      Gaussian        */
template<typename T>
float canny<T>::gauss(float x, float sigma)
{
    float xx;

    if (sigma == 0) return 0.0;
    xx = (float)exp((double) ((-x*x)/(2*sigma*sigma)));
    return xx;
}

template<typename T>
float canny<T>::meanGauss (float x, float sigma)
{
    float z;

    z = (gauss(x,sigma)+gauss(x+0.5,sigma)+gauss(x-0.5,sigma))/3.0;
    z = z/(M_PI*2.0*sigma*sigma);
    return z;
}

/*      First derivative of Gaussian    */
template<typename T>
float canny<T>::dGauss (float x, float sigma)
{
    return -x/(sigma*sigma) * gauss(x, sigma);
}



/*      HYSTERESIS thersholding of edge pixels. Starting at pixels with a
  value greater than the HIGH threshold, trace a connected sequence
  of pixels that have a value greater than the LOW threhsold.        */

template<typename T>
void canny<T>::hysteresis ( const image< float > &mag, image< int > &im,
                            int fringe )
{

    //cout << "Beginning hysteresis thresholding..." << std::endl;
    for ( int i=0; i<im.getHeight(); i++)
        for ( int j=0; j<im.getWidth(); j++)
          im.pixel(j,i)=0;


    if ((autothreshold) || (highthres<lowthres))
    {
      estimate_thresh( mag, fringe );
        //cout << "Estimated hysteresis thresholds (from image): highthres=" << highthres << " lowthres=" << lowthres << std::endl;
    }


    /* For each edge with a magnitude above the high threshold, begin
      tracing edge pixels that are above the low threshold.                */

    for ( int i=0; i<im.getHeight(); i++)
        for ( int j=0; j<im.getWidth(); j++)
            if (mag.pixel(j,i) * MAG_SCALE >= highthres)
              trace ( mag, im, i, j);
    return; /*************************************************/
}


template<typename T>
float canny<T>::distancePointToPoint(int ax, int ay, int bx, int by)
{
    return(sqrt( ((bx-ax)*(bx-ax)) + ((by-ay)*(by-ay)) ));
}

/* Formula from comp.graphics.algorithms FAQ */
template<typename T>
float canny<T>::distancePointToLine(int ax, int ay, int bx, int by, int px, int py)
{
    return(fabs( ((ay-py)*(bx-ax) - (ax-px)*(by-ay)) / distancePointToPoint (ax,ay,bx,by)  ));
}


/*******************************/

/*      TRACE - recursively trace edge pixels that have a threshold > the low
  edge threshold, continuing from the pixel at (i,j).                     */



template<typename T>
int canny<T>::trace ( const image< float > &mag, image< int > &im, int i, int j)
{
    char flag = 0;

    if (im.pixel(j,i) == 0)
    {
      im.pixel(j,i) = 255;
        flag=0;
        for ( int n= -1; n<=1; n++)
        {
            for ( int m= -1; m<=1; m++)
            {
                if (n==0 && m==0) continue;
                if (((i+n)>=0) && ((i+n)<mag.getHeight()) &&
                    ((j+m)>=0) && ((j+m)<mag.getWidth())) {
                  if (mag.pixel(j+m, i+n) * MAG_SCALE >=lowthres)
                    {
                      if (trace( mag, im, i+n, j+m ))
                        {
                          flag=1;
                          break;
                        }
                    }
                };
            }
            if (flag) break;
        }
        return(1);
    }
    return(0);
}






template<typename T>
void canny<T>::separable_convolution( const image< int > &im,
                                      const std::vector< float > &gau,
                                      image< float > &smx,
                                      image< float > &smy )
{
    int
      nr = im.getHeight(),
      nc = im.getWidth();
    smx.init( nc, nr );
    smy.init( nc, nr );

    for ( int i=0; i<nr; i++)
        for ( int j=0; j<nc; j++)
        {
          float x = gau[0] * im.pixel(j,i);
          float y = gau[0] * im.pixel(j,i);
          for ( int k=1; k<(signed)gau.size(); k++)
            {
              int
                I1 = (i+k)%nr,
                I2 = (i-k+nr)%nr;
              y += gau[k]*im.pixel(j,I1) + gau[k]*im.pixel(j,I2);
              int
                J1 = (j+k)%nc,
                J2 = (j-k+nc)%nc;
              x += gau[k]*im.pixel(J1,i) + gau[k]*im.pixel(J2,i);
            }
            smx.pixel(j,i) = x; smy.pixel(j,i) = y;
        }
}

template<typename T>
  void canny<T>::dxy_separable_convolution( const std::vector< float > &dgau,
                                            const image< float > &smx,
                                            const image< float > &smy,
                                            image< float > &dx,
                                            image< float > &dy )
  {
    int
      nr = smx.getHeight(),
      nc = smx.getWidth();
    
    dx.init( nc, nr );
    dy.init( nc, nr );

    for ( int i=0; i<nr; i++)
        for ( int j=0; j<nc; j++)
        {
          float
            x = 0.0,
            y = 0.0;
          for ( int k=1; k<(signed)dgau.size(); k++)
            {
              {
                int
                  I1 = (i+k)%nr,
                  I2 = (i-k+nr)%nr;
                y += -dgau[k]*smy.pixel(j,I1) + dgau[k]*smy.pixel(j,I2);
              };
              {
                int
                  I1 = (j+k)%nc,
                  I2 = (j-k+nc)%nc;
                x += -dgau[k]*smx.pixel(I1,i) + dgau[k]*smx.pixel(I2,i);
              }
            }
          dx.pixel(j,i) = x;
          dy.pixel(j,i) = y;
        }
}

template<typename T>
void canny<T>::nonmax_suppress( const image< float > &dx,
                                const image< float > &dy,
                                image< float > &mag )
{

  // image< int > ori; ori.init( mag.getWidth(), mag.getHeight() );

    for ( int i=1; i<mag.getHeight()-1; i++)
    {
        for ( int j=1; j<mag.getWidth()-1; j++)
        {
          mag.pixel(j,i)=0;

            /* Treat the x and y derivatives as components of a vector */
          float
            xc = dx.pixel(j,i),
            yc = dy.pixel(j,i);
            if ( ( fabs( xc ) + fabs( yc ) ) * MAG_SCALE < lowthres )
              continue;

            float g  = norm (xc, yc);

            float xx, yy, g2, g1, g3, g4;

            /* Follow the gradient direction, as indicated by the direction of
               the vector (xc, yc); retain pixels that are a local maximum. */

            if (fabs(yc) > fabs(xc))
            {

                /* The Y component is biggest, so gradient direction is basically UP/DOWN */
                xx = fabs(xc)/fabs(yc);
                yy = 1.0;

                g2 = norm (dx.pixel(j,i-1), dy.pixel(j,i-1));
                g4 = norm (dx.pixel(j,i+1), dy.pixel(j,i+1));
                if (xc*yc > 0.0)
                {
                    g3 = norm (dx.pixel(j+1,i+1), dy.pixel(j+1,i+1));
                    g1 = norm (dx.pixel(j-1,i-1), dy.pixel(j-1,i-1));
                }
                else
                {
                    g3 = norm (dx.pixel(j-1,i+1), dy.pixel(j-1,i+1));
                    g1 = norm (dx.pixel(j+1,i-1), dy.pixel(j+1,i-1));
                }

            } else
            {

                /* The X component is biggest, so gradient direction is basically LEFT/RIGHT */
                xx = fabs(yc)/fabs(xc);
                yy = 1.0;

                g2 = norm (dx.pixel(j+1,i), dy.pixel(j+1,i));
                g4 = norm (dx.pixel(j-1,i), dy.pixel(j-1,i));
                if (xc*yc > 0.0)
                {
                    g3 = norm (dx.pixel(j-1,i-1), dy.pixel(j-1,i-1));
                    g1 = norm (dx.pixel(j+1,i+1), dy.pixel(j+1,i+1));
                }
                else
                {
                    g1 = norm (dx.pixel(j+1,i-1), dy.pixel(j+1,i-1));
                    g3 = norm (dx.pixel(j-1,i+1), dy.pixel(j-1,i+1));
                }
            }

            /* Compute the interpolated value of the gradient magnitude */
            if ( (g > (xx*g1 + (yy-xx)*g2)) &&
                    (g > (xx*g3 + (yy-xx)*g4)) )
            {
              if (g*MAG_SCALE <= 255)
                mag.pixel(j,i) = g;
              else
                mag.pixel(j,i) = 255 / MAG_SCALE;
              // ori.setPixel(j,i,(T)(atan2 (yc, xc) * ORI_SCALE));
            } else
            {
              mag.pixel(j,i)= 0;
              // ori.setPixel(j,i,0);
            }

        }
        //  cerr << i << " " << j << std::endl;

    }
}

template<typename T>
void canny<T>::estimate_thresh( const image< float > &mag, int fringe )
{
    int i,j,k, hist[256], count;

    /* Build a histogram of the magnitude image. */
    for (k=0; k<256; k++) hist[k] = 0;

    for (i=fringe; i<mag.getHeight()-fringe; i++)
        for (j=fringe; j<mag.getWidth()-fringe; j++)
          hist[(int)(mag.pixel(j,i)*MAG_SCALE)]++;

    /* The high threshold should be > 80 or 90% of the pixels
      j = (int)(ratio*mag->info->nr*mag->info->nc);
    */
    j = mag.getHeight();
    if (j<mag.getWidth()) j = mag.getWidth();
    j = (int)(0.9*j);
    k = 255;

    count = hist[255];
    while (count < j)
    {
        k--;
        if (k<0) break;
        count += hist[k];
    }
    highthres = k;

    i=0;
    while (hist[i]==0) i++;

    lowthres = (int)((highthres+i)/2.0);
}

}

#endif /* CANNY_H */

---