/***************************************************************************
 * The Greyscale Compass Operator -- an algorithm to find edges in         *
 *   greyscale images by representing regions with distributions           *
 *                                                                         *
 * Mark A. Ruzon                                                           *
 * Stanford Vision Laboratory                                              *
 * October, 1997                                                           *
 * Copyright 1997-1999.  All rights reserved.                              *
 *                                                                         *
 * Details -- The compass operator uses a circular window centered at a    *
 *   junction where 4 pixel squares meet.  The needle is a diameter at a   *
 *   given orientation.  The intensity distributions of the two            *
 *   semicircles are computed, and the distance between them is measured   *
 *   using the Earth Mover's Distance (EMD).  The maximum EMD over all     *
 *   orientations gives the edge strength and orientation at that point.   *
 *   Uncertainty and abnormality are also measured.                        *
 *                                                                         *
 * This file contains the Unix wrapper for the compass operator.           *
 *                                                                         *
 * 2 August 1999 -- extended for use with greyscale images                 *
 *                                                                         *
 * 7 October 1999 -- creation of the Unix wrapper.  Most of the            *
 *   functionality of the MATLAB version (e.g. subimages, spacing, plot,   *
 *   maxclusters, multiple sigmas, abnormality, uncertainty, etc.)         *
 *   is being removed.                                                     *
 *                                                                         *
 * 21 October 1999 -- allowed only compass operator and NMS routine, or    *
 *                    only thresholding, or both to be computed.           *
 ***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <opencv2/highgui/highgui_c.h>

#include "greycompass.h"

#define MAXVAL 255
#define MAXLINE 2000
#define PI 3.14159265358979323846
#define MODEC 0x2
#define MODET 0x1

#warning this should generate an error since it should be exported
static void Error(char *text)
{
  fprintf(stderr,text);
  exit(1);
}


static void PrintHelp(char *programname)
{
  fprintf(stderr,"Function: %s <infile> <outfile> [options]\n", programname);
  fprintf(stderr,"   <infile>/<outfile> are input/output images\n");
  fprintf(stderr,"   options include:\n");
  fprintf(stderr,"   -s <sigma>, the standard deviation of the Gaussian (default 1.0)\n");
  fprintf(stderr,"   -m C|T, the mode, compass operator or thresholding (default: both)\n"); 
  fprintf(stderr,"   -l <low>, the low threshold for edge detection (default 0.5)\n");
  fprintf(stderr,"   -h <high>, the high threshold for edge detection (default 0.7)\n");
  exit(1);
}


static CvMat* readIntoMat( const char* filename, int type )
{
  CvMat* img = cvLoadImageM( filename, CV_LOAD_IMAGE_GRAYSCALE );
  if( !img ) return NULL;

  if( type == img->type )
    return img;

  CvMat* img_wantedtype = cvCreateMat( img->rows, img->cols, type );
  if( !img_wantedtype ) return NULL;

  cvConvert( img, img_wantedtype );
  cvRelease( &img );

  return img_wantedtype;
}

static void ParseArguments(int argc, char *argv[], CvMat** input,
                           double sigma[], double *low, double *high, int *mode)
{
  int curarg = 3;
  char modechar;

  if (argc < 3)
    PrintHelp(argv[0]);

  while (curarg < argc) {
    if (argv[curarg][0] == '-') {  /* New argument has appeared */
      switch(argv[curarg][1]) {
      case 's': /* Sigma */
	if (curarg < argc) {
	  sigma[0] = strtod(argv[curarg+1],NULL);
	  if (sigma[0] <= 0.0)
	    Error("Sigma must be positive\n");
	} else
	  Error("No sigma value specified\n");
	break;
      case 'l': /* Low */
	if (curarg < argc) {
	  *low = strtod(argv[curarg+1],NULL);
	  if (*low < 0.0 || *low > 1.0)
	    Error("Low threshold must lie in [0,1]\n");
	} else
	  Error("No low threshold specified\n");
	break;
      case 'h': /* High */
	if (curarg < argc) {
	  *high = strtod(argv[curarg+1],NULL);
	  if (*high < 0.0 || *high > 1.0)
	    Error("High threshold must lie in [0,1]\n");
	} else
	  Error("No high threshold specified\n");
	break;
      case 'm': /* Mode */
	if (curarg < argc) {
	  modechar = argv[curarg+1][0];
	  if (modechar == 'C' || modechar == 'c')
	    *mode = MODEC;
	  else if (modechar == 'T' || modechar == 't')
	    *mode = MODET;
	  else
	    Error("Mode must be 'C' or 'T'\n");
	} else
	  Error("No mode specified\n");
	break;	
      default: /* Help */
	fprintf(stderr,"Unrecognized command option\n");
	PrintHelp(argv[0]);
	break;
      }
      curarg += 2;
    }
  }
  if (*high < *low)
    Error("High threshold cannot be below low threshold\n");

  *input = readIntoMat( argv[1], (*mode & MODEC) ? CV_8UC1 : CV_64FC1 );
  if( ! *input )
    Error("Error reading input\n");
}


/*  NonMaximalSuppression
 *
 *  The standard Canny non-maximal suppression algorithm; looks in a 3x3
 *  neighborhood to determine if the pixel is a maximum in the direction
 *  perpendicular to that of the edge orientation.
 */
static Matrix *NonMaximalSuppression(Matrix *strength, Matrix *orientation)
{
  int x, y, i, rows = strength->rows, cols = strength->cols;
  int pixels = strength->rows * strength->cols, maximum;
  double str1, str2;              /* interpolated edge strength */
  double a1, a2, b1, b2, c1, c2;  /* nearest pixels' edge strength */
  float ux, uy;                   /* weights of a, b, c, and str */
  double ori, str;       /* strength and orientation at center */
  Matrix *newstrength;

  /* Newstrength holds the NMS'ed strength values */
  newstrength = (Matrix *)malloc(sizeof(Matrix));
  newstrength->rows = strength->rows;
  newstrength->sheets = 1;
  newstrength->cols = strength->cols;
  newstrength->ptr = (double *)calloc(pixels, sizeof(double));

  return(newstrength);
}