# include "data.h"

enum {SUM, MEAN, MEDIAN};

int cut (int argc, char **argv) {
  
  int i, j, N, Nx, Ny, Mode;
  float *Vin, *Vbuf, value;
  int sx, sy, nx, ny;
  Vector *xvec, *yvec;
  Buffer *buf;

  Mode = SUM;
  if ((N = get_argument (argc, argv, "-median"))) {
    remove_argument (N, &argc, argv);
    Mode = MEDIAN;
  }
  if ((N = get_argument (argc, argv, "-mean"))) {
    remove_argument (N, &argc, argv);
    Mode = MEAN;
  }

  if (argc != 9) {
    gprint (GP_ERR, "USAGE: cut <buffer> <X vector> <Y vector> <X|Y> sx sy nx ny\n");
    return (FALSE);
  }

  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
 
  sx = atof (argv[5]);
  sy = atof (argv[6]);
  nx = atof (argv[7]);
  ny = atof (argv[8]);

  Nx = buf[0].matrix.Naxis[0];
  Ny = buf[0].matrix.Naxis[1];

  if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
    gprint (GP_ERR, "region out of range\n");
    return (FALSE);
  }

  if ((xvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((yvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);

  switch (argv[4][0]) {
  case 'x':
  case 'X':
    /* create output vectors */
    ResetVector (xvec, OPIHI_FLT, nx);
    ResetVector (yvec, OPIHI_FLT, nx);
    bzero (yvec[0].elements.Flt, nx*sizeof(opihi_flt));
    for (i = 0; i < nx; i++) {
      xvec[0].elements.Flt[i] = i + sx; 
    }
    ALLOCATE (Vbuf, float, MAX (ny, 1));

    for (i = 0; i < nx; i++) {
      /* accumulate values */
      Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i; 
      for (j = 0; j < ny; j++, Vin += Nx) {
	Vbuf[j] = *Vin;
      }
      /* apply stat of choice */
      if (Mode == MEDIAN) {
	fsort (Vbuf, ny);
	value = Vbuf[(int)(0.5*ny)];
      } else {
	value = 0;
	for (j = 0; j < ny; j++) {
	  value += Vbuf[j];
	}
	if (Mode == MEAN) { value /= ny; }
      }
      yvec[0].elements.Flt[i] = value;
    }
    free (Vbuf);
    break;
    
  case 'y':
  case 'Y':
    ResetVector (xvec, OPIHI_FLT, ny);
    ResetVector (yvec, OPIHI_FLT, ny);
    bzero (yvec[0].elements.Flt, ny*sizeof(opihi_flt));
    for (i = 0; i < ny; i++) {
      xvec[0].elements.Flt[i] = i + sy; 
    }
    ALLOCATE (Vbuf, float, MAX (nx, 1));

    for (i = 0; i < ny; i++) {
      /* accumulate values */
      Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx; 
      for (j = 0; j < nx; j++, Vin ++) {
	Vbuf[j] = *Vin;
      }
      /* apply stat of choice */
      if (Mode == MEDIAN) {
	fsort (Vbuf, nx);
	value = Vbuf[(int)(0.5*nx)];
      } else {
	value = 0;
	for (j = 0; j < nx; j++) {
	  value += Vbuf[j];
	}
	if (Mode == MEAN) { value /= nx; }
      }
      yvec[0].elements.Flt[i] = value;
    }
    free (Vbuf);
    break;

  default:
    gprint (GP_ERR, "<dir> can be X or Y\n");
    return (FALSE);
    break;
  }

  return (TRUE);
}

