# include "psastroInternal.h"

// look for and exclude objects in clumps in the input list
psArray *psastroRemoveClumps (psArray *input, int scale) {

    // determine data range
    pmAstromObj *obj = input->data[0];
    float Xmin = obj->FP->x;
    float Xmax = obj->FP->x;
    float Ymin = obj->FP->y;
    float Ymax = obj->FP->y;
    for (int i = 0; i < input->n; i++) {
	obj = (pmAstromObj *)input->data[i];
	if (!isfinite(obj->FP->x)) continue;
	if (!isfinite(obj->FP->y)) continue;
	Xmin = PS_MIN (Xmin, obj->FP->x);
	Xmax = PS_MAX (Xmax, obj->FP->x);
	Ymin = PS_MIN (Ymin, obj->FP->y);
	Ymax = PS_MAX (Ymax, obj->FP->y);
    }

    int nX = (Xmax - Xmin) / scale + 10;
    int nY = (Ymax - Ymin) / scale + 10;
    psImage *count = psImageAlloc (nX, nY, PS_TYPE_U32);
    psImageInit (count, 0);

    // accumulate 2D histogram in image
    for (int i = 0; i < input->n; i++) {
	obj = (pmAstromObj *)input->data[i];
	if (!isfinite(obj->FP->x)) continue;
	if (!isfinite(obj->FP->y)) continue;
	int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols);
	int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows);
	count->data.U32[Yi][Xi] ++;
    }

    // determine image statistics
    psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MAX | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    if (!psImageStats(stats, count, NULL, 0)) {
	psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
	psFree(stats);
	psFree(count);
	return NULL;
    }

    if (stats->max < 1) {
	psError(PS_ERR_UNKNOWN, false, "no valid sources in image\n");
	psFree(stats);
	psFree(count);
	return NULL;
    }

    // XXX make this a user option
    float limit = PS_MAX (5.0*stats->sampleStdev, 5.0);
    psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit);

    // find and exclude objects in bad pixels
    psArray *output = psArrayAllocEmpty (input->n);
    for (int i = 0; i < input->n; i++) {
	obj = (pmAstromObj *)input->data[i];
	if (!isfinite(obj->FP->x)) continue;
	if (!isfinite(obj->FP->y)) continue;
	int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols);
	int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows);
	if (count->data.U32[Yi][Xi] > limit) continue;
	psArrayAdd (output, 16, obj);
    }

    psFree(stats);
    psFree(count);
    return output;
}

# if (0)
// make a list of the outlier pixels
psArray *badpix = psArrayAllocEmpty (16);
for (int iy = 0; iy < count->numRows; iy++) {
    for (int ix = 0; ix < count->numCols; ix++) {
	if (count->data.U32[iy][ix] > limit) {
	    psPlane *pixel = psPlaneAlloc();
	    pixel->x = ix;
	    pixel->y = iy;
	    psArrayAdd (badpix, 16, pixel);
	    psFree (pixel);
	}
    }
}

if (badpix->n == 0) {
    psArray *output = psMemIncrRefCounter (input);
    psFree (stats);
    psFree (count);
    psFree (badpix);
    return output;
}
# endif

