# include "relastro.h"

# define USE_IMAGE_ID 1

static off_t       **bin;     // link from catalog,measure to image
static int         **clist;   // catalog which supplied measurement on image
static off_t       **mlist;   // measure reference for measurement on image
static off_t        *Nlist;   // number of measurements on image
static off_t        *NLIST;   // allocated number of measurements on image   

static Image        *image;   // list of available images
static off_t        Nimage;   // number of available images

static int         *Ncatlist;  // catalogs associated with each image
static int         *NCATLIST;  // catalogs associated with each image
static int         **catlist;  // catalogs associated with each image

// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
// index

# if USE_IMAGE_ID
static off_t        *imageIDs; // list of all image IDs
static off_t        *imageIdx; // list of index for image IDs 
# else
static unsigned int *start;
static unsigned int *stop;
# endif

Image *getimages (off_t *N) {
  *N = Nimage;
  return (image);
}

Image *getimage (off_t N) {
  return (&image[N]);
}

int *getCatlist (int *N, off_t im) {

  *N = Ncatlist[im];
  return (catlist[im]);
}

void initImages (Image *input, off_t N) {

  off_t i;

  image = input;
  Nimage = N;

# if USE_IMAGE_ID
  ALLOCATE (imageIDs, off_t, Nimage);
  ALLOCATE (imageIdx, off_t, Nimage);

  for (i = 0; i < Nimage; i++) {
    imageIdx[i] = i;
    imageIDs[i] = image[i].imageID;
  }
  llsortpair (imageIDs, imageIdx, Nimage);
# else
  ALLOCATE (start, unsigned, Nimage);
  ALLOCATE (stop, unsigned, Nimage);

  for (i = 0; i < Nimage; i++) {
    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
  }
# endif
}

off_t getImageByID (off_t ID) {

  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
  // use bisection to find the specified image ID

# if USE_IMAGE_ID
  off_t Nlo, Nhi, N;

  Nlo = 0; Nhi = Nimage;
  while (Nhi - Nlo > 10) {
    N = 0.5*(Nlo + Nhi);
    if (imageIDs[N] < ID) {
      Nlo = MAX(N, 0);
    } else {
      Nhi = MIN(N + 1, Nimage);
    }
  }

  for (N = Nlo; N < Nhi; N++) {
    if (imageIDs[N] == ID)
      return (imageIdx[N]);
  }
# endif

  return (-1);
}

// these are really image & catalog indexes
void initImageBins (Catalog *catalog, int Ncatalog, int FULLINIT) {

  off_t i, j;

  ALLOCATE (bin, off_t *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
  }

  ALLOCATE (Nlist, off_t,   Nimage);
  ALLOCATE (NLIST, off_t,   Nimage);
  ALLOCATE (clist, int *,   Nimage);
  ALLOCATE (mlist, off_t *, Nimage);

  for (i = 0; i < Nimage; i++) {
    Nlist[i] =   0;
    NLIST[i] = 100;
    ALLOCATE (clist[i], int, NLIST[i]);
    ALLOCATE (mlist[i], off_t, NLIST[i]);
  }

  if (FULLINIT) {
    ALLOCATE (Ncatlist, int,  Nimage);
    ALLOCATE (NCATLIST, int,  Nimage);
    ALLOCATE (catlist, int *, Nimage);

    for (i = 0; i < Nimage; i++) {
      Ncatlist[i] =  0;
      NCATLIST[i] = 32;
      ALLOCATE (catlist[i], int, NCATLIST[i]);
    }
  }
}

void freeImageBins (int Ncatalog) {

  off_t i;

  for (i = 0; i < Ncatalog; i++) {
    free (bin[i]);
  }
  free (bin);
  for (i = 0; i < Nimage; i++) {
    free (clist[i]);
    free (mlist[i]);
  }
  free (clist);
  free (mlist);
}

/* match measurements to images */
void findImages (Catalog *catalog, int Ncatalog, int MATCHCAT) {

  off_t i, j;
  char *name;

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      // XXX : update to optionally restrict by photcode equiv a la relphot
      // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
      // if (photcode[0].code != ecode) continue;
      matchImage (catalog, j, i, MATCHCAT);
    }
  }

  for (i = 0; VERBOSE2 && (i < Nimage); i++) {
    name = GetPhotcodeNamebyCode (image[i].photcode);
    fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i],
             ohana_sec_to_date(image[i].tzero), name);
  }
}

# if USE_IMAGE_ID
// this is the imageID-based match
void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {

  off_t idx, ID;
  Measure *measure;
  int i, found;

  measure = &catalog[cat].measure[meas];

  ID = measure[0].imageID;
  idx = getImageByID (ID);
  if (idx == -1) {
    if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
    return;
  }

  // index for (catalog, measure) -> image
  bin[cat][meas] = idx;

  // index for image, Nentry -> catalog
  clist[idx][Nlist[idx]] = cat;

  // index for image, Nentry -> measure
  mlist[idx][Nlist[idx]] = meas;
  Nlist[idx] ++;

  if (Nlist[idx] == NLIST[idx]) {
    NLIST[idx] += 100;
    REALLOCATE (clist[idx], int,   NLIST[idx]);
    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
  }

  if (MATCHCAT) {
    // index for image -> catalog list
    found = FALSE;
    for (i = 0; !found && (i < Ncatlist[idx]); i++) {
      if (catlist[idx][i] == cat) found = TRUE;
    }
    if (!found) {
      catlist[idx][Ncatlist[idx]] = cat;
      Ncatlist[idx] ++;
      if (Ncatlist[idx] == NCATLIST[idx]) {
	NCATLIST[idx] += 32;
	REALLOCATE (catlist[idx], int, NCATLIST[idx]);
      }
    }
  }

  return;
}

# else
// this is the time-based match
void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {

  off_t i;
  Measure *measure;

  measure = &catalog[cat].measure[meas];

  /* find the image that supplied this measurement */
  for (i = 0; i < Nimage; i++) {
    // let's try the very slow method first before adding a bisection search
    // if (image[i].imageID != measure[0].imageID) continue;
    if (image[0].photcode == -1) continue;
    if (measure[0].photcode != image[i].photcode) continue;
    if (measure[0].t < start[i]) continue;
    if (measure[0].t > stop[i]) continue;

    // index for (catalog, measure) -> image
    bin[cat][meas] = i;

    // index for image, Nentry -> catalog
    clist[i][Nlist[i]] = cat;

    // index for image, Nentry -> measure
    mlist[i][Nlist[i]] = meas;
    Nlist[i] ++;

    if (Nlist[i] == NLIST[i]) {
      NLIST[i] += 100;
      REALLOCATE (clist[i], int,   NLIST[i]);
      REALLOCATE (mlist[i], off_t, NLIST[i]);
    }
    return;
  }
  if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
  return;
}
# endif

/*
Coords *getCoords (off_t meas, int cat) {

  off_t i;

  i = bin[cat][meas];
  if (i == -1) return (NULL);
  return (&image[i].coords);
}
*/

void plot_images () {

  off_t i, bin;
  double *xlist, *Mlist, *dlist;
  Graphdata graphdata;

  ALLOCATE (xlist, double, Nimage);
  ALLOCATE (dlist, double, Nimage);
  ALLOCATE (Mlist, double, Nimage);

  /**** dMcal vs airmass ****/
  for (i = 0; i < Nimage; i++) {
    Mlist[i] = image[i].Mcal;
    dlist[i] = image[i].dMcal;
    xlist[i] = image[i].secz;
  }

  plot_defaults (&graphdata);
  graphdata.ymin = PlotdMmin;
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png");
  plot_defaults (&graphdata);
  plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL);

# define NBIN 20
  REALLOCATE (xlist, double, NBIN);
  REALLOCATE (Mlist, double, NBIN);

  /**** dMcal histgram ****/
  for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Nimage; i++) {
    bin = image[i].dMcal / 0.00025;
    bin = MAX (0, MIN (NBIN - 1, bin));
    Mlist[bin] += 1.0;
  }

  plot_defaults (&graphdata);
  graphdata.style = 1;
  plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png");

  free (dlist);
  free (xlist);
  free (Mlist);
}

// return StarData values for detections in the specified image, converting coordinates from the
// chip positions: X,Y -> L,M -> P,Q -> R,D
void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {

  off_t i, m, c, n, nPos;
  double X, Y, L, M, P, Q, R, D, dR, dD;
  double dPos;

  Mosaic *mosaic;
  Coords *moscoords, *imcoords;

  // check if this image is bad and should be skipped
  if (badCoords(im)) return;

  // WRP images need to have an associated mosaic
  moscoords = NULL;
  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
    mosaic = getMosaicForImage (im);
    if (mosaic == NULL) return;  // if we cannot find the associated image, skip it
    moscoords = &mosaic[0].coords;
  }
  imcoords = &image[im].coords;

  // accumulate the rms position offsets.  if this value, or any specific entry, is too
  // large, we will reset the image to the original coords at the end of the analysis

  dPos = 0.0;
  nPos = 0;

  for (i = 0; i < Nlist[im]; i++) {
    m = mlist[im][i];
    c = clist[im][i];

    X = catalog[c].measure[m].Xccd;
    Y = catalog[c].measure[m].Yccd;
    n = catalog[c].measure[m].averef;

    if (moscoords == NULL) {
      // this is a Simple image (not a mosaic)
      // note that for a Simple image, L,M = P,Q
      XY_to_LM (&L, &M, X, Y, imcoords);
      LM_to_RD (&R, &D, L, M, imcoords);
    } else {
      XY_to_LM (&L, &M, X, Y, imcoords);
      XY_to_LM (&P, &Q, L, M, moscoords);
      LM_to_RD (&R, &D, P, Q, moscoords);
    }

    // new dR, dD : test
    dR = 3600.0*(catalog[c].average[n].R - R);
    dD = 3600.0*(catalog[c].average[n].D - D);

    // complain if the new location is far from the average location
    if (fabs(dR) > 2.0) {
      fprintf (stderr, "!");
      setBadCoords (im); // report a failure for this image
      return;
    }
    if (fabs(dD) > 2.0) {
      fprintf (stderr, "*");
      setBadCoords (im); // report a failure for this image
      return;
    }

    // complain if the new location is far from the old location
    if (fabs(catalog[c].measure[m].dR - dR) > 2.0) {
      fprintf (stderr, "@");
      setBadCoords (im); // report a failure for this image
      return;
    }
    if (fabs(catalog[c].measure[m].dD - dD) > 2.0) {
      fprintf (stderr, "#");
      setBadCoords (im); // report a failure for this image
      return;
    }

    dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD);
    nPos ++;

    catalog[c].measure[m].dR = dR;
    catalog[c].measure[m].dD = dD;

    if (catalog[c].measure[m].dR > +180.0*3600.0) {
      // average on high end of boundary, move star up
      R += 360.0;
      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    }
    if (catalog[c].measure[m].dR < -180.0*3600.0) {
      // average on low end of boundary, move star down
      R -= 360.0;
      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    }
  }

  saveOffsets (dPos, nPos, im);

  return;
}

// return StarData values for detections in the specified image, converting coordinates from the
// chip positions: X,Y -> L,M -> P,Q -> R,D
void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im) {

  off_t i, m, c, n;
  double X, Y, L, M, P, Q, R, D, dR, dD;

  Mosaic *mosaic;
  Coords *moscoords, *imcoords, *oldcoords;

  // check if this image is bad and should be skipped
  if (!badCoords(im)) {
    fprintf (stderr, "ERROR: inconsistent result?");
    abort();
  }

  // replace the current coords with the old coords:
  oldcoords = getCoords (im);
  memcpy (&image[im].coords, oldcoords, sizeof(Coords));

  // WRP images need to have an associated mosaic
  moscoords = NULL;
  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
    mosaic = getMosaicForImage (im);
    if (mosaic == NULL) return;  // if we cannot find the associated image, skip it
    moscoords = &mosaic[0].coords;
  }
  imcoords = &image[im].coords;

  for (i = 0; i < Nlist[im]; i++) {
    m = mlist[im][i];
    c = clist[im][i];

    X = catalog[c].measure[m].Xccd;
    Y = catalog[c].measure[m].Yccd;
    n = catalog[c].measure[m].averef;

    dR = dD = 0.0;
    if (moscoords == NULL) {
      // this is a Simple image (not a mosaic)
      // note that for a Simple image, L,M = P,Q
      XY_to_LM (&L, &M, X, Y, imcoords);
      LM_to_RD (&R, &D, L, M, imcoords);
    } else {
      XY_to_LM (&L, &M, X, Y, imcoords);
      XY_to_LM (&P, &Q, L, M, moscoords);
      LM_to_RD (&R, &D, P, Q, moscoords);
    }

    catalog[c].measure[m].dR = dR;
    catalog[c].measure[m].dD = dD;

    if (catalog[c].measure[m].dR > +180.0*3600.0) {
      // average on high end of boundary, move star up
      R += 360.0;
      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    }
    if (catalog[c].measure[m].dR < -180.0*3600.0) {
      // average on low end of boundary, move star down
      R -= 360.0;
      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
    }
  }

  return;
}

// return StarData values for detections in the specified image, converting coordinates from the
// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
// which the detections have already been filtered when they were loaded (bcatalog)
StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {

  off_t i, m, c, n;

  Mosaic *mosaic;
  Coords *moscoords;
  StarData *raw;

  mosaic = NULL;
  moscoords = NULL;
  if (mode == MODE_MOSAIC) {
    mosaic = getMosaicForImage (im);
    if (mosaic == NULL) {
      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
      return NULL;
    }
    moscoords = &mosaic[0].coords;
  }

  ALLOCATE (raw, StarData, Nlist[im]);

  for (i = 0; i < Nlist[im]; i++) {
    m = mlist[im][i];
    c = clist[im][i];

    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    raw[i].X = catalog[c].measure[m].Xccd;
    raw[i].Y = catalog[c].measure[m].Yccd;

    raw[i].Mag  = catalog[c].measure[m].M;
    raw[i].dMag = catalog[c].measure[m].dM;
    raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);

    n = catalog[c].measure[m].averef;

    // an object with only one detection provides no information about the image calibration
    //XXX this is already taken care of in bcatalog
    raw[i].mask = FALSE;
    int mask = FALSE;
    if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
      mask = TRUE;
    }
    if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
      mask = TRUE;
    }

    raw[i].mask = mask;


    switch (mode) {
      case MODE_SIMPLE:
        /* note that for a Simple image, L,M = P,Q */
        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
        raw[i].P = raw[i].L;
        raw[i].Q = raw[i].M;
        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
        break;
      case MODE_MOSAIC:
        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
        XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
        break;
      default:
	fprintf (stderr, "error: invalid mode in getImageRaw");
	abort ();
    }
  }

  *Nstars = Nlist[im];
  return (raw);
}

// return StarData values for averages positions in the specified image, converting coordinates from
// the sky positions: R,D -> P,Q -> L,M -> X,Y

StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {

  off_t i, m, c, n;

  Mosaic *mosaic;
  Coords *moscoords;
  StarData *ref;

  mosaic = NULL;
  moscoords = NULL;
  if (mode == MODE_MOSAIC) {
    mosaic = getMosaicForImage (im);
    if (mosaic == NULL) {
      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
      return NULL;
    }
    moscoords = &mosaic[0].coords;
  }

  ALLOCATE (ref, StarData, Nlist[im]);

  for (i = 0; i < Nlist[im]; i++) {
    m = mlist[im][i];
    c = clist[im][i];
    n = catalog[c].measure[m].averef;

    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    ref[i].R = catalog[c].average[n].R;
    ref[i].D = catalog[c].average[n].D;

    ref[i].Mag  = catalog[c].measure[m].M;
    ref[i].dMag = catalog[c].measure[m].dM;
    ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);

    ref[i].mask = FALSE;

    /* note that for a Simple image, L,M = P,Q */
    switch (mode) {
      case MODE_SIMPLE:
	RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
	ref[i].L = ref[i].P;
	ref[i].M = ref[i].Q;
	LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
	break;
      case MODE_MOSAIC:
	RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
	LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
	LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
	break;
      default:
        fprintf (stderr, "invalid case");
        abort();
    }
  }

  *Nstars = Nlist[im];
  return (ref);
}

/** lifted from relphot/StarOps.clean_measures */
void FlagOutliers2D(Catalog *catalog);

void FlagOutliers (Catalog *catalog) {
  // XXX FlagOutliers is just using FlagOutliers2D
  FlagOutliers2D(catalog);
  return;

  int Ndel, Nave;
  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
  double Ns, theta, x, y;
  double *R, *D, *dR, *dD;
  StatType statsR, statsD;

  Nsecfilt = GetPhotcodeNsecfilt();
  assert(catalog[0].Nsecfilt == Nsecfilt);

  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
  Nmax = 0;
  for (i = 0; i < catalog[0].Naverage; i++) {
    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
  }

  ALLOCATE (R, double, Nmax);
  ALLOCATE (D, double, Nmax);
  ALLOCATE (dR, double, Nmax);
  ALLOCATE (dD, double, Nmax);

  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);

  Ns = CLIP_THRESH;
  Ndel = Nave = 0;
      
  /* loop over each object in the catalog */
  for (j = 0; j < catalog[0].Naverage; j++) {
    
    /* accumulate list of valid measurements */
    m = catalog[0].average[j].measureOffset;
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
      
      // skip measurements based on user selected criteria
      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
      R[N] = catalog[0].measure[m].dR;
      D[N] = catalog[0].measure[m].dD;
      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
      if (isnan(R[N]) || isnan(D[N])) continue;
      N++;
    }
    if (N <= TOOFEW) continue;
    
    /* 3-sigma clip based on stats of inner 50% */
    initstats ("MEAN");
    liststats (R, dR, N, &statsR);
    liststats (D, dD, N, &statsD);
    
    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
    
    /* compare per-object distance to this standard deviation, and flag outliers*/
    m = catalog[0].average[j].measureOffset;
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
      //reset flag on each invocation
      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;

      // skip measurements based on user selected criteria
      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
      
      x = catalog[0].measure[m].dR - statsR.median;
      y = catalog[0].measure[m].dD - statsD.median;
      theta = atan2(y,x);
      if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + 
			 SQR(statsD.sigma * Ns * sin(theta)))) {   
	catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
	Ndel++;
      }
      N++;
      Nave ++;
    }

    //examine results
    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
			       catalog[0].average[j].Nmeasure, 
			       statsR, statsD, Ns);
  }
  
  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
  free (R);
  free(dR);
  free(D);
  free(dD); 
}


/** an alternative outlier rejection scheme */
void FlagOutliers2D (Catalog *catalog) {

  int Ndel, Nave;
  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
  double *index;
  double Ns, theta, x, y;
  double *R, *D, *dR, *dD, *d2;
  StatType statsR, statsD;

  Nsecfilt = GetPhotcodeNsecfilt();
  assert(catalog[0].Nsecfilt == Nsecfilt);

  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
  Nmax = 0;
  for (i = 0; i < catalog[0].Naverage; i++) {
    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
  }

  ALLOCATE (R, double, Nmax);
  ALLOCATE (D, double, Nmax);
  ALLOCATE (dR, double, Nmax);
  ALLOCATE (dD, double, Nmax);
  ALLOCATE (d2, double, Nmax);
  ALLOCATE (index, double, Nmax);

  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);

  Ns = CLIP_THRESH;
  Ndel = Nave = 0;
      
  /* loop over each object in the catalog */
  for (j = 0; j < catalog[0].Naverage; j++) {
    
    /* accumulate list of valid measurements */
    m = catalog[0].average[j].measureOffset;
    N = 0;
    
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
      
      //reset flag on each invocation
      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
      
      // skip measurements based on user selected criteria
      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
      R[N] = catalog[0].measure[m].dR;
      D[N] = catalog[0].measure[m].dD;
      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
      if (isnan(R[N]) || isnan(D[N])) continue;
      N++;
    }
    if (N <= TOOFEW) continue;
    
    /* calculate mean of all points*/
    initstats ("MEAN");
    liststats (R, dR, N, &statsR);
    liststats (D, dD, N, &statsD);
    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
    
    /* calculate deviations of all points*/
    m = catalog[0].average[j].measureOffset;
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
      //skip bad measurements
      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
      x = catalog[0].measure[m].dR - statsR.median;
      y = catalog[0].measure[m].dD - statsD.median;
      theta = atan2(y,x);
      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
			     SQR(statsD.sigma * Ns * sin(theta)));      
      index[N] = m;
      N++;
    }
    
    // sort d2
    dsortpair(d2, index, N);
    N = (N/2 > (N-1)) ? N/2 : N-1;

    // recalculate image center, sigma based on closest 50% of points
    for(k = 0;  k < N; k++) {
      off_t ind = (off_t) index[k];
      R[k] = catalog[0].measure[ind].dR;
      D[k] = catalog[0].measure[ind].dD;
      dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA);
      dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC);
    }
    liststats (R, dR, N, &statsR);
    liststats (D, dD, N, &statsD);
    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
    
    // use these new statistics to flag outliers 
    m = catalog[0].average[j].measureOffset;
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
      //skip bad measurements
      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
      x = catalog[0].measure[m].dR - statsR.median;
      y = catalog[0].measure[m].dD - statsD.median;
      theta = atan2(y,x);
      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
			     SQR(statsD.sigma * Ns * sin(theta)));      
      if ((d2[N]) > 1) {
	catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
	Ndel ++;
      }
      N++;
      Nave++;
    }  //done rejecting outliers

    //examine results
    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
			       catalog[0].average[j].Nmeasure, 
			       statsR, statsD, Ns);
    
  } //done looping over objects
  
  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
  free (R);
  free(dR);
  free(D);
  free(dD); 
  free(d2);
  free(index);
}



/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
int MeasFilterTest(Measure *measure) {
  int found, k;
  long mask;
  PhotCode *code;
  float mag;

  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
  if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
  
  /* select measurements by photcode, or equiv photcode, if specified */
  if (NphotcodesKeep > 0) {
    found = FALSE;
    for (k = 0; (k < NphotcodesKeep) && !found; k++) {
      if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
      if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    }
    if (!found) return FALSE;
  }
  
  if (NphotcodesSkip > 0) {
    found = FALSE;
    for (k = 0; (k < NphotcodesSkip) && !found; k++) {
      if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
      if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    }
    if (found) return FALSE;
  }  
  
  /* select measurements by time */
  if (TimeSelect) {
    if (measure[0].t < TSTART) return FALSE;
    if (measure[0].t > TSTOP) return FALSE;
  }
  
  /* select measurements by quality */
  if (PhotFlagSelect) {
    if (PhotFlagBad) {
      mask = PhotFlagBad;
    } else {
      code = GetPhotcodebyCode (measure[0].photcode);
      mask = code[0].astromBadMask;
    }
    if (mask & measure[0].photFlags) return FALSE;
  }

  /* select measurements by measurement error */
  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
  
  /* select measurements by mag limit */
  if (ImagSelect) {
    mag = PhotInst (measure);
    if (mag < ImagMin || mag > ImagMax) return FALSE;
  }
  
  return TRUE;
}
