#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <assert.h>
#include "pslib.h"
#include "stac.h"

bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
              float *xMapDiff, float *yMapDiff, // Difference applied to maps
              const psVector *xSizes,   // Sizes of images in x
              const psVector *ySizes,   // Sizes of images in y
              psArray *maps             // Transformation maps
    )
{
    int nImages = xSizes->n;            // Number of input images
    assert(ySizes->n == nImages);
    assert(xSizes->type.type == PS_TYPE_S32);
    assert(ySizes->type.type == PS_TYPE_S32);
    assert(maps->n == nImages);
    assert(outnx && outny);             // Must be able to write to these

    psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
    psPlane *outCoord = psAlloc(sizeof(psPlane)); // Output coordinates

    // Initial "guess" at limits
    float xMin = INFINITY;
    float xMax = 0.0;
    float yMin = INFINITY;
    float yMax = 0.0;

    for (int i = 0; i < nImages; i++) {
        int xSize = xSizes->data.S32[i];// Size of image in x
        int ySize = ySizes->data.S32[i];// Size of image in y
        psPlaneTransform *map = maps->data[i]; // The map

        psTrace("stac.size", 4, "Image %d:\n", i);

        // Lower left corner
        inCoord->x = 0.0;
        inCoord->y = 0.0;
        (void)psPlaneTransformApply(outCoord, map, inCoord);
        if (outCoord->x < xMin) {
            xMin = outCoord->x;
        } else if (outCoord->x > xMax) {
            xMax = outCoord->x;
        }
        if (outCoord->y < yMin) {
            yMin = outCoord->y;
        } else if (outCoord->y > yMax) {
            yMax = outCoord->y;
        }
        psTrace("stac.size", 4, "Lower left: %f %f\n", outCoord->x, outCoord->y);

        // Lower right corner
        inCoord->x = xSize;
        (void)psPlaneTransformApply(outCoord, map, inCoord);
        if (outCoord->x < xMin) {
            xMin = outCoord->x;
        } else if (outCoord->x > xMax) {
            xMax = outCoord->x;
        }
        if (outCoord->y < yMin) {
            yMin = outCoord->y;
        } else if (outCoord->y > yMax) {
            yMax = outCoord->y;
        }
        psTrace("stac.size", 4, "Lower right: %f %f\n", outCoord->x, outCoord->y);

        // Upper right corner
        inCoord->x = 0.0;
        inCoord->y = ySize;
        (void)psPlaneTransformApply(outCoord, map, inCoord);
        if (outCoord->x < xMin) {
            xMin = outCoord->x;
        } else if (outCoord->x > xMax) {
            xMax = outCoord->x;
        }
        if (outCoord->y < yMin) {
            yMin = outCoord->y;
        } else if (outCoord->y > yMax) {
            yMax = outCoord->y;
        }
        psTrace("stac.size", 4, "Upper right: %f %f\n", outCoord->x, outCoord->y);

        // Upper left corner
        inCoord->x = xSize;
        (void)psPlaneTransformApply(outCoord, map, inCoord);
        if (outCoord->x < xMin) {
            xMin = outCoord->x;
        } else if (outCoord->x > xMax) {
            xMax = outCoord->x;
        }
        if (outCoord->y < yMin) {
            yMin = outCoord->y;
        } else if (outCoord->y > yMax) {
            yMax = outCoord->y;
        }
        psTrace("stac.size", 4, "Upper left: %f %f\n", outCoord->x, outCoord->y);
    }

    // Tweak the maps to account for the offset
    if (xMapDiff != NULL) {
        *xMapDiff = floorf(xMin);
    }
    if (yMapDiff != NULL) {
        *yMapDiff = floorf(yMin);
    }
    for (int i = 0; i < nImages; i++) {
        psPlaneTransform *map = maps->data[i]; // The map of interest
        map->x->coeff[0][0] -= floorf(xMin);
        map->y->coeff[0][0] -= floorf(yMin);
    }

    *outnx = (int)ceilf(xMax - xMin);
    *outny = (int)ceilf(yMax - yMin);

    psTrace("stac.size", 1, "Output size is to be %dx%d\n", *outnx, *outny);

    psFree(inCoord);
    psFree(outCoord);

    return true;
}
