/* 
   "sobel.cc" - Example Sobel code

   Last modified: Simon Baker, 4/16/99.
*/

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#define PI 3.1415927

class ucImage 
{
    int srow, scol;                               /* Starting row, column */
    int nrows, ncols;                          /* Number of rows, columns */
    int nbands;                /* Number of bands: 1 = greyscale, 3 = rgb */
    unsigned char *data;                              /* Holds image data */
    
public:
    /* Constructor and Destructor */
    ucImage();                                     /* Creates empty image */
    ~ucImage();
    
    /* Access Routines */
    int NBands() { return nbands; }
    int SRow()   { return srow; }
    int SCol()   { return scol; }
    int NRows()  { return nrows; }
    int NCols()  { return ncols; }
    int InImage(double r, double c)
	{
	    return (r >= srow && r <= srow + nrows - 1 &&
		    c >= scol && c <= scol + ncols - 1);
	}

    /* Returns exact pixel value for integer pixel coordinates */
    unsigned char *Pixel(int r, int c)
	{
	    assert(r >= srow && r < srow + nrows &&
		   c >= scol && r < scol + ncols);
	    return data+((r-srow)*ncols+(c-scol))*nbands;
	}

    /* Returns interpolated pixel value in array pv[nbands] */
    int Pixel(double* pv, double r, double c, int interp = 1);
        
    /* Reallocation of space */
    void ReAllocate(int sr, int sc, int nr, int nc, int nb = 1);
    void ReAllocate(int nr, int nc, int nb = 1)
	{ ReAllocate(0,0,nr,nc,nb); }

    /* Adds a border to the image */
    void AddBorder(int width, unsigned char intensity = 0);
    
    /* IO Routines */
    int ReadPPM(char *ifile);
    friend istream &operator>>(istream &input, ucImage &im); 
    int WritePPM(char *ofile) const;
    friend ostream &operator<<(ostream &output, const ucImage &im);
};

void usage(char** argv);
int comp(const void* n, const void* m);


int
main(int argc,char** argv){

  ofstream output;
  ifstream input;
  ucImage image;
  
  /*Process command line arguments*/
  if (argc != 3) {
      usage(argv);
  }
  char* ifile2=argv[1];
  int number = atoi(argv[2]);

  /* Read in the Image */
  input.open(ifile2);
  input >> image;
  input.close();

  /* Do the detection */
  int m = image.NRows();
  int n = image.NCols();
  int pixels = n*m;
  double *dat = new double[4*pixels];
  double delta1, delta2;
  for(int i=0; i<n; i++){
      for(int j=0; j<m; j++){
	  dat[4*(i+j*n)+1] = i;
	  dat[4*(i+j*n)+2] = m-1-j;
	  if (i < 1 || j < 1 || i > n - 2 || j > m - 2) {
	      dat[4*(i+n*j)+3] = 0.0;
	      dat[4*(i+n*j)] = 0.0;
	  }
	  else {
	      delta1 = 2*image.Pixel(j,i+1)[0] + image.Pixel(j+1,i+1)[0]
                       + image.Pixel(j-1,i+1)[0] - 2*image.Pixel(j,i-1)[0]
		       - image.Pixel(j+1,i-1)[0] - image.Pixel(j-1,i-1)[0];
	      delta2 = -(2*image.Pixel(j+1,i)[0] + image.Pixel(j+1,i+1)[0]
                       + image.Pixel(j+1,i-1)[0] - 2*image.Pixel(j-1,i)[0]
		       - image.Pixel(j-1,i+1)[0] - image.Pixel(j-1,i-1)[0]);
	      dat[4*(i+j*n)] = sqrt(delta1*delta1+delta2*delta2);
	      dat[4*(i+j*n)+3] = 90+180.0*atan2(delta2,delta1)/PI;
	  }
      }
  }
  
  /* Sort by distance */
  qsort((void*) dat, pixels, 4*sizeof(double), comp);
  for(int i=0; i< ((number < pixels) ? number : pixels); i++) {
      cout << dat[4*i+1] << " " << dat[4*i+2] << " " << dat[4*i+3]
	   << " " << dat[4*i] << "\n";
  }
      
  return 0;
}


void usage(char** argv)
{
  cerr << "Usage: " << argv[0] << " image.pgm number\n";
  exit(1);
}

int comp(const void* n, const void* m)
{
    double *x = (double*) n;
    double *y = (double*) m;

    if (*x > *y)
	return -1;
    if (*x < *y)
	return 1;
    return 0;
}


ucImage::ucImage(void)
{
    nbands = srow = scol = nrows = ncols = 0;
    data = 0;
}

ucImage::~ucImage(void)
{
    delete [] data;
}

void ucImage::ReAllocate(int sr, int sc, int nr, int nc, int nb)
{
    assert(nr > 0 && nc > 0 && nb > 0);
    
    delete [] data;

    srow   = sr;
    scol   = sc;
    nrows  = nr;
    ncols  = nc;
    nbands = nb;

    data = new unsigned char[nb*nr*nc];
}

int ucImage::ReadPPM(char *ifile)
{
    assert(ifile);

    /* Open file */
    ifstream input(ifile);
    if (input.bad()) {
	cerr << "Error opening file: " << ifile << "\n";
	return -1;
    }

    /* Read image */
    input >> *this;

    /* Check status */
    if (input.bad()) {
	cerr << "Error reading file: " << ifile << "\n";
	return -1;
    }
    else {
	return 0;
    }
}

istream &operator>>(istream &input, ucImage &im)
{
    int grey_levels;
    unsigned char dummy=0;
    
    /* Skip any comments */
    input.get(dummy);
    while(dummy == '#' || dummy == '\n') {
	do {
	    input.get(dummy);
	}
	while(dummy != '\n');
	input.get(dummy);
    }

    /* Read in the P5 or P6*/
    if (dummy != 'P') {
	cerr << "Image file format error trying to read P\n";
	return input;
    }
    input.get(dummy);
    if (dummy == '5') {
	im.nbands = 1;
    }
    else if (dummy == '6') {
	im.nbands = 3;
    }
    else {
	cerr << "Image file format error trying to read 5 or 6\n";
	return input;
    }
    input.get(dummy);
    if (dummy != '\n') {
	cerr << "Image file format error trying to read P5 or P6\n";
	return input;
    }

    /* Skip any comments */
    input.get(dummy);
    while(dummy == '#' || dummy == '\n') {
	do {
	    input.get(dummy);
	}
	while(dummy != '\n');
	input.get(dummy);
    }
    input.putback(dummy);
    
    /* Read in the size and grey levels */
    input >> im.ncols >> im.nrows;
    if (im.ncols <= 0 || im.nrows <= 0) {
	cerr << "Negative number of rows or columns\n";
	return input;
    }
    input >> grey_levels;
    input.get(dummy);
    
    /* Read in the image */
    im.ReAllocate(im.nrows, im.ncols, im.nbands);
    input.read(im.data, im.ncols*im.nrows*im.nbands);

    /* Check for errors */
    if (input.bad()) {
	cerr << "File Error reading image\n";
	return input;
    }

    return input;
}

int ucImage::WritePPM(char *ofile) const
{
    assert(ofile);

    /* Open file */
    ofstream output(ofile);
    if (output.bad()) {
	cerr << "Error opening file for writing: " << ofile << "\n";
	return -1;
    }

    /* Write image */
    output << *this;

    /* Check status */
    if (output.bad()) {
	cerr << "Error writing file: " << ofile << "\n";
	return -1;
    }
    else {
	return 0;
    }
}

ostream &operator<<(ostream &output, const ucImage &im)
{
    /* Write the header */
    output << 'P';
    if (im.nbands == 1) {
	output << '5';
    }
    else if (im.nbands == 3) {
	output << '6';
    }
    else {
	cerr << "Wrong number of bands for output to ppm\n";
	return output;
    }
    output << "\n" << im.ncols << " " << im.nrows << "\n255\n";

    /* Write the data */
    output.write(im.data, im.ncols*im.nrows*im.nbands);

    /* Check for errors */
    if (output.bad()) {
	cerr << "File error writing image\n";
	return output;
    }

    return output;
}

int ucImage::Pixel(double* pv, double r, double c, int interp)
{
    switch (interp) {
    case 0:                               /* Nearest pixel interpolation */
    {
	/* Find nearest pixel */
	int row = int(r+0.5);
	int col = int(c+0.5);
	int old_row = row;
	int old_col = col;
	
	/* Check in range */
	row = (row < srow) ? srow : row;
	row = (row >= srow+nrows) ? srow+nrows-1 : row;
	col = (col < scol) ? scol : col;
	col = (col >= scol+ncols) ? scol+ncols-1 : col;

	/* Get nearest pixel value and copy into pv[nbands] */
	unsigned char *pixel =  Pixel(row, col);
	for(int i = 0; i < nbands; i++)
	    pv[i] = pixel[i];

	if (row == old_row && col == old_col) {
	    return 1;
	}
	else {
	    return 0;
	}
    }
    
    case 1:                                    /* Bilinear interpolation */
    {
	/* Find bottom left pixel */
	int row = int(r);
	int col = int(c);

	/* If "grid square" lies entirely in image */ 
	if (row >= srow && row + 1 <= srow + nrows - 1 &&
	    col >= scol && col + 1 <= scol + ncols - 1) {

	    /* Do bilinear interpolation */
	    double t = r - row;
	    double u = c - col;
	    for(int i = 0; i < nbands; i++) {
		pv[i] =
		(1-t)*(1-u)*data[((row-srow)*ncols+(col-scol))*nbands+i] +
		t*(1-u)*data[((row+1-srow)*ncols+(col-scol))*nbands+i] +
		t*u*data[((row+1-srow)*ncols+(col+1-scol))*nbands+i] +
	        (1-t)*u*data[((row-srow)*ncols+(col+1-scol))*nbands+i];
	    }
	    return 1;
	}
	else {
	    return Pixel(pv, r, c, 0);
	}

	break;
    }
    
    case 2:                                     /* Bicubic interpolation */
	cerr << "Bicubic interpolation not implemented yet.\n";
	exit(-1);
	break;
    
    default:
	cerr << "Error: Invalid interpolation option.\n";
	exit(-1);
    }
}

void ucImage::AddBorder(int width, unsigned char intensity)
{
    if (width > 0) {
	/* Allocate new image */
	unsigned char *old_data = data;
	data = new unsigned char[(ncols+2*width)*(nrows+2*width)*nbands];

	/* Copy image */
	for(int i = -width; i < nrows + width; i++) {
	    for(int j = -width; j < ncols + width; j++) {
		for(int k = 0; k < nbands; k++) {
		    if (i < 0 || i >= nrows || j < 0 || j >= ncols) {
			data[(((i+width)*(ncols+2*width))+(j+width))*nbands+k]
			    = intensity;
		    }
		    else {
			data[(((i+width)*(ncols+2*width))+(j+width))*nbands+k]
			    = old_data[((i*ncols)+j)*nbands+k];
		    }
		}
	    }
	}

	/* Update image size */
	nrows += 2*width;
	ncols += 2*width;
	srow -= width;
	scol -= width;

	/* Clean up */
	delete [] old_data;
    }
}