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

/* given coordinate, find images in list that contain the point */
off_t MatchCoords (Image *dbImages, off_t NdbImages, Point *points, int Npoints) {
  
  off_t i;
  int j;
  int totalMatches = 0;
  double Xi[4], Yi[4];  /* image and original corners */
  double xmin, xmax, ymin, ymax;


  /* 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]);
    }

    for (j = 0; j < Npoints; j++) {
        Point *pt = points + j;

        // transform input point to image coords
        double x, y;
        int status = RD_to_XY(&x, &y, pt->ra, pt->dec, &dbImages[i].coords);

        if (!status) {
            // avoid matching antipodal skycells
            continue;
        }

#ifdef notdef
        printf("%ld %8.1f %8.1f %s\n", i, x, y, dbImages[i].name);
#endif

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

            pt->matches[pt->Nmatches] = i;
            pt->Nmatches ++;

            if (pt->Nmatches == pt->arrayLength) {
                pt->arrayLength += 20;
                REALLOCATE (pt->matches, int, pt->arrayLength);
            }
        }
    }
 }
 if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", totalMatches);
  
 return (totalMatches);
}
  
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;
  }
}
