#include <stdio.h>
#include <stdlib.h>
#define _INCLUDE_XOPEN_SOURCE /* for DBL_MAX on HP */
#include <math.h>
#include <limits.h>
#include <float.h>
#include "ppm.h"
#include "paraboloid.h"

#ifndef DBL_MAX /* for DBL_MAX in Linux */
	#define DBL_MAX = 0
#endif

#define IDX(i, j, im) (i + j * im->width)
#define SQR(a) (a * a)
#define Normalize(x, y) \
  {double n = sqrt(x * x + y * y); x /= n, y /= n;}

#define DUPLICATE_THRESHOLD 2.0
#define N_RESERVE_PCT 1.20 /* keep 120% of requested number of points */

typedef struct 
{
  double u, v, val, dx, dy, refineDist;
} 
Response;

/* these used to be constants, but are now user-settable variables */
int WINDOW = 1;
double SIGMA = 1.0;

/* for computing derivatives and weighting */
double Gauss(double u)
{
    return exp(-0.5 * SQR(u) / SQR(SIGMA));
}

double DerivOfGauss(double u)
{
    return u * Gauss(u) / -SQR(SIGMA);
}

int InterpVerbosity = 0;
#define ONESIXTH  0.16666666666666666666
#define TWOTHIRDS 0.66666666666666666666
double CubicUniformBMatrix[4][4] = 
{
  {-ONESIXTH, 0.5, -0.5, ONESIXTH},
  { 0.5, -1.0, 0.5, 0.0},
  {-0.5,  0.0, 0.5, 0.0},
  {ONESIXTH, TWOTHIRDS, ONESIXTH, 0.0}
};

double CatmullRomMatrix[4][4] = /* Initially, letting c = 0.5 */
{ 
  {-0.5,  1.5, -1.5,  0.5},  /* -c 2-c  c-2  c */
  { 1.0, -2.5,  2.0, -0.5},  /* 2c c-3 3-2c -c */
  {-0.5,  0.0,  0.5,  0.0},  /* -c  0    c   0 */
  { 0.0,  1.0,  0.0,  0.0}   /*  0  1    0   0 */
};

/* Using (x,fx) pairs, find the interpolant at v (x[1] <= v <= x[2]) */
/* Currently using Catmull-Rom curve for interpolation */
#ifndef CUBE
#define CUBE(a) (a * a * a)
#endif

double Interp(double* x, double* fx, double v)
{
    int i, j;
    double p[4], vv;

    if( v < x[ 1 ] || v > x[ 2 ] ) {
      printf( "illegal parameters to Interp!\n" );
      return 0.0;
    }

    for( i = 0; i < 4; i++ ) {
      for( p[ i ] = 0.0, j = 0; j < 3; j++ )
	p[ i ] += CatmullRomMatrix[ i ][ j ] * fx[ j ];
    }

    vv = v - x[ 1 ];
#ifdef DEBUG
    if( InterpVerbosity == 1 ) {
      printf( "Interpolating samples:\n" );
      for( i = 0; i < 4; i++ )
	printf( "%8.4lf %8.4lf\n", x[ i ], fx[ i ] );
      printf( "P( t ) = %8.4lf * ( x - %d )**3 + %8.4lf * ( x - %d )**2 + %8.4lf * ( x - %d ) + %8.4lf",
	      p[ 0 ], (int) x[ 1 ], p[ 1 ], (int) x[ 1 ], p[ 2 ], (int) x[ 1 ], p[ 3 ] );
      printf( ": %8.4lf   %8.4lf\n", v,
	      CUBE( vv ) * p[ 0 ] + SQR( vv ) * p[ 1 ] + vv * p[ 2 ] + p[ 3 ] );
    }
#endif /* DEBUG */

    return CUBE( vv ) * p[ 0 ] + SQR( vv ) * p[ 1 ] + ( vv ) * p[ 2 ] + p[ 3 ];
}
#undef DEBUG

double
Harris( MyImage *im, int r, int c, double k,
	double *h_im, double *xGrad, double *yGrad )
{
    static int printIntens = 1;
    int i, j;
    double dx, dy, A, B, C, H, H_minus, **I;
    double det, tr;
    unsigned char graylevel;
    PixelType p;
    double *g, *DoG, sg, sdg;

    /* set up */
    g = (double *) malloc(( WINDOW + 1 ) * sizeof( double ));
    DoG = (double *) malloc(( WINDOW + 1 ) * sizeof( double ));
    sg = g[ 0 ] = 1.0; sdg = DoG[ 0 ] = 0.0;
    for( i = 1; i <= WINDOW; i++ ) {
      sg += 2.0 * ( g[ i ] = Gauss( (double) i ));
	sdg += 2.0 * ( DoG[ i ] = DerivOfGauss( (double) i ));
    }
    for( i = 0; i <= WINDOW; i++ ) {
      g[ i ] /= sg;
      DoG[ i ] /= sdg;
    }

    I = (double **) malloc(( 2 * WINDOW + 1 ) * sizeof( double * ));
    for( i = 0; i < 2 * WINDOW + 1; i++ ) {
      I[ i ] = (double *) malloc(( 2 * WINDOW + 1 ) * sizeof( double ));
      for( j = 0; j < 2 * WINDOW + 1; j++ )
	I[ i ][ j ] = 0.0;
    }

    /* make sure we don't go beyond bounds */
    if( r < WINDOW || r > im->width - ( WINDOW + 1 ) ||
	c < WINDOW || c > im->height - ( WINDOW + 1 ))
      return 0.0;

    /* compute the intensity values needed */
    for( i = -WINDOW; i <= WINDOW; i++ )
      for( j = -WINDOW; j <= WINDOW; j++ ) {
	p = im->data[ IDX( r + i, c + j, im )];
	I[ WINDOW + j ][ WINDOW + i ] = ( 0.2990 * (double) GET_RED( p ) +
					  0.5866 * (double) GET_GREEN( p ) +
					  0.1144 * (double) GET_BLUE( p ))
	  / (double) im->maxvalue;
      } /* end for( j, i ) */

#ifdef DEBUG
    if(( r >= rBEGIN && r <= rEND ) && ( c >= cBEGIN && c <= cEND )) {
      printf( "dx g  =" );
      for( i = 0; i < WINDOW + 1; i++ )
	printf( " %13.10lf", g[ i ] );
      printf( "\n" );
      printf( "dx g' =" );
      for( i = 0; i < WINDOW + 1; i++ )
	printf( " %13.10lf", DoG[ i ] );
      printf( "\n" );
    }
#endif /* DEBUG */

    /* compute the derivatives in x and y */
    for( dx = 0.0, dy = 0.0, i = -WINDOW; i <= WINDOW; i++ ) {
      for( j = 1; j <= WINDOW; j++ ) {
        dx += ( I[ WINDOW + i ][ WINDOW + j ] -
                I[ WINDOW + i ][ WINDOW - j ] ) * g[ abs( i )] * DoG[ j ];
        dy += ( I[ WINDOW + j ][ WINDOW + i ] -
                I[ WINDOW - j ][ WINDOW + i ] ) * g[ abs( i )] * DoG[ j ];
#ifdef DEBUG
        if(( r >= rBEGIN && r <= rEND ) && ( c >= cBEGIN && c <= cEND )) {
          printf( "dx %2d %2d = %13.10lf %13.10lf\n", WINDOW + i, WINDOW + j,
                  g[ abs( i )], DoG[ j ] );
	  printf( "dx %2d %2d = %13.10lf %13.10lf\n", WINDOW + i, WINDOW - j,
                  g[ abs( i )], -DoG[ j ] );
	  printf( "dy %2d %2d = %13.10lf %13.10lf\n", WINDOW + j, WINDOW + i,
                  DoG[ j ], g[ abs( i )] );
	  printf( "dy %2d %2d = %13.10lf %13.10lf\n", WINDOW - j, WINDOW + i,
                  -DoG[ j ], g[ abs( i )] );
	}
#endif /* DEBUG */
      } /* end for( j ) */
      if( xGrad != NULL )
	xGrad[ IDX( r, c, im )] = dx;
      if( yGrad != NULL )
	yGrad[ IDX( r, c, im )] = dy;
#define rBEGIN 414
#define rEND 422
#define cBEGIN 428
#define cEND 436
#ifdef DEBUG
      if(( r >= rBEGIN && r <= rEND ) && ( c >= cBEGIN && c <= cEND )) {
	printf( "dx %2d  2 = %13.10lf %13.10lf\n",
		WINDOW + i, g[ abs( i )], 0.0 );
	printf( "dy  2 %2d = %13.10lf %13.10lf\n",
		WINDOW + i, 0.0, g[ abs( i )] );
      }
#endif /* DEBUG */
    }

    /* compute the Harris function */
    A = SQR( dx );
    B = SQR( dy );
    C = dx * dy;
    det = A * B + SQR( C ); /* NOT the determinant anymore!! */
    tr = A + B;
    H_minus = 1000.0 * (( A * B - SQR( C )) - k * SQR( tr ));
    H = 1000.0 * ( det - k * SQR( tr ));
    if( h_im != NULL )
      h_im[ IDX( r, c, im )] = H;

/*     if( r > 156 && r < 162 && c > 64 && c < 70 ); */ /* for rotated image */
#ifdef DEBUG
    if(( r >= rBEGIN && r <= rEND ) && ( c >= cBEGIN && c <= cEND )) {
#if 1
      if( printIntens &&
	  r == ( rBEGIN + rEND ) / 2 && c == ( cBEGIN + cEND ) / 2 ) {
	printf( "( %3d, %3d )-------------------------------------\n", r, c );
	for( i = rBEGIN; i <= rEND; i++ ) {
	  for( j = cBEGIN; j <= cEND; j++ ) {
	    p = im->data[ IDX( i, j, im )];
	    printf( "%8d", (int) ( 0.2990 * (double) GET_RED( p ) +
				   0.5866 * (double) GET_GREEN( p ) +
				   0.1144 * (double) GET_BLUE( p ) ));
	  }
	  printf( "\n" );
	}
	printIntens = 0;
      }
      /*
      printf( "( %3d, %3d ): A = %g\tB = %g\tC = %g\tdet = %g\ttr = %g\n",
	      r, c, A, B, C, det, tr );
      */
#endif
      /*
      printf( "dx( %3d, %3d ) = %12g\t\tdy = %12g\n", r, c, dx, dy );
      printf( "Harris-( %3d, %3d ) = %8.4lf, I = %12g\n",
	      r, c, H_minus, I[ WINDOW ][ WINDOW ] * 255.0 );
      printf( "Harris+( %3d, %3d ) = %8.4lf, I = %12g\n",
	      r, c, H, I[ WINDOW ][ WINDOW ] * 255.0 );
      */
      if( c == cBEGIN )
	printf( "\n" );
      printf( "H( %3d, %3d ) = %8.4lf", r, c, H );
    }
#endif /* DEBUG */

    /* clean up */
    free( g );
    free( DoG );
    for( i = 0; i < 2 * WINDOW + 1; i++ )
      free( I[ i ] );
    free( I );

    return H;
}
#undef DEBUG

/* is val[i][j] greater than all eight neighbors that are not border pix? */
int
MaxAmongNeighbors( double *val, int i, int j, MyImage *im, int border )
{
    int m = 1;
    double my;

    my = val[ IDX( i, j, im )];

    if( i > border ) {
      if( j > border && val[ IDX( i - 1, j - 1, im )] > my )
	m = 0;
      if( val[ IDX( i - 1, j, im )] > my )
	m = 0;
      if( j < im->height - border - 1 && val[ IDX( i - 1, j + 1, im )] > my )
	m = 0;
    } else
      m = 0;

    if( j > border ) {
      if( val[ IDX( i, j - 1, im )] > my )
	m = 0;
    } else
      m = 0;
    if( j < im->height - border - 1 ) {
      if( val[ IDX( i, j + 1, im )] > my )
	m = 0;
    } else
      m = 0;

    if( i < im->width - border - 1 ) {
      if( j > border && val[ IDX( i + 1, j - 1, im )] > my )
	m = 0;
      if( val[ IDX( i + 1, j, im )] > my )
	m = 0;
      if( j < im->height - border - 1 && val[ IDX( i + 1, j + 1, im )] > my )
	m = 0;
    } else
      m = 0;

    return m;
}

/* add v to sorted list if it is greater than the N values there now */
void
Insert( Response *list, int N, double newval, int r, int c )
{
    int i, j;

    for( i = 0; i < N && newval < list[ i ].val; i++ )
      /* EMPTY */;

    if( i < N - 1 ) {
      /* scoot rest of list down */
      for( j = N - 2; j >= i; j-- ) {
	list[ j + 1 ].val = list[ j ].val;
	list[ j + 1 ].u   = list[ j ].u;
	list[ j + 1 ].v   = list[ j ].v;
	/* ignoring dx, dy for now ... */
      }
    }

    if( i < N ) { /* puts in middle or at very end, as appropriate */
      list[ i ].val = newval;
      list[ i ].u   = (double) r;
      list[ i ].v   = (double) c;
      /* ... here too */
    }

    return;
}

void
Refine( Response *H, double *imgVals, double *dxIm, double *dyIm, MyImage *im )
{
    int uI, vI, i, j;
    double x[9], y[9], z[9];
#ifdef CRderiv /* see below */
    double gSamp[4], dxValA[4], dxValB[4], dyValA[4], dyValB[4];
#endif
    double deriv, ux, lx, uy, ly;

    uI = (int) H->u, vI = (int) H->v;

    for( i = 0; i < 3; i++ )
      for( j = 0; j < 3; j++ ) {
	x[ i + j * 3 ] = (double) ( uI + i - 1 );
	y[ i + j * 3 ] = (double) ( vI + j - 1 );
	z[ i + j * 3 ] = imgVals[ IDX( uI + i - 1, vI + j - 1, im )];
      }

#ifdef DEBUG
    H->val = paraboloid( x, y, z, 9, &H->u, &H->v, 0 );
/* 			 ( uI == 135 && vI == 77 ) ? 1 : 0 ); */
    if( uI == 135 && ( vI == 77 || vI == 90 )) {
      printf( "refine [ %3d, %3d ]:\npoints\n", uI, vI );
      for( i = 0; i < 9; i++ )
	printf( "%8.4lf %8.4lf %8.4lf\n", x[ i ], y[ i ], z[ i ] );
      printf( "max = %8.4lf %8.4lf %8.4lf\n", H->u, H->v, H->val );
    }
#else
    H->val = paraboloid( x, y, z, 9, &H->u, &H->v, 0 );
#endif
    H->refineDist = sqrt( SQR( H->u - (double) uI ) +
			  SQR( H->v - (double) vI ));

#ifdef DEBUG
    if( uI == 135 && ( vI == 77 || vI == 90 ))
      InterpVerbosity = 1;
    else
      InterpVerbosity = 0;
#endif
    if( dxIm != NULL && dyIm != NULL ) {
      uI = (int) H->u, vI = (int) H->v;
#ifdef CRderiv /* Catmull-Rom Interpolation for derivatives -- why?!? */
      for( i = -1; i < 3; i++ )
	gSamp[ 1 + i ] = (double) ( uI + i );
      for( j = -1; j < 3; j++ ) {
	for( i = -1; i < 3; i++ ) {
	  dxValA[ 1 + i ] = dxIm[ IDX( uI + 1 + i, vI + 1 + j, im )];
	  dyValA[ 1 + i ] = dyIm[ IDX( uI + 1 + i, vI + 1 + j, im )];
	}
	dxValB[ 1 + j ] = Interp( gSamp, dxValA, H->u );
	dyValB[ 1 + j ] = Interp( gSamp, dyValA, H->u );
      } /* end for( j ) */
      for( j = -1; j < 3; j++ )
	gSamp[ 1 + j ] = (double) ( vI + j );
      H->dx = Interp( gSamp, dxValB, H->v );
      H->dy = Interp( gSamp, dyValB, H->v );
#elif defined( GaussDeriv )
      /* Gaussian interpolation of derivatives; use same WINDOW and SIGMA */
      /* compute weighted sum of X derivative values */
      sum_weight = 0.0, deriv = 0.0;
      for( i = -WINDOW; i <= WINDOW; i++ )
	for( j = -WINDOW; j <= WINDOW; j++ ) {
	  weight = Gauss( (double) i ) * Gauss( (double) j );
	  deriv += dxIm[ IDX( uI + i, vI + j, im ) ] * weight;
	  sum_weight += weight;
	}
      H->dx = deriv / sum_weight;
      /* compute weighted sum of Y derivative values */
      sum_weight = 0.0, deriv = 0.0;
      for( i = -WINDOW; i <= WINDOW; i++ )
	for( j = -WINDOW; j <= WINDOW; j++ ) {
	  weight = Gauss( H->u - (double) ( uI + i )) *
	           Gauss( H->v - (double) ( vI + j ));
	  deriv += dyIm[ IDX( uI + i, vI + j, im ) ] * weight;
	  sum_weight += weight;
	}
      H->dy = deriv / sum_weight;
      if(( uI > 133 && uI < 139 &&
	   ( vI > 74 && vI < 80 || vI > 88 && vI < 94 )) ||
	 ( uI > 146 && uI < 152 &&
	   ( vI > 74 && vI < 80 || vI > 88 && vI < 94 ))) {
	printf( "derivative at ( %lf, %lf ) = %lf, %lf\n",
		H->u, H->v, H->dx, H->dy );
      }
#else
      /* Bilerp derivatives */
      lx = H->u - (double) uI;  ux = 1.0 - lx;
      ly = H->v - (double) vI;  uy = 1.0 - ly;
      H->dx  = dxIm[ IDX( uI    , vI    , im ) ] * ux * uy;
      H->dx += dxIm[ IDX( uI + 1, vI    , im  ) ] * lx * uy;
      H->dx += dxIm[ IDX( uI    , vI + 1, im  ) ] * ux * ly;
      H->dx += dxIm[ IDX( uI + 1, vI + 1, im  ) ] * lx * ly;
      H->dy  = dyIm[ IDX( uI    , vI    , im  ) ] * ux * uy;
      H->dy += dyIm[ IDX( uI + 1, vI    , im  ) ] * lx * uy;
      H->dy += dyIm[ IDX( uI    , vI + 1, im  ) ] * ux * ly;
      H->dy += dyIm[ IDX( uI + 1, vI + 1, im  ) ] * lx * ly;

#ifdef DEBUG
      if(( uI > 133 && uI < 139 &&
	   ( vI > 74 && vI < 80 || vI > 88 && vI < 94 )) ||
	 ( uI > 146 && uI < 152 &&
	   ( vI > 74 && vI < 80 || vI > 88 && vI < 94 ))) {
	printf( "derivative at ( %lf, %lf ) = %lf, %lf, atan = %lf\n",
		H->u, H->v, H->dx, H->dy,
		atan2( H->dy, H->dx ) / M_PI * 180.0 );
      }
#endif /* DEBUG */
#endif /* choice of derivative interpolation method */
#if 0 /* use paraboloid interpolation for derivatives ??
	 functions don't look very paraboloidic.... */
      for( i = 0; i < 3; i++ )
	for( j = 0; j < 3; j++ ) {
	  x[ i + j * 3 ] = (double) ( uI + i - 1 );
	  y[ i + j * 3 ] = (double) ( vI + j - 1 );
	  z[ i + j * 3 ] = dxIm[ IDX( uI + i - 1, vI + j - 1, im )];
	}
      /* this is wrong, see above for a correct call; we want something
	 different here; not "find minimum" but "find paraboloid and eval
	 at these coords" */
      H->dx = paraboloid( x, y, z, 9, H->u, H->v, 0 );
#endif
      Normalize( H->dx, H->dy );
    }

    return;
}
#undef DEBUG

PixelType *
ScaleData( double *d, int w, int h, double min, double max )
{
    int i;
    double *s, high, low;
    PixelType *r;

    if( max > 255.0 || min < 0.0 || max <= min ) {
      printf( "error in data scaling parameters: min = %lf, max = %lf\n",
	      min, max );
      exit( -2 );
    }

    for( low = high = d[ 0 ], i = 1; i < w * h; i++ ) {
      if( low > d[ i ] )
	low = d[ i ];
      if( high < d[ i ] )
	high = d[ i ];
    }
#ifdef DEBUG
    printf( "scale data: low = %lf, high = %lf\n", low, high );
#endif

    s = (double *) malloc( w * h * sizeof( double ));
    r = (PixelType *) malloc( w * h * sizeof( PixelType ));

    for( i = 0; i < w * h; i++ ) {
      s[ i ] = min + ( max - min ) * ( d[ i ] - low ) / ( high - low );
      SET_RED( r[ i ], (unsigned char) ( s[ i ] + 0.5 ));
      SET_GREEN( r[ i ], (unsigned char) ( s[ i ] + 0.5 ));
      SET_BLUE( r[ i ], (unsigned char) ( s[ i ] + 0.5 ));
    }

    free( s );
    return r;
}

void
Sort( Response *list, int N )
{
    int i, j;
    double d;
    Response tmp;

    /*
    ** First, pass through the list and punish any point that is close to
    ** another point that didn't move as much during refinement.  This is
    ** a heuristic to eliminate points that are created by aliasing near
    ** a true corner.
    */
    for( i = 0; i < N ; i++ ) {
      for( j = 0; j < N; j++ ) {
	if( j != i ) {
	  d = sqrt( SQR( list[ i ].u - list[ j ].u ) +
		    SQR( list[ i ].v - list[ j ].v ));
	  if( d < DUPLICATE_THRESHOLD &&
	      list[ i ].refineDist < list[ j ].refineDist ) {
	    list[ i ].val = 0.0;
	  } /* if( close to other point and moved more ) */
	} /* if not comparing point to itself */
      } /* for( j ) */
    } /* for( i ) */

    /*
    ** Second, perform bubble sort on list.  This is the best choice, since
    ** the list was created in order and can only have changed a little with
    ** the changes in responses during the refinement procedure.
    */
    for( i = 0; i < N - 1; i++ ) {
      for( j = N; j > i; j-- ) {
	if( list[ j ].val > list[ j - 1 ].val ) {
	  /* swap */
	  tmp = list[ j ];
	  list[ j ] = list[ j - 1 ];
	  list[ j - 1 ] = tmp;
	} /* if */
      } /* for( j ) */
    } /* for( i ) */

    return;
}

void
PrintList( Response *list, int N )
{
    int i;

    for( i = 0; i < N && list[ i ].val > 0.0; i++ )
      printf( "%12.6lf %4d %4d\n", list[ i ].val, list[ i ].u, list[ i ].v );

    return;
}

int
main( int argc, char *argv[] )
{
    int i, j, w, h;
    double *response, *xGrad, *yGrad, *harrisIm;
    char gradFile[80];
    FILE *fp;
    MyImage *im, dxIm, dyIm, hIm;
    Response *strongest;
    /* user-settable parameters */
    int border, N, N_print;
    char *filenameIN, *gradFileBase;
    double harrisK;

    /* set defaults */
    border = 5;
    N_print = 40, N = (int)((double)N_print * N_RESERVE_PCT);
    filenameIN = NULL;
    harrisK = 0.04; /* apparently from Harris, see Pollefeys p. 33 */
    gradFileBase = NULL;
    xGrad = yGrad = NULL;

    /* parse arguments */
    for(i = 1; i < argc; i++)
	 {
    	if(strcmp(argv[i], "-f") == 0)
			filenameIN = argv[++i];
      else if(strcmp(argv[i], "-g") == 0)
			gradFileBase = argv[++i];
      else if(strncmp(argv[ i ], "-b", 2) == 0) 
		{
			if(sscanf(argv[i], "-b%d", &border) != 1)
	  			if(sscanf(argv[++i], "%d", &border) != 1) 
				{
	    			printf("No argument for '-b' option\n");
	    			exit(-1);
	  			}
      } 
		else if(strncmp(argv[i], "-w", 2) == 0) 
		{
			if(sscanf(argv[i], "-w%d", &WINDOW) != 1)
	  			if(sscanf(argv[++i], "%d", &WINDOW) != 1) 
				{
	    			printf("No argument for '-w' option\n");
	    			exit(-1);
	  			}
      } 
		else if(strncmp(argv[i], "-n", 2) == 0) 
		{
			if(sscanf(argv[i], "-n%d", &N) != 1) 
			{
	  			if(sscanf(argv[++i], "%d", &N) != 1) 
				{
	    			printf("No argument for '-n' option\n");
	    			exit(-1);
	  			}
			} 
			else 
			{
	  			N_print = N;
	  			N = (int)((double)N_print * N_RESERVE_PCT);
			}
      } 
		else if(strncmp(argv[i], "-k", 2) == 0) 
		{
			if( sscanf( argv[ i ], "-k%lf", &harrisK ) != 1 )
	  			if( sscanf( argv[ ++i ], "%lf", &harrisK ) != 1 ) 
				{
	    			printf( "No argument for '-k' option\n" );
	    			exit( -1 );
	  			}
      } 
		else if( strncmp( argv[ i ], "-s", 2 ) == 0 ) 
		{
			if( sscanf( argv[ i ], "-s%lf", &SIGMA ) != 1 )
	  			if( sscanf( argv[ ++i ], "%lf", &SIGMA ) != 1 ) 
				{
	    			printf( "No argument for '-s' option\n" );
	    			exit( -1 );
	  			}
      } 
		else printf( "Unknown argument '%s' ignored\n", argv[ i ] );
    } /* end for( i ) */

    /* read the named file */
    if( filenameIN != NULL ) {
      fp = fopen( filenameIN, "r" );
      if( fp == NULL ) {
	printf( "unable to open image '%s'\n", filenameIN );
	exit( -1 );
      } /* else */
      im = readppm( fp );
      fclose( fp );
    } else {
      fp = stdin;
      im = readppm( fp );
    }

    /* allocate an image with values of type double for the filter response */
    response = (double *) malloc( im->width * im->height * sizeof( double ));
    if( response == NULL ) {
      printf( "unable to allocate space\n" );
      exit( -1 );
    }
    memset( (void *) response, 0, im->width * im->height * sizeof( double ));

    /* if user requested it, allocate gradient images */
    if( gradFileBase != NULL ) {
      w = im->width, h = im->height;
      hIm.width = w, hIm.height = h, hIm.maxvalue = 255;
      harrisIm = (double *) malloc( w * h * sizeof( double ));
      memset( (void *) harrisIm, 0, w * h * sizeof( double ));
      dxIm.width = w, dxIm.height = h, dxIm.maxvalue = 255;
      xGrad = (double *) malloc( w * h * sizeof( double ));
      memset( (void *) xGrad, 0, w * h * sizeof( double ));
      dyIm.width = w, dyIm.height = h, dyIm.maxvalue = 255;
      yGrad = (double *) malloc( w * h * sizeof( double ));
      memset( (void *) yGrad, 0, w * h * sizeof( double ));
    }

    /* compute the Harris operator at each pixel, except for border pixels */
    for( i = border; i < im->width - border; i++ )
      for( j = border; j < im->height - border; j++ ) {
	response[ IDX( i, j, im )] = Harris( im, i, j, harrisK,
					     harrisIm, xGrad, yGrad );
      }

    if( gradFileBase != NULL ) {
      hIm.data = ScaleData( harrisIm, w, h, 0.0, 255.0 );
      sprintf( gradFile, "%s_h.ppm", gradFileBase );
      fp = fopen( gradFile, "w" );
      writeppm( fp, &hIm );
      fclose( fp );

      dxIm.data = ScaleData( xGrad, w, h, 0.0, 255.0 );
      sprintf( gradFile, "%s_x.ppm", gradFileBase );
      fp = fopen( gradFile, "w" );
      writeppm( fp, &dxIm );
      fclose( fp );

      dyIm.data = ScaleData( yGrad, w, h, 0.0, 255.0 );
      sprintf( gradFile, "%s_y.ppm", gradFileBase );
      fp = fopen( gradFile, "w" );
      writeppm( fp, &dyIm );
      fclose( fp );
    }

    /* initialize array of strongest responses */
    strongest = (Response *) malloc( N * sizeof( Response ));
    if( strongest == NULL ) {
      printf( "unable to allocate space\n" );
      exit( -1 );
    }
    for( i = 0; i < N; i++ )
      strongest[ i ].val = -DBL_MAX;

    /* find the 110% of N largest responses among pixels that have a stronger
       response than all eight of their neighbors */
    for( i = border; i < im->width - border; i++ )
      for( j = border; j < im->height - border; j++ )
	if( MaxAmongNeighbors( response, i, j, im, border ))
	  /* insertion will be no-op if value is less than last on full list */
	  Insert( strongest, N, response[ IDX( i, j, im )], i, j );

#if 0
    /* report those N locations */
    printf( "%d\n", N_print );
    for( i = 0; i < N_print; i++ ) {
      /*
      ** This +1 for the V coord corrects for a difference in indexing scheme
      ** between the ray trace addressing (base-1) and the image processing
      ** (base-0).  Sigh.
      */
      printf( "%8.4lf  %8.4lf  %8.4lf\n",
	      strongest[ i ].u, strongest[ i ].v + 1.0, strongest[ i ].val );
    }
#endif

    /* for each of those N pixels, compute a subpixel accurate location
       for the local maximum strength response to the filter and the gradient
       at that location */
    for( i = 0; i < N; i++ )
      Refine( &strongest[ i ], response, xGrad, yGrad, im );

    Sort( strongest, N );

    /* report the top N_print locations */
    printf( "%d\n", N_print );
    for( i = 0; i < N_print; i++ ) {
      /*
      ** This +1 for the V coord corrects for a difference in indexing scheme
      ** between the ray trace addressing (base-1) and the image processing
      ** (base-0).  Sigh.
      */
      printf( "%8.4lf  %8.4lf  %8.4lf    %9.5lf %9.5lf\n",
	      strongest[ i ].u, strongest[ i ].v + 1.0, strongest[ i ].val,
	      strongest[ i ].dx, strongest[ i ].dy );
    }

    return 0;
}