hysteresis.c

#include <stdio.h>
#include <stdlib.h>
#include "delogo.h"

/*
 * Recursive routine that traces edgs along all paths whose magnitude values 
 * remain above some specifyable lower threshhold. 
 */
void
follow_edges( u8 * edgemapptr, short *edgemagptr, short lowval, int cols )
{
  short        *tempmagptr;
  u8           *tempmapptr;
  int           i;
  int           x[8] = { 1, 1, 0, -1, -1, -1, 0, 1 }, y[8] =
  {
  0, 1, 1, 1, 0, -1, -1, -1};

  for( i = 0; i < 8; i++ )
  {
    tempmapptr = edgemapptr - y[i] * cols + x[i];
    tempmagptr = edgemagptr - y[i] * cols + x[i];

    if( ( *tempmapptr == POSSIBLE_EDGE ) && ( *tempmagptr > lowval ) )
    {
      *tempmapptr = ( u8 ) EDGE;
      follow_edges( tempmapptr, tempmagptr, lowval, cols );
    }
  }
}

/*
 * Finds edges that are above some high threshhold or are connected to a
 * high pixel by a path of pixels greater than a low threshold.  
 */
void
apply_hysteresis( short int *mag, u8 * nms, int rows, int cols,
    float tlow, float thigh, u8 * edge )
{
  int           r, c, pos, numedges, highcount, lowthreshold,
      highthreshold, hist[32768];
  short int     maximum_mag;

  maximum_mag = 0;
  /* 
   * Initialize the edge map to possible edges everywhere the non-maximal
   * suppression suggested there could be an edge except for the border. At
   * the border we say there can not be an edge because it makes the
   * follow_edges algorithm more efficient to not worry about tracking an
   * edge off the side of the image. 
   */
  for( r = 0, pos = 0; r < rows; r++ )
  {
    for( c = 0; c < cols; c++, pos++ )
    {
      if( nms[pos] == POSSIBLE_EDGE )
	edge[pos] = POSSIBLE_EDGE;
      else
	edge[pos] = NOEDGE;
    }
  }

  for( r = 0, pos = 0; r < rows; r++, pos += cols )
  {
    edge[pos] = NOEDGE;
    edge[pos + cols - 1] = NOEDGE;
  }
  pos = ( rows - 1 ) * cols;
  for( c = 0; c < cols; c++, pos++ )
  {
    edge[c] = NOEDGE;
    edge[pos] = NOEDGE;
  }

  /* 
   * Compute the histogram of the magnitude image. Then use the histogram to
   * compute hysteresis thresholds.
   */
  for( r = 0; r < 32768; r++ )
    hist[r] = 0;
  for( r = 0, pos = 0; r < rows; r++ )
  {
    for( c = 0; c < cols; c++, pos++ )
    {
      if( edge[pos] == POSSIBLE_EDGE )
	hist[mag[pos]]++;
    }
  }

  // Compute the number of pixels that passed the nonmaximal suppression.
  for( r = 1, numedges = 0; r < 32768; r++ )
  {
    if( hist[r] != 0 )
      maximum_mag = r;
    numedges += hist[r];
  }

  highcount = ( int ) ( numedges * thigh + 0.5 );

  /* 
   * Compute the high threshold value as the (100 * thigh) percentage point
   * in the magnitude of the gradient histogram of all the pixels that passes
   * non-maximal suppression. Then calculate the low threshold as a fraction
   * of the computed high threshold value. John Canny said in his paper
   * "A Computational Approach to Edge Detection" that "The ratio of the
   * high to low threshold in the implementation is in the range two or three
   * to one." That means that in terms of this implementation, we should
   * choose tlow ~= 0.5 or 0.33333.
   */
  r = 1;
  numedges = hist[1];
  while( ( r < ( maximum_mag - 1 ) ) && ( numedges < highcount ) )
  {
    r++;
    numedges += hist[r];
  }
  highthreshold = r;
  lowthreshold = ( int ) ( highthreshold * tlow + 0.5 );

  /* 
   * This loop looks for pixels above the highthreshold to locate edges and
   * then calls follow_edges to continue the edge.
   */
  for( r = 0, pos = 0; r < rows; r++ )
  {
    for( c = 0; c < cols; c++, pos++ )
    {
      if( ( edge[pos] == POSSIBLE_EDGE ) && ( mag[pos] >= highthreshold ) )
      {
	edge[pos] = EDGE;
	follow_edges( ( edge + pos ), ( mag + pos ), lowthreshold, cols );
      }
    }
  }

  // Set all the remaining possible edges to non-edges.
  for( r = 0, pos = 0; r < rows; r++ )
  {
    for( c = 0; c < cols; c++, pos++ )
      if( edge[pos] != EDGE )
	edge[pos] = NOEDGE;
  }
}

// Applies non-maximal suppression to the magnitude of the gradient image.
void
non_max_supp( short *mag, short *gradx, short *grady, int nrows, int ncols,
    u8 * result )
{
  int           rowcount, colcount, count;
  short        *magrowptr, *magptr;
  short        *gxrowptr, *gxptr;
  short        *gyrowptr, *gyptr, z1, z2;
  short         m00, gx, gy;
  float         mag1, mag2, xperp, yperp;
  u8           *resultrowptr, *resultptr;

  gx = gy = 0;
  xperp = yperp = 0;

  // Zero the edges of the result image.
  for( count = 0, resultrowptr = result, resultptr =
      result + ncols * ( nrows - 1 ); count < ncols;
      resultptr++, resultrowptr++, count++ )
  {
    *resultrowptr = *resultptr = ( u8 ) 0;
  }

  for( count = 0, resultptr = result, resultrowptr = result + ncols - 1;
      count < nrows; count++, resultptr += ncols, resultrowptr += ncols )
  {
    *resultptr = *resultrowptr = ( u8 ) 0;
  }

  // Suppress non-maximum points.
  for( rowcount = 1, magrowptr = mag + ncols + 1, gxrowptr =
      gradx + ncols + 1, gyrowptr = grady + ncols + 1, resultrowptr =
      result + ncols + 1; rowcount < nrows - 2;
      rowcount++, magrowptr += ncols, gyrowptr += ncols, gxrowptr +=
      ncols, resultrowptr += ncols )
  {
    for( colcount = 1, magptr = magrowptr, gxptr = gxrowptr, gyptr =
	gyrowptr, resultptr = resultrowptr; colcount < ncols - 2;
	colcount++, magptr++, gxptr++, gyptr++, resultptr++ )
    {
      m00 = *magptr;
      if( m00 == 0 )
      {
	*resultptr = ( u8 ) NOEDGE;
      }
      else
      {
	xperp = -( gx = *gxptr ) / ( ( float ) m00 );
	yperp = ( gy = *gyptr ) / ( ( float ) m00 );
      }

      if( gx >= 0 )
      {
	if( gy >= 0 )
	{
	  if( gx >= gy )
	  {
	    // 111 
	    // Left point 
	    z1 = *( magptr - 1 );
	    z2 = *( magptr - ncols - 1 );

	    mag1 = ( m00 - z1 ) * xperp + ( z2 - z1 ) * yperp;

	    // Right point 
	    z1 = *( magptr + 1 );
	    z2 = *( magptr + ncols + 1 );

	    mag2 = ( m00 - z1 ) * xperp + ( z2 - z1 ) * yperp;
	  }
	  else
	  {
	    // 110 
	    // Left point 
	    z1 = *( magptr - ncols );
	    z2 = *( magptr - ncols - 1 );

	    mag1 = ( z1 - z2 ) * xperp + ( z1 - m00 ) * yperp;

	    // Right point 
	    z1 = *( magptr + ncols );
	    z2 = *( magptr + ncols + 1 );

	    mag2 = ( z1 - z2 ) * xperp + ( z1 - m00 ) * yperp;
	  }
	}
	else
	{
	  if( gx >= -gy )
	  {
	    // 101 
	    // Left point 
	    z1 = *( magptr - 1 );
	    z2 = *( magptr + ncols - 1 );

	    mag1 = ( m00 - z1 ) * xperp + ( z1 - z2 ) * yperp;

	    // Right point 
	    z1 = *( magptr + 1 );
	    z2 = *( magptr - ncols + 1 );

	    mag2 = ( m00 - z1 ) * xperp + ( z1 - z2 ) * yperp;
	  }
	  else
	  {
	    // 100 
	    // Left point 
	    z1 = *( magptr + ncols );
	    z2 = *( magptr + ncols - 1 );

	    mag1 = ( z1 - z2 ) * xperp + ( m00 - z1 ) * yperp;

	    // Right point 
	    z1 = *( magptr - ncols );
	    z2 = *( magptr - ncols + 1 );

	    mag2 = ( z1 - z2 ) * xperp + ( m00 - z1 ) * yperp;
	  }
	}
      }
      else
      {
	if( ( gy = *gyptr ) >= 0 )
	{
	  if( -gx >= gy )
	  {
	    // 011 
	    // Left point 
	    z1 = *( magptr + 1 );
	    z2 = *( magptr - ncols + 1 );

	    mag1 = ( z1 - m00 ) * xperp + ( z2 - z1 ) * yperp;

	    // Right point 
	    z1 = *( magptr - 1 );
	    z2 = *( magptr + ncols - 1 );

	    mag2 = ( z1 - m00 ) * xperp + ( z2 - z1 ) * yperp;
	  }
	  else
	  {
	    // 010 
	    // Left point 
	    z1 = *( magptr - ncols );
	    z2 = *( magptr - ncols + 1 );

	    mag1 = ( z2 - z1 ) * xperp + ( z1 - m00 ) * yperp;

	    // Right point 
	    z1 = *( magptr + ncols );
	    z2 = *( magptr + ncols - 1 );

	    mag2 = ( z2 - z1 ) * xperp + ( z1 - m00 ) * yperp;
	  }
	}
	else
	{
	  if( -gx > -gy )
	  {
	    // 001 
	    // Left point 
	    z1 = *( magptr + 1 );
	    z2 = *( magptr + ncols + 1 );

	    mag1 = ( z1 - m00 ) * xperp + ( z1 - z2 ) * yperp;

	    // Right point 
	    z1 = *( magptr - 1 );
	    z2 = *( magptr - ncols - 1 );

	    mag2 = ( z1 - m00 ) * xperp + ( z1 - z2 ) * yperp;
	  }
	  else
	  {
	    // 000 
	    // Left point 
	    z1 = *( magptr + ncols );
	    z2 = *( magptr + ncols + 1 );

	    mag1 = ( z2 - z1 ) * xperp + ( m00 - z1 ) * yperp;

	    // Right point 
	    z1 = *( magptr - ncols );
	    z2 = *( magptr - ncols - 1 );

	    mag2 = ( z2 - z1 ) * xperp + ( m00 - z1 ) * yperp;
	  }
	}
      }

      // Now determine if the current point is a maximum point 
      if( ( mag1 > 0.0 ) || ( mag2 > 0.0 ) )
      {
	*resultptr = ( u8 ) NOEDGE;
      }
      else
      {
	if( mag2 == 0.0 )
	  *resultptr = ( u8 ) NOEDGE;
	else
	  *resultptr = ( u8 ) POSSIBLE_EDGE;
      }
    }
  }
}