# include "relastro.h"

/* this function returns a list of all images which overlap the given SkyList (set of
   SkyRegions).  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 */

typedef struct {
  double Xc[5];
  double Yc[5];
  double Rc;
  double Dc;
} SkyRegionCoords;

Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) {
  
  Image *image;
  int i, j, k, m, found, nimage, NIMAGE;
  int InRange;
  double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
  int *line_number;
  Coords tcoords;
  SkyRegionCoords *skycoords;
  
  if (skylist[0].Nregions < 1) {
    *Nimage = 0;
    *LineNumber = NULL;
    if (VERBOSE) fprintf (stderr, "no matching sky regions\n");
    return NULL;
  }

  // the comparison is made in the catalog local projection. below we set crval1,2
  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");

  ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);

  /* compare with each region file */
  for (i = 0; i < skylist[0].Nregions; i++) { 

    /* we make positional comparisons in the projection of catalog */
    skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin);
    skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin);
    tcoords.crval1 = skycoords[i].Rc;
    tcoords.crval2 = skycoords[i].Dc;

    /* define catalog corners */
    RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords);
    RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords);
    RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords);
    RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords);
    skycoords[i].Xc[4] = skycoords[i].Xc[0];    
    skycoords[i].Yc[4] = skycoords[i].Yc[0];    

    dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
    dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
    skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy;
    skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy;
    skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy;
    skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
  }

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

  nimage = 0;
  NIMAGE = 100;
  ALLOCATE (image, Image, NIMAGE);
  ALLOCATE (line_number, int, NIMAGE);
  
  // go through the complete list of images, selecting ones which overlap any region
  for (i = 0; i < Ntimage; i++) {
      
    /* exclude images by photcode, if selected */
    if (Nphotcodes > 0) {
      found = FALSE;
      for (j = 0; (j < Nphotcodes) && !found; j++) {
	if (photcodes[j][0].code == timage[i].photcode) found = TRUE;
      }
      if (!found) continue;
    }

    /* exclude images by time */
    if (TimeSelect) {
      if (timage[i].tzero < TSTART) continue;
      if (timage[i].tzero > TSTOP) continue;
    }
    
    if (!FindMosaicForImage (timage, Ntimage, i)) continue;

    /* 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 ra,dec */
    for (j = 0; j < 5; j++) {
      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
    }

    /* compare with each region file */
    for (m = 0; (m < skylist[0].Nregions) && !found; m++) { 

      /* we make positional comparisons in the projection of catalog */
      tcoords.crval1 = skycoords[m].Rc;
      tcoords.crval2 = skycoords[m].Dc;

      /* transform to ra,dec */
      InRange = TRUE;
      for (j = 0; (j < 5) && InRange; j++) {
	InRange = RD_to_XY (&Xi[j], &Yi[j], Ri[j], Di[j], &tcoords);
      }
      if (!InRange) continue;

      /* check if image corner inside catalog */
      for (j = 0; (j < 4) && !found; j++) {
	found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]);
      }
      /* check if catalog corner inside image */
      for (j = 0; (j < 4) && !found; j++) {
	found = corner_check (&skycoords[m].Xc[j], &skycoords[m].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], &skycoords[m].Xc[k], &skycoords[m].Yc[k]);
	}
      }
      if (!found) continue;

      image[nimage] = timage[i]; 
      /* always allow 'few' images to succeed, if possible */
      if (image[nimage].code & ID_IMAGE_FEW) { 
	image[nimage].code &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR);
      }
      if (RESET) {
	// XXX do we need this in relastro?
	assignMcal (&image[nimage], (double *) NULL, -1);
	image[nimage].dMcal = NAN;
	image[nimage].code &= ~ID_IMAGE_POOR;
      }
      line_number[nimage] = i;
      nimage ++;
      if (nimage == NIMAGE) {
	NIMAGE += 100;
	REALLOCATE (image, Image, NIMAGE);
	REALLOCATE (line_number, int, NIMAGE);
      }
    }
  }
      
  if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);

  REALLOCATE (image, Image, MAX (nimage, 1));
  REALLOCATE (line_number, int, MAX (nimage, 1));
  free (skycoords);

  *Nimage  = nimage;
  *LineNumber = line_number;
  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);

}

