# include "data.h"
# define LL { \
 dx =  d00 / (*v01 - *v00); \
 dy = -d00 / (*v10 - *v00); \
 x = i - 0.5; \
 y = j - dy - 0.5; }

# define UL { \
 tmp = d00 / (*v10 - *v00); \
 dy = 1 - tmp; \
 dx =  d10 / (*v11 - *v10); \
 x = i - 0.5; \
 y = j + tmp - 0.5; }
      
# define LR { \
 tmp = d00 / (*v01 - *v00); \
 dx = 1 - tmp; \
 dy = d01 / (*v11 - *v01); \
 x = i + tmp - 0.5; \
 y = j - 0.5; }

# define UR { \
 tmp = d10 / (*v11 - *v10); \
 dx = 1 - tmp; \
 dy = -d11 / (*v01 - *v11); \
 x = i + tmp - 0.5; \
 y = j + 1 - 0.5; }
      
# define HZ { \
 tmp = d00 / (*v10 - *v00); \
 dy = d01 / (*v11 - *v01) - tmp; \
 dx = 1; \
 x = i - 0.5; \
 y = j + tmp - 0.5; }

# define VT { \
 tmp = d00 / (*v01 - *v00); \
 dx = d10 / (*v11 - *v10) - tmp; \
 x = i + tmp - 0.5; \
 dy = 1; \
 y = j - 0.5; }

# define DUMP { \
sprintf (&buffer[Nline*128], "%14s %20.10f %20.10f %20.10f %20.10f ", "LINE", Npix*x, Npix*y, Npix*dx, Npix*dy); \
Nline ++; \
if (Nline == 512) { \
  sprintf (line, "NLINES  %7d ", Nline); \
  write (Ximage, line, 16); \
  write (Ximage, buffer, Nline*128); \
  bzero (buffer, 65536); \
  Nline = 0; \
} \
continue; }

# define NUMP { \
sprintf (&buffer[Nline*128], "%14s %20.10f %20.10f %20.10f %20.10f ", "LINE", Npix*x, Npix*y, Npix*dx, Npix*dy); \
Nline ++; \
if (Nline == 512) { \
  sprintf (line, "NLINES  %7d ", Nline); \
  write (Ximage, line, 16); \
  write (Ximage, buffer, Nline*128); \
  bzero (buffer, 65536); \
  Nline = 0; \
} }

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

  int i, j, ii, jj, n, Npix, Nx, Ny, Nline;
  float level, d00, d01, d10, d11, tmp;
  float x, y, dx, dy;
  float *v00, *v01, *v10, *v11;
  float *Vout, *Vin, *matrix;
  char *buffer, line[17];
  int Ximage, Nimage, N;
  Buffer *buf;
  
  Nimage = -1;
  if (N = get_argument (argc, argv, "-n")) {
    remove_argument (N, &argc, argv);
    Nimage = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }
  if (!GetImage (&Ximage, &Nimage)) return (FALSE);

  if ((argc != 4) && (argc != 5)) {
    fprintf (stderr, "USAGE: contour <buffer> (overlay) level [Npix]\n");
    return (FALSE);
  }

  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
  if (!SelectOverlay (argv[2], &n)) return (FALSE);

  level = atof (argv[3]);
  if (argc == 5) 
    Npix = (int) MAX (atof (argv[4]), 1); 
  else 
    Npix = 1;
  level *= Npix*Npix;

  /*** make rebinned image ***/
  Nx = buf[0].header.Naxis[0]/Npix;
  Ny = buf[0].header.Naxis[1]/Npix;
  if (Npix != 1) {
    fprintf (GetOutfile(), "rebin by a factor of %d (%d,%d)\n", Npix, Nx, Ny);
    ALLOCATE (matrix, float, Nx*Ny);
    bzero (matrix, Nx*Ny*sizeof(float));
	
    for (j = 0; j < Ny; j++) {
      for (jj = 0; jj < Npix; jj++) {
	Vout = matrix + j*Nx;
	Vin  = (float *)(buf[0].matrix.buffer) + (j*Npix + jj)*buf[0].header.Naxis[0];
	for (i = 0; i < Nx; i++, Vout++) {
	  for (ii = 0; ii < Npix; ii++, Vin++) {
	    *Vout += *Vin;
	  }
	}
      }
    }
  }
  else {
    fprintf (stderr, "using scale of 1\n");
    matrix = (float *)(buf[0].matrix.buffer);
  }

  write (Ximage, "LOAD", 4); /* force Ximage to look for the incoming image */
  sprintf (line, "OVERLAY %7d ", n);
  write (Ximage, line, 16);
  ALLOCATE (buffer, char, 65536);  /* space for 512 lines of 128 bytes */
  bzero (buffer, 65536);
  Nline = 0;

  v00 = matrix;
  v01 = matrix + 1;
  v10 = matrix + Nx;
  v11 = matrix + Nx + 1;
  d01 = (level - *v00);
  d11 = (level - *v10);
  for (j = 1; j < Ny; j++) {
    if (!(j%10)) fprintf (stderr, ".");
    for (i = 1; i < Nx; i++, v00++, v01++, v10++, v11++) {

      d00 = d01;
      d10 = d11;
      d01 = (level - *v01);
      d11 = (level - *v11);

      if (((d00 > 0) && (d01 > 0) && (d10 > 0) && (d11 > 0)) ||
	  ((d00 < 0) && (d01 < 0) && (d10 < 0) && (d11 < 0)))
	continue;

      if ((d00 > 0) && (d10 <= 0)) { 
	if ((d01 <= 0) && (d11 <= 0)) { /* \  */
	  LL;
	  DUMP;
	}
	if ((d01 > 0) && (d11 <= 0)) { /* -  */
	  HZ;
	  DUMP;
	}
	if ((d01 > 0) && (d11 > 0)) { /* /  */
	  UL;
	  DUMP;
	}
	if ((d01 <= 0) && (d11 > 0)) { /* \\  */
	  LL;
	  NUMP;
	  UR;
	  DUMP;
	}
      }

      if ((d00 <= 0) && (d10 > 0)) {
	if ((d01 > 0) && (d11 > 0)) { /* \  */
	  LL;
	  DUMP;
	}
	if ((d01 <= 0) && (d11 > 0)) { /* -  */
	  HZ;
	  DUMP;
	}
	if ((d01 <= 0) && (d11 <= 0)) { /* /  */
	  UL;
	  DUMP;
	}
	if ((d01 > 0) && (d11 <= 0)) { /* //  */
	  UL;
	  NUMP;
	  LR;
	  DUMP;
	}
      }
      

      if ((d00 <= 0) && (d10 <= 0)) { 
	if ((d01 > 0) && (d11 <= 0)) {
	  LR;
	  DUMP;
	}
	if ((d01 <= 0) && (d11 > 0)) {
	  UR;
	  DUMP;
	}
	if ((d01 > 0) && (d11 > 0)) {
	  VT;
	  DUMP;
	}
      }

      if ((d00 > 0) && (d10 > 0)) { 
	if ((d01 <= 0) && (d11 > 0)) {
	  LR;
	  DUMP;
	}
	if ((d01 > 0) && (d11 <= 0)) {
	  UR;
	  DUMP;
	}
	if ((d01 <= 0) && (d11 <= 0)) {
	  VT;
	  DUMP;
	}
      }

    }
    /* skip left-hand edge */
    v00++; v01++; v10++; v11++;
    d01 = (level - *v00);
    d11 = (level - *v10);
  }

  /****** bottom line *******/
  v00 = matrix;  
  v01 = matrix + 1;  
  y = 0;
  dx = 0;
  dy = -0.5;
  for (i = 1; i < Nx; i++, v00++, v01++) { /* do the edges */
    if (((*v00 > level) && (*v01 <= level)) || ((*v00 <= level) && (*v01 > level))) {
      x = i + (level - *v01)/(*v01 - *v00);
      NUMP;
    }
  }

  /********** top line *******/
  v00 = matrix + Nx*(Ny - 1);  
  v01 = v00 + 1;
  y = Ny - 1;
  dx = 0;
  dy = 0.5;
  for (i = 1; i < Nx; i++, v00++, v01++) { /* do the edges */
    if (((*v00 > level) && (*v01 <= level)) || ((*v00 <= level) && (*v01 > level))) {
      x = i + (level - *v01)/(*v01 - *v00);
      NUMP;
    }
  }

  /******** left line *********/
  v00 = matrix; 
  v01 = matrix + Nx;
  x = 0;
  dx = -0.5;
  dy = 0;
  for (j = 1; j < Ny; j++, v00+=Nx, v01+=Nx) { /* do the edges */
    if (((*v00 > level) && (*v01 <= level)) || ((*v00 <= level) && (*v01 > level))) {
      y = j + (level - *v01)/(*v01 - *v00);
      NUMP;
    }
  }

  /******** right line *********/
  v00 = matrix + Nx - 1; 
  v01 = v00 + Nx;
  x = Nx - 1;
  dx = 0.5;
  dy = 0;
  for (j = 1; j < Ny; j++, v00+=Nx, v01+=Nx) { /* do the edges */
    if (((*v00 > level) && (*v01 <= level)) || ((*v00 <= level) && (*v01 > level))) {
      y = j + (level - *v01)/(*v01 - *v00);
      NUMP;
    }
  }
  
  if (Nline > 0) {
    sprintf (line, "NLINES  %7d ", Nline); 
    write (Ximage, line, 16); 
    write (Ximage, buffer, Nline*128); 
  }

  free (buffer);
  sprintf (line, "DONE ");
  write (Ximage, line, 16);
  if (Npix != 1) {
    free (matrix);
  }
  fprintf (stderr, "\n");
  return (TRUE);
}
