vircam_utils.c

/* $Id: vircam_utils.c,v 1.56 2008/07/10 13:04:03 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2005 Cambridge Astronomy Survey Unit
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: jim $
 * $Date: 2008/07/10 13:04:03 $
 * $Revision: 1.56 $
 * $Name:  $
 */

/* Includes */

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

#include <sys/time.h>
#include <time.h>
#include <libgen.h>
#include <string.h>
#include <unistd.h>

#include <cpl.h>
#include <math.h>

#include "vircam_utils.h"
#include "vircam_stats.h"
#include "vircam_fits.h"
#include "catalogue/imcore.h"

#define SZKEY 32
#define SZVAL 64

/* Define some columns for illumination correction tables */

#define NI_COLS 5
static const char *illcor_cols[NI_COLS] = {"xmin","xmax","ymin","ymax",
                                           "illcor"};

/* Static subroutine prototypes */

static float madfunc(int npts, float *xt, float *yt, float b);

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern const char *vircam_get_license(void) {
    const char  *vircam_license = 
        "This file is part of the VIRCAM Instrument Pipeline\n"
        "Copyright (C) 2006 Cambridge Astronomy Survey Unit\n"
        "\n"
        "This program is free software; you can redistribute it and/or modify\n"
        "it under the terms of the GNU General Public License as published by\n"
        "the Free Software Foundation; either version 2 of the License, or\n"
        "(at your option) any later version.\n"
        "\n"
        "This program is distributed in the hope that it will be useful,\n"
        "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
        "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n"
        "GNU General Public License for more details.\n"
        "\n"
        "You should have received a copy of the GNU General Public License\n"
        "along with this program; if not, write to the Free Software\n"
        "Foundation, Inc., 59 Temple Place, Suite 330, Boston, \n"
        "MA  02111-1307  USA";
    return(vircam_license);
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_compare_tags(const cpl_frame *frame1, 
                               const cpl_frame *frame2) {
    char *v1,*v2;

    /* Test entries */

    if (frame1 == NULL || frame2 == NULL) 
        return(-1);

    /* Get the tags */

    if ((v1 = (char *)cpl_frame_get_tag(frame1)) == NULL) 
        return(-1);
    if ((v2 = (char *)cpl_frame_get_tag(frame2)) == NULL) 
        return(-1);

    /* Compare the tags */

    if (strcmp(v1,v2)) 
        return(0);
    else 
        return(1);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern cpl_frameset *vircam_frameset_subgroup(cpl_frameset *frameset, 
                                              int *labels, int nlab, 
                                              const char *tag) {
    int i;
    cpl_frameset *cur_set,*ret_set;
    cpl_frame *cur_frame;
    char *cur_tag;

    ret_set = NULL;
    for (i = 0; i < nlab; i++) {
        cur_set = cpl_frameset_extract(frameset,labels,i);
        if (cur_set == NULL)
            break;
        cur_frame = cpl_frameset_get_frame(cur_set,0);
        cur_tag = (char *)cpl_frame_get_tag(cur_frame);
        if (!strcmp(cur_tag,tag)) {
            ret_set = cur_set;
            break;
        }
        cpl_frameset_delete(cur_set);
    }
    return(ret_set);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern cpl_frame *vircam_frameset_subgroup_1(cpl_frameset *frameset, 
                                             int *labels, int nlab, 
                                             const char *tag) {
    cpl_frameset *cur_set;
    cpl_frame *cur_frame,*new_frame;

    if ((cur_set = vircam_frameset_subgroup(frameset,labels,nlab,tag)) == NULL) {
        return(NULL);
    } else {
        cur_frame = cpl_frameset_get_frame(cur_set,0);
        new_frame = cpl_frame_duplicate(cur_frame);
        cpl_frameset_delete(cur_set);
        return(new_frame);
    }
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_exten_range(int inexten, const cpl_frame *fr, int *out1, 
                               int *out2) {
    int nvircam = 16,n,nmax;
    char *fctid = "vircam_exten_range";

    /* Right, how many frames actually exist? */

    n = cpl_frame_get_nextensions(fr);

    /* If there are less than the cannonical number, then issue a warning
       message here */

    if (n < nvircam) {
        cpl_msg_warning(fctid,"Only %d extensions out of %d are present",
                        n,nvircam);
        nmax = n;
    } else {
        nmax = nvircam;
    }

    /* If the extension number requested is zero, then do all the extensions
       that are available */
    
    if (inexten == 0) {
        *out1 = 1;
        *out2 = nmax;

    /* Otherwise, check that the requested extension is actually present */

    } else {

        if (inexten > nmax) {
            cpl_msg_error(fctid,"Requested extension %d is not present",
                          inexten);
            *out1 = -1;
            *out2 = -1;
        } else {
            *out1 = inexten;
            *out2 = inexten;
        }
    }
    return;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_madfit(int npts, float *xdata, float *ydata, 
                          float *intercept, float *slope) {
    int j;
    float sx,sy,sxx,sxy,det,aa,bb,temp,chisq,sigb,b1,f1,b2,f2,f;

    /* Do a linear least squares for a first estimate */

    sx = 0.0;
    sy = 0.0;
    sxx = 0.0;
    sxy = 0.0;
    for (j = 0; j < npts; j++) {
        sx += xdata[j];
        sy += ydata[j];
        sxx += xdata[j]*xdata[j];
        sxy += xdata[j]*ydata[j];
    }
    det = (float)npts*sxx - sx*sx;
    if (det == 0.0) {
        *slope = 0.0;
        *intercept = 0.0;
        return;
    }
    aa = (sxx*sy - sx*sxy)/det;
    bb = ((float)npts*sxy - sx*sy)/det;
    chisq = 0.0;
    for (j = 0; j < npts; j++) {
        temp = ydata[j] - (aa + bb*xdata[j]);
        chisq += temp*temp;
    }
    sigb = sqrt(chisq/det);
    if (sigb == 0.0) {
        *slope = bb;
        *intercept = aa;
        return;
    }

    /* Now bracket the solution */

    b1 = bb;
    f1 = madfunc(npts,xdata,ydata,b1);
    b2 = bb + ((f1 > 0.0) ? fabs(3.0*sigb) : -fabs(3.0*sigb));
    f2 = madfunc(npts,xdata,ydata,b2);
    while (f1*f2 > 0.0) {
        bb = 2.0*b2 - b1;
        b1 = b2;
        f1 = f2;
        b2 = bb;
        f2 = madfunc(npts,xdata,ydata,b2);
    }

    /* Iterate to find the minimum value */

    sigb = 0.01*sigb;
    while (fabs(b2 - b1) > sigb) {
        bb = 0.5*(b1 + b2);
        if (bb == b1 || bb == b2) 
            break;
        f = madfunc(npts,xdata,ydata,bb);
        if (f*f1 >= 0.0) {
            f1 = f;
            b1 = bb;
        } else {
            f2 = f;
            b2 = bb;
        }
    }
    *intercept = aa;
    *slope = bb;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/


static float madfunc(int npts, float *xt, float *yt, float b) {
    float *arr,aa,d,sum;
    int j;

    arr = cpl_malloc(npts*sizeof(*arr));
    for (j = 0; j < npts; j++) 
        arr[j] = yt[j] - b*xt[j];
    aa = vircam_med(arr,NULL,(long)npts);
    sum = 0.0;
    for (j = 0; j < npts; j++) {
        d = yt[j] - (b*xt[j] + aa);
        sum += d > 0.0 ? xt[j] : -xt[j];
    }
    cpl_free(arr);
    return(sum);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_linfit(int npts, double *xdata, double *ydata, 
                          double *intercept, double *slope, double *sig) {
    int j;
    double sx,sy,sxx,sxy,det,aa,bb,temp,sum,sumsq;

    /* Do a linear least squares */

    sx = 0.0;
    sy = 0.0;
    sxx = 0.0;
    sxy = 0.0;
    for (j = 0; j < npts; j++) {
        sx += xdata[j];
        sy += ydata[j];
        sxx += xdata[j]*xdata[j];
        sxy += xdata[j]*ydata[j];
    }
    det = (double)npts*sxx - sx*sx;
    if (det == 0.0) {
        *slope = 0.0;
        *intercept = 0.0;
        *sig = 0.0;
        return;
    }

    /* Intercept and slope */

    aa = (sxx*sy - sx*sxy)/det;
    bb = ((double)npts*sxy - sx*sy)/det;

    /* Work out RMS of fit */

    sum = 0.0;
    sumsq = 0.0;
    for (j = 0; j < npts; j++) {
        temp = ydata[j] - (aa + bb*xdata[j]);
        sum += temp;
        sumsq += temp*temp;
    }
    sum /= (double)npts;

    /* Set return values */

    *sig = sqrt(sumsq/(double)npts - sum*sum);
    *slope = bb;
    *intercept = aa;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_solve_gauss(double **a, double *b, int m) {
    double temp,big,pivot,rmax;
    int i,iu,j,k,jl,ib,ir;
    int l = 0;

    iu = m - 1;
    for (i = 0; i < iu; i++) {
        big = 0.0;

        /* find largest remaining term in ith column for pivot */

        for (k = i; k < m; k++) {
            rmax = fabs(a[i][k]);
            if (rmax > big) {
                big = rmax;
                l = k;
            }
        }

        /* check for non-zero term */

        if (big == 0.0) {
            for (ib = 0; ib < m; ib++)
                b[ib] = 0.0;
            cpl_msg_error("vircam_solve_gauss","Zero Determinant\n");
            return(VIR_FATAL);
        }

        if (i != l) {

            /* switch rows */

            for (j = 0; j < m; j++) {
                temp = a[j][i];
                a[j][i] = a[j][l];
                a[j][l] = temp;
            }
            temp = b[i];
            b[i] = b[l];
            b[l] = temp;
        }


        /* pivotal reduction */

        pivot = a[i][i];
        jl = i+1;

        for (j = jl; j < m; j++) {
            temp = a[i][j]/pivot;
            b[j] -= temp*b[i];
            for (k = i; k < m; k++)
                a[k][j] -= temp*a[k][i];
        }
    }

    /* back substitution for solution */

    for (i = 0; i < m; i++) {
        ir = m - 1 - i;
        if (a[ir][ir] != 0.0) {
            temp = b[ir];
            if (ir != m - 1) {
                for (j = 1; j <= i; j++) {
                    k = m - j;
                    temp -= a[k][ir]*b[k];
                }
            }
            b[ir] = temp/a[ir][ir];
        } else
            b[ir] = 0.0;
    }
    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_polyfit(const cpl_array *xarray, const cpl_array *yarray, 
                          int ncoefs, int ilim, int niter, float lclip, 
                          float hclip, cpl_array **polycf, double *sigfit) {
    const char *fctid = "vircam_polyfit";
    int npts,iter,i,j,nnew,k,retval,n;
    double *xdata,*ydata,*pdata,*res,**a,*b,temp,sum,sumsq,val;
    double lcut,hcut;
    unsigned char *pm;

    /* Initialise a few things */

    *polycf = NULL;
    *sigfit = -1.0;

    /* How many data points do we have? */

    npts = cpl_array_get_size(xarray);

    /* Do we have enough points for the required order of the fit? */

    if (npts < ncoefs) {
        cpl_msg_warning(fctid,"Not data for fit, Npts = %d, Ncoefs = %d",
                        npts,ncoefs);
        return(VIR_FATAL);
    }

    /* Create some arrays */

    a = cpl_malloc(ncoefs*sizeof(double *));
    b = cpl_calloc(ncoefs,sizeof(double));
    for (i = 0; i < ncoefs; i++) 
        a[i] = cpl_calloc(ncoefs,sizeof(double));
    pm = cpl_calloc(npts,sizeof(unsigned char));
    res = cpl_malloc(npts*sizeof(double));

    /* Get pointers to the input arrays */

    xdata = (double *)cpl_array_get_data_double_const(xarray);
    ydata = (double *)cpl_array_get_data_double_const(yarray);

    /* Get an output array */

    *polycf = cpl_array_new(ncoefs,CPL_TYPE_DOUBLE);
    pdata = cpl_array_get_data_double(*polycf);

    /* Iteration loop */

    for (iter = 0; iter < niter; iter++) {

        /* Zero some accumulators */

        for (i = 0; i < ncoefs; i++)  {
            for (j = 0; j < ncoefs; j++) 
                a[i][j] = 0.0;
            b[i] = 0.0;
        }
        nnew = 0;

        /* Cumulate sums */
 
        for (i = 0; i < npts; i++) {
            if (pm[i] == 1)
                continue;
            for (k = 0; k < ncoefs; k++) {
                temp = 1.0;
                if (k + ilim != 0)
                    temp = pow(xdata[i],(double)(k+ilim));
                b[k] += ydata[i]*temp;
                for (j = 0; j <= k; j++) {
                    temp = 1.0;
                    if (k + j + 2*ilim != 0)
                        temp = pow(xdata[i],(double)(k+j+2*ilim));
                    a[j][k] += temp;
                }
            }
        }
        for (k = 1; k < ncoefs; k++) 
            for (j = 0; j < k; j++) 
                a[k][j] = a[j][k];

        /* Solve linear equations */

        retval = vircam_solve_gauss(a,b,ncoefs);
        if (retval != VIR_OK) {
            cpl_msg_warning(fctid,"Fit failed");
            freearray(*polycf);
            freespace2(a,ncoefs);
            freespace(b);
            freespace(pm);
            freespace(res);
            return(VIR_FATAL);
        }

        /* Ok, assuming this is OK, then fill the polynomial coefficients */

        for (i = 0; i < ncoefs; i++)
            pdata[i] = b[i];

        /* Calculate the fit quality */

        sum = 0.0;
        sumsq = 0.0;
        n = 0;
        for (i = 0; i < npts; i++) {
            if (pm[i] == 1)
                continue;
            val = 0.0;
            for (j = ilim; j <= ncoefs; j++)
                val += pdata[j-ilim]*pow(xdata[i],(double)j);
            res[i] = val - ydata[i];
            sum += res[i];
            sumsq += pow(res[i],2.0);
            n++;
        }
        sum /= (double)n;
        *sigfit = sqrt(sumsq/(double)n - sum*sum);

        /* If this is not the last iteration, then do some clipping */

        lcut = sum - lcut*(*sigfit);
        hcut = sum + hcut*(*sigfit);
        if (iter < niter - 1) {
            for (i = 0; i < npts; i++) {
                if (pm[i] == 1) 
                    continue;
                if (res[i] > hcut || res[i] < lcut) {
                    nnew++;
                    pm[i] = 1;
                }
            }
        }
        
        /* If no new points have been clipped, then get out of here now... */
        
        if (nnew == 0) 
            break;
    }

    /* Tidy up and get out of here */

    freespace2(a,ncoefs);
    freespace(b);
    freespace(pm);
    freespace(res);
    return(VIR_OK);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_difference_image(cpl_image *master, cpl_image *prog, 
                                    unsigned char *bpm, cpl_table *chantab, 
                                    int ncells, int oper, float *global_diff, 
                                    float *global_rms, cpl_image **diffim, 
                                    cpl_table **diffimstats) {
    float *ddata,*work,mean,sig,med,mad;
    long nx,ny,npts;
    int nrows,i,nc1,nc2,nr,ixmin,ixmax,iymin,iymax,cnum,cx,cy,idx,idy;
    int icx,icy,indy1,indy2,indx1,indx2,jp,jcx,jj,jcy,ii,ncx,ncy;
    const char *fctid = "vircam_difference_image";

    /* Initialise the output */

    *diffim = NULL;
    *diffimstats = NULL;
    *global_diff = 0.0;
    *global_rms = 0.0;

    /* Is there a programme frame or a master? */

    if (prog == NULL || master == NULL)
        return;

    /* Start by subtracting the master image from the programme image */

    switch (oper) {
    case 1:
        *diffim = cpl_image_subtract_create(prog,master);
        break;
    case 2:
        *diffim = cpl_image_divide_create(prog,master);
        break;
    default:
        *diffim = NULL;
        cpl_msg_error(fctid,"Invalid operation requested %d",oper);
        break;
    }      
    if (*diffim == NULL)
        return;

    /* Work out a median difference */

    ddata = cpl_image_get_data_float(*diffim);
    nx = cpl_image_get_size_x(*diffim);
    ny = cpl_image_get_size_y(*diffim);
    npts = nx*ny;
    vircam_medmad(ddata,bpm,npts,global_diff,global_rms);

    /* Is there a channel table? */

    if (chantab == NULL)
        return;

    /* Before going any further, check that the channel table has all of 
       the columns we need */

    if (! cpl_table_has_column(chantab,"ixmin") ||
        ! cpl_table_has_column(chantab,"ixmax") ||
        ! cpl_table_has_column(chantab,"iymin") ||
        ! cpl_table_has_column(chantab,"iymax") ||
        ! cpl_table_has_column(chantab,"channum")) {
            cpl_msg_error(fctid,"Channel table is missing one of the required columns");

            return;
    }

    /* Work out how to divide the channels */

    switch (ncells) {
    case 1:
        nc1 = 1;
        nc2 = 1;
        break;
    case 2:
        nc1 = 2;
        nc2 = 1;
        break;
    case 4:
        nc1 = 4;
        nc2 = 1;
        break;
    case 8:
        nc1 = 8;
        nc2 = 1;
        break;
    case 16:
        nc1 = 16;
        nc2 = 1;
        break;
    case 32:
        nc1 = 16;
        nc2 = 2;
        break;
    case 64:
        nc1 = 32;
        nc2 = 2;
        break;
    default:
        nc1 = 32;
        nc2 = 2;
        break;
    }

    /* Create a difference image stats table */

    nrows = cpl_table_count_selected(chantab);
    *diffimstats = vircam_create_diffimg_stats(nrows*nc1*nc2);

    /* Loop for each data channel now */

    nr = 0;
    for (i = 0; i < nrows; i++) {
        ixmin = cpl_table_get_int(chantab,"ixmin",i,NULL);
        ixmax = cpl_table_get_int(chantab,"ixmax",i,NULL);
        iymin = cpl_table_get_int(chantab,"iymin",i,NULL);
        iymax = cpl_table_get_int(chantab,"iymax",i,NULL);
        cnum = cpl_table_get_int(chantab,"channum",i,NULL);

        /* Which is the long axis? If the channels are rectangular then
           divide the long axis by the greater number of cells. If the number
           of cells is less than 8, then just divide the long axis. If the
           channel is square, then divide the X axis more finely */

        cx = ixmax - ixmin + 1;
        cy = iymax - iymin + 1;
        if (cx > cy) {
            ncx = max(nc1,nc2);
            ncy = min(nc1,nc2);
        } else if (cx < cy) {
            ncy = max(nc1,nc2);
            ncx = min(nc1,nc2);
        } else {
            ncx = max(nc1,nc2);
            ncy = min(nc1,nc2);
        }
        
        /* How big is the cell? */

        idy = cy/ncy;
        idx = cx/ncx;
        work = cpl_malloc(idx*idy*sizeof(*work));

        /* Now loop for aach cell */

        for (icy = 0; icy < ncy; icy++) {
            indy1 = idy*icy;
            indy2 = min(iymax,indy1+idy-1);
            for (icx = 0; icx < ncx; icx++) {
                indx1 = idx*icx;
                indx2 = min(ixmax,indx1+idx-1);
                jp = 0;
                for (jcy = indy1; jcy < indy2; jcy++) {
                    jj = jcy*nx;
                    for (jcx = indx1; jcx < indx2; jcx++) {
                        ii = jj + jcx;
                        if (bpm != NULL && bpm[ii] == 0)
                            work[jp++] = ddata[ii];
                    }
                }
                (void)vircam_meansig(work,NULL,(long)jp,&mean,&sig);
                (void)vircam_medmad(work,NULL,(long)jp,&med,&mad);
                cpl_table_set_int(*diffimstats,"xmin",nr,indx1+1);
                cpl_table_set_int(*diffimstats,"xmax",nr,indx2+1);
                cpl_table_set_int(*diffimstats,"ymin",nr,indy1+1);
                cpl_table_set_int(*diffimstats,"ymax",nr,indy2+1);
                cpl_table_set_int(*diffimstats,"chan",nr,cnum);
                cpl_table_set_float(*diffimstats,"mean",nr,mean);
                cpl_table_set_float(*diffimstats,"median",nr,med);
                cpl_table_set_float(*diffimstats,"variance",nr,(sig*sig));
                cpl_table_set_float(*diffimstats,"mad",nr++,mad);
            }
        }
        cpl_free(work);
    }
}               

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

cpl_table *vircam_create_diffimg_stats(int nrows) {
    cpl_table *diffimstats;

    diffimstats = cpl_table_new(nrows);
    cpl_table_new_column(diffimstats,"xmin",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"xmin","pixels");
    cpl_table_new_column(diffimstats,"xmax",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"xmax","pixels");
    cpl_table_new_column(diffimstats,"ymin",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"ymin","pixels");
    cpl_table_new_column(diffimstats,"ymax",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"ymax","pixels");
    cpl_table_new_column(diffimstats,"chan",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"chan","pixels");
    cpl_table_new_column(diffimstats,"mean",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"mean","ADU");
    cpl_table_new_column(diffimstats,"median",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"median","ADU");
    cpl_table_new_column(diffimstats,"variance",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"variance","ADU**2");
    cpl_table_new_column(diffimstats,"mad",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"mad","ADU");
    return(diffimstats);
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_sort(float **a, int n, int m) {
    int iii,ii,i,ifin,j,it;
    float *t;

    t = cpl_malloc(m*sizeof(*t));

    iii = 2;
    while (iii < n)
        iii *= 2;
    iii = min(n,(3*iii)/4 - 1);

    while (iii > 1) {
        iii /= 2;
        ifin = n - iii;
        for (ii = 0; ii < ifin; ii++) {
            i = ii;
            j = i + iii;
            if (a[0][i] > a[0][j]) {
                for (it = 0; it < m; it++)
                   t[it] = a[it][j];
                while (1) {
                    for (it = 0; it < m; it++) 
                        a[it][j] = a[it][i];
                    j = i;
                    i = i - iii;
                    if (i < 0 || a[0][i] <= t[0]) 
                        break;
                }
                for (it = 0; it < m; it++) 
                    a[it][j] = t[it];
            }
        }
    }
    cpl_free(t);
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern long vircam_getnpts(cpl_image *in) {
    int nx,ny;
    long npts;
    const char *fctid = "vircam_getnpts";

    if ((nx = cpl_image_get_size_x(in)) == -1) {
        cpl_msg_error(fctid,"NULL image input");
        return(0);
    }
    if ((ny = cpl_image_get_size_y(in)) == -1) {
        cpl_msg_error(fctid,"NULL image input");
        return(0);
    }
    npts = (long)nx*ny;
    return(npts);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_fndmatch(float x, float y, float *xlist, float *ylist, 
                           int nlist, float err) {
    int isp,ifp,index,i;
    float errsq,errmin,dx,dy,poserr;

    /* Find lower limit index */

    isp = 0;
    ifp = nlist - 1;
    errsq = err*err;
    index = (isp + ifp)/2;
    while (ifp-isp >= 2) {
        if (ylist[index] < y - err) {
            isp = index;
            index = (index+ifp)/2;
        } else if (ylist[index] > y - err) {
            ifp = index;
            index = (index+isp)/2;
        } else {
            isp = index;
            break;
        }
    }

    /* Now find nearest one within limit */

    index = -1;
    errmin = errsq;
    for (i = isp; i < nlist; i++) {
        if (ylist[i] > y+err)
            break;
        dx = x - xlist[i];
        dy = y - ylist[i];
        poserr = dx*dx + dy*dy;
        if (poserr < errsq) {
            if (poserr <= errmin) {
                index = i;
                errmin = poserr;
            }
        }
    }
    return(index);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int *vircam_dummy_confidence(long n) {
    int *cdata,i;
 
    cdata = cpl_malloc(n*sizeof(*cdata));
    for (i = 0; i < n; i++)
        cdata[i] = 100;
    return(cdata);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_compare_dims(cpl_image *im1, cpl_image *im2) {
    
    if ((cpl_image_get_size_x(im1) != cpl_image_get_size_x(im2)) ||
        cpl_image_get_size_y(im1) != cpl_image_get_size_y(im2))
        return(VIR_FATAL);
    else
        return(VIR_OK);
}
   
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_prov(cpl_propertylist *p, vir_fits **inlist, int n) {
    int i;
    char keyword[SZKEY],value[SZVAL],*fn,*base;

    /* Delete all the provenance keywords that might already exist */

    cpl_propertylist_erase_regexp(p,"ESO DRS PROV*",0);

    /* Add the new provenance keywords */

    for (i = 0; i < n; i++) {
        (void)snprintf(keyword,SZKEY,"ESO DRS PROV%04d",i+1);
        fn = cpl_strdup(vircam_fits_get_fullname(inlist[i]));
        base = basename(fn);
        (void)snprintf(value,SZVAL,"%s",base);
        cpl_free(fn);
        cpl_propertylist_update_string(p,keyword,value);
        (void)snprintf(value,SZVAL,"Input file # %d",i+1);
        cpl_propertylist_set_comment(p,keyword,value);
    }
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_merge_propertylists(cpl_propertylist *p1, 
                                       cpl_propertylist *p2) {
    int i;
    
    /* Return if either propertylist is NULL */

    if (p1 == NULL || p2 == NULL)
        return;

    /* Copy each property from the second list into the first one */

    for (i = 0; i < cpl_propertylist_get_size(p2); i++)
        cpl_propertylist_copy_property(p1,p2,cpl_property_get_name(cpl_propertylist_get(p2,i)));
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_dummy_property(cpl_propertylist *p) {

    /* Check for silly input */

    if (p == NULL)
        return;

    /* Add the property now */

    cpl_propertylist_update_bool(p,"ESO DRS IMADUMMY",TRUE);
    cpl_propertylist_set_comment(p,"ESO DRS IMADUMMY",
                                 "This is a dummy product");
    return;
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_is_dummy(cpl_propertylist *p) {

    /* Check for silly input */

    if (p == NULL)
        return(0);

    /* Check the propertylist and return the result */

    return(cpl_propertylist_has(p,"ESO DRS IMADUMMY"));
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_overexp(vir_fits **fitslist, int *n, float lthr, 
                           float hthr, int ditch) {
    int i,m;
    cpl_image *im;
    double val;

    /* Loop for each of the fits items */

    m = 0;
    for (i = 0; i < *n; i++) {
        im = vircam_fits_get_image(fitslist[i]);
        val = cpl_image_get_mean_window(im,100,100,200,200);
        if (val > lthr && val < hthr)
            fitslist[m++] = fitslist[i];
        else
            if (ditch) 
                vircam_fits_delete(fitslist[i]);
    }
    for (i = m; i < *n; i++)
        fitslist[i] = NULL;
    *n = m;
    
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern cpl_image *vircam_dummy_image(vir_fits *model) {
    cpl_image *im;

    /* Copy the input model image */

    im = cpl_image_duplicate(vircam_fits_get_image(model));

    /* Now change it all to zeros */

    cpl_image_multiply_scalar(im,0.0);

    /* Return the result */

    return(im);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern cpl_table *vircam_dummy_catalogue(int type) {

    cattype = type;
    tabinit(NULL);
    return(tab);
}


/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern cpl_table *vircam_illcor_newtab(int nrows) {
    cpl_table *illcor;
    int i;

    /* Create the table structure and fill in the columns */

    illcor = cpl_table_new(nrows);
    for (i = 0; i < NI_COLS; i++) 
        cpl_table_new_column(illcor,illcor_cols[i],CPL_TYPE_FLOAT);
    return(illcor);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_timestamp(char *out, int n) {
    struct timeval tv;
    struct tm *tm;
    float sec;
    
    /* Get the Greenwich Mean Time */

    (void)gettimeofday(&tv,NULL);
    tm = gmtime(&(tv.tv_sec));
    sec = (float)tm->tm_sec + 1.0e-6*(float)tv.tv_usec;

    /* Now format it */

    (void)snprintf(out,n,"%04d-%02d-%02dT%02d:%02d:%07.4f",1900+tm->tm_year,
                   tm->tm_mon+1,tm->tm_mday,tm->tm_hour,tm->tm_min,sec);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_backmap(vir_fits *in, vir_mask *mask, int nbsize, 
                           cpl_image **out, float *med) {
    int i,j,nx,ny,ifracx,ifracy,nbsizx,nbsizy,nbx,nby,*nps,jx,jy;
    int jy1,jy2,np,nout,jyp1,jxp1,jz1,jz2,jz3,jz4,nbsize2;
    float fracx,fracy,*bmap,**rowpoints,*data,*ptr,dely,delx,t1,t2;
    unsigned char *bpm;

    /* Check to see if nbsize is close to exact divisor */

    nx = cpl_image_get_size_x(vircam_fits_get_image(in));
    ny = cpl_image_get_size_y(vircam_fits_get_image(in));
    fracx = ((float)nx)/((float)nbsize);
    fracy = ((float)ny)/((float)nbsize);
    ifracx = (int)(fracx + 0.1);
    ifracy = (int)(fracy + 0.1);
    nbsizx = nx/ifracx;
    nbsizy = ny/ifracy;
    nbsize = max(vircam_nint(0.9*nbsize),min(nbsize,min(nbsizx,nbsizy)));
    nbsize = min(nx,min(ny,nbsize)); /* trap for small maps */
    nbsize2 = nbsize/2;

    /* Divide the map into partitions */

    nbx = nx/nbsize;
    nby = ny/nbsize;

    /* Get some space for to use for accumulating stats */

    bmap = cpl_malloc(nbx*nby*sizeof(float));
    rowpoints = cpl_malloc(nbx*sizeof(float *));
    nps = cpl_malloc(nbx*sizeof(int));

    /* Get the data from the input file and the mask */

    data = cpl_image_get_data_float(vircam_fits_get_image(in));
    bpm = vircam_mask_get_data(mask);

    /* Work out the global median */

    *med = vircam_med(data,bpm,(long)(nx*ny));

    /* For each cell allocate some nbsize*nbsize pixels in the row pointers.
       It's ok to do this here since this will be the maximum that each cell
       will have... */

    for (i = 0; i < nbx; i++) 
        rowpoints[i] = cpl_malloc(nbsize*nbsize*sizeof(float));
    
    /* Loop for each row of cells and work out which rows in the input
       map these relate to */

    for (jy = 0; jy < nby; jy++) {
        jy1 = jy*nbsize;
        jy2 = min((jy+1)*nbsize - 1,ny-1);

        /* Zero the counter for each cell */

        for (jx = 0; jx < nbx; jx++)
            nps[jx] = 0;

        /* Loop for the strip of the input map covered by this row of cells.
           Work out which cell the current pixel is in and then put it into
           the correct accumulator assuming the bad pixel mask says it's ok */

        for (j = jy1; j <= jy2; j++) {
            for (i = 0; i < nx; i++) {
                jx = min(nbx-1,i/nbsize);
                if (bpm[j*nx + i] == 0) {
                    ptr = rowpoints[jx];
                    np = nps[jx];
                    ptr[np++] = data[j*nx + i];
                    nps[jx] = np;
                }
            }
        }

        /* Now do the medians for each of the cells in this row */

        for (jx = 0; jx < nbx; jx++) 
            bmap[jy*nbx+jx] = vircam_med(rowpoints[jx],NULL,(long)(nps[jx])) - *med;
    }

    /* Free up some workspace */

    for (jx = 0; jx < nbx; jx++) 
        freespace(rowpoints[jx]);
    freespace(rowpoints);
    freespace(nps);

    /* Create a new image for the output array */

    *out = cpl_image_new(nx,ny,CPL_TYPE_FLOAT);
    ptr = cpl_image_get_data_float(*out);

    /* Ok, we now have the block averaged background map. Do a two point
       interpolation to create the full map now */

    nout = 0;
    for (j = 1; j <= ny; j++) {
        jy = (j + nbsize2)/nbsize;
        jyp1 = jy + 1;
        jy = min(nby,max(1,jy));
        jyp1 = min(nby,jyp1);
        dely = (float)(j - nbsize*jy + nbsize2)/(float)nbsize;
        dely = max(0.0,min(1.0,dely));
        for (i = 1; i <= nx; i++) {
            jx = (i + nbsize2)/nbsize;
            jxp1 = jx + 1;
            jx = min(nbx,max(1,jx));
            jxp1 = min(nbx,jxp1);
            delx = (float)(i - nbsize*jx + nbsize2)/(float)nbsize;
            delx = max(0.0,min(1.0,delx));
            jz1 = (jy - 1)*nbx + jx;
            jz2 = (jyp1 - 1)*nbx + jx;
            if (jx == nbx) {
                jz3 = jz1;
                jz4 = jz2;
            } else {
                jz3 = jz1 + 1;
                jz4 = jz2 + 1;
            }
            t1 = (1.0 - delx)*bmap[jz1-1] + delx*bmap[jz3-1];
            t2 = (1.0 - delx)*bmap[jz2-1] + delx*bmap[jz4-1];
            ptr[nout++] = (1.0 - dely)*t1 + dely*t2;
        }
    }

    /* Tidy up before leaving... */

    freespace(bmap)
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_findcol(cpl_propertylist *p, char *col) {

    if (!strcmp(col,"X")) {
        if (cpl_propertylist_has(p,"ESO DRS XCOL"))
            return(cpl_propertylist_get_int(p,"ESO DRS XCOL"));
        else 
            return(-1);
    } else if (!strcmp(col,"Y")) {
        if (cpl_propertylist_has(p,"ESO DRS YCOL"))
            return(cpl_propertylist_get_int(p,"ESO DRS YCOL"));
        else 
            return(-1);
    }
    return(-1);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern void vircam_rename_property(cpl_propertylist *p, char *oldname,
                                   char *newname) {
    cpl_propertylist *temp;
    cpl_property *property;

    /* First get the old property. Note that we have to do this in a 
       particularly silly way since you cannot reference an individual property
       in a propertylist by its name. You can only do it by its index 
       number. Remeber to change this when CPL comes to it's senses... */

    if (! cpl_propertylist_has(p,oldname))
        return;
    temp = cpl_propertylist_new();
    cpl_propertylist_copy_property(temp,p,oldname);
    property = cpl_propertylist_get(temp,0);

    /* Now change its name */

    cpl_property_set_name(property,newname);

    /* Now insert this into the propertylist and delete the old one */

    cpl_propertylist_append(p,temp);
    cpl_propertylist_erase(p,oldname);
    cpl_propertylist_delete(temp);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_catpars(cpl_frame *index, char **catpath, char **catname) {
    cpl_propertylist *p;
    char *unk = "unknown",*fname;
    int status;
    const char *fctid = "vircam_catpars";

    /* Initialise a few things */

    *catpath = NULL;
    *catname = NULL;

    /* First get the full path to the index file and make sure it exists */

    fname = cpl_strdup(cpl_frame_get_filename(index));
    if (access((const char *)fname,R_OK) != 0) {
        cpl_msg_error(fctid,"Can't access index file %s",fname);
        cpl_free(fname);
        return(VIR_FATAL);
    }
    *catpath = cpl_strdup(dirname(fname));

    /* Load the propertylist if you can. If you can't then signal a fatal
       error since this probably means the whole file is messed up */

    if ((p = cpl_propertylist_load(cpl_frame_get_filename(index),0)) == NULL) {
        freespace(*catpath);
        cpl_msg_error(fctid,"Can't load index file header %s",fname);
        cpl_free(fname);
        return(VIR_FATAL);
    }
        
    /* If there is a catalogue name in the header then send it back. If there
       isn't then give a default name and send a warning */

    if (cpl_propertylist_has(p,"CATNAME")) {
        *catname = cpl_strdup(cpl_propertylist_get_string(p,"CATNAME"));
        status = VIR_OK;
    } else {
        *catname = cpl_strdup(unk);
        cpl_msg_warn(fctid,"Property CATNAME not in index file header %s",
                     fname);
        status = VIR_WARN;
    }
    cpl_free(fname);
    freepropertylist(p);

    /* Get out of here */

    return(status);
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_gaincor_calc(cpl_frame *frame, int *n, float **cors,
                               int *status) {
    float sum,val;
    int i,ngood;
    unsigned char *iflag;
    cpl_propertylist *p;

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);

    /* Find the number of extensions in the file and allocate some workspace
       to hold the results. Allocate a small workspace to flag dummy 
       extensions */

    *n = cpl_frame_get_nextensions(frame);
    *cors = cpl_malloc(*n*sizeof(float));
    iflag = cpl_calloc(*n,sizeof(iflag));

    /* Ok, loop through the extensions and read the propertylists */

    sum = 0.0;
    ngood = 0;
    for (i = 0; i < *n; i++) {
        p = cpl_propertylist_load(cpl_frame_get_filename(frame),i+1);
        if (cpl_propertylist_has(p,"ESO DRS IMADUMMY")) {
            iflag[i] = 1;
        } else if (! cpl_propertylist_has(p,"ESO DRS MEDFLAT")) {
            iflag[i] = 1;
        } else {
            val = cpl_propertylist_get_double(p,"ESO DRS MEDFLAT");
            if (val == 0.0) {
                iflag[i] = 1;
            } else {
                sum += val;
                (*cors)[i] = val;
                ngood++;
            }
        }
        cpl_propertylist_delete(p);
    }

    /* If any of them are good, then work out what the average is and
       create the correction factors. If the image was a dummy or
       there was no MEDFLAT keyword in the header, then just give it
       a factor of 1 */

    if (ngood > 0)
        sum /= (float)ngood;
    for (i = 0; i < *n; i++) {
        if (iflag[i] == 0) {
            (*cors)[i] = sum/(*cors)[i];
        } else {
            (*cors)[i] = 1.0;
        }
    }
    cpl_free(iflag);

    /* Get out of here */

    GOOD_STATUS
}

/*

$Log: vircam_utils.c,v $
Revision 1.56  2008/07/10 13:04:03  jim
Fixed some comments

Revision 1.55  2008/05/06 12:15:02  jim
Changed vircam_catpars to check that we can read the index file and to send
out error messages if there are problems

Revision 1.54  2008/05/06 08:38:29  jim
Modified call to cpl_propertylist_get_ from float to double in
vircam_gaincor_calc.

Revision 1.53  2007/11/14 14:47:53  jim
vircam_linfit now works only with doubles

Revision 1.52  2007/11/14 10:46:07  jim
Fixed bugs in vircam_polyfit

Revision 1.51  2007/10/25 17:33:31  jim
Modified to add vircam_polyfit and to remove lint warnings

Revision 1.50  2007/10/19 09:25:10  jim
Fixed problems with missing includes

Revision 1.49  2007/10/19 06:55:06  jim
Modifications made to use new method for directing the recipes to the
standard catalogues using the sof

Revision 1.48  2007/10/15 12:50:28  jim
Modified for compatibility with cpl_4.0

Revision 1.47  2007/07/09 13:21:06  jim
Changed argument list to vircam_exten_range to include a cpl_frame

Revision 1.46  2007/05/08 10:41:01  jim
Added vircam_gaincor_calc

Revision 1.45  2007/05/03 11:15:33  jim
Fixed little problem with table wcs

Revision 1.44  2007/05/02 09:14:48  jim
Added vircam_findcol and vircam_rename_property

Revision 1.43  2007/03/01 12:42:42  jim
Modified slightly after code checking

Revision 1.42  2007/02/14 12:53:34  jim
Removed vircam_paf_print_header

Revision 1.41  2007/02/07 10:12:15  jim
Removed vircam_ndit_correct as it is now no longer necessary

Revision 1.40  2007/01/15 12:36:27  jim
Fixed bug in nditcor where factor was being divided by instead of multiplied

Revision 1.39  2006/12/06 12:58:41  jim
Fixed a small bug in vircam_backmap

Revision 1.38  2006/12/01 14:08:33  jim
added vircam_backmap

Revision 1.37  2006/11/28 20:52:17  jim
Added vircam_illcor_newtab

Revision 1.36  2006/11/27 12:00:17  jim
cosmetic changes

Revision 1.35  2006/10/05 09:23:00  jim
Small modifications to a couple of cpl calls to bring them into line with
cpl v3.0

Revision 1.34  2006/08/27 20:28:15  jim
added vircam_is_dummy

Revision 1.33  2006/08/21 09:09:29  jim
Added routines vircam_create_diffimg_stats and vircam_dummy_property

Revision 1.32  2006/07/11 14:53:58  jim
Fixed vircam_ndit_correct so that images with bad status are not corrected

Revision 1.31  2006/07/04 09:19:06  jim
replaced all sprintf statements with snprintf

Revision 1.30  2006/06/30 21:32:09  jim
Fixed bug in vircam_prov so that sprintf is replaced by snprintf

Revision 1.29  2006/06/20 18:59:51  jim
Fixed small bug in vircam_ndit_correct

Revision 1.28  2006/06/09 22:24:47  jim
Added vircam_ndit_correct

Revision 1.27  2006/06/09 11:26:26  jim
Small changes to keep lint happy

Revision 1.26  2006/06/06 13:07:54  jim
Change the stats routine used in vircam_difference_image

Revision 1.25  2006/05/30 13:31:47  jim
Added vircam_timestamp

Revision 1.24  2006/05/24 13:35:49  jim
Added vircam_dummy_image and vircam_dummy_catalogue

Revision 1.23  2006/05/04 15:16:33  jim
Fixed memory bug in vircam_overexp

Revision 1.22  2006/04/20 11:29:41  jim
Added vircam_overexp

Revision 1.21  2006/03/22 11:42:32  jim
Moved some routines into vircam_mask

Revision 1.20  2006/03/15 10:43:42  jim
Fixed a few things

Revision 1.19  2006/03/08 14:32:22  jim
Lots of little modifications

Revision 1.18  2006/03/01 10:31:29  jim
Now uses new vir_fits objects

Revision 1.17  2006/02/18 11:48:55  jim
*** empty log message ***

Revision 1.16  2006/01/23 10:30:50  jim
Mainly documentation mods

Revision 1.15  2005/12/14 22:17:33  jim
Updated docs

Revision 1.14  2005/12/01 16:26:03  jim
Added vircam_dummy_confidence

Revision 1.13  2005/11/25 15:33:22  jim
Some code fixes to keep splint happy

Revision 1.12  2005/11/25 09:56:15  jim
Tidied up some more documentation

Revision 1.11  2005/11/25 09:36:23  jim
Added vircam_linfit

Revision 1.10  2005/11/07 13:15:16  jim
Fixed lots of bugs and added some error checking

Revision 1.9  2005/11/03 15:16:28  jim
Lots of changes mainly to strengthen error reporting

Revision 1.8  2005/11/03 13:28:50  jim
All sorts of changes to tighten up error handling

Revision 1.7  2005/10/14 13:21:04  jim
Added some new routines

Revision 1.6  2005/09/20 15:07:46  jim
Fixed a few bugs and added a few things

Revision 1.5  2005/09/07 12:47:22  jim
renamed a malloc and free to the cpl equivallent

Revision 1.4  2005/08/10 09:55:05  jim
Modified vircam_madfit so that if determinant is zero, then it sends a
zero slope and intercept back

Revision 1.3  2005/08/10 08:42:00  jim
Modified vircam_madfit so that if the initial least squares fit gets a
perfect result, then it just returns with those values rather than to try
and do the median fit

Revision 1.2  2005/08/09 15:30:00  jim
vircam_exten_range had number of chips wrong!

Revision 1.1.1.1  2005/08/05 08:29:09  jim
Initial import


*/

---