/** @file pswarpTransformReadout.c
 *
 *  @brief
 *
 *  @ingroup pswarp
 *
 *  @author IfA
 *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-18 02:58:59 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

#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;             ///< Region to transform

    psTimerStart("warp");

    // Get warp parameters
    bool mdok;                          ///< Status of MD lookup
    int nGridX = psMetadataLookupS32(NULL, config->arguments, "GRID.NX"); ///< Number of grid points in x
    int nGridY = psMetadataLookupS32(NULL, config->arguments, "GRID.NY"); ///< Number of grid points in y
    psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
                                                                   "INTERPOLATION.MODE"); ///< Mode for interp
    int numKernels = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.NUM"); ///< Number of kernels

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

    // output mask bits
    psImageMaskType maskIn   = psMetadataLookupImageMask(&mdok, recipe, "MASK.INPUT");
    psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config);
    if (!maskPoor) {
        maskPoor = pmConfigMaskGet("POOR.WARP", config);
    }
    psImageMaskType maskBad  = pmConfigMaskGet("CONV.BAD", config);
    if (!maskBad) {
        maskBad  = pmConfigMaskGet("BAD.WARP", config);
    }
    psAssert(mdok, "MASK.INPUT was not defined");

    if (!input->covariance) {
        input->covariance = psImageCovarianceNone();
    }

    int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads
    if (!mdok) {
        nThreads = 0;
    }
    float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); ///< Flux fraction for "poor"

    // 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); // Maximum (positional) error from using grid
    psLogMsg("pswarp", 3, "maximum error using this grid sampling: %lf\n", maxError);

    // Get range of interest
    int xOutMin, xOutMax, yOutMin, yOutMax; ///< Output pixel range
    pswarpMatchRange(&xOutMin, &yOutMin, &xOutMax, &yOutMax, input, output);

    // Check the range of the grid coordinates
#define CHECK_GRID_RANGE() { \
        xGridMin = PS_MIN(xGridMin, xGrid); \
        xGridMax = PS_MAX(xGridMax, xGrid); \
        yGridMin = PS_MIN(yGridMin, yGrid); \
        yGridMax = PS_MAX(yGridMax, yGrid); \
    }

    int xGridMin = nGridX, xGridMax = 0, yGridMin = nGridY, yGridMax = 0; ///< Grid range
    int xGrid, yGrid;                   ///< Grid coordinates
    pswarpMapGridSetGrid(grid, xOutMin, yOutMin, &xGrid, &yGrid);
    CHECK_GRID_RANGE();
    pswarpMapGridSetGrid(grid, xOutMin, yOutMax, &xGrid, &yGrid);
    CHECK_GRID_RANGE();
    pswarpMapGridSetGrid(grid, xOutMax, yOutMin, &xGrid, &yGrid);
    CHECK_GRID_RANGE();
    pswarpMapGridSetGrid(grid, xOutMax, yOutMax, &xGrid, &yGrid);
    CHECK_GRID_RANGE();

    if (!pmReadoutMaskNonfinite(input, pmConfigMaskGet("SAT", config))) {
        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in input.");
        return false;
    }

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

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

    // Ensure threading is off for the covariance calculation, since we are threading on a different level.
    psImageCovarianceSetThreads(false);

    psAssert (xGridMin >= 0, "xGridMin too small\n");
    psAssert (yGridMin >= 0, "yGridMin too small\n");
    psAssert (xGridMax < grid->nXpts, "xGridMax too big\n");
    psAssert (yGridMax < grid->nYpts, "yGridMax too big\n");

    // create jobs and supply them to the threads
    for (int gridY = yGridMin; gridY <= yGridMax; gridY++) {
        for (int gridX = xGridMin; gridX <= xGridMax; gridX++) {
            pswarpTransformTileArgs *args = pswarpTransformTileArgsAlloc();
            args->input = psMemIncrRefCounter(input);
            args->output = psMemIncrRefCounter(output);
            args->grid = psMemIncrRefCounter(grid);
            args->interp = psMemIncrRefCounter(interp);
            args->region = psMemIncrRefCounter(region);

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

            // allocate a job
            psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
            psArrayAdd(job->args, 1, args);
            if (!psThreadJobAddPending(job)) {
                psError(psErrorCodeLast(), 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 (false, true)) {
        psError(psErrorCodeLast(), false, "Unable to interpolate image.");
        return false;
    }

    // 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;
    int xMin = output->image->numCols, xMax = 0, yMin = output->image->numRows, yMax = 0; // Bounds
    int goodPixels = psMetadataLookupS32(&mdok, output->analysis, PSWARP_ANALYSIS_GOODPIX); // Number of pixels
    psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
                                              PSWARP_ANALYSIS_COVARIANCES); // Collection of covar. matrices
    if (!covariances) {
        covariances = psListAlloc(NULL);
        psMetadataAddList(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_COVARIANCES, 0,
                          "Collection of covariance matrices", covariances);
        psFree(covariances);            // Drop reference; still have the copy on the analysis metadata
    }
    double jacobian = psMetadataLookupF64(&mdok, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
    if (!isfinite(jacobian)) {
        jacobian = 0.0;
    }

    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;
            xMin = PS_MIN(args->xMin, xMin);
            xMax = PS_MAX(args->xMax, xMax);
            yMin = PS_MIN(args->yMin, yMin);
            yMax = PS_MAX(args->yMax, yMax);
            if (args->covariance) {
                psListAdd(covariances, PS_LIST_TAIL, args->covariance);
            }
            if (args->goodPixels > 0 && isfinite(args->jacobian)) {
                jacobian += args->jacobian * args->goodPixels;
            }
        }
        psFree(job);
    }
    psFree(grid);

    psMetadataAddS32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, PS_META_REPLACE,
                     "Number of good pixels", goodPixels);
    psMetadataAddF64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_JACOBIAN, PS_META_REPLACE,
                     "Jacobian of transformation", jacobian);

    if (xMin < xMax && yMin < yMax) {
        psTrace("pswarp.transform", 1, "Bounds [%d:%d,%d:%d]\n", xMin, xMax, yMin, yMax);
    } else {
        psTrace("pswarp.transform", 1, "No overlap\n");
    }
    psFree(interp);

    if (goodPixels > 0 && psMetadataLookupBool(&mdok, recipe, "SOURCES")) {
        if (!pswarpTransformSources(output, input, config)) {
            psError(psErrorCodeLast(), false, "Unable to interpolate image.");
            return false;
        }
    }

    if (goodPixels > 0) {
        // Data is only written out if there are 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;
}
