# include "photdbc.h"

/* this function returns a list of all images which overlap the given
   set of region files.  All images in the image catalog are tested
   once, so there is no check that an image already has been included.
   LineNum stores the locations in the Image database of the list of
   images */

Image *find_images (GSCRegion *region, int Nregion, int *Nimages, int **LineNum) {
  
  Header header;
  Image *timage, *image;
  int i, j, k, m, found, nimage, Nimage, NIMAGE, NTIMAGE, Nloop, Nlast;
  int n, Nim, status, InRange, np, Np;
  FILE *f;
  double Xc[5], Yc[5], Xi[5], Yi[5], r, d, dx, dy;
  char line[256];
  int *line_number;
  Coords tcoords;

  if (VERBOSE) fprintf (stderr, "finding images\n");

  /* we make positional comparisons in the projection of catalog */
  tcoords.crval1 = 0.5*(region[0].RA[0]  + region[0].RA[1]);
  tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]);
  tcoords.crpix1 = tcoords.crpix2 = 0.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;
  strcpy (tcoords.ctype, "RA---TAN");

  /* check if image datafile exists, get header, number of images */
  if (!fits_read_header (ImageCat, &header)) {
    fprintf (stderr, "ERROR: No images in catalog %s (1)\n", ImageCat);
    exit (0);
  }
  Nimage = 0;
  fits_scan (&header, "NIMAGES", "%d", 1, &Nimage);
  if (Nimage == 0) {
    fprintf (stderr, "ERROR: No images in catalog %s (1)\n", ImageCat);
    exit (0);
  }

  /* get ready to read data on images */ 
  f = fopen (ImageCat, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "ERROR: trouble opening Image catalog: %s (2)\n", ImageCat);
    exit (0);
  }
  fseek (f, header.size, SEEK_SET); 

  /* load all images */
  ALLOCATE (timage, Image, Nimage);
  status = Fread (timage, sizeof(Image), Nimage, f, "image");
  if (status != Nimage) {
    fprintf (stderr, "ERROR: couldn't read images from image catalog: %s\n", ImageCat);
    exit (0);
  }

  nimage = 0;
  NIMAGE = 100;
  ALLOCATE (image, Image, NIMAGE);
  ALLOCATE (line_number, int, NIMAGE);
  
  for (i = 0; i < Nimage; i++) {
      
# if (0)
    /* select images by photcode */
    np = timage[i].source;
    Np = photcodes.hashcode[np];
    if (Np < 0) continue;
    if (photcodes.code[Np].equiv != PHOTCODE) continue;

    /* select images by time */
    if (TimeSelect) {
      if (timage[i].tzero < TSTART) continue;
      if (timage[i].tzero > TSTOP) continue;
    }
# endif
    
    /* define image corners */
    Xi[0] = 0;            Yi[0] = 0;
    Xi[1] = timage[i].NX; Yi[1] = 0;
    Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
    Xi[3] = 0;            Yi[3] = timage[i].NY;
    Xi[4] = 0;            Yi[4] = 0;
    found = FALSE;
    /* transform to tcoords */
    for (j = 0; j < 5; j++) {
      XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords);
      InRange = RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords);
      if (!InRange) goto imskip;
    }
    /* compare with each region file */
    for (m = 0; (m < Nregion) && !found; m++) { 
      /* define catalog corners */
      Xc[0] = region[m].RA[0]; Yc[0] = region[m].DEC[0];
      Xc[1] = region[m].RA[1]; Yc[1] = region[m].DEC[0];
      Xc[2] = region[m].RA[1]; Yc[2] = region[m].DEC[1];
      Xc[3] = region[m].RA[0]; Yc[3] = region[m].DEC[1];
      Xc[4] = region[m].RA[0]; Yc[4] = region[m].DEC[0];
      for (j = 0; j < 5; j++) {
	r = Xc[j]; d = Yc[j];
	RD_to_XY (&Xc[j], &Yc[j], r, d, &tcoords);
      }
      dx = 0.02*(Xc[2] - Xc[0]);
      dy = 0.02*(Yc[2] - Yc[0]);
      Xc[0] -= dx; Yc[0] -= dy;
      Xc[1] += dx; Yc[1] -= dy;
      Xc[2] += dx; Yc[2] += dy;
      Xc[3] -= dx; Yc[3] += dy;
      Xc[4] -= dx; Yc[4] -= dy;
      
      /* check if image corner inside catalog */
      for (j = 0; (j < 4) && !found; j++) {
	found |= corner_check (&Xi[j], &Yi[j], &Xc[0], &Yc[0]);
      }
      /* check if catalog corner inside image */
      for (j = 0; (j < 4) && !found; j++) {
	found |= corner_check (&Xc[j], &Yc[j], &Xi[0], &Yi[0]);
      }
      /* check if edges cross */
      for (j = 0; (j < 4) && !found; j++) {
	for (k = 0; (k < 4) && !found; k++) {
	  found |= edge_check (&Xi[j], &Yi[j], &Xc[k], &Yc[k]);
	}
      }
      if (found) {
	image[nimage] = timage[i]; 
	if (image[nimage].code == ID_IMAGE_NOCAL) {
	  image[nimage].code &= ~ID_IMAGE_NOCAL;
	}	    
	line_number[nimage] = i;
	nimage ++;
	if (nimage == NIMAGE) {
	  NIMAGE += 100;
	  REALLOCATE (image, Image, NIMAGE);
	  REALLOCATE (line_number, int, NIMAGE);
	}
      }
    }
  imskip:
  }
      
  if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);

  /* we are going to use image[nimage] to represent images not in the database:
     these are measurements from external sources, like USNO */

  assignMcal (&image[nimage], (double *) NULL, -1);
  image[nimage].code = ID_IMAGE_NEW;
  
  REALLOCATE (image, Image, MAX (nimage + 1, 1));
  REALLOCATE (line_number, int, MAX (nimage + 1, 1));
  free (timage);
  *Nimages = nimage;
  *LineNum = line_number;
  fclose (f);
  return (image);
}

/* check if line between points 0 and 1 of x1
   crosses line between points 0 and 1 of x2 */
int edge_check (double *x1, double *y1, double *x2, double *y2) {

  double theta1, theta2;
  double Theta1, Theta2;

  theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); 
  theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); 

  if (theta1*theta2 < 0.0) {
    return (FALSE);
  }

  if (fabs(theta1) < fabs(theta2)) {
    return (FALSE);
  }

  Theta1 = theta1;
  Theta2 = theta2;
  theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); 
  theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); 
  
 
  if (theta1*theta2 < 0.0) {
    return (FALSE);
  }

  if (fabs(theta1) < fabs(theta2)) {
    return (FALSE);
  }

  return (TRUE);

}

/* check if point x1,y1 is in box formed by x2[0-4] */
int corner_check (double *x1, double *y1, double *x2, double *y2) {

  int i;
  double theta;

  theta = 0;

  for (i = 0; i < 4; i++) {
    theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]); 
  }
  if (fabs(theta) > 6) {
    return (TRUE);
  } else {
    return (FALSE);
  }
}

/* returns the opening angle between the three points (2 is in middle) 
   in range -pi to pi */

double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) {

  double dx1, dy1, dx2, dy2, ct, st, theta;

  dx1 = x1 - x2;
  dy1 = y1 - y2;
  
  dx2 = x3 - x2;
  dy2 = y3 - y2;
  
  ct = (dx1*dx2 + dy1*dy2);
  st = (dx1*dy2 - dx2*dy1);

  theta = atan2 (st, ct);

  return (theta);

}

