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

#define SIGMA 3.0                       // Gaussian width
#define NUM_SIGMA 3                     // Size of convolution

static psImage *readImage(const char *name)
{
    psFits *fits = psFitsOpen(name, "r");
    psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0);
    psFitsClose(fits);
    return image;
}

void writeImage(const psImage *image,   // Image to write
                const char *name        // Name of FITS file
    )
{
    psFits *fits = psFitsOpen(name, "w");
    psFitsWriteImage(fits, NULL, image, 0, NULL);
    psFitsClose(fits);
}

int main(int argc, char *argv[])
{
    psImage *bright = readImage("bright.conv.fits");
    psImage *faint = readImage("faint.conv.fits");
    psImage *weight = readImage("weight.conv.fits");
    psImage *mask = readImage("mask.conv.fits");

    // smooth the image, applying the mask as we go
    psImageSmoothMaskF32(bright, mask, 0xFF, SIGMA, NUM_SIGMA);
    psImageSmoothMaskF32(faint, mask, 0xFF, SIGMA, NUM_SIGMA);
    psImageSmoothMaskF32(weight, mask, 0xFF, SIGMA/M_SQRT2, NUM_SIGMA);

    float factor = 4.0 * M_PI * PS_SQR(SIGMA);

    for (int j = 0; j < bright->numRows; j++) {
        for (int i = 0; i < bright->numCols; i++) {
            if (mask->data.PS_TYPE_MASK_DATA[j][i]) {
                bright->data.F32[j][i] = 0.0;
                faint->data.F32[j][i] = 0.0;
            } else {
                bright->data.F32[j][i] = factor * PS_SQR(bright->data.F32[j][i]) /
                    weight->data.F32[j][i];
                faint->data.F32[j][i] = factor * PS_SQR(faint->data.F32[j][i]) /
                    weight->data.F32[j][i];
            }
        }
    }

    psFree(weight);
    psFree(mask);

    writeImage(bright, "bright.sig.fits");
    writeImage(faint, "faint.sig.fits");

    psFree(bright);
    psFree(faint);

    return PS_EXIT_SUCCESS;
}
