# include "dvoImagesAtCoords.h"
# ifndef FLT_MAX
# define FLT_MAX 1e32
# endif

/* given coordinate, find images in list that contain the point */
int *MatchCoords (Image *dbImages, int NdbImages, double ra, double dec, int *Nmatch) {
  
  int i, j, N, addtolist, status;
  int NMATCH, nmatch, *match;
  Coords tcoords;
  double r, d;
  double Xi[4], Yi[4], Xo[4], Yo[4];  /* image and original corners */
  double Xmin, Xmax, Ymin, Ymax;
  double xmin, xmax, ymin, ymax;

  *Nmatch = 0;

  /* match represents the subset of overlapping images */
  nmatch = 0;
  NMATCH = 20;
  ALLOCATE (match, int, NMATCH);

  /* setup links for mosaic WRP and DIS entries */
  BuildChipMatch (dbImages, NdbImages);

  /* run through image table and search for overlaps
     also define the vtable entries for the images we keep  */

  for (i = 0; i < NdbImages; i++) {

    if (!WITH_PHU && !strcmp (&dbImages[i].coords.ctype[4], "-DIS")) continue;
    if ( SOLO_PHU &&  strcmp (&dbImages[i].coords.ctype[4], "-DIS")) continue;

    /* if any of these images are WRP images, need to find matching DIS */
    if (!FindMosaicForImage (dbImages, NdbImages, i)) continue;

    /* define image corners */
    SetImageCorners (Xi, Yi, &dbImages[i]);
    // Xi[4] = Xi[0]; Yi[4] = Yi[0];


    ymin = xmin = +FLT_MAX;
    ymax = xmax = -FLT_MAX;
    for (j = 0; j < 4; j++) {
        xmin = MIN(xmin, Xi[j]);
        xmax = MAX(xmax, Xi[j]);
        ymin = MIN(ymin, Yi[j]);
        ymax = MAX(ymax, Yi[j]);
    }

    // transform input point to image coords
    double x, y;
    RD_to_XY(&x, &y, ra, dec, &dbImages[i].coords);

    if ((x >= xmin && x <= xmax) && (y >= ymin && y <= ymax)) {

        match[nmatch] = i;
        nmatch ++;
        if (nmatch == NMATCH) {
          NMATCH += 20;
          REALLOCATE (match, int, NMATCH);
        }
    }
      
  }
  if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", nmatch);
  
  *Nmatch = nmatch;
  return (match);
}
  
void SetImageCorners (double *X, double *Y, Image *image) {

  if (!strcmp(&image[0].coords.ctype[4], "-DIS")) {
    X[0] = -0.5*image[0].NX; Y[0] = -0.5*image[0].NY;
    X[1] = +0.5*image[0].NX; Y[1] = -0.5*image[0].NY;
    X[2] = +0.5*image[0].NX; Y[2] = +0.5*image[0].NY;
    X[3] = -0.5*image[0].NX; Y[3] = +0.5*image[0].NY;
  } else {
    X[0] = 0;           Y[0] = 0;
    X[1] = image[0].NX; Y[1] = 0;
    X[2] = image[0].NX; Y[2] = image[0].NY;
    X[3] = 0;           Y[3] = image[0].NY;
  }
}
