#include <stdio.h>
#include <iostream>
#include <octave/oct.h>
#include <octave/file-ops.h>
#include <octave/pt-plot.h>
#include <octave/parse.h>
#include <octave/toplev.h>

static double mmin( const Matrix &m ) {
  double d;
  int i,j,nr,nc;

  nr = m.rows();
  nc = m.columns();

  if( nr==0 || nc==0 ) return 0.0;

  d = m(0,0);

  for(i = 0; i < nr; i++ ) {
    for( j = 0; j < nc; j++ ) {
      d = (m(i,j)<d) ? m(i,j) : d;
    }
  }
  return d;
}

static double mmax( const Matrix &m ) {
  double d;
  int i,j,nr,nc;

  nr = m.rows();
  nc = m.columns();

  if( nr==0 || nc==0 ) return 0.0;

  d = m(0,0);

  for(i = 0; i < nr; i++ ) {
    for( j = 0; j < nc; j++ ) {
      d = (m(i,j)>d) ? m(i,j) : d;
    }
  }
  return d;
}


DEFUN_DLD( __cntr__, args, , 
    "__cntr__ (x,y,z,levels): plot contour lines of a function\n\n"
    "Plot contour lines by interpolating the values of a function\n"
    "on a grid.  If x,y,z are matrices with the same dimensions\n"
    "the the corresponding elements represent points on the function.\n"
    "levels contains a vector of heights at which to draw level curves.\n" )
{
  octave_value_list retval;
  int nargin = args.length();
  int nr, nc, nl, i, j, k, count;
  double x11, x21, x12, x22,
  y11, y21, y12, y22,
  z11, z21, z12, z22;
  double c, x[8], y[8];
  double minx, miny, maxx, maxy;
  std::string fname;
  FILE *f;
  char cmd[256];
  int status;

  if( nargin != 4 ) {
    print_usage( "cnt" );
    return retval;
  }

  Matrix xx = args(0).matrix_value();
  Matrix yy = args(1).matrix_value();
  Matrix zz = args(2).matrix_value();
  ColumnVector level = ColumnVector( args(3).vector_value());

  nl = level.length();
  nr = xx.rows();
  nc = xx.columns();

  minx = mmin(xx);
  maxx = mmax(xx);
  miny = mmin(yy);
  maxy = mmax(yy);

  for( k = 0; k < nl; k++ ) {
    c = level(k);
    fname = tempnam ("", "oct-");
    f = fopen( fname.c_str(), "w" );
    mark_for_deletion( fname );

    for( i = 0; i < nr-1; i++ ) {
      for( j = 0; j < nc-1; j++ ) {
	x11 = xx(i,j );
	x21 = xx(i,j+1);
	x12 = xx(i+1,j);
	x22 = xx(i+1,j+1);

	y11 = yy(i,j );
	y21 = yy(i,j+1);
	y12 = yy(i+1,j);
	y22 = yy(i+1,j+1);

	z11 = zz(i,j );
	z21 = zz(i,j+1);
	z12 = zz(i+1,j);
	z22 = zz(i+1,j+1);

	count = 0;

	// find bottom intersection with contour
	if(( z11>c && z21<c ) || ( z11<c && z21>c )) {
	  x[count] = (c - z11)*(x21-x11)/(z21-z11) + x11;
	  y[count] = (c - z11)*(y21-y11)/(z21-z11) + y11;
	  count++;
	}

	// find top intersection with contour line
	if(( z12>c && z22<c ) || ( z12<c && z22>c )) {
	  x[count] = (c - z12)*(x22-x12)/(z22-z12) + x12;
	  y[count] = (c - z12)*(y22-y12)/(z22-z12) + y12;
	  count++;
	}

	// find left intersection with contour line
	if(( z11>c && z12<c ) || ( z11<c && z12>c )) {
	  x[count] = (c - z11)*(x12-x11)/(z12-z11) + x11;
	  y[count] = (c - z11)*(y12-y11)/(z12-z11) + y11;
	  count++;
	}

	// find right intersection with contour line
	if(( z21>c && z22<c ) || ( z21<c && z22>c )) {
	  x[count] = (c - z21)*(x22-x21)/(z22-z21) + x21;
	  y[count] = (c - z21)*(y22-y21)/(z22-z21) + y21;
	  count++;
	}

	// bottom left corner
	if( z11 == c ) {
	    x[count] = x11;
	    y[count] = y11;
	    count++;
	}

	// bottom right corner
	if( z21 == c ) {
	    x[count] = x21;
	    y[count] = y21;
	    count++;
	}

	// top left corner
	if( z12 == c ) {
	    x[count] = x12;
	    y[count] = y12;
	    count++;
	}

	// top right corner
	if( z22 == c ) {
	    x[count] = x22;
	    y[count] = y22;
	    count++;
	}

	// Only draw the contour if it crosses exactly two points on the boundary
	if( count == 2 ) {
	  fprintf( f, "%lf %lf\n", x[0], y[0] );
	  fprintf( f, "%lf %lf\n", x[1], y[1] );
	  fprintf( f, "\n" );
	}
      }
    }
    fclose( f );

    snprintf( cmd, sizeof(cmd), "gplot \"%s\" title \"%g\"", fname.c_str(), c );
    eval_string( std::string(cmd), false, status );

    if( k == 0 ) eval_string( std::string("hold on"), false, status );
  }

  return retval;
}