// Use a mask description (which can be generated from a ds9 region file using ds9_regions_mask.pl)
// to mask extra pixels.

// Compilation: gcc mask_extra.c -o mask_extra -std=gnu99 -g -Wall `pslib-config --cflags --libs`
// Usage: mask_extra OLD_MASK.fits MASK_DESCRIPTION.dat NEW_MASK.fits

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

static bool readImage(psImage **image, psMetadata **header, const char *name)
{
    psFits *fits = psFitsOpen(name, "r");
    if (header) {
        *header = psFitsReadHeader(NULL, fits);
    }
    if (image) {
        *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0);
    }
    psFitsClose(fits);
    return true;
}

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

int main(int argc, char *argv[])
{
    if (argc != 4) {
        fprintf(stderr, "Not enough arguments\n\n");
        exit(EXIT_FAILURE);
    }

    const char *origName = argv[1];     // Original mask name
    const char *descName = argv[2];     // Mask description file name
    const char *outName = argv[3];      // Output file name

    psImage *mask;                      // Mask image
    psMetadata *header;                 // Header for image

    readImage(&mask, &header, origName);
    int numCols = mask->numCols, numRows = mask->numRows; // Size of image

    psArray *desc = psVectorsReadFromFile(descName, "%d %f %f %f"); // Mask description
    psVector *type = desc->data[0];     // Type of mask
    psVector *x = desc->data[1];        // x coordinates
    psVector *y = desc->data[2];        // y coordinates
    psVector *rad = desc->data[3];      // Radius
    int num = type->n;                  // Number of masks

    int maxType = 0;                    // Maximum number of types
    for (int i = 0; i < num; i++) {
        if (type->data.S32[i] > maxType) {
            maxType = type->data.S32[i];
        }
    }
    psVector *count = psVectorAlloc(maxType + 1, PS_TYPE_S32); // Count of each type
    psVectorInit(count, 0);

    for (int i = 0; i < num; i++) {
        float xCen = x->data.F32[i] - 1, yCen = y->data.F32[i] - 1; // Coordinates of centre
        float radius = rad->data.F32[i];// Radius of circle
        float radius2 = PS_SQR(radius); // Radius of circle, squared

        // Bounds of circle
        int xMin = PS_MAX(xCen - radius, 0), xMax = PS_MIN(xCen + radius, numCols - 1);
        int yMin = PS_MAX(yCen - radius, 0), yMax = PS_MIN(yCen + radius, numRows - 1);

        int numMask = 0;                // Number of pixels masked
        for (int v = yMin; v <= yMax; v++) {
            int v2 = PS_SQR(v - yCen);  // Distance
            for (int u = xMin; u <= xMax; u++) {
                if (PS_SQR(u - xCen) + v2 > radius2) {
                    continue;
                }
                mask->data.PS_TYPE_MASK_DATA[v][u] = 0xFF;
                numMask++;
            }
        }
        count->data.S32[type->data.S32[i]] += numMask;
    }

    psFree(desc);

    long numPix = numCols * numRows;    // Number of pixels
    for (int i = 0; i <= maxType; i++) {
        printf("Type %d: %d (%f%%)\n", i, count->data.S32[i], count->data.S32[i] / (float)numPix * 100.0);
    }

    psFree(count);

    writeImage(mask, header, outName);

    psFree(mask);
    psFree(header);

    return PS_EXIT_SUCCESS;
}
