# include "pswarp.h"

// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config)
{
    // XXX this implementation currently ignores the use of the region
    psImage *region = NULL;

    psTimerStart ("warp");

    // Get warp parameters
    bool mdok;
    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX");
    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY");
    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE");

    // load the recipe
    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
    psAssert (recipe, "missing recipe %s", PSWARP_RECIPE);

    // output mask bits
    psMaskType maskIn   = psMetadataLookupU8(&mdok, recipe, "MASK.INPUT"); 
    psMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); 
    psMaskType maskBad  = pmConfigMaskGet("BAD.WARP", config); 
    psAssert (mdok, "MASK.INPUT was not defined");

    int nThreads = psMetadataLookupS32 (&mdok, config->arguments, "NTHREADS");
    if (!mdok) nThreads = 0;

    // Flux fraction for "poor"
    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); 

    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
    // output coordinates to input coordinates
    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, nGridX, nGridY);

    // XXX optionally modify the grid based on this result and force the maxError < XXX
    double maxError = pswarpMapGridMaxError (grid);
    psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError);

    // Interpolation options : move these from the arguments to explicit assignments
    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(interpolationMode, input->image, input->weight, input->mask, maskIn, NAN, NAN, maskBad, maskPoor, poorFrac);

    if (input->weight && !output->weight) {
	output->weight = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_F32);
	psImageInit(output->weight, NAN);
    }
    if ((input->mask || maskPoor || maskBad) && !output->mask) {
	output->mask = psImageAlloc(output->image->numCols, output->image->numRows, PS_TYPE_MASK);
	psImageInit(output->mask, maskBad);
    }

    // total number of good pixels across all tiles (summed below)
    int goodPixels = 0;

    // create jobs and supply them to the threads
    for (int gridY = 0; gridY < grid->nYpts; gridY++) {
	for (int gridX = 0; gridX < grid->nXpts; gridX++) {

	    pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();

	    // these items are just views to the data; they are not freed with args
	    args->input = input;
	    args->output = output;
	    args->grid = grid;
	    args->interp = interp;
	    args->region = region;

	    args->gridX = gridX;
	    args->gridY = gridY;
	    args->goodPixels = 0;

	    // allocate a job
	    psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");

	    // construct the arguments for this job
	    // job is pswarpTransformTile (gridX, gridY);
	    psArrayAdd (job->args, 1, args);
	    // fprintf (stderr, "adding job %d,%d, Nargs: %ld\n", gridX, gridY, job->args->n);

	    // call: pswarpTransformTile (args);
	    if (!psThreadJobAddPending (job)) {
		psError(PS_ERR_UNKNOWN, false, "Unable to warp image.");
		return false;
	    }
	    psFree (args);
	}
    }

    // wait for the threads to finish and manage results
    // wait here for the threaded jobs to finish
    if (!psThreadPoolWait ()) {
	psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
	return false;
    }
    fprintf (stderr, "success for threaded jobs\n");

    // each job records its own goodPixel values; sum them here
    // we have only supplied one type of job, so we can assume the types here
    psThreadJob *job = NULL;
    while ((job = psThreadJobGetDone()) != NULL) {
	if (job->args->n < 1) {
	    fprintf (stderr, "error with job\n");
	} else {
	    pswarpTransformTileArgs *args = job->args->data[0];
	    // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
	    goodPixels += args->goodPixels;
	}
	psFree (job);
    }
    psFree(interp);
    psFree(grid);


    // Store the variance factor and number of good pixels
    if (goodPixels > 0) {
	// Variance factor: large factor --> small scale
        float varFactor = psImageInterpolateVarianceFactor(input->image->numCols / 2.0 + input->image->col0, input->image->numRows / 2.0 + input->image->row0, interp); 
        psMetadataItem *vfItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_VARFACTOR);
        if (vfItem) {
            psMetadataItem *goodpixItem = psMetadataLookup(output->analysis, PSWARP_ANALYSIS_GOODPIX);
            psAssert(goodpixItem, "It should be where we left it!");
            psAssert(vfItem->type == PS_TYPE_F32 && goodpixItem->type == PS_TYPE_S64,
                     "Should be the type we said.");

            vfItem->data.F32 += varFactor * goodPixels;
            goodpixItem->data.S64 += goodPixels;
        } else {
            psMetadataAddF32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_VARFACTOR, 0,
                             "Variance factor weighted by the good pixels", varFactor * goodPixels);
            psMetadataAddS64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, 0,
                             "Number of good pixels", goodPixels);
        }
    }

    if (goodPixels > 0) {
	if (!pswarpTransformSources (output, input, config)) {
	    psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
	    return false;
	}
    }

    // XXX should we not write anything out if there are no good pixels?
    output->data_exists = true;
    output->parent->data_exists = true;
    output->parent->parent->data_exists = true;

    psLogMsg ("pswarp", 3, "warping analysis: %f sec\n", psTimerMark ("warp"));

    return true;
}
