# include "relastro.h"

# define USE_IMAGE_ID 1

// we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage)
static Image        *image;   // list of available images
static off_t        Nimage;   // number of available images

// if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row
// each image belongs
static off_t       *LineNumber; // match of subset to full image table

static off_t        *N_onImage;   // number of measurements on image
static off_t        *N_ONIMAGE;   // allocated number of measurements on image   

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

static IDX_T       **MeasureToImage;     // link from catalog,measure to image
static IDX_T       **ImageToCatalog;   // catalog which supplied measurement on image
static IDX_T       **ImageToMeasure;   // measure reference for measurement on 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

// MeasureToImage was 'bin'
// ImageToCatalog was 'clist'
// ImageToMeasure was 'mlist'

// N_onImage was 'Nlist'
// N_ONIMAGE was 'NLIST'

Image *getimages (off_t *N, off_t **line_number) {

  *N = Nimage;
  if (line_number) *line_number = LineNumber;
  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 *line_number, off_t N) {

  off_t i;

  image = input;
  LineNumber = line_number;
  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 (MeasureToImage, IDX_T *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (MeasureToImage[i], IDX_T, MAX (catalog[i].Nmeasure, 1));
    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToImage[i][j] = -1;
  }

  ALLOCATE (N_onImage, off_t,   Nimage);
  ALLOCATE (N_ONIMAGE, off_t,   Nimage);
  ALLOCATE (ImageToCatalog, IDX_T *,   Nimage);
  ALLOCATE (ImageToMeasure, IDX_T *, Nimage);

  for (i = 0; i < Nimage; i++) {
    N_onImage[i] =   0;
    N_ONIMAGE[i] =  30;
    ImageToCatalog[i] = NULL;  // we allocate these iff they are needed in matchImage
    ImageToMeasure[i] = NULL;  // we allocate these iff they are needed in matchImage
  }

  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 (MeasureToImage[i]);
  }
  free (MeasureToImage);
  for (i = 0; i < Nimage; i++) {
    if (ImageToCatalog[i]) { free (ImageToCatalog[i]); }
    if (ImageToMeasure[i]) { free (ImageToMeasure[i]); }
  }
  free (ImageToCatalog);
  free (ImageToMeasure);
  free (N_onImage);
  free (N_ONIMAGE);
}

/* 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++) {
      if (NphotcodesSkip > 0) {
	int k;
	int found = FALSE;
	for (k = 0; (k < NphotcodesSkip) && !found; k++) {
	  if (photcodesSkip[k][0].code == catalog[i].measureT[j].photcode) found = TRUE;
	  if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measureT[j].photcode)) found = TRUE;
	}
	if (found) continue;
      }  
      matchImage (catalog, j, i, MATCHCAT);
    }
  }

  // watch for chips with few stars
  int Nfew = 0;
  int Nbad = 0;
  for (i = 0; i < Nimage; i++) {
    // ignore the PHU entries
    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
    if (VERBOSE2 || (MATCHCAT && (N_onImage[i] < 20))) {
      name = GetPhotcodeNamebyCode (image[i].photcode);
      fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)",  i,  image[i].imageID, image[i].name, N_onImage[i], image[i].nstar,
	       ohana_sec_to_date(image[i].tzero), name);
      if (N_onImage[i] < 20) {
	fprintf (stderr, "*");
	Nfew ++;
      } 
      if (N_onImage[i] < 15) {
	fprintf (stderr, "*** warning ***");
	Nbad ++;
      } 
      fprintf (stderr, "\n");
    }
  }
  fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad);
}

# 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;
  MeasureTiny *measure;
  int i, found;

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

  ID = measure[0].imageID;

  if (catalog[cat].measure) {
    Measure *measureBig = &catalog[cat].measure[meas];
    int TESTPT = FALSE;
    TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measureBig->imageID == CAT_ID_SRC) && (measureBig->detID == OBJ_ID_SRC);
    TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measureBig->imageID == CAT_ID_DST) && (measureBig->detID == OBJ_ID_DST);
    if (TESTPT) {
      fprintf (stderr, "got test det\n");
    }
  }

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

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

  // if we need to allocate an image index table, do so here
  if (!ImageToCatalog[idx]) {
      ALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]);
  }
  if (!ImageToMeasure[idx]) {
      ALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]);
  }

  // index for image, Nentry -> catalog
  ImageToCatalog[idx][N_onImage[idx]] = cat;

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

  if (N_onImage[idx] == N_ONIMAGE[idx]) {
    N_ONIMAGE[idx] += 30;
    REALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]);
    REALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[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;
  MeasureTiny *measure;

  measure = &catalog[cat].measureT[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
    MeasureToImage[cat][meas] = i;

    // index for image, Nentry -> catalog
    ImageToCatalog[i][N_onImage[i]] = cat;

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

    if (N_onImage[i] == N_ONIMAGE[i]) {
      N_ONIMAGE[i] += 30;
      REALLOCATE (ImageToCatalog[i], IDX_T, N_ONIMAGE[i]);
      REALLOCATE (ImageToMeasure[i], IDX_T, N_ONIMAGE[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 = MeasureToImage[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);
}

void dump_measures (Average *average, Measure *measure) {

  off_t j, off;

  for (j = 0; j < average[0].Nmeasure; j++) {
    off = average[0].measureOffset + j;
    fprintf (stderr, "dR, dD, mag, dMag: %f, %f, %f, %f\n", measure[off].dR, measure[off].dD, measure[off].M, measure[off].dM);
  }
  return;
}

static int NcatTotal = 0;

// 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, dPosSys, DPOS_MAX_ASEC;

  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;

  if (moscoords) {
    DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(moscoords[0].cdelt1, moscoords[0].cdelt2);
  } else {
    DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(imcoords[0].cdelt1, imcoords[0].cdelt2);
  }

  // these are used to 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;

  // convert the image systematic error in pixels to a value in arcsec
  { 
    double dLsig, dMsig;
    double Ro, Do, Rx, Dx, dP0, dP1;
    Coords *coords;

    // these values are in pixels, but we to convert to arcsec
    dLsig = image[0].dXpixSys;
    dMsig = image[0].dYpixSys;

    if (moscoords == NULL) {
      coords = imcoords;
    } else {
      coords = moscoords;
    }
    XY_to_LM (&Ro, &Do, 0.0, 0.0, coords);
    XY_to_LM (&Rx, &Dx, dLsig, 0.0, coords);
    dP0 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec
    XY_to_LM (&Rx, &Dx, 0.0, dMsig, coords);
    dP1 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec
    dPosSys = 0.5 * (dP0 + dP1);
  }      

  int NoffRAave = 0;  int NoffRAori = 0; 
  int NoffDECave = 0; int NoffDECori = 0;

  for (i = 0; i < N_onImage[im]; i++) {
    m = ImageToMeasure[im][i];
    c = ImageToCatalog[im][i];

    Measure *measure = &catalog[c].measure[m];
    MeasureTiny *measureT = &catalog[c].measureT[m];

    int TESTPT = FALSE;

    TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC);
    TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST);
    if (TESTPT) {
      fprintf (stderr, "got test det\n");
    }

    if (USE_FIXED_PIXCOORDS) {
      X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
      Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
    } else {
      X = measure[0].Xccd;
      Y = measure[0].Yccd;
    }
    n = measure[0].averef;
    Average *average = &catalog[c].average[n];

    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*(average[0].R - R);
    dD = 3600.0*(average[0].D - D);

    // make sure detection is on the same side of the 0,360 boundary as the average
    // this will give some funny results withing ~1 arcsec of the pol
    if (dR > +180.0*3600.0) {
      // average on high end of boundary, move star up
      R += 360.0;
      dR = 3600.0*(average[0].R - R);
    }
    if (dR < -180.0*3600.0) {
      // average on low end of boundary, move star down
      R -= 360.0;
      dR = 3600.0*(average[0].R - R);
    }

    float csdec = cos(average[0].D * RAD_DEG);

    // complain if the new location is far from the average location
    // NOTE: This should never happen, or our StarMap tests are not working
    if (fabs(dR*csdec) > 3.0*ADDSTAR_RADIUS) {
      NoffRAave ++;
      if (VERBOSE2) {
	fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, dR*csdec, dD);
	dump_measures (&average[0], catalog[c].measure);
      }
      // abort ();
    }
    if (fabs(dD) > 3.0*ADDSTAR_RADIUS) {
      NoffDECave ++;
      if (VERBOSE2) {
	fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", average[0].R, average[0].D, dR, dD);
	dump_measures (&average[0], catalog[c].measure);
      }
      // abort ();
    }

    // complain if the new location is far from the old location
    if (fabs(csdec*(measure[0].dR - dR)) > DPOS_MAX_ASEC) {
      NoffRAori ++;
      if (VERBOSE2) {
	fprintf (stderr, "measurement moves far from original location (R): %f %f (%f %f %f %f)\n", average[0].R, average[0].D, measure[0].dR, dR, csdec*(measure[0].dR - dR), dD);
	dump_measures (&average[0], catalog[c].measure);
      }
      // abort();
    }
    if (fabs(measure[0].dD - dD) > DPOS_MAX_ASEC) {
      NoffDECori ++;
      if (VERBOSE2) {
	fprintf (stderr, "measurement moves far from original location (D): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, measure[0].dD, dD);
	dump_measures (&average[0], catalog[c].measure);
      }      // abort();
    }

    dPos += SQ(measure[0].dR - dR) + SQ(measure[0].dD - dD);
    nPos ++;

    measure[0].dR = dR;
    measure[0].dD = dD;
    measureT[0].dR = dR;
    measureT[0].dD = dD;
    
    // set the systematic error for this image:
    measure[0].dRsys = ToShortPixels(dPosSys);
    measureT[0].dRsys = ToShortPixels(dPosSys);
  }

  NcatTotal += nPos;

  int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori;
  if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);
  saveOffsets (dPos, nPos, im);

  return;
}

void printNcatTotal () {
  fprintf (stderr, "NcatTotal: %d\n", NcatTotal);
}

// 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;

  fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n");
  abort();

  // 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 < N_onImage[im]; i++) {
    m = ImageToMeasure[im][i];
    c = ImageToCatalog[im][i];

    // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated
    Measure *measure = &catalog[c].measure[m];

    if (USE_FIXED_PIXCOORDS) {
      X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
      Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
    } else {
      X = measure[0].Xccd;
      Y = measure[0].Yccd;
    }
    n = measure[0].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);
    }

    measure[0].dR = dR;
    measure[0].dD = dD;

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

  for (i = 0; i < N_onImage[im]; i++) {
    m = ImageToMeasure[im][i];
    c = ImageToCatalog[im][i];

    MeasureTiny *measure = &catalog[c].measureT[m];

    /* apply the current image transformation or use the current value of R+dR, D+dD? */
    if (USE_FIXED_PIXCOORDS) {
      raw[i].X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
      raw[i].Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
    } else {
      raw[i].X = measure[0].Xccd;
      raw[i].Y = measure[0].Yccd;
    }

    raw[i].Mag  = measure[0].M;
    raw[i].dMag = measure[0].dM;
    raw[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS);

    n = measure[0].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 = MARK_MEAS_DEFAULT;
    if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
      raw[i].mask |= MARK_TOO_FEW_MEAS;
    }
    if (!finite(measure[0].dR) || !finite(measure[0].dD)) {
      raw[i].mask |= MARK_NAN_POS_ERROR;
    }
    raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is

    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 = N_onImage[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, N_onImage[im]);

  for (i = 0; i < N_onImage[im]; i++) {
    m = ImageToMeasure[im][i];
    c = ImageToCatalog[im][i];

    MeasureTiny *measure = &catalog[c].measureT[m];
    n = measure[0].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  = measure[0].M;
    ref[i].dMag = measure[0].dM;
    ref[i].dPos = GetAstromErrorTiny (&measure[0], 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 = N_onImage[im];
  return (ref);
}

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

// operates on Full values (not tiny)
void FlagOutliers (Catalog *catalog) {

  // XXX FlagOutliers is now 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 (VERBOSE2) 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++) {
    
    // pointer to this set of measurements
    m = catalog[0].average[j].measureOffset;
    Measure *measure = &catalog[0].measure[m];

    /* accumulate list of valid measurements */
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
      // skip measurements based on user selected criteria
      if (!MeasFilterTest(&measure[k], FALSE)) continue;
      R[N] = measure[k].dR;
      D[N] = measure[k].dD;
      dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA);
      dD[N] = GetAstromError (&measure[k], 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*/
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
      //reset flag on each invocation
      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;

      // skip measurements based on user selected criteria
      if (!MeasFilterTest(&measure[k], FALSE)) continue;
      
      x = measure[k].dR - statsR.median;
      y = measure[k].dD - statsD.median;
      theta = atan2(y,x);
      if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + 
			 SQR(statsD.sigma * Ns * sin(theta)))) {   
	measure[k].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 (VERBOSE2) 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++) {
    
    // pointer to this set of measurements
    m = catalog[0].average[j].measureOffset;
    Measure *measure = &catalog[0].measure[m];

    /* accumulate list of valid measurements */
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {

      // reset flag on each invocation
      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
      
      // skip measurements based on user selected criteria
      if (!MeasFilterTest(&measure[k], FALSE)) continue;
      R[N] = measure[k].dR;
      D[N] = measure[k].dD;
      dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA);
      dD[N] = GetAstromError(&measure[k], 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*/
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
      // skip bad measurements
      if (!MeasFilterTest(&measure[k], FALSE)) continue;  
      x = measure[k].dR - statsR.median;
      y = measure[k].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] = k;
      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] = measure[ind].dR;
      D[k] = measure[ind].dD;
      dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA);
      dD[k] = GetAstromError(&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 
    N = 0;
    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
      //skip bad measurements
      if (!MeasFilterTest(&measure[k], FALSE)) continue;  
      x = measure[k].dR - statsR.median;
      y = measure[k].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) {
	measure[k].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 */ 
// we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) {
  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?
  if (!finite(measure[0].dM)) 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;
  }  
  
  // if (MinBadQF > 0.0) {
  //   if (measure[0].psfQF < MinBadQF) 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 (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
    return FALSE;
  }
  
  /* select measurements by mag limit */
  if (ImagSelect) {
    mag = PhotInstTiny (measure);
    if (mag < ImagMin || mag > ImagMax) return FALSE;
  }
  
  return TRUE;
}

/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
// we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
int MeasFilterTest(Measure *measure, int applySigmaLim) {
  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?
  if (!finite(measure[0].dM)) 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;
  }  
  
  if (MinBadQF > 0.0) {
    if (measure[0].psfQF < MinBadQF) 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 (applySigmaLim && (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;
}
