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

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

#define BUFFER 100                      // Size of buffer for incrementally reading coordinates

psArray *stacReadImages(psArray **headers, // The image headers, to be returned
    psArray *filenames // The file names of the images
    )
{
    int nFiles = filenames->n;          // The number of input files
    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
    assert(!headers || ! *headers || (*headers)->n == nFiles);
    if (headers && ! *headers) {
        *headers = psArrayAlloc(nFiles);
    }

    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read

    psTrace("stac.read.images", 1, "Reading input images....\n");
    for (int i = 0; i < nFiles; i++) {
        psTrace("stac.read.images", 2, "Reading input image %s....\n", (char *)filenames->data[i]);
        psFits *imageFile = psFitsOpen(filenames->data[i], "r");
        if (!imageFile) {
            psError(PS_ERR_IO, false, "Unable to open FITS file %s.\n", (char *)filenames->data[i]);
            exit(EXIT_FAILURE);
        }
        // We only read PHUs --- not mucking around with extensions for now
        if (headers) {
            (*headers)->data[i] = psFitsReadHeader(NULL, imageFile);
        }
        psImage *image = psFitsReadImage(imageFile, imageRegion, 0);
        if (image == NULL) {
            psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", (char *)filenames->data[i]);
            exit(EXIT_FAILURE);
        }
        psFitsClose(imageFile);
        psTrace("stac.read.images", 4, "Image %s is %dx%d\n", (char *)filenames->data[i], image->numCols,
                image->numRows);
        // Convert to 32-bit floating point, in necessary
        if (image->type.type != PS_TYPE_F32) {
            psTrace("stac.read.images", 3, "Converting %s to floating point in memory....\n",
                    (char *)filenames->data[i]);
            psImage *copy = psImageCopy(NULL, image, PS_TYPE_F32);
            psFree(image);
            image = copy;
        }
        int numNaN = psImageClipNaN(image, 0.0);
        if (numNaN) {
            psTrace("stac.read.images", 5, "Clipped %d NaN pixels to zero.\n", numNaN);
        }
        images->data[i] = image;
    }
    psTrace("stac.read.images",1,"%d input images read.\n",nFiles);

    return images;
}

psArray *stacReadCoords(const char *filename)
{
    FILE *file = fopen(filename, "r");
    if (file == NULL) {
        psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
        return NULL;
    }

    psTrace("stac.read.coords", 5, "Reading coordinate file, %s\n", filename);

    psArray *coords = psArrayAllocEmpty(BUFFER); // The array of coordinates to be returned
    float x, y;                         // Coordinates to read
    while (fscanf(file, "%f %f\n", &x, &y) == 2) {
        psPlane *coord = psPlaneAlloc();// A coordinate
        coord->x = x;
        coord->y = y;
        coords->data[coords->n] = coord;
        coords->n++;
        if (coords->n % BUFFER == 0) {
            coords = psArrayRealloc(coords, coords->n + BUFFER);
        }
    }

    psTrace("stac.read.coords", 5, "%ld coordinates read.\n", coords->n);

    fclose(file);
    return coords;
}


psPlaneTransform *stacReadMap(const char *filename)
{
    FILE *mapfp = fopen(filename, "r");
    if (mapfp == NULL) {
        psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
        return NULL;
    }
    // Read the file
    psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);

    // Format is now:
    // order
    // x coefficients
    // y coefficients
    //
    // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
    // and the coefficients are read by the following pseudo-code:
    //
    //  for (int k = 0; k < order + 1; k++)
    //      for (int j = 0; j < k; j++)
    //          int i = k - j;
    //              read coefficient of x^i y^j
    //
    // This produces the ordering:
    // x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
    //
    // This is, of course, for third order polynomials.
    // For lower orders, the list is truncated at the appropriate level.

    int order = 0;                      // Polynomial order
    if (fscanf(mapfp, "%d", &order) != 1) {
        psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
        fclose(mapfp);
        return NULL;
    }

    psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);

    psPlaneTransform *map = psPlaneTransformAlloc(order, order); // The transformation
    // Set coefficient masks
    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order + 1) {
                map->x->mask[i][j] = 1;
                map->y->mask[i][j] = 1;
            } else {
                map->x->mask[i][j] = 0;
                map->y->mask[i][j] = 0;
            }
        }
    }

    // Read x coefficients
    psTrace("stac.read.map", 7, "x' = \n");
    for (int k = 0; k <= order; k++) {
        for (int j = 0; j <= k; j++) {
            int i = k - j;
            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
                fclose(mapfp);
                psFree(map);
                return NULL;
            }
            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
        }
    }
    // Read y coefficients
    psTrace("stac.read.maps", 7, "y' = \n");
    for (int k = 0; k <= order; k++) {
        for (int j = 0; j <= k; j++) {
            int i = k - j;
            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
                fclose(mapfp);
                psFree(map);
                return NULL;
            }
            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
        }
    }

    fclose(mapfp);

    return map;
}



psArray *stacReadMaps(psArray *filenames // The file names of the images whose maps are to be read
    )
{
    int nFiles = filenames->n;          // The number of input files
    psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
    char mapfile[MAXCHAR];              // Filename of map

    psTrace("stac.read.maps",1,"Reading maps....\n");
    for (int i = 0; i < nFiles; i++) {
        if (strlen(filenames->data[i]) > MAXCHAR - 4) {
            psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",(char *)filenames->data[i]);
            exit(EXIT_FAILURE);
        }
        // Read the file
        sprintf(mapfile, "%s.map", (const char*)filenames->data[i]);
        maps->data[i] = stacReadMap(mapfile);
        if (maps->data[i] == NULL) {
            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
        }

    }
    psTrace("stac.read.maps",1,"%d maps read.\n",nFiles);

    return maps;
}
