# include "psphotInternal.h" 

bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects);
 
psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view) 
{
    bool status = true;

    psArray *objects = psArrayAllocEmpty(100);

    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");

    // loop over the available readouts
    for (int i = 0; i < num; i++) {
        if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) {
            psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
	    psFree (objects);
            return NULL;
        }
    }

    psphotMatchSourcesGenerate (config, view, objects);

    return objects;
}

bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { 
 
    bool status = false;

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    psAssert (recipe, "missing recipe?");

    int RADIUS = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.MATCH.RADIUS");
    psAssert (status, "programming error: must define PSPHOT.STACK.MATCH.RADIUS");

    // find the currently selected readout
    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    psAssert (file, "missing file?");

    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    psAssert (readout, "missing readout?");

    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    psAssert (detections, "missing detections?");
    psAssert (detections->newSources, "new sources not defined?");
    psAssert (!detections->allSources, "all sources already defined?");

    // XXX TEST:
    if (detections->newSources) {
        psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS);
    }

    return true;
}

# define NEXT1 { i++; continue; } 
# define NEXT2 { j++; continue; } 
bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { 
 
    float dx, dy; 
 
    float RADIUS2 = RADIUS*RADIUS;

    // sort the source list by X 
    sources = psArraySort (sources, pmSourceSortByX); 
    objects = psArraySort (objects, pmPhotObjSortByX); 
 
    psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8);
    psVectorInit (found, 0);

    // match sources to existing objects

    psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld)", sources->n, objects->n);

    int i, j; 
    for (i = j = 0; (i < sources->n) && (j < objects->n); ) { 
 
        pmSource  *src = sources->data[i]; 
        pmPhotObj *obj = objects->data[j]; 
 
        if (!src) NEXT1; 
        if (!src->peak) NEXT1; 
        if (!finite(src->peak->xf)) NEXT1; 
        if (!finite(src->peak->yf)) NEXT1; 
 
        if (!obj) NEXT2; 
        if (!finite(obj->x)) NEXT2; 
        if (!finite(obj->y)) NEXT2; 
 
        dx = src->peak->xf - obj->x; 
        if (dx < -1.02*RADIUS) NEXT1; 
        if (dx > +1.02*RADIUS) NEXT2; 
 
        // we are within match range, look for matches: 
	int Jmin = -1;
	float Rmin = RADIUS2;
        for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { 
 
	    obj = objects->data[J]; 
	    
	    dx = src->peak->xf - obj->x; 
            dy = src->peak->yf - obj->y; 
 
            float dr = dx*dx + dy*dy; 
            if (dr > RADIUS2) continue; 
	    if (dr < Rmin) {
		Rmin = dr;
		Jmin  = J;
	    }
	}

	// no match, try next source
	if (Jmin == -1) {
	    i++;
	    continue;
	}
	obj = objects->data[Jmin]; 

	// add to object
	pmPhotObjAddSource (obj, src);
	found->data.U8[i] = 1;
        i++; 
    } 

    // add missed sources to new objects

    for (i = 0; i < sources->n; i++) {

	if (found->data.U8[i]) continue;

        pmSource *src = sources->data[i]; 

	pmPhotObj *obj = pmPhotObjAlloc();
	pmPhotObjAddSource(obj, src);
	psArrayAdd (objects, 100, obj);
	psFree (obj);
    }
    psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n);

    psFree (found);
    return true;
} 

bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) {

    bool status = false;

    // select the appropriate recipe information
    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    psAssert (recipe, "missing recipe?");

    // determine properties (sky, moments) of initial sources
    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    psAssert (status, "missing SKY_OUTER_RADIUS in recipe?");

    int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");

    // generate look-up arrays for detections and readouts
    psArray *detArrays = psArrayAlloc(nImages);
    psArray *readouts = psArrayAlloc(nImages);

    for (int i = 0; i < nImages; i++) {

	// find the currently selected readout
	pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
	psAssert (file, "missing file?");

	pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
	psAssert (readout, "missing readout?");

	pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
	psAssert (detections, "missing detections?");

	detArrays->data[i] = psMemIncrRefCounter(detections);
	readouts->data[i] = psMemIncrRefCounter(readout);
    }

    // vector to track if source for an image is found
    psVector *found = psVectorAlloc(nImages, PS_TYPE_U8);

    for (int i = 0; i < objects->n; i++) { 
        pmPhotObj *obj = objects->data[i]; 

	// mark the images for which sources have been found
	psVectorInit (found, 0);
	for (int j = 0; j < obj->sources->n; j++) {

	    pmSource *src = obj->sources->data[j]; 
	    int index = src->imageID;
	    psAssert (index >= 0, "invalid index");
	    psAssert (index < found->n, "invalid index");

	    found->data.U8[index] = 1;
	}

	// generate new sources for the image that are missing
	for (int index = 0; index < found->n; index++) {
	    if (found->data.U8[index]) continue;

	    pmDetections *detections = detArrays->data[index];
	    pmReadout *readout = readouts->data[index];
	    int row0 = readout->image->row0;
	    int col0 = readout->image->col0;

	    // XXX the peak type is not really used in psphot
	    // PM_PEAK_LONE is certainly not true, but irrelevant
	    float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)];
	    pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE);
	    peak->flux = peakFlux;
	    peak->SN = 1.0;
	    peak->xf = obj->x;
	    peak->yf = obj->y;
	    peak->dx = NAN;
	    peak->dy = NAN;
	    
	    // XXX assign to a footprint?

	    // create a new source
	    pmSource *source = pmSourceAlloc();
	    source->imageID = index;

	    // add the peak
	    source->peak = psMemIncrRefCounter(peak);

	    // allocate space for moments
	    source->moments = pmMomentsAlloc();

	    // allocate image, weight, mask arrays for each peak (square of radius OUTER)
	    pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);

	    peak->assigned = true;
	    pmPhotObjAddSource(obj, source);
	    psArrayAdd (detections->newSources, 100, source);
	    psFree (source);
	    psFree (peak);
	}
    }

    // how many sources do we have now?
    int nSources = 0;
    for (int i = 0; i < objects->n; i++) { 
        pmPhotObj *obj = objects->data[i]; 
	nSources += obj->sources->n;
    }
    psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages);


    psFree (found);
    psFree (detArrays);
    psFree (readouts);
    return true;
}
