susan边缘及角点检测vc代码.txt

SUSAN边缘检测及角点检测的vc实现代码

/******************************/
/* }}} */
/* {{{ Readme First */
/******************************/  
See following section for different machine information. Please report any bugs (and fixes). There are a few optional changes that  can be made in the "defines" section which follows shortly.
Usage: type "susan" to get usage. Only PGM format files can be input and output. Utilities such as the netpbm package and XV can be used to convert to and from other formats. Any size of image can be  processed.
THRESHOLDS:  There are two thresholds which can be set at run-time. These are the brightness threshold (t) and the distance threshold (d).
SPATIAL CONTROL: d
    In SUSAN smoothing d controls the size of the Gaussian mask; its  default is 4.0. Increasing d gives more smoothing. In edge finding,  a fixed flat mask is used, either 37 pixels arranged in a "circle"  (default), or a 3 by 3 mask which gives finer detail. In corner  finding, only the larger 37 pixel mask is used; d is not  variable. In smoothing, the flat 3 by 3 mask can be used instead of  a larger Gaussian mask; this gives low smoothing and fast operation.
BRIGHTNESS CONTROL: t
  In all three algorithms, t can be varied (default=20); this is the  main threshold to be varied. It determines the maximum difference in  greylevels between two pixels which allows them to be considered  part of the same "region" in the image. Thus it can be reduced to  give more edges or corners, i.e. to be more sensitive, and vice  versa. In smoothing, reducing t gives less smoothing, and vice  versa. Set t=10 for the test image available from the SUSAN web  page.
ITERATIONS:
  With SUSAN smoothing, more smoothing can also be obtained by iterating the algorithm several times. This has a different effect  from varying d or t.
FIXED MASKS: 
37 pixel mask: 
              




  
CORNER ATTRIBUTES dx, dy and I  (Only read this if you are interested in the C implementation or in  using corner attributes, e.g., for corner matching)  Corners reported in the corner list have attributes associated with  them as well as positions. This is useful, for example, when attempting to match corners from one image to another, as these  attributes can often be fairly unchanged between images. The  attributes are dx, dy and I. I is the value of image brightness at  the position of the corner. In the case of susan_corners_quick, dx  and dy are the first order derivatives (differentials) of the image  brightness in the x and y directions respectively, at the position  of the corner. In the case of normal susan corner finding, dx and dy  are scaled versions of the position of the centre of gravity of the  USAN with respect to the centre pixel (nucleus).
BRIGHTNESS FUNCTION LUT IMPLEMENTATION:  (Only read this if you are interested in the C implementation)
The SUSAN brightness function is implemented as a LUT  (Look-Up-Table) for speed. The resulting pointer-based code is a  little hard to follow, so here is a brief explanation. In  setup_brightness_lut() the LUT is setup. This mallocs enough space  for *bp and then repositions the pointer to the centre of the  malloced space. The SUSAN function e^-(x^6) or e^-(x^2) is  calculated and converted to a uchar in the range 0-100, for all  possible image brightness differences (including negative  ones). Thus bp[23] is the output for a brightness difference of 23  greylevels. In the SUSAN algorithms this LUT is used as follows:
  p=in + (i-3)*x_size + j - 1;
  p points to the first image pixel in the circular mask surrounding  point (x,y).
  cp=bp + in[i*x_size+j];
  cp points to a position in the LUT corresponding to the brightness  of the centre pixel (x,y).  now for every pixel within the mask surrounding (x,y),
  n+=*(cp-*p++);
  the brightness difference function is found by moving the cp pointer  down by an amount equal to the value of the pixel pointed to by p,  thus subtracting the two brightness values and performing the  exponential function. This value is added to n, the running USAN  area. in SUSAN smoothing, the variable height mask is implemented by  multiplying the above by the moving mask pointer, reset for each new  centre pixel.
  tmp = *dpt++ * *(cp-brightness);
/******************************/
/* }}} */
/* {{{ defines, includes and typedefs */
/* ********** Optional settings */
#ifdef PPC
typedef int        TOTAL_TYPE; /* this is faster for "int" but should be "float" for large d masks */
#else
typedef float      TOTAL_TYPE; /* for my PowerPC accelerator only */
#endif
/*#define FOPENB*/           /* uncomment if using djgpp gnu C for DOS or certain Win95 compilers */
#define SEVEN_SUPP           /* size for non-max corner suppression; SEVEN_SUPP or FIVE_SUPP */
#define MAX_CORNERS   15000  /* max corners per frame */
/* Leave the rest - but you may need to remove one or both of sys/file.h and malloc.h lines */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <iostream.h>
#include <malloc.h>      /* may want to remove this line */
#define  exit_error(IFB,IFC) { fprintf(stderr,IFB,IFC); exit(0); }
#define  FTOI(a) ( (a) < 0 ? ((int)(a-0.5)) : ((int)(a+0.5)) )
typedef  unsigned char uchar;
typedef  struct {int x,y,info, dx, dy, I;} CORNER_LIST[MAX_CORNERS];

/* {{{ get_image(filename,in,x_size,y_size) */
/* {{{ int getint(fp) derived from XV */
int getint(FILE *fd)
  {
  int c, i;
  char dummy[10000];

  c = getc(fd);
  while (1) /* find next integer */
  {
    if (c=='#')    /* if we're at a comment, read to end of line */
      fgets(dummy,9000,fd);
    if (c==EOF)
      exit_error("Image %s not binary PGM.\n","is");
    if (c>='0' && c<='9')
      break;   /* found what we were looking for */
    c = getc(fd);
  }
/* we're at the start of a number, continue until we hit a non-number */
  i = 0;
  while (1) {
    i = (i*10) + (c - '0');
    c = getc(fd);
    if (c==EOF) return (i);
    if (c<'0' || c>'9') break;
  }
  return (i);
}
/* }}} */
void get_image(char filename[200],unsigned char  **in,int *x_size,int *y_size)
 {
FILE  *fd;
char header [100];
int  tmp;

#ifdef FOPENB
  if ((fd=fopen(filename,"rb")) == NULL)
#else
  if ((fd=fopen(filename,"r")) == NULL)
#endif
    exit_error("Can't input image %s.\n",filename);
  /* {{{ read header */
  header[0]=fgetc(fd);
  header[1]=fgetc(fd);
  if(!(header[0]=='P' && header[1]=='5'))
    exit_error("Image %s does not have binary PGM header.\n",filename);
  *x_size = getint(fd);
  *y_size = getint(fd);
  tmp = getint(fd);
/* }}} */
  *in = (uchar *) malloc((*x_size) * (*y_size));
  if (fread(*in,1,(*x_size) * (*y_size),fd) == 0)
    exit_error("Image %s is wrong size.\n",filename);
  fclose(fd);
}
/* }}} */
/* {{{ put_image(filename,in,x_size,y_size) */
void put_image(char filename[100],char *in,int x_size,int y_size)
  {
FILE  *fd;
#ifdef FOPENB
  if ((fd=fopen(filename,"wb")) == NULL) 
#else
  if ((fd=fopen(filename,"w")) == NULL) 
#endif
    exit_error("Can't output image%s.\n",filename);
  fprintf(fd,"P5\n");
  fprintf(fd,"%d %d\n",x_size,y_size);
  fprintf(fd,"255\n");
    if (fwrite(in,x_size*y_size,1,fd) != 1)
    exit_error("Can't write image %s.\n",filename);
  fclose(fd);
}
/* }}} */
/* {{{ int_to_uchar(r,in,size) */
void int_to_uchar(int *r,uchar *in,int size)
 {
int i,
    max_r=r[0],
    min_r=r[0];
  for (i=0; i<size; i++)
    {
      if ( r[i] > max_r )
        max_r=r[i];
      if ( r[i] < min_r )
        min_r=r[i];
    }
   max_r-=min_r;
  for (i=0; i<size; i++)
    in[i] = (uchar)((int)((int)(r[i]-min_r)*255)/max_r);
}
void setup_brightness_lut(uchar **bp,int thresh,int form)
 {
int   k;
float temp;
  *bp=(unsigned char *)malloc(516);
  *bp=*bp+258;
  for(k=-256;k<257;k++)
  {
    temp=((float)k)/((float)thresh);
    temp=temp*temp;
    if (form==6)
      temp=temp*temp*temp;
    temp=(float)exp(-temp)*(float)100.0;
    *(*bp+k)= (uchar)temp;
  }
}
/* }}} */
/* {{{ susan principle */

/* {{{ susan_principle(in,r,bp,max_no,x_size,y_size) */
void susan_principle(uchar *in,int *r,uchar *bp,int max_no,int x_size,int y_size)
  {
int   i, j, n;
uchar *p,*cp;
  memset (r,0,x_size * y_size * sizeof(int));
  for (i=3;i<y_size-3;i++)
    for (j=3;j<x_size-3;j++)
    {
      n=100;
      p=in + (i-3)*x_size + j - 1;
      cp=bp + in[i*x_size+j];

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-3; 

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-5;

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-6;

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=2;
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-6;

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-5;

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-3;

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);

      if (n<=max_no)
        r[i*x_size+j] = max_no - n;
    }
}
/* }}} */
/* {{{ susan_principle_small(in,r,bp,max_no,x_size,y_size) */
void susan_principle_small(uchar *in,int *r,uchar *bp,int max_no,int x_size,int y_size)
  {
int   i, j, n;
uchar *p,*cp;

  memset (r,0,x_size * y_size * sizeof(int));
  max_no = 730; /* ho hum ;) */
  for (i=1;i<y_size-1;i++)
    for (j=1;j<x_size-1;j++)
    {
      n=100;
      p=in + (i-1)*x_size + j - 1;
      cp=bp + in[i*x_size+j];

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);
      p+=x_size-2; 

      n+=*(cp-*p);
      p+=2;
      n+=*(cp-*p);
      p+=x_size-2;

      n+=*(cp-*p++);
      n+=*(cp-*p++);
      n+=*(cp-*p);

      if (n<=max_no)
        r[i*x_size+j] = max_no - n;
    }
}
/* }}} */
/* {{{ smoothing */
/* {{{ median(in,i,j,x_size) */
uchar median(uchar *in,int i,int j,int x_size)
  {
int p[8],k,l,tmp;
  p[0]=in[(i-1)*x_size+j-1];
  p[1]=in[(i-1)*x_size+j  ];
  p[2]=in[(i-1)*x_size+j+1];
  p[3]=in[(i  )*x_size+j-1];
  p[4]=in[(i  )*x_size+j+1];
  p[5]=in[(i+1)*x_size+j-1];
  p[6]=in[(i+1)*x_size+j  ];
  p[7]=in[(i+1)*x_size+j+1];

  for(k=0; k<7; k++)
    for(l=0; l<(7-k); l++)
      if (p[l]>p[l+1])
      {
        tmp=p[l]; p[l]=p[l+1]; p[l+1]=tmp;
      }
  return( (p[3]+p[4]) / 2 );
}
/* }}} */
/* {{{ enlarge(in,tmp_image,x_size,y_size,border) */
/* this enlarges "in" so that borders can be dealt with easily */
void enlarge(uchar **in,uchar *tmp_image,int *x_size,int *y_size,int border)
  {
int   i, j;
  for(i=0; i<*y_size; i++)   /* copy *in into tmp_image */
memcpy(tmp_image+(i+border)*(*x_size+2*border)+border, *in+i* *x_size, *x_size);
  for(i=0; i<border; i++) /* copy top and bottom rows; invert as many as necessary */
  {
 memcpy(tmp_image+(border-1-i)*(*x_size+2*border)+border,*in+i* *x_size,*x_size);
    memcpy(tmp_image+(*y_size+border+i)*(*x_size+2*border)+border,*in+(*y_size-i-1)* *x_size,*x_size);
  }
  for(i=0; i<border; i++) /* copy left and right columns */
    for(j=0; j<*y_size+2*border; j++)
    {
      tmp_image[j*(*x_size+2*border)+border-1-i]=tmp_image[j*(*x_size+2*border)+border+i];
      tmp_image[j*(*x_size+2*border)+ *x_size+border+i]=tmp_image[j*(*x_size+2*border)+ *x_size+border-1-i];
    }
  *x_size+=2*border;  /* alter image size */
  *y_size+=2*border;
  *in=tmp_image;      /* repoint in */
}

/* }}} */
/* {{{ void susan_smoothing(three_by_three,in,dt,x_size,y_size,bp) */

void susan_smoothing(int three_by_three,uchar *in,float dt,int x_size,int y_size,uchar *bp)
  {
/* {{{ vars */
float temp;
int   n_max, increment, mask_size,
      i,j,x,y,area,brightness,tmp,centre;
uchar *ip, *dp, *dpt, *cp, *out=in,
      *tmp_image;
TOTAL_TYPE total;
/* }}} */
  /* {{{ setup larger image and border sizes */
  if (three_by_three==0)
    mask_size = ((int)(1.5 * dt)) + 1;
  else
    mask_size = 1;
  total=0.1; /* test for total's type */
  if ( (dt>15) && (total==0) )
  {
    printf("Distance_thresh (%f) too big for integer arithmetic.\n",dt);
    printf("Either reduce it to <=15 or recompile with variable \"total\"\n");
    printf("as a float: see top \"defines\" section.\n");
    exit(0);
  }
  if ( (2*mask_size+1>x_size) || (2*mask_size+1>y_size) )
  {
    printf("Mask size (1.5*distance_thresh+1=%d) too big for image (%dx%d).\n",mask_size,x_size,y_size);
    exit(0);
  }
  tmp_image = (uchar *) malloc( (x_size+mask_size*2) * (y_size+mask_size*2) );
  enlarge(&in,tmp_image,&x_size,&y_size,mask_size);
/* }}} */
  if (three_by_three==0)
  {     /* large Gaussian masks */
    /* {{{ setup distance lut */
  n_max = (mask_size*2) + 1;
  increment = x_size - n_max;
  dp     = (unsigned char *)malloc(n_max*n_max);
  dpt    = dp;
  temp   = -(dt*dt);

  for(i=-mask_size; i<=mask_size; i++)
    for(j=-mask_size; j<=mask_size; j++)
    {
      x = (int) (100.0 * exp( ((float)((i*i)+(j*j))) / temp ));
      *dpt++ = (unsigned char)x;
    }
/* }}} */
    /* {{{ main section */
  for (i=mask_size;i<y_size-mask_size;i++)
  {
    for (j=mask_size;j<x_size-mask_size;j++)
    {
      area = 0;
      total = 0;
      dpt = dp;
      ip = in + ((i-mask_size)*x_size) + j - mask_size;
      centre = in[i*x_size+j];
      cp = bp + centre;
      for(y=-mask_size; y<=mask_size; y++)
      {
        for(x=-mask_size; x<=mask_size; x++)
	{
          brightness = *ip++;
          tmp = *dpt++ * *(cp-brightness);
          area += tmp;
          total += tmp * brightness;
        }
        ip += increment;
      }
      tmp = area-10000;
      if (tmp==0)
        *out++=median(in,i,j,x_size);
      else
        *out++=((total-(centre*10000))/tmp);
    }
  }

/* }}} */
  }
  else
  {     /* 3x3 constant mask */
    /* {{{ main section */
  for (i=1;i<y_size-1;i++)
  {
    for (j=1;j<x_size-1;j++)
    {
      area = 0;
      total = 0;
      ip = in + ((i-1)*x_size) + j - 1;
      centre = in[i*x_size+j];
      cp = bp + centre;

      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
      brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
      ip += x_size-2;
      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
      brightness=*ip++; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
      brightness=*ip; tmp=*(cp-brightness); area += tmp; total += tmp * brightness;
      ip += x_size-2;