# include "relastro.h"
# include "relastroVisual.h"

int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {

  int i, NstatFull, Nstat, Niter;
  float dLsig, dMsig, dRsig;
  float dLsigFull, dMsigFull, dRsigFull;
  float dL, dM, dR, dRmax;

  CoordFit *fit = NULL;

  dRmax = 0.0;

  for (Niter = 0; Niter < IMFIT_CLIP_NITER; Niter ++) {
    
    // measure the scatter for the unmarked (good) measurements with dM < limit
    // SIGMA_LIM here is redundant with ImageOps.c:915
    GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, SIGMA_LIM);
    dRmax = IMFIT_CLIP_NSIGMA*dRsig;

    // (a) skip unfittable points
    // (b) mark good and bad points for fit (in and outliers)

    // measure the scatter distribution (use only the bright end detections)
    for (i = 0; i < Nmatch; i++) {
      if (raw[i].mask) continue;
      if (isnan(raw[i].dMag)) { 
	raw[i].mask |= MARK_NAN_MAG_ERROR; 
	continue; 
      }
      if (raw[i].dMag > SIGMA_LIM) { 
	raw[i].mask |= MARK_BIG_MAG_ERROR; 
	continue; 
      } // is this redundant with ImageOps.c:915?

      dL = raw[i].L - ref[i].L;
      dM = raw[i].M - ref[i].M;
      dR = hypot (dL, dM);
      if (dR > dRmax) {
	raw[i].mask |= MARK_BIG_OFFSET;
	continue;
      }
    }

    // figures to assess the fitting process:
    // x vs dx, x vs dy, y vs dx, y vs dy : 
    // residual vector field 
    // dx vs mag, dy vs mag
    relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);

    // fit the requested order polynomial
    if (CHIPORDER > 0) {
      int Ndof_min = 3;
      int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3);
      int order_use = MIN (CHIPORDER, order_max);
      if (order_use < 1) {
	if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch);
	image[0].flags |= ID_IMAGE_ASTROM_FEW;
	return FALSE;
      }
      image[0].coords.Npolyterms = order_use;
    }

    if (fit) fit_free (fit);
    fit = fit_init (image[0].coords.Npolyterms);

    // generate the fit matches
    for (i = 0; i < Nmatch; i++) {
      if (raw[i].mask) continue;
      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
    }

# if (0)
    // check if the fit has enough data points for the polynomial order
    skip = FALSE;
    switch (image[0].coords.Npolyterms) {
      case 0:
      case 1:
        skip = (fit[0].Npts < 8);
        break;
      case 2:
        skip = (fit[0].Npts < 11);
        break;
      case 3:
        skip = (fit[0].Npts < 15);
        break;
      default:
        fprintf (stderr, "invalid chip order %d\n", image[0].coords.Npolyterms);
	skip = TRUE;
    }
    if (skip) {
      if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms);
      fit_free (fit);
      image[0].flags |= ID_IMAGE_ASTROM_FEW;
      return FALSE;
    }
# endif

    // measure the fit, update the coords & object coordinates
    if (!fit_eval (fit)) {
      fprintf (stderr, "failed to fit new model\n");
      image[0].flags |= ID_IMAGE_ASTROM_FAIL;
      if (fit) fit_free (fit);
      return FALSE;
    }

    if (!fit_apply_coords (fit, &image[0].coords, FALSE)) {
      fprintf (stderr, "failed to fit new model\n");
      image[0].flags |= ID_IMAGE_ASTROM_FAIL;
      if (fit) fit_free (fit);
      return FALSE;
    }

    for (i = 0; i < Nmatch; i++) {
      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords);
    }

  }

  // count mask classes
  int nMask1 = 0;
  int nMask2 = 0;
  int nMask3 = 0;
  int nMask4 = 0;
  int nMask5 = 0;
  for (i = 0; i < Nmatch; i++) {
    if (!raw[i].mask) continue;
    if (raw[i].mask & MARK_TOO_FEW_MEAS ) nMask1 ++;
    if (raw[i].mask & MARK_NAN_POS_ERROR) nMask2 ++;
    if (raw[i].mask & MARK_NAN_MAG_ERROR) nMask3 ++;
    if (raw[i].mask & MARK_BIG_MAG_ERROR) nMask4 ++;
    if (raw[i].mask & MARK_BIG_OFFSET   ) nMask5 ++;
  }

  GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM);
  GetScatterRawRef(&dLsig,     &dMsig,     &dRsig,     &Nstat,     raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);

  int Nm = 0;
  int Ns = 0;
  for (i = 0; i < Nmatch; i++) {
    if (raw[i].mask) continue;
    Nm += raw[i].Nmeas;
    Ns++;
  }
  image[0].nLinkAstrom = (Nm / Ns);

  if (VERBOSE2) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5);

  image[0].dXpixSys = dLsig;
  image[0].dYpixSys = dMsig;
  image[0].nFitAstrom = fit[0].Npts;

  if (fit) fit_free (fit);

  return TRUE;
}

// measure the scatter distribution (limit selected objects by dMag)
int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit) {

  int i, Ns;
  float dL, dM, dLsum, dLsum2, dMsum, dMsum2, *dR;

  ALLOCATE (dR, float, Nstars);

  Ns = 0;
  dLsum = dLsum2 = dMsum = dMsum2 = 0.0;
  for (i = 0; i < Nstars; i++) {
    if (raw[i].mask) continue;
    if (isnan(raw[i].dMag)) continue;
    if ((SigmaLimit > 0.0) && (raw[i].dMag > SigmaLimit)) continue;
    
    dL = raw[i].L - ref[i].L;
    dM = raw[i].M - ref[i].M;

    dLsum  += dL;
    dLsum2 += dL*dL;
    dMsum  += dM;
    dMsum2 += dM*dM;

    dR[Ns] = hypot (dL, dM);
    Ns++;
  }

  *dLsig = sqrt((dLsum2 / Ns) - SQ(dLsum/Ns));
  *dMsig = sqrt((dMsum2 / Ns) - SQ(dMsum/Ns));
  *nKeep = Ns;

  if (Ns < 5) {
    *dRsig = NAN;
    free  (dR);
    return (FALSE);
  }

  // for a 2D Gaussian, 40% of the points are within R = 1 sigma
  fsort (dR, Ns);
  *dRsig = dR[(int)(0.40*Ns)];
  
  free  (dR);
  return (TRUE);
}

/* in the mosaic case, we have four coord systems of interest:
   R,D : the sky
   P,Q : the tangent plane
   L,M : the focal plane
   X,Y : the chip

   R,D -> P,Q (projection)
   P,Q -> L,M (polynomial transformation : DIS)
   L,M -> X,Y (polynomial transformation : WRP)

Some details about this function:

- the initial value of the clipping radius is set based on the distribution of the dR
  values for the bright sources.

- NITER clipping passes are performed

- the per-star astrometric errors are included in the fit.  The photcode table controls
  whether only the reported positional errors are used, or if the magnitudes errors are
  scaled to determine the error, and if there is a systematic floor for all sources.

- input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog

- outlying detections are clipped in the iterative passes, but the clipped detections are
  not recorded

*/

/* example using fit_apply() :

  f = fopen ("test3.dat", "w");

  // apply new coords to raw (X,Y -> L,M)
  for (i = 0; i < Nmatch; i++) {
    fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
    fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
    fprintf (f, "%f %f\n", L1, M1);
  }
  fclose (f);
*/

