# include "getstar.h"

int clean_stars (Stars *stars, int *nstars, Image *image) {

  int i, j, Nsubset, *subset, *skip, *found, Nstars, check, Nskip;
  double Xo, Yo, Xt, Yt, radius;
  double Mo, Mt, dM;

  Nstars = *nstars;
  ALLOCATE (found,  int, Nstars);
  ALLOCATE (subset, int, Nstars);
  ALLOCATE (skip,   int, Nstars);
  bzero (skip, Nstars*sizeof(int));

  Nskip = 0;
  Nsubset = 0;
  for (i = 0; i < Nstars; i++) {
    RD_to_XY (&stars[i].X, &stars[i].Y, stars[i].R, stars[i].D, &image[0].coords);
    check = (stars[i].found & BLEND_IMAGE) || (stars[i].found & BLEND_IMAGE_NEIGHBOR) || (stars[i].found & BLEND_CATALOG);
    if (check) {
      subset[Nsubset] = i;
      found[Nsubset] = FALSE;
      Nsubset ++;
    }
  }

  if (Nsubset > 1) {
    for (i = 0; i < Nsubset; i++) {
      if (found[i]) continue;
      Xo = stars[subset[i]].X;
      Yo = stars[subset[i]].Y;
      Mo = stars[subset[i]].M;
      for (j = 0; j < Nsubset; j++) {
	if (j == i) continue;
	if (found[j]) continue;
	Xt = stars[subset[j]].X;
	Yt = stars[subset[j]].Y;
	Mt = stars[subset[j]].M;
	radius = hypot ((Xo-Xt), (Yo-Yt));
	dM = fabs (Mo-Mt);
	if ((radius < 5.0) && (dM < 0.005)) {
	  skip[subset[j]] = TRUE;
	  found[j] = TRUE;
	  Nskip ++;
	}
      }
      found[i] = TRUE;
    }
 
    for (i = j = 0; i < Nstars; i++) {
      if (skip[i]) continue;
      if (i == j) { 
	j++;
	continue;
      }
      stars[j] = stars[i];
      j++;
    }
    *nstars = j;
  }
  fprintf (stderr, "%d skipped\n", Nskip);
}

