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

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

int main(int argc, char *argv[])
{
    if (argc != 3) {
        printf("Usage: %s FITS_IN FITS_OUT\n", argv[0]);
        exit(EXIT_FAILURE);
    }
    const char *inName = argv[1];       // Input filename
    const char *outName = argv[2];      // Output filename

    psFits *inFile = psFitsOpen(inName, "r"); // File to read
    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read

    // We only read PHUs --- not mucking around with extensions for now
    psImage *image = psFitsReadImage(inFile, imageRegion, 0);
    if (image == NULL) {
        psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", inName);
        exit(EXIT_FAILURE);
    }
    psTrace("calcGradients.read", 4, "Image %s is %dx%d\n", inName, image->numCols, image->numRows);
    // Convert to 32-bit floating point, in necessary
    if (image->type.type != PS_TYPE_F32) {
        psTrace("calcGradients.read", 3, "Converting %s to floating point in memory....", inName);
        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
        psFree(image);
        image = temp;
    } else {
        int numNaN = psImageClipNaN(image, FLT_MAX);
        if (numNaN) {
            psTrace("calcGradients.read", 5, "Clipped %d NaN pixels.\n", numNaN);
        }
    }
    psFitsClose(inFile);

    psImage *grad = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
    for (int y = 0; y < image->numRows; y++) {
        for (int x = 0; x < image->numCols; x++) {
            grad->data.F32[y][x] = stacGradient(image, x, y);
        }
    }

    psFree(image);

    psFits *outFile = psFitsOpen(outName, "w");
    if (!psFitsWriteImage(outFile, NULL, grad, 0, NULL)) {
        psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
    }
    psFitsClose(outFile);

    return EXIT_SUCCESS;
}

