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

// Brain-dead simple stack of multiple frames
// Makes no effort to zap bad pixels, or rescale images.
// Assumes that images are already registered

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

int main(int argc, char *argv[])
{
    const char *outputName = argv[1];   // Output file name
    psArray *inputNames = psArrayAlloc(argc-2);
    for (int i = 2; i < argc; i++) {
        inputNames->data[i-2] = psStringCopy(argv[i]);
    }

    int numCols = 0;                    // Number of columns
    int numRows = 0;                    // Number of rows
    psImage *output = NULL;             // Output image
    psImage *scale = NULL;              // Scale image
    psMetadata *header = NULL;          // FITS header for output image
    for (int i = 0; i < inputNames->n; i++) {
        psFits *fits = psFitsOpen(inputNames->data[i], "r"); // FITS file
        psRegion readRegion = {0, 0, 0, 0}; // Region to read
        psImage *image = psFitsReadImage(fits, readRegion, 0);
        if (numCols == 0 && numRows == 0) {
            // First run through --- set the size, initialise the output, read the header
            numCols = image->numCols;
            numRows = image->numRows;
            output = psImageAlloc(numCols, numRows, PS_TYPE_F32);
            scale = psImageAlloc(numCols, numRows, PS_TYPE_U8);
            for (int y = 0; y < numCols; y++) {
                for (int x = 0; x < numRows; x++) {
                    output->data.F32[y][x] = 0.0;
                    scale->data.U8[y][x] = 0;
                }
            }
            header = psFitsReadHeader(NULL, fits);
        } else if (image->numCols != numCols || image->numRows != numRows) {
            psError(PS_ERR_IO, true, "Input images have different dimensions.  %s is %dx%d, but should be "
                    "%dx%d\n", (char *)inputNames->data[i], image->numCols, image->numRows, numCols, numRows);
            return EXIT_FAILURE;
        }
        psFitsClose(fits);

        // Convert the type, if necessary
        if (image->type.type != PS_TYPE_F32) {
            psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
            psFree(image);
            image = temp;
        }

        // Add in
        for (int y = 0; y < numCols; y++) {
            for (int x = 0; x < numRows; x++) {
                if (image->data.F32[y][x] > 0) {
                    output->data.F32[y][x] += image->data.F32[y][x];
                    scale->data.U8[y][x]++;
                }
            }
        }

        psFree(image);
    }
    psFree(inputNames);

    for (int y = 0; y < numCols; y++) {
        for (int x = 0; x < numRows; x++) {
            if (scale->data.U8[y][x] > 0) {
                output->data.F32[y][x] /= (float)scale->data.U8[y][x];
            } else {
                output->data.F32[y][x] = 0.0;
            }
        }
    }

    psFits *fits = psFitsOpen(outputName, "w");
    (void)psFitsWriteImage(fits, header, output, 0, NULL);
    psFitsClose(fits);
    psFree(output);
    // Pau.
    return EXIT_SUCCESS;
}
