# include "data.h"

int extract (int argc, char **argv) {
  
  int i, j;
  float *Vin, *Vout;
  int sx, sy, nx, ny, NX, NY;
  int Sx, Sy, Nx, Ny;
  Buffer *in, *out;

  if (argc != 11) {
    fprintf (stderr, "USAGE: extract <from> <to> sx sy nx ny Sx Sy Nx Ny\n");
    return (FALSE);
  }

  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
  NX = in[0].matrix.Naxis[0];
  NY = in[0].matrix.Naxis[1];

  sx = atof (argv[3]);
  sy = atof (argv[4]);
  nx = atof (argv[5]);
  ny = atof (argv[6]);

  Sx = atof (argv[7]);
  Sy = atof (argv[8]);
  Nx = atof (argv[9]);
  Ny = atof (argv[10]);

  if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
    fprintf (stderr, "mismatch between source and dest regions\n");
    fprintf (stderr, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx);
    return (FALSE);
  }

  /* region is not on first image */
  if ((sx + nx < 0) || (sy + ny < 0) || 
      (sx > in[0].matrix.Naxis[0]) || 
      (sy > in[0].matrix.Naxis[1])) {
    fprintf (stderr, "region outside of source image\n");
    return (FALSE);
  }

  if ((Sx + nx > Nx) || (Sy + ny > Ny)) {
    fprintf (stderr, "source region larger than dest region\n");
    return (FALSE);
  }
  if ((Sx < 0) || (Sy < 0)) {
    fprintf (stderr, "dest region out of range\n");
    return (FALSE);
  }

  if ((out = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) {
    if ((out = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    fits_free_matrix (&out[0].matrix);
    fits_free_header (&out[0].header);

    out[0].bitpix = in[0].bitpix;
    out[0].unsign = in[0].unsign;
    out[0].bscale = in[0].bscale;
    out[0].bzero  = in[0].bzero;
    fits_copy_header (&in[0].header, &out[0].header);
    fits_modify (&out[0].header, "NAXIS1", "%d", 1, Nx);
    fits_modify (&out[0].header, "NAXIS2", "%d", 1, Ny);
    out[0].header.Naxis[0] = Nx;
    out[0].header.Naxis[1] = Ny;
    fits_create_matrix (&out[0].header, &out[0].matrix);
  } else {
    if ((out[0].header.Naxis[1] != Ny) || (out[0].header.Naxis[0] != Nx)) {
      fprintf (stderr, "matrix sizes mis-matched\n");
      fprintf (stderr, "%d x %d  vs  %d x %d\n", Nx, Ny, 
	       out[0].header.Naxis[0], out[0].header.Naxis[1]);
      return (FALSE);
    }
  }

  for (j = 0; j < ny; j++) {
    if (j + sy < 0) continue;
    if (j + sy >= NY) continue;
    Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx;  
    Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
    for (i = 0; i < nx; i++, Vin++, Vout++) {
      if (i + sx < 0) continue;
      if (i + sx >= NX) continue;
      *Vout = *Vin;
    }
  }

  return (TRUE);

}

