Index: /trunk/tools/mask_extra.c
===================================================================
--- /trunk/tools/mask_extra.c	(revision 20601)
+++ /trunk/tools/mask_extra.c	(revision 20601)
@@ -0,0 +1,103 @@
+// 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;
+                num++;
+            }
+        }
+        count->data.S32[type->data.S32[i]] += numMask;
+    }
+
+    psFree(desc);
+
+    for (int i = 0; i <= maxType; i++) {
+        printf("Type %d: %d\n", i, count->data.S32[i]);
+    }
+
+    psFree(count);
+
+    writeImage(mask, header, outName);
+
+    psFree(mask);
+    psFree(header);
+
+    return PS_EXIT_SUCCESS;
+}