filter_nonmaxsup.h

Go to the documentation of this file.

//
// more filters. this time non maxima suppression
// (static)
// stumeikle Sun May  5 20:48:34 2002
//
//

#ifndef FILTER_NONMAXSUP_H
#define FILTER_NONMAXSUP_H

#include <boost/shared_array.hpp>
#include "mimasexception.h"
#include "filter.h"
#include "image.h"
#include "property_image.h"
#include "edgel.h"
#include <cmath>

namespace mimas  {
  /*
    General image filter routines... continued. NonMax suppression filter is
    placed in this file. Chose to break up the original filter.h file to
    keep the file sizes down.
 $Header: /cvs/mimas2/include/filter_nonmaxsup.h,v 1.1.1.1 2005/08/09 15:37:45 engmb Exp $
   */
template<typename T>
void filter<T>::nonMaximaSuppression( const image<T>& gradx,
                                         const image<T>& grady,
                                         const image<T>& gradsq,
                                         double threshold,
                                         property_image& output )
  throw (mimasexception)
{
    double   sqthres = (double)(threshold * threshold);
    output.init( gradsq.getWidth(), gradsq.getHeight());

    MMERROR( gradx.getWidth() == grady.getWidth() &&
             gradx.getWidth() == gradsq.getWidth() &&
             gradx.getHeight() == grady.getHeight() &&
             gradx.getHeight() == gradsq.getHeight(), mimasexception, ,
             "Trying to do non maxima suppression with images of different "
             "sizes" );
    
    int      w = gradx.getWidth();
    int      h = gradx.getHeight();
    boost::shared_array< double >
      rowx( new double[w] ),
      rowy( new double[w] ),
      grad( new double[w] ),
      grad_m1( new double[w] ),
      grad_p1( new double[w] );

    for(int y=1;y<(h-1);++y)
    {
  int x;

  //extract lots of rows and convert to double (@ WHY TINA?)
  for(x=0;x<w;++x)
  {
      rowx[x] = (double)gradx.pixel(x,y);
      rowy[x] = (double)grady.pixel(x,y);
      grad[x] = (double)gradsq.pixel(x,y);
      grad_p1[x]=(double)gradsq.pixel(x,y+1);
      grad_m1[x]=(double)gradsq.pixel(x,y-1);
  }
    
  for(x=1;x<(w-1);++x)
  {
      // now the messy guts ...
      double   a, b;
      double   st;
      double   del;
      double   dr, dc;
      double   gradratio, fabsgradratio;

      if (grad[x] < sqthres)
    continue;

      st = grad[x];
      gradratio = rowx[x] / rowy[x];
      fabsgradratio = (float)fabs(gradratio);
      if (fabsgradratio < 1.5 && fabsgradratio > 0.67)  // diagonal
      {
    if ((st < grad[x - 1] || st <= grad[x + 1]) &&
        (st < grad_m1[x]  || st <= grad_p1[x]))
        continue;
        
    if (gradratio > 0)  // right diagonal
    {
        if (st < grad_m1[x - 1] || st <= grad_p1[x + 1])
      continue;
        st = sqrt(st);
        a  = sqrt(grad_m1[x - 1]);
        b  = sqrt(grad_p1[x + 1]);
        del = (a - b) / ((a + b - st * 2) * 2);
        dr = (del + 0.5);
        dc = (del + 0.5);
    } else
    {
        if (st < grad_p1[x - 1] || st <= grad_m1[x + 1])
      continue;
        st = sqrt(st);
        a =  sqrt(grad_p1[x - 1]);
        b =  sqrt(grad_m1[x + 1]);
        del = (a - b) / ((a + b - st * 2) * 2);
        dr = (0.5 - del);
        dc = (del + 0.5);
    }
      } else if (fabs(rowx[x]) > fabs(rowy[x]))
      {
    if (st < grad[x - 1] || st <= grad[x + 1])
        continue;
    st = sqrt(st);
    a = sqrt(grad[x - 1]);
    b = sqrt(grad[x + 1]);
    del = (a - b) / ((a + b - st * 2) * 2);
    dr = 0.5;
    dc = (del + 0.5);
      } else
      {
    if (st < grad_m1[x] || st <= grad_p1[x])
        continue;
    st = sqrt(st);
    a = sqrt(grad_m1[x]);
    b = sqrt(grad_p1[x]);
    del = (a - b) / ((a + b - st * 2) * 2);
    dr = (del + 0.5);
    dc = 0.5; // pixel centre 
      }

      {
    int     ii, jj; // actual storage location 

    edgel_ptr ej( new edgel );

    //positions are measured relative to the image plane
    //(IMPORTANT)
    double  posx = (double)x + dc;
    ii = (int)floor(posx);
    double  posy = (double)y + dr;
    jj = (int)floor(posy);
    
    if (output.getPixel(ii,jj)) {
        // delete ej; //dont store >1 edgel per pixel
    } else
    {
        ej->setPosition( posx,posy );
        ej->setOrientation( angle( rowx[x], rowy[x]) );
        ej->setContrast( st );
        output.setPixel( ii,jj, ej );
    }
      }
  }
    }
}


}

#endif

---