corner_tool.h

Go to the documentation of this file.

//
// harris and stevens corner detection
// stumeikle Tue May  7 21:16:45 2002
//

#ifndef HSCORNER_FILTER_H
#define HSCORNER_FILTER_H

#include "filter.h"
#include "image.h"
#include "image_op.h"
#include "gauss.h"
#include "filter_grad.h"
#include "property_image.h"
#include "corner.h"

namespace mimas {
  
template<typename T>
class corner_tool
{
    private:
  //actually most of these aren't used!
  double   sigma, maxError;   
  double   thres;       
  double   lowd, upd, width;    
  double   disp ;       
  double   uniqueness, correlate, asymetry ;
  double   tuniqueness, tcorrelate, tasymetry;
  double   theight, twidth ;    
  double   patch_sigma;     
  int      stereo_images;     
  int      fixed_matches;     
  int      good_matches;      

    public:
  corner_tool()
  {
      sigma = 1.0, maxError = 0.5 / 256.0, thres = 6.0;
      lowd = -0.5, upd = 0.5, width = 3.0;
      disp = 0;
      uniqueness = 0.004, correlate = 0.985, asymetry= 0.85;
      tuniqueness = 0.001, tcorrelate = 0.980, tasymetry = 0.85;
      theight = 20.0, twidth = 60.0;
      patch_sigma = 3.0;
      stereo_images = 1;
      fixed_matches = 0;
      good_matches = 0;
  }

  corner_tool(double t)
  {
      sigma = 1.0, maxError = 0.5 / 256.0, thres = t;
      lowd = -0.5, upd = 0.5, width = 3.0;
      disp = 0;
      uniqueness = 0.004, correlate = 0.985, asymetry= 0.85;
      tuniqueness = 0.001, tcorrelate = 0.980, tasymetry = 0.85;
      theight = 20.0, twidth = 60.0;
      patch_sigma = 3.0;
      stereo_images = 1;
      fixed_matches = 0;
      good_matches = 0;
  }

  void  setThreshold( double t ) 
  {
      thres = t;
  }

  void filterImage( image<T> &input, corner_image &corner_im )
  {
      image<T>    gradx,grady,outsq;
      image<T>    gim;
      image<T>    peak_im;
      //corner_image  corner_im;

      doCornerExtract( input, sigma, maxError, 0.05, peak_im, gradx, grady );
      gim = gaussBlur< float >( peak_im, sigma, maxError );
      peak_im.clear();
      findMaxima( gim, 1024.0*thres,0,corner_im);
      cornerLocate( gim, gradx, grady, corner_im );
  }


private:
  void  doCornerExtract( image<T>& im, double sigma, double maxError,
         double edge_sup, image<T>& output_im, image<T>& gradx,
         image<T>& grady )
  {
     if ( &output_im != &im && &output_im != &gradx && &output_im != &grady )
       output_im.init( im.getWidth(), im.getHeight() );

      filter<T>::gradX( im, gradx );
      filter<T>::gradY( im, grady );

      image<T>
         gradsqx ( gradx * gradx ),
         gradsqy ( grady * grady ),
         gradsqxy( gradx * grady ),
         gradsqxsqy( gradsqx + gradsqy );

      image<T>  gimx, gimy, gimxy, gimedge;

      gimx = gaussBlur< float >( gradsqx, sigma, maxError );
      gradsqx.clear();
      gimy = gaussBlur< float >( gradsqy, sigma, maxError );
      gradsqy.clear();
      gimxy = gaussBlur< float >( gradsqxy, sigma, maxError );
      gradsqxy.clear();

      gimedge = gaussBlur< float >( gradsqxsqy, 1.1 * sigma, maxError );
      gradsqxsqy.clear();

      //now do corner im 2!!!!!! argggg
      doCornerExtract2(gimx,gimy,gimxy,gimedge,edge_sup, output_im);

  }


  void  doCornerExtract2( const image<T>& im1, const image<T>& im2,
                           const image<T>& im3, const image<T>& im4,
                           double edge_supp, image<T>&  output ) const
  {
     MMERROR( im1.getSize() == im2.getSize(),
              mimasexception, , "corner_tool::doCornerExtract2 - "
              "Images 1 and 2 are of different size." );
     MMERROR( im1.getSize() == im3.getSize(),
              mimasexception, , "corner_tool::doCornerExtract2 - "
              "Images 1 and 3 are of different size." );
     MMERROR( im1.getSize() == im4.getSize(),
              mimasexception, , "corner_tool::doCornerExtract2 - "
              "Images 1 and 4 are of different size." );
     MMERROR( im1.getSize() == output.getSize(),
              mimasexception, , "corner_tool::doCornerExtract2 - "
              "Image 1 and output image are of different size." );
     const T
       *p1 = im1.rawData(),
       *p2 = im2.rawData(),
       *p3 = im3.rawData(),
       *p4 = im4.rawData();
     T *q = output.rawData();

     for ( int i=0; i<(signed)output.getSize(); i++ ) {
       *q = (T)( (double)*p1 * (double)*p2 -
                 (double)*p3 * (double)*p3 -
                 edge_supp * (double)*p4 * (double)*p4 );
       p1++; p2++; p3++; p4++; q++;
     };
  }

  void  findMaxima(const image<T> &input, double thres, int num,
      corner_image& output) 
  {
       output.init( input.getWidth(), input.getHeight() );

      int wm1 = input.getWidth()-1;
      int hm1 = input.getHeight()-1;
      int nmax;
      double pix;

      for(int y=1;y<hm1;++y)
      {
         const T
           *max = input.rawData() + y * input.getWidth(),
           *max_m1 = max - input.getWidth(),
           *max_p1 = max + input.getWidth();
    for(int x=1;x<wm1;++x)
    {

        pix = (double) max[x];
        if (pix<thres)
          continue;

        nmax = 0;

        if (pix<max_m1[x-1] && ++nmax > num) continue;
        if (pix<max_m1[x]   && ++nmax > num) continue;
        if (pix<max_m1[x+1] && ++nmax > num) continue;
        if (pix<max[x-1]    && ++nmax > num) continue;
        if (pix<max[x+1]    && ++nmax > num) continue;
        if (pix<max_p1[x-1] && ++nmax > num) continue;
        if (pix<max_p1[x]   && ++nmax > num) continue;
        if (pix<max_p1[x+1] && ++nmax > num) continue;

        corner_ptr c( new corner );

        c->setPosition((double)x+0.5,(double)y+0.5);
        c->setOrientation(angle(0));
        c->setContrast( sqrt(pix) );

        if (!output.pixel(x,y))
          output.pixel( x,y ) = c;
    }
      }
  }


  void  cornerLocate( image<T>& corner_im,
            image<T>& gradx, image<T>& grady,
            corner_image& im ) 
  {
      int x,y;
      double  rx,ry;

      for(y=0;y<im.getHeight();++y)
      {
    for(x=0;x<im.getWidth();x++)
    {
        // we ASSUME here that the propery image contains only corners
          corner_ptr c( boost::dynamic_pointer_cast< corner >
                           ( im.pixel(x,y) ) );

        if (!c) continue;

        //process the corner and provide estimates of its orientation and
        //actual sub pixel position
        c->setContrast( sqrt(quadFit( corner_im, x,y, &rx, &ry )) );
        c->setPosition( rx,ry );
        c->setOrientation( cornerOrient(rx,ry, 10.0, gradx, grady ) );
    }
      }
  }

    double  quadFit( image<T>& image, int xin, int yin, double * px, double *py )
    {
  double  x,y;
  x = (double)xin;
  y = (double)yin;

  int n,m;
  double  pixval[3][3];
  for(n=0;n<3;++n)
  {
      for(m=0;m<3;++m)
      {
          pixval[n][m] = image.getPixel( xin+m-1, yin+n-1 );
      }
  }

  double           a, b, c, d, e, f;
  double          inter;
        double           temp;
  double    xs,ys;

        a = pixval[1][1];
        b = (pixval[0][2] - pixval[0][0]
    + pixval[1][2] - pixval[1][0]
    + pixval[2][2] - pixval[2][0]) / 6.0;
        c = (pixval[2][0] - pixval[0][0]
    + pixval[2][1] - pixval[0][1]
    + pixval[2][2] - pixval[0][2]) / 6.0;
        d = (pixval[0][0] - 2.0 * pixval[0][1] + pixval[0][2]
    + 3.0 * pixval[1][0] - 6.0 * pixval[1][1] + 3.0 * pixval[1][2]
    + pixval[2][0] - 2.0 * pixval[2][1] + pixval[2][2]) / 10.0;
        e = (pixval[0][0] - pixval[2][0]
    + pixval[2][2] - pixval[0][2]) / 4.0;
        f = (pixval[0][0] + 3.0 * pixval[0][1] + pixval[0][2]
    - 2.0 * pixval[1][0] - 6.0 * pixval[1][1] - 2.0 * pixval[1][2]
    + pixval[2][0] + 3.0 * pixval[2][1] + pixval[2][2]) / 10.0;

        temp = 4.0 * d * f - e * e;
        *px = 1.5 + x + (e * c - 2 * f * b) / temp;
        *py = 1.5 + y + (e * b - 2 * d * c) / temp;

        xs = (e * c - 2 * f * b) / temp;
        ys = (e * b - 2 * d * c) / temp;

        inter = a + b * xs + c * ys + d * xs * xs + e * xs * ys + f * ys * ys;

        return (inter); 
    }


    angle 
    cornerOrient(double rx,double ry, double range,
     image<T>& gradx, image<T>& grady )
    {

        int lx,ux,ly,uy,x,y;
        double meanx=0.0,meany=0.0;
        double dx,dy, delx,dely;
        double twosigma2 = range*range/2.0; // set the scale of radial weighting
  double weight;
        double thresh=0.0;
        double pix,v;
  angle orient;

        lx = (int)(rx-range-1.0);
  ux = (int)(rx+range+2);
        ly = (int)(ry-range-1.0);
        uy = (int)(ry+range+2);

        image<double> grad;
        grad.init( ux-lx, uy-ly );

  // calculate a robust mean with which to threshold the likely edge pixels 
  for (y=ly;y<uy;++y)
  {
          for (x=lx;x<ux;++x)
          {
          dx = (double)gradx.getPixel( x,y );
          dy = (double)grady.getPixel( x,y );
    v  = dx*dx+dy*dy;
    grad.setPixel( x-lx, y-ly, v );
          thresh += v;
          }
  }

  thresh /= (double)((ux-lx)*(uy-ly));

  // calculate the orientation from the moments of the thresholded gradient image
  for (y=ly+1;y<uy-1;++y)
  {
          for (x=lx+1;x<ux-1;++x)
          {
          if ((pix=grad.getPixel(x-lx,y-ly)) > thresh)
          {
                  delx = (double)x + 0.5 - rx;
                  dely = (double)y + 0.5 - ry;

                  if (delx*delx + dely*dely < twosigma2*2.0) weight = 1.0;
                  else weight = 0.0;
                  meanx += delx*weight;
                  meany += dely*weight;
          }
          }
  }

        if (meanx ==0) orient = angle(0); /* ### should be undefined */
  else
          orient = angle(meany,meanx);

        return(orient);
    }

};

};

#endif

---