# include "addstar.h"
void SetImageCorners (double *X, double *Y, Image *image);

/* given image, find catalog images which overlap it */
Image *gimages (FITS_DB *db, Image *image, Coords *mosaic, int *Npimage) {
  
  int i, j, k, addtolist;
  int NPIMAGE, npimage, Ntimage;
  Image *timage, *pimage;
  Coords tcoords;
  double r, d, x, y;
  double Xb[5], Yb[5], Xi[5], Yi[5], Xo[5], Yo[5];  /* image and original corners */
  double zeropt;

  /* if image.db is opened successfully, f will not be NULL.
     if image.db is empty, dbstate will be LCK_EMPTY 
     if image.db is not locked & opened successfully, we will not reach here!
  */
  if (db[0].dbstate == LCK_EMPTY) {
    *Npimage = 0;
    return (NULL);
  }

  /* read entire db table (use block read eventually?) */
  if (!dvo_image_load (db, VERBOSE, FORCE_READ)) {
    Shutdown ("can't read image catalog %s", db[0].filename);
  }
  timage = gfits_table_get_Image (&db[0].ftable, &Ntimage, &db[0].swapped);

  gfits_scan (&db[0].header, "ZERO_PT", "%lf", 1, &zeropt);
  if (fabs (ZeroPt - zeropt) > 1e-4) {
    Shutdown ("zero point in image (%f:%s) inconsistent with zero point in catalog (%f)",
	     ZeroPt, image[0].name, zeropt);
  }

  /* project onto rectilinear grid with 1 arcsec pixels */
  /* we keep the original crpix1,2 and crref1,2 */
  /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */
  if (!strcmp (&image[0].coords.ctype[4], "-WRP")) {
    tcoords = mosaic[0];
    tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
    tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
    tcoords.Npolyterms = 1;
    strcpy (tcoords.ctype, "RA---TAN");
    RegisterMosaic (mosaic);
    /* register so image->sky conversions below have correct mosaic */
  } else {
    tcoords = image[0].coords;
    tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
    tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
    tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
    strcpy (tcoords.ctype, "RA---TAN");
  }

  /* define original corners */
  SetImageCorners (Xo, Yo, &image[0]);
  Xo[4] = Xo[0]; Yo[4] = Yo[0];  /* so we can make a loop easily */

  for (j = 0; j < 5; j++) {
    /* XY-to_RD is two-level if ctype == WRP */
    XY_to_RD (&r, &d, Xo[j], Yo[j], &image[0].coords);
    RD_to_XY (&Xo[j], &Yo[j], r, d, &tcoords);
  }
  
  /* pimage represents the subset of overlapping images */
  npimage = 0;
  NPIMAGE = 20;
  ALLOCATE (pimage, Image, NPIMAGE);

  /* setup links for mosaic WRP and DIS entries */
  BuildChipMatch (timage, Ntimage);

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

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

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

    addtolist = FALSE;

    /* define image corners */
    SetImageCorners (Xi, Yi, &timage[i]);
    Xi[4] = Xi[0]; Yi[4] = Yi[0];
    /* transform to tcoords */
    for (j = 0; j < 5; j++) {
      XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords);
      RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords);
    }
    /* check if any edges cross */
    for (j = 0; (j < 4) && !addtolist; j++) {
      for (k = 0; (k < 4) && !addtolist; k++) {
	addtolist |= edge_check (&Xi[j], &Yi[j], &Xo[k], &Yo[k]);
      }
    }
    
    /* note that 0,0 is a valid pixel for both DIS and !DIS */
    /* if no edges cross, check if 0,0 pixel of timage[i] is inside image[0] */
    SetImageCorners (Xb, Yb, &image[0]);
    if (!addtolist) {
      XY_to_RD (&r, &d, 0.0, 0.0, &timage[i].coords);
      RD_to_XY (&x, &y, r, d, &image[0].coords);
      addtolist = TRUE;
      addtolist &= (x < Xb[0]);
      addtolist &= (y < Yb[0]);
      addtolist &= (x > Xb[2]);
      addtolist &= (y > Yb[2]);
    }

    /* if no edges cross, check if one corner of image[0] is inside timage[i] */
    SetImageCorners (Xb, Yb, &timage[i]);
    if (!addtolist) {
      XY_to_RD (&r, &d, 0.0, 0.0, &image[0].coords);
      RD_to_XY (&x, &y, r, d, &timage[i].coords);
      addtolist = TRUE;
      addtolist &= (x < Xb[0]);
      addtolist &= (y < Yb[0]);
      addtolist &= (x > Xb[2]);
      addtolist &= (y > Yb[2]);
    }

    if (addtolist) {
      pimage[npimage] = timage[i];
      npimage ++;
      if (npimage == NPIMAGE) {
	NPIMAGE += 20;
	REALLOCATE (pimage, Image, NPIMAGE);
      }
    }
  }

  BuildChipMatch (pimage, npimage);
  
  if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", npimage);

  *Npimage = npimage;
  return (pimage);
}
  
/* this function tries to find overlapping images with the current image.
   it loads the image table in blocks (up to 1000 at a time) and looks
   for overlaps within the set.  the search is not very robust. in 
   addition, it needs to look for mosaic DIS images matching any WRP images.
   if the WRP image is loaded near the beginning of a block, then the corresponding DIS 
   may be in the previous block, and thus missed 

   we are not adding the matching DIS image to the table.

   we are now loading in the entire image db table into memory.  

*/

# if (0)
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;
  }
}
# endif
