add Tom's masking code

makemask/cookbook.c

+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/*
+  Every program which uses the CFITSIO interface must include the
+  the fitsio.h header file.  This contains the prototypes for all
+  the routines and defines the error status values and other symbolic
+  constants used in the interface.  
+*/
+#include <fitsio.h>
+
+int main( void );
+void writeimage( void );
+void writeascii( void );
+void writebintable( void );
+void copyhdu( void );
+void selectrows( void );
+void readheader( void );
+void readimage( void );
+void readtable( void );
+void printerror( int status);
+
+int main()
+{
+/*************************************************************************
+   This is a simple main program that calls the following routines:
+
+    writeimage    - write a FITS primary array image
+    writeascii    - write a FITS ASCII table extension
+    writebintable - write a FITS binary table extension
+    copyhdu       - copy a header/data unit from one FITS file to another
+    selectrows    - copy selected row from one HDU to another
+    readheader    - read and print the header keywords in every extension
+    readimage     - read a FITS image and compute the min and max value
+    readtable     - read columns of data from ASCII and binary tables
+
+**************************************************************************/
+
+    writeimage();
+    writeascii();
+    writebintable();
+    copyhdu();
+    selectrows();
+    readheader();
+    readimage();
+    readtable();
+
+    printf("\nAll the cfitsio cookbook routines ran successfully.\n");
+    return(0);
+}
+/*--------------------------------------------------------------------------*/
+void writeimage( void )
+
+    /******************************************************/
+    /* Create a FITS primary array containing a 2-D image */
+    /******************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status, ii, jj;
+    long  fpixel, nelements, exposure;
+    unsigned short *array[200];
+
+    /* initialize FITS image parameters */
+    char filename[] = "atestfil.fit";             /* name for new FITS file */
+    int bitpix   =  USHORT_IMG; /* 16-bit unsigned short pixel values       */
+    long naxis    =   2;  /* 2-dimensional image                            */    
+    long naxes[2] = { 300, 200 };   /* image is 300 pixels wide by 200 rows */
+
+    /* allocate memory for the whole image */ 
+    array[0] = (unsigned short *)malloc( naxes[0] * naxes[1]
+                                        * sizeof( unsigned short ) );
+
+    /* initialize pointers to the start of each row of the image */
+    for( ii=1; ii<naxes[1]; ii++ )
+      array[ii] = array[ii-1] + naxes[0];
+
+    remove(filename);               /* Delete old file if it already exists */
+
+    status = 0;         /* initialize status before calling fitsio routines */
+
+    if (fits_create_file(&fptr, filename, &status)) /* create new FITS file */
+         printerror( status );           /* call printerror if error occurs */
+
+    /* write the required keywords for the primary array image.     */
+    /* Since bitpix = USHORT_IMG, this will cause cfitsio to create */
+    /* a FITS image with BITPIX = 16 (signed short integers) with   */
+    /* BSCALE = 1.0 and BZERO = 32768.  This is the convention that */
+    /* FITS uses to store unsigned integers.  Note that the BSCALE  */
+    /* and BZERO keywords will be automatically written by cfitsio  */
+    /* in this case.                                                */
+
+    if ( fits_create_img(fptr,  bitpix, naxis, naxes, &status) )
+         printerror( status );          
+
+    /* initialize the values in the image with a linear ramp function */
+    for (jj = 0; jj < naxes[1]; jj++)
+    {   for (ii = 0; ii < naxes[0]; ii++)
+        {
+            array[jj][ii] = ii + jj;
+        }
+    }
+
+    fpixel = 1;                               /* first pixel to write      */
+    nelements = naxes[0] * naxes[1];          /* number of pixels to write */
+
+    /* write the array of unsigned integers to the FITS file */
+    if ( fits_write_img(fptr, TUSHORT, fpixel, nelements, array[0], &status) )
+        printerror( status );
+      
+    free( array[0] );  /* free previously allocated memory */
+
+    /* write another optional keyword to the header */
+    /* Note that the ADDRESS of the value is passed in the routine */
+    exposure = 1500;
+    if ( fits_update_key(fptr, TLONG, "EXPOSURE", &exposure,
+         "Total Exposure Time", &status) )
+         printerror( status );           
+
+    if ( fits_close_file(fptr, &status) )                /* close the file */
+         printerror( status );           
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void writeascii ( void )
+
+    /*******************************************************************/
+    /* Create an ASCII table extension containing 3 columns and 6 rows */
+    /*******************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status;
+    long firstrow, firstelem;
+
+    int tfields = 3;       /* table will have 3 columns */
+    long nrows  = 6;       /* table will have 6 rows    */
+
+    char filename[] = "atestfil.fit";           /* name for new FITS file */
+    char extname[] = "PLANETS_ASCII";             /* extension name */
+
+    /* define the name, datatype, and physical units for the 3 columns */
+    char *ttype[] = { "Planet", "Diameter", "Density" };
+    char *tform[] = { "a8",     "I6",       "F4.2"    };
+    char *tunit[] = { "\0",      "km",       "g/cm"    };
+
+    /* define the name diameter, and density of each planet */
+    char *planet[] = {"Mercury", "Venus", "Earth", "Mars","Jupiter","Saturn"};
+    long  diameter[] = {4880,    12112,    12742,   6800,  143000,   121000};
+    float density[]  = { 5.1f,     5.3f,      5.52f,   3.94f,    1.33f,    0.69f};
+
+    status=0;
+
+    /* open with write access the FITS file containing a primary array */
+    if ( fits_open_file(&fptr, filename, READWRITE, &status) ) 
+         printerror( status );
+
+    /* append a new empty ASCII table onto the FITS file */
+    if ( fits_create_tbl( fptr, ASCII_TBL, nrows, tfields, ttype, tform,
+                tunit, extname, &status) )
+         printerror( status );
+
+    firstrow  = 1;  /* first row in table to write   */
+    firstelem = 1;  /* first element in row  (ignored in ASCII tables) */
+
+    /* write names to the first column (character strings) */
+    /* write diameters to the second column (longs) */
+    /* write density to the third column (floats) */
+
+    fits_write_col(fptr, TSTRING, 1, firstrow, firstelem, nrows, planet,
+                   &status);
+    fits_write_col(fptr, TLONG, 2, firstrow, firstelem, nrows, diameter,
+                   &status);
+    fits_write_col(fptr, TFLOAT, 3, firstrow, firstelem, nrows, density,
+                   &status);
+
+    if ( fits_close_file(fptr, &status) )       /* close the FITS file */
+         printerror( status );
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void writebintable ( void )
+
+    /*******************************************************************/
+    /* Create a binary table extension containing 3 columns and 6 rows */
+    /*******************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status, hdutype;
+    long firstrow, firstelem;
+
+    int tfields   = 3;       /* table will have 3 columns */
+    long nrows    = 6;       /* table will have 6 rows    */
+
+    char filename[] = "atestfil.fit";           /* name for new FITS file */
+    char extname[] = "PLANETS_Binary";           /* extension name */
+
+    /* define the name, datatype, and physical units for the 3 columns */
+    char *ttype[] = { "Planet", "Diameter", "Density" };
+    char *tform[] = { "8a",     "1J",       "1E"    };
+    char *tunit[] = { "\0",      "km",       "g/cm"    };
+
+    /* define the name diameter, and density of each planet */
+    char *planet[] = {"Mercury", "Venus", "Earth", "Mars","Jupiter","Saturn"};
+    long  diameter[] = {4880,     12112,   12742,   6800,  143000,   121000};
+    float density[]  = { 5.1f,      5.3f,     5.52f,   3.94f,   1.33f,     0.69f};
+
+    status=0;
+
+    /* open the FITS file containing a primary array and an ASCII table */
+    if ( fits_open_file(&fptr, filename, READWRITE, &status) ) 
+         printerror( status );
+
+    if ( fits_movabs_hdu(fptr, 2, &hdutype, &status) ) /* move to 2nd HDU */
+         printerror( status );
+
+    /* append a new empty binary table onto the FITS file */
+    if ( fits_create_tbl( fptr, BINARY_TBL, nrows, tfields, ttype, tform,
+                tunit, extname, &status) )
+         printerror( status );
+
+    firstrow  = 1;  /* first row in table to write   */
+    firstelem = 1;  /* first element in row  (ignored in ASCII tables) */
+
+    /* write names to the first column (character strings) */
+    /* write diameters to the second column (longs) */
+    /* write density to the third column (floats) */
+
+    fits_write_col(fptr, TSTRING, 1, firstrow, firstelem, nrows, planet,
+                   &status);
+    fits_write_col(fptr, TLONG, 2, firstrow, firstelem, nrows, diameter,
+                   &status);
+    fits_write_col(fptr, TFLOAT, 3, firstrow, firstelem, nrows, density,
+                   &status);
+
+    if ( fits_close_file(fptr, &status) )       /* close the FITS file */
+         printerror( status );
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void copyhdu( void)
+{
+    /********************************************************************/
+    /* copy the 1st and 3rd HDUs from the input file to a new FITS file */
+    /********************************************************************/
+    fitsfile *infptr;      /* pointer to the FITS file, defined in fitsio.h */
+    fitsfile *outfptr;                 /* pointer to the new FITS file      */
+
+    char infilename[]  = "atestfil.fit";  /* name for existing FITS file   */
+    char outfilename[] = "btestfil.fit";  /* name for new FITS file        */
+
+    int status, morekeys, hdutype;
+
+    status = 0;
+
+    remove(outfilename);            /* Delete old file if it already exists */
+
+    /* open the existing FITS file */
+    if ( fits_open_file(&infptr, infilename, READONLY, &status) ) 
+         printerror( status );
+
+    if (fits_create_file(&outfptr, outfilename, &status)) /*create FITS file*/
+         printerror( status );           /* call printerror if error occurs */
+
+    /* copy the primary array from the input file to the output file */
+    morekeys = 0;     /* don't reserve space for additional keywords */
+    if ( fits_copy_hdu(infptr, outfptr, morekeys, &status) )
+         printerror( status );
+
+    /* move to the 3rd HDU in the input file */
+    if ( fits_movabs_hdu(infptr, 3, &hdutype, &status) )
+         printerror( status );
+
+    /* copy 3rd HDU from the input file to the output file (to 2nd HDU) */
+    if ( fits_copy_hdu(infptr, outfptr, morekeys, &status) )
+         printerror( status );
+
+    if (fits_close_file(outfptr, &status) ||
+        fits_close_file(infptr, &status)) /* close files */
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void selectrows( void )
+
+    /*********************************************************************/
+    /* select rows from an input table and copy them to the output table */
+    /*********************************************************************/
+{
+    fitsfile *infptr, *outfptr;  /* pointer to input and output FITS files */
+    unsigned char *buffer;
+    char card[FLEN_CARD];
+    int status, hdutype, nkeys, keypos, nfound, colnum, anynulls, ii;
+    long naxes[2], frow, felem, noutrows, irow;
+    float nullval, density[6];
+
+    char infilename[]  = "atestfil.fit";  /* name for existing FITS file   */
+    char outfilename[] = "btestfil.fit";  /* name for new FITS file        */
+
+    status = 0;
+
+    /* open the existing FITS files */
+    if ( fits_open_file(&infptr,  infilename,  READONLY,  &status) ||
+         fits_open_file(&outfptr, outfilename, READWRITE, &status) ) 
+         printerror( status );
+
+    /* move to the 3rd HDU in the input file (a binary table in this case) */
+    if ( fits_movabs_hdu(infptr, 3, &hdutype, &status) )
+         printerror( status );
+
+    if (hdutype != BINARY_TBL)  {
+        printf("Error: expected to find a binary table in this HDU\n");
+        return;
+    }
+    /* move to the last (2rd) HDU in the output file */
+    if ( fits_movabs_hdu(outfptr, 2, &hdutype, &status) )
+         printerror( status );
+
+    /* create new extension in the output file */
+    if ( fits_create_hdu(outfptr, &status) )
+         printerror( status );
+
+    /* get number of keywords */
+    if ( fits_get_hdrpos(infptr, &nkeys, &keypos, &status) ) 
+         printerror( status );
+
+    /* copy all the keywords from the input to the output extension */
+    for (ii = 1; ii <= nkeys; ii++)  {
+        fits_read_record (infptr, ii, card, &status); 
+        fits_write_record(outfptr,    card, &status); 
+    }
+
+    /* read the NAXIS1 and NAXIS2 keyword to get table size */
+    if (fits_read_keys_lng(infptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
+         printerror( status );
+
+    /* find which column contains the DENSITY values */
+    if ( fits_get_colnum(infptr, CASEINSEN, "density", &colnum, &status) )
+         printerror( status );
+
+    /* read the DENSITY column values */
+    frow = 1;
+    felem = 1;
+    nullval = -99.;
+    if (fits_read_col(infptr, TFLOAT, colnum, frow, felem, naxes[1], 
+        &nullval, density, &anynulls, &status) )
+        printerror( status );
+
+    /* allocate buffer large enough for 1 row of the table */
+    buffer = (unsigned char *) malloc(naxes[0]);
+
+    /* If the density is less than 3.0, copy the row to the output table */
+    for (noutrows = 0, irow = 1; irow <= naxes[1]; irow++)  {
+      if (density[irow - 1] < 3.0)  {
+        noutrows++;
+        fits_read_tblbytes( infptr, irow,      1, naxes[0], buffer, &status); 
+        fits_write_tblbytes(outfptr, noutrows, 1, naxes[0], buffer, &status); 
+    } }
+
+    /* update the NAXIS2 keyword with the correct number of rows */
+    if ( fits_update_key(outfptr, TLONG, "NAXIS2", &noutrows, 0, &status) )
+         printerror( status );
+
+    if (fits_close_file(outfptr, &status) || fits_close_file(infptr, &status))
+        printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void readheader ( void )
+
+    /**********************************************************************/
+    /* Print out all the header keywords in all extensions of a FITS file */
+    /**********************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+
+    int status, nkeys, keypos, hdutype, ii, jj;
+    char filename[]  = "atestfil.fit";     /* name of existing FITS file   */
+    char card[FLEN_CARD];   /* standard string lengths defined in fitsioc.h */
+
+    status = 0;
+
+    if ( fits_open_file(&fptr, filename, READONLY, &status) ) 
+         printerror( status );
+
+    /* attempt to move to next HDU, until we get an EOF error */
+    for (ii = 1; !(fits_movabs_hdu(fptr, ii, &hdutype, &status) ); ii++) 
+    {
+        /* get no. of keywords */
+        if (fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
+            printerror( status );
+
+        printf("Header listing for HDU #%d:\n", ii);
+        for (jj = 1; jj <= nkeys; jj++)  {
+            if ( fits_read_record(fptr, jj, card, &status) )
+                 printerror( status );
+
+            printf("%s\n", card); /* print the keyword card */
+        }
+        printf("END\n\n");  /* terminate listing with END */
+    }
+
+    if (status == END_OF_FILE)   /* status values are defined in fitsioc.h */
+        status = 0;              /* got the expected EOF error; reset = 0  */
+    else
+       printerror( status );     /* got an unexpected error                */
+
+    if ( fits_close_file(fptr, &status) )
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void readimage( void )
+
+    /************************************************************************/
+    /* Read a FITS image and determine the minimum and maximum pixel values */
+    /************************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status,  nfound, anynull;
+    long naxes[2], fpixel, nbuffer, npixels, ii;
+
+#define buffsize 1000
+    float datamin, datamax, nullval, buffer[buffsize];
+    char filename[]  = "atestfil.fit";     /* name of existing FITS file   */
+
+    status = 0;
+
+    if ( fits_open_file(&fptr, filename, READONLY, &status) )
+         printerror( status );
+
+    /* read the NAXIS1 and NAXIS2 keyword to get image size */
+    if ( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
+         printerror( status );
+
+    npixels  = naxes[0] * naxes[1];         /* number of pixels in the image */
+    fpixel   = 1;
+    nullval  = 0;                /* don't check for null values in the image */
+    datamin  = 1.0E30f;
+    datamax  = -1.0E30f;
+
+    while (npixels > 0)
+    {
+      nbuffer = npixels;
+      if (npixels > buffsize)
+        nbuffer = buffsize;     /* read as many pixels as will fit in buffer */
+
+      /* Note that even though the FITS images contains unsigned integer */
+      /* pixel values (or more accurately, signed integer pixels with    */
+      /* a bias of 32768),  this routine is reading the values into a    */
+      /* float array.   Cfitsio automatically performs the datatype      */
+      /* conversion in cases like this.                                  */
+
+      if ( fits_read_img(fptr, TFLOAT, fpixel, nbuffer, &nullval,
+                  buffer, &anynull, &status) )
+           printerror( status );
+
+      for (ii = 0; ii < nbuffer; ii++)  {
+        if ( buffer[ii] < datamin )
+            datamin = buffer[ii];
+
+        if ( buffer[ii] > datamax )
+            datamax = buffer[ii];
+      }
+      npixels -= nbuffer;    /* increment remaining number of pixels */
+      fpixel  += nbuffer;    /* next pixel to be read in image */
+    }
+
+    printf("\nMin and max image pixels =  %.0f, %.0f\n", datamin, datamax);
+
+    if ( fits_close_file(fptr, &status) )
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void readtable( void )
+
+    /************************************************************/
+    /* read and print data values from an ASCII or binary table */
+    /************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status, hdunum, hdutype,  nfound, anynull, ii;
+    long frow, felem, nelem, longnull, dia[6];
+    float floatnull, den[6];
+    char strnull[10], *name[6], *ttype[3]; 
+
+    char filename[]  = "atestfil.fit";     /* name of existing FITS file   */
+
+    status = 0;
+
+    if ( fits_open_file(&fptr, filename, READONLY, &status) )
+         printerror( status );
+
+    for (ii = 0; ii < 3; ii++)      /* allocate space for the column labels */
+        ttype[ii] = (char *) malloc(FLEN_VALUE);  /* max label length = 69 */
+
+    for (ii = 0; ii < 6; ii++)    /* allocate space for string column value */
+        name[ii] = (char *) malloc(10);   
+
+    for (hdunum = 2; hdunum <= 3; hdunum++) /*read ASCII, then binary table */
+    {
+      /* move to the HDU */
+      if ( fits_movabs_hdu(fptr, hdunum, &hdutype, &status) ) 
+           printerror( status );
+
+      if (hdutype == ASCII_TBL)
+          printf("\nReading ASCII table in HDU %d:\n",  hdunum);
+      else if (hdutype == BINARY_TBL)
+          printf("\nReading binary table in HDU %d:\n", hdunum);
+      else
+      {
+          printf("Error: this HDU is not an ASCII or binary table\n");
+          printerror( status );
+      }
+
+      /* read the column names from the TTYPEn keywords */
+      fits_read_keys_str(fptr, "TTYPE", 1, 3, ttype, &nfound, &status);
+
+      printf(" Row  %10s %10s %10s\n", ttype[0], ttype[1], ttype[2]);
+
+      frow      = 1;
+      felem     = 1;
+      nelem     = 6;
+      strcpy(strnull, " ");
+      longnull  = 0;
+      floatnull = 0.;
+
+      /*  read the columns */  
+      fits_read_col(fptr, TSTRING, 1, frow, felem, nelem, strnull,  name,
+                    &anynull, &status);
+      fits_read_col(fptr, TLONG, 2, frow, felem, nelem, &longnull,  dia,
+                    &anynull, &status);
+      fits_read_col(fptr, TFLOAT, 3, frow, felem, nelem, &floatnull, den,
+                    &anynull, &status);
+
+      for (ii = 0; ii < 6; ii++)
+        printf("%5d %10s %10ld %10.2f\n", ii + 1, name[ii], dia[ii], den[ii]);
+    }
+
+    for (ii = 0; ii < 3; ii++)      /* free the memory for the column labels */
+        free( ttype[ii] );
+
+    for (ii = 0; ii < 6; ii++)      /* free the memory for the string column */
+        free( name[ii] );
+
+    if ( fits_close_file(fptr, &status) ) 
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void printerror( int status)
+{
+    /*****************************************************/
+    /* Print out cfitsio error messages and exit program */
+    /*****************************************************/
+
+
+    if (status)
+    {
+       fits_report_error(stderr, status); /* print error report */
+
+       exit( status );    /* terminate the program, returning error status */
+    }
+    return;
+}

makemask/getMaxFits.c

+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*
+  protgram to compute the 
+*/
+#include <fitsio.h>
+
+
+
+int main(int argc, char **argv) {
+
+  fitsfile *inImage;      /* pointer to input image */
+
+
+  int status;
+  int ii;
+    
+  int maxdim;
+  long  nelements;
+  long naxes[10];
+  int bitpix = FLOAT_IMG ; /* output image is floats; BITPIX = -32*/
+    
+  float *imArray;
+    
+  int nullval ;
+  int anynull;
+  float sig;
+  
+  status = 0;
+  maxdim = 10;
+
+        
+  /* open the existing image */
+  fits_open_file(&inImage, argv[1], READONLY, &status) ;
+    
+  /* get size of image. assume x and y are axes 0 and 1 */  
+  fits_get_img_size(inImage, maxdim, naxes, &status);
+
+  /* number of pixels in image */
+  nelements = naxes[0]*naxes[1];
+  /* allocate memory for image */
+  imArray = (float *)malloc(nelements*sizeof(float));
+  /* set maimum to something low */
+  sig = -1100100101001.0;
+
+  /* don't check for null values */
+  nullval = 0;
+
+  /* read the image */
+  fits_read_img(inImage, TFLOAT, 1, nelements, &nullval,
+                imArray, &anynull, &status) ;
+
+  /* find maximum value */
+  for (ii = 0; ii < nelements; ii++) {
+    if (imArray[ii] > sig) {
+      sig = imArray[ii];
+    }
+  }
+
+  /* print maximum */
+  printf("%12.8f\n",sig);
+    
+  /* close all files */
+  fits_close_file(inImage,&status);
+
+
+
+
+}

makemask/locNoiseMed.c

+/*
+$Id: mirRenzo.c,v 1.6 2009/04/06 10:58:50 oosterloo Exp oosterloo $ 
+*/
+
+
+#include <port.h>
+#include <mirutil.h>
+#include <miriad.h>
+
+
+
+
+int main (int argc, char *argv[]) {
+
+  char  inpFilename[80];
+  char  inpFilename2[80];
+  char  outFilename[80];
+  char  outFilename2[80];
+  char  charBuff[80];
+  
+  int   xlen, ylen, zlen, n, m, i, j, k;
+  int   boxSize,  ind;
+  long int numP;
+  
+  float   maskVa, mean, sig;
+  float  tmpAr[1024];
+  
+  
+  float sia, value;  
+  long  index;  
+  float *sigs;
+  char *timeString;
+
+  
+  DataCube  *inpCube, *inpCube2, *outCube, *outCube2;
+  
+  
+  /* init command line */
+  keyini_c(argc, argv);
+  
+  /* get name of input file */
+  keya_c("in", inpFilename, "bla");
+  /* get name of input file */
+  keya_c("inima", inpFilename2, "bla");
+
+  /* get name of output file */
+  keya_c("out", outFilename, "bla");
+  /* get name of output file */
+  keya_c("outrat", outFilename2, "bla");
+
+  /* we have read all input */
+  keyfin_c();
+
+  inpCube = newDataCube();
+  
+  /* open output file */
+  openDataCube(inpFilename, inpCube, "old");
+  readDataCube(inpCube, FALSE);
+  xlen = inpCube->axes[0];
+  ylen = inpCube->axes[1];
+  zlen = inpCube->axes[2];
+
+  inpCube2 = newDataCube();
+  
+  /* open output file */
+  openDataCube(inpFilename2, inpCube2, "old");
+  readDataCube(inpCube2, FALSE);
+
+  boxSize=10*4/3;
+  
+  outCube = cloneDataCube(inpCube);
+  openDataCube(outFilename, outCube, "new");
+  outCube2 = cloneDataCube(inpCube);
+  openDataCube(outFilename2, outCube2, "new");
+  
+  outCube->mask = NULL;
+  outCube2->mask = NULL;
+  /* copy complete (?) header */
+  mirHeadCopy(inpCube->fileHandle, outCube->fileHandle);
+  mirHeadCopy(inpCube->fileHandle, outCube2->fileHandle);
+
+  numP = xlen*ylen*zlen;
+  
+  for (i = 0; i < numP; i++) {
+    outCube->data[i]= 0.0;
+    outCube2->data[i]= 0.0;
+  }
+
+
+  numP = (2*boxSize+1)*(2*boxSize+1);
+  for (i = 0; i < zlen; i++) {
+    for (j = boxSize; j < xlen-boxSize; j=j+5) {
+      for (k = boxSize; k < ylen-boxSize; k= k+5) {
+
+        ind = 0;
+        for (m = -boxSize; m <= boxSize; m++) {
+          for (n = -boxSize; n <= boxSize; n++) {
+            index = j+n + xlen*(k+m) + xlen*ylen*i;
+            tmpAr[ind] = fabs(inpCube->data[index]);
+            ind++;
+          }
+        }
+        mirRealSort(numP, tmpAr);
+        
+        sig = tmpAr[numP/2]*1.4826;
+	for (m = -2; m < 3; m++) {
+          for (n = -2; n < 3; n++) {
+            index = j+m + xlen*(k+n) + xlen*ylen*(i);
+            outCube->data[index]= sig;
+            if (sig > 0.0) {
+              outCube2->data[index]= inpCube2->data[index]/sig;
+            }
+	  }
+	}
+      }
+    }
+  }
+
+
+
+  writeDataCube(outCube);
+  writeDataCube(outCube2);
+  
+  /* sign off */
+  closeDataCube(inpCube);
+  closeDataCube(inpCube2);
+  closeDataCube(outCube);
+  closeDataCube(outCube2);
+  
+  return 0;
+  
+}
+

makemask/makeMaskFits.c

+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*
+  protgram to compute the 
+*/
+#include <fitsio.h>
+
+
+
+int main(int argc, char **argv) {
+
+  fitsfile *inImage;      /* pointer to input image */
+
+  fitsfile *outMask;      /* pointer to output mask */
+
+  int status, morekeys, hdutype;
+  int  i, j, ind;
+  int naxis;
+  int xlen, ylen;
+    
+  int maxdim;
+  long  nelements, index;
+  long naxes[10];
+  int bitpix = SHORT_IMG ; /* output image is ints; BITPIX = -16*/
+    
+  float *imArray;
+  int   *maskArray;
+     
+  float clipVal;
+  int nullval ;
+  int anynull;
+   
+  status = 0;
+  maxdim = 10;
+  clipVal = atof(argv[3]);
+ 
+          
+  /* open the existing image */
+  fits_open_file(&inImage, argv[1], READONLY, &status) ;
+
+    
+  /* get number of dimensions */
+  fits_get_img_dim(inImage, &naxis, &status);
+  /* get size of image. assume x and y are axes 0 and 1 */  
+  fits_get_img_size(inImage, maxdim, naxes, &status);
+
+  /* number of pixels in image */
+  nelements = naxes[0]*naxes[1];
+  xlen = naxes[0];
+  ylen = naxes[1];
+
+  /* allocare memeory for the arrays */
+  imArray = (float *)malloc(nelements*sizeof(float));
+  maskArray = (int *)malloc(nelements*sizeof(int));
+
+  /* don't check for null values */
+  nullval = 0;
+
+  /* read the image */
+  fits_read_img(inImage, TFLOAT, 1, nelements, &nullval,
+                imArray, &anynull, &status) ;
+
+  /* loop over image */
+  for (i = 0; i < xlen; i++) {
+    for (j = 0; j < ylen;  j++) {
+      index = i + xlen*(j) ;
+      if (imArray[index] > clipVal) {
+         maskArray[index]= 1;
+      } else {
+        maskArray[index]= 0;
+      }
+    }
+  }
+
+  /* delete output mask file */
+  remove(argv[2]);
+  /* create output mask file, copy header and write maskfits
+     file */
+  fits_create_file(&outMask, argv[2], &status);
+  fits_copy_hdu(inImage, outMask, 0, &status);
+  fits_write_img(outMask,TINT,1,nelements,maskArray,&status);
+
+  /* close all files */
+  fits_close_file(inImage,&status);
+  fits_close_file(outMask,&status);
+
+
+
+
+}

makemask/makeNoiseMapFits.c

+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+/*
+  protgram to compute the 
+*/
+#include <fitsio.h>
+
+int comp(const void *a,const void *b) {
+float *x = (float *) a;
+float *y = (float *) b;
+// return *x - *y; // this is WRONG...
+if (*x < *y) return -1;
+else if (*x > *y) return 1; return 0;
+}
+
+int main(int argc, char **argv) {
+
+  fitsfile *inImage;      /* pointer to input image */
+  fitsfile *inRes;        /* pointer to input residual */
+
+  fitsfile *outNoise;     /* pointer to output noise image */
+  fitsfile *outSN;        /* pointer to output S/N image */
+
+  int status, morekeys, hdutype;
+  int ii, jj,  numP, i, j, n, m, ind;
+  int naxis;
+  int xlen, ylen;
+    
+  int maxdim;
+  long  nelements, index;;
+  long naxes[10];
+  int bitpix = FLOAT_IMG ; /* output image is floats; BITPIX = -32*/
+    
+  float *imArray, *resArray, *noiseArray, *snArray;
+  float tmpAr[1024];
+    
+  int boxsize;
+  int nullval ;
+  int anynull;
+  float sig;
+  
+  status = 0;
+  maxdim = 10;
+  boxsize = 13;
+        
+  /* open the existing image */
+  fits_open_file(&inImage, argv[1], READONLY, &status) ;
+  /* open the residual image */
+  fits_open_file(&inRes, argv[2], READONLY, &status) ;
+    
+  /* get number of dimensions */
+  fits_get_img_dim(inImage, &naxis, &status);
+  /* get size of image. assume x and y are axes 0 and 1 */  
+  fits_get_img_size(inImage, maxdim, naxes, &status);
+
+  /* number of pixels in image */
+  nelements = naxes[0]*naxes[1];
+  xlen = naxes[0];
+  ylen = naxes[1];
+
+  /* allocare memeory for the arrays */
+  imArray = (float *)malloc(nelements*sizeof(float));
+  resArray = (float *)malloc(nelements*sizeof(float));
+  noiseArray = (float *)malloc(nelements*sizeof(float));
+  snArray = (float *)malloc(nelements*sizeof(float));
+
+  /* don't check for null values */
+  nullval = 0;
+
+  /* read the image */
+  fits_read_img(inImage, TFLOAT, 1, nelements, &nullval,
+                imArray, &anynull, &status) ;
+  /* read the residual */
+  fits_read_img(inRes, TFLOAT, 1, nelements, &nullval,
+               resArray, &anynull, &status) ;
+
+  /* set outpit noise array to zero */
+  for (ii = 0; ii < nelements; ii++) {
+    noiseArray[ii] = 0.0;
+  }
+
+  /* number of pixels in the box for which we compute the median */
+  numP = (2*boxsize+1)*(2*boxsize+1);
+  /* compute index that will point to the median */
+  numP = numP/2;
+
+  /* loop over image in steps of 5 pixels */
+  for (i = boxsize; i < xlen-boxsize; i=i+5) {
+    for (j = boxsize; j < ylen-boxsize; j= j+5) {
+
+      /* get the median for the current box */
+
+      /* reset index */
+      ind = 0;
+      /* loop over box, putting absolute value of the image value in temp array */
+      for (m = -boxsize; m <= boxsize; m++) {
+        for (n = -boxsize; n <= boxsize; n++) {
+          index = i+n + xlen*(j+m);
+          tmpAr[ind] = fabs(resArray[index]);
+          ind++;
+        }
+      }
+
+      /* sort the temp array and take the central value */
+      qsort(tmpAr,ind,sizeof(float),comp);
+      /* multiply by 1.4826 to turn the median absolute value into sigma */
+      sig = tmpAr[numP]*1.4826;
+      /* and put this value in the output noise image */
+      for (m = -2; m < 3; m++) {
+        for (n = -2; n < 3; n++) {
+          index = i+m + xlen*(j+n) ;
+          noiseArray[index]= sig;
+          /* if not zero, compute the S/N, i.e. the image value / noise value */
+          if (sig > 0.0) {
+            snArray[index]= imArray[index]/sig;
+          }
+
+        }
+      }
+    }
+  }
+
+  /* delete output noise file */
+  remove(argv[3]);
+  /* create output noise file, copy header and write noise values into fits
+     file */
+  fits_create_file(&outNoise, argv[3], &status);
+  fits_copy_hdu(inRes, outNoise, 0, &status);
+  fits_write_img(outNoise,TFLOAT,1,nelements,noiseArray,&status);
+
+  /* same as above, but now for S/N image */
+  remove(argv[4]);
+  
+  fits_create_file(&outSN, argv[4], &status);
+  fits_copy_hdu(inRes, outSN, 0, &status);
+  fits_write_img(outSN,TFLOAT,1,nelements,snArray,&status);
+
+  /* close all files */
+  fits_close_file(inImage,&status);
+  fits_close_file(inRes,&status);
+  fits_close_file(outNoise,&status);
+  fits_close_file(outSN,&status);
+
+
+
+}