vircam_stats.c

/* $Id: vircam_stats.c,v 1.16 2007/10/25 17:34:01 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: 2007/10/25 17:34:01 $
 * $Revision: 1.16 $
 * $Name:  $
 */

/* Includes */

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

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

#include "vircam_stats.h"
#include "vircam_utils.h"

/* Static subroutine prototypes */

static float kselect(float *a, int n, int k);
static double dkselect(double *a, int n, int k);
static float histexam(int *histo, int nhist, int level);

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

extern float vircam_med(float *data, unsigned char *bpm, long npts) {
    int i,j,is_even,ilevel;
    float *buf,value;

    /* If there is not BPM, then just do a straight forward median */

    buf = cpl_malloc(npts*sizeof(*buf));
    if (bpm == NULL) {
        is_even = !(npts & 1);
        memmove((char *)buf,(char *)data,npts*sizeof(float));
        if (is_even) {
            ilevel = npts/2 - 1;
            value = kselect(buf,npts,ilevel);
            ilevel = npts/2;
            value = 0.5*(value + kselect(buf,npts,ilevel));
        } else {
            ilevel = npts/2;
            value = kselect(buf,npts,ilevel);
        }

    /* Otherwise get rid of the dodgy values and then do the median */

    } else {
        j = 0;
        for (i = 0; i < npts; i++) {
            if (bpm[i] == 0)
                buf[j++] = data[i];
        }
        if (j == 0) {
            cpl_free(buf);
            value = CX_MAXFLOAT;
            return(value);
        }
        is_even = !(j & 1);
        if (is_even) {
            ilevel = j/2 - 1;
            value = kselect(buf,j,ilevel);
            ilevel = j/2;
            value = 0.5*(value + kselect(buf,j,ilevel));
        } else {
            ilevel = j/2;
            value = kselect(buf,j,ilevel);
        }
    }
    cpl_free(buf);
    return(value);
}

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

extern double vircam_dmed(double *data, unsigned char *bpm, long npts) {
    int i,j,is_even,ilevel;
    double *buf,value;

    /* If there is not BPM, then just do a straight forward median */

    buf = cpl_malloc(npts*sizeof(*buf));
    if (bpm == NULL) {
        is_even = !(npts & 1);
        memmove((char *)buf,(char *)data,npts*sizeof(double));
        if (is_even) {
            ilevel = npts/2 - 1;
            value = dkselect(buf,npts,ilevel);
            ilevel = npts/2;
            value = 0.5*(value + dkselect(buf,npts,ilevel));
        } else {
            ilevel = npts/2;
            value = dkselect(buf,npts,ilevel);
        }

    /* Otherwise get rid of the dodgy values and then do the median */

    } else {
        j = 0;
        for (i = 0; i < npts; i++) {
            if (bpm[i] == 0)
                buf[j++] = data[i];
        }
        if (j == 0) {
            cpl_free(buf);
            value = CX_MAXDOUBLE;
            return(value);
        }
        is_even = !(j & 1);
        if (is_even) {
            ilevel = j/2 - 1;
            value = dkselect(buf,j,ilevel);
            ilevel = j/2;
            value = 0.5*(value + dkselect(buf,j,ilevel));
        } else {
            ilevel = j/2;
            value = dkselect(buf,j,ilevel);
        }
    }
    cpl_free(buf);
    return(value);
}

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

extern float vircam_mean(float *data, unsigned char *bpm, long npts) {
    int i,n;
    float sum,value;

    /* Separate sections depending on whether there is a BPM or not */

    sum = 0.0;
    if (bpm == NULL) {
        n = npts;
        for (i = 0; i < npts; i++) 
            sum += data[i];
    } else {
        n = 0;
        for (i = 0; i < npts; i++) {
            if (bpm[i] == 0) {
                sum += data[i];
                n++;
            }
        }
    }
    if (n > 0)
        value = sum/(float)n;
    else
        value = CX_MAXFLOAT;
    return(value);
}

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

extern double vircam_dmean(double *data, unsigned char *bpm, long npts) {
    int i,n;
    double sum,value;

    /* Separate sections depending on whether there is a BPM or not */

    sum = 0.0;
    if (bpm == NULL) {
        n = npts;
        for (i = 0; i < npts; i++) 
            sum += data[i];
    } else {
        n = 0;
        for (i = 0; i < npts; i++) {
            if (bpm[i] == 0) {
                sum += data[i];
                n++;
            }
        }
    }
    if (n > 0)
        value = sum/(float)n;
    else
        value = CX_MAXDOUBLE;
    return(value);
}

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

extern int vircam_meansig(float *data, unsigned char *bpm, long npts, 
                          float *mean, float *sig) {
    int i,n;
    double sum,sum2,d;
    const char *fctid = "vircam_meansig";

    /* Separate sections depending on whether there is a BPM or not */

    sum = 0.0;
    sum2 = 0.0;
    if (bpm == NULL) {
        n = npts;
        for (i = 0; i < npts; i++) {
            d = (double)(data[i]);
            sum += d;
            sum2 += d*d;
        }
    } else {
        n = 0;
        for (i = 0; i < npts; i++) {
            if (bpm[i] == 0) {
                d = (double)(data[i]);
                sum += d;
                sum2 += d*d;
                n++;
            }
        }
    }

    /* Check whether we can do the mean and sigma calculations */

    switch (n) {
    case 0:
        *mean = CX_MAXFLOAT;
        *sig = CX_MAXFLOAT;
        cpl_msg_warning(fctid,"All values flagged as bad\n");
        return(VIR_WARN);
    case 1:
        *mean = (float)sum;
        *sig = 0.0;
        return(VIR_OK);
    default:
        sum /= (double)n;
        *mean = (float)sum;
        sum2 = sum2/(double)n - sum*sum;
        *sig = (float)sqrt(max(1.0e-12,sum2));
        return(VIR_OK);
    }
}

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

extern int vircam_meansigcut(float *data, unsigned char *bpm, long npts,
                             float lcut, float hcut, float *mean, float *sig) {
    int i,n;
    double sum,sum2;
    const char *fctid = "vircam_meansigcut";

    /* Separate sections depending on whether there is a BPM or not */

    sum = 0.0;
    sum2 = 0.0;
    if (bpm == NULL) {
        n = 0;
        for (i = 0; i < npts; i++) {
            if (data[i] > lcut && data[i] < hcut) {
                sum += data[i];
                sum2 += data[i]*data[i];
                n++;
            }
        }
    } else {
        n = 0;
        for (i = 0; i < npts; i++) {
            if (bpm[i] == 0 && data[i] > lcut && data[i] < hcut) {
                sum += data[i];
                sum2 += data[i]*data[i];
                n++;
            }
        }
    }

    /* Check whether we can do the mean and sigma calculations */

    switch (n) {
    case 0:
        *mean = CX_MAXFLOAT;
        *sig = CX_MAXFLOAT;
        cpl_msg_warning(fctid,"All values flagged as bad\n");
        return(VIR_WARN);
    case 1:
        *mean = (float)sum;
        *sig = 0.0;
        return(VIR_OK);
    default:
        sum /= (double)n;
        *mean = (float)sum;
        sum2 = sum2/(double)n - sum*sum;
        *sig = (float)sqrt(max(1.0e-12,sum2));
        return(VIR_OK);
    }
}

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

extern void vircam_qmedsig(float *data, unsigned char *bpm, long npts,
                           float thresh, int niter, float lowv, float highv,
                           float *median, float *sigma) {
    int *histo,nbins,nhist,ilev,iclip,nhist2,halflev,quartlev;
    int irej,jst,j;
    long i;
    float mlev,qlev;
    unsigned char *b;

    /* Right, first thing is to histogram the data.  Cut values below
       and above the 'ceiling' values */

    if (bpm == NULL) 
        b = cpl_calloc(npts,sizeof(unsigned char));
    else
        b = bpm;
    nbins = vircam_nint(highv - lowv + 1.0);
    histo = cpl_calloc(nbins,sizeof(*histo));
    nhist = 0;
    for (i = 0; i < npts; i++) {
        if (b[i] || data[i] < lowv || data[i] > highv)
            continue;
        ilev = vircam_nint(data[i] - lowv);
        ilev = max(0,min(nbins-1,ilev));
        histo[ilev] += 1;
        nhist += 1;
    }
    if (bpm == NULL)
        freespace(b);

    /* Right, find the median value and the first quartile. */

    iclip = nbins - 1;
    nhist2 = nhist;
    for (i = 0; i <= niter; i++) {
        halflev = (nhist2 + 1)/2;
        quartlev = (nhist2 + 3)/4;
        mlev = histexam(histo,nbins,halflev);
        *median = mlev + lowv;
        qlev = histexam(histo,nbins,quartlev);
        *sigma = (mlev - qlev)*1.48;
        if (i == niter)
            break;
        irej = 0;
        jst = vircam_nint(mlev + thresh*(*sigma));
        for (j = jst; j <= iclip; j++)
            irej += histo[j];
        if (irej == 0)
            break;
        iclip = jst - 1;
        nhist2 -= irej;
    }
    cpl_free(histo);
}

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

extern void vircam_medmad(float *data, unsigned char *bpm, long np, float *med,
                          float *mad) {
    int i;
    float *work;

    /* First find the median value */

    *med = vircam_med(data,bpm,np);

    /* Now work out the MAD. Start by getting a bit of workspace and filling
       it with absolute residuals */
    
    work = cpl_malloc(np*sizeof(*work));
    for (i = 0; i < np; i++)
        work[i] = (float)fabs((double)(data[i] - *med));

    /* Now get the median value of the absolute deviations */

    *mad = vircam_med(work,bpm,np);

    /* Tidy and exit */

    cpl_free(work);
}

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

extern void vircam_medmadcut(float *data, unsigned char *bpm, long np, 
                             float lcut, float hcut, float *med, float *mad) {
    int i;
    float *work;
    unsigned char *bad;

    /* Get a workspace for the pseudo-bad pixel mask... */

    bad = cpl_calloc(np,sizeof(*bad));
    if (bpm != NULL) {
        for (i = 0; i < np; i++)
            if (bpm[i] != 0 || data[i] < lcut || data[i] > hcut)
                bad[i] = 1;
    } else {
        for (i = 0; i < np; i++) 
            if (data[i] < lcut || data[i] > hcut)
                bad[i] = 1;
    }

    /* First find the median value */

    *med = vircam_med(data,bad,np);

    /* Now work out the MAD. Start by getting a bit of workspace and filling
       it with absolute residuals */
    
    work = cpl_malloc(np*sizeof(*work));
    for (i = 0; i < np; i++)
        work[i] = (float)fabs((double)(data[i] - *med));

    /* Now get the median value of the absolute deviations */

    *mad = vircam_med(work,bad,np);

    /* Tidy and exit */

    cpl_free(work);
    cpl_free(bad);
}

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

extern void vircam_medsig(float *data, unsigned char *bpm, long np, float *med,
                          float *sig) {
    int i,n;
    float sum,resid;

    /* First find the median value */

    *med = vircam_med(data,bpm,np);

    /* If no bpm is present the just use them all */

    if (bpm == NULL) {
        sum = 0.0;
        for (i = 0; i < np; i++) {
            resid = data[i] - *med;
            sum += resid*resid;
        }
        *sig = sqrt(sum/(float)np);

    /* Otherwise test the bpm */

    } else {
        sum = 0.0;
        n = 0;
        for (i = 0; i < np; i++) {
            if (bpm[i] == 0) {
                n++;
                resid = data[i] - *med;
                sum += resid*resid;
            }
        }
        if (n > 0) 
            *sig = sqrt(sum/(float)n);
        else
            *sig = 0.0;
    }
}

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

extern int vircam_sumbpm(unsigned char *bpm, int npts, int *sumb) {
    int j;

    *sumb = 0;
    for (j = 0; j < npts; j++)
        *sumb += bpm[j];
    return(VIR_OK);
}

static float histexam(int *histo, int nhist, int level) {
    int ilev,ii;
    float value;

    ii = 0;
    ilev = -1;
    while (ii < level && ilev < nhist-1)
        ii += histo[++ilev];
    value = (float)ilev - (float)(ii - level)/(float)histo[ilev] + 0.5;
    return(value);
}

/*
 * given an array a of n elements, return the element that would be at
 * position k, (0 <= k < n), if the array were sorted.  from Algorithms
 * and Data Structures in C++ by Leendert Ammeraal, pg. 82.  O(n).
 *
 * NB: partially reorders data in array
 */

/* Stolen from C. Sabbey */


static float kselect(float *a, int n, int k) {
    while (n > 1) {
        int i = 0, j = n - 1;
        float x = a[j/2], w;

        do {
            while (a[i] < x) i++;
            while (a[j] > x) j--;
            if (i < j) {
                w = a[i]; a[i] = a[j]; a[j] = w;
            } else {
                if (i == j) i++;
                break;
            }
        } while (++i <= --j);

        if (k < i)
            n = i;
        else {
            a += i; n -= i; k -= i;
        }
    }

    return a[0];
}

static double dkselect(double *a, int n, int k) {
    while (n > 1) {
        int i = 0, j = n - 1;
        double x = a[j/2], w;

        do {
            while (a[i] < x) i++;
            while (a[j] > x) j--;
            if (i < j) {
                w = a[i]; a[i] = a[j]; a[j] = w;
            } else {
                if (i == j) i++;
                break;
            }
        } while (++i <= --j);

        if (k < i)
            n = i;
        else {
            a += i; n -= i; k -= i;
        }
    }

    return a[0];
}

/*

$Log: vircam_stats.c,v $
Revision 1.16  2007/10/25 17:34:01  jim
Modified to remove lint warnings

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

Revision 1.14  2007/02/06 11:59:17  jim
Added vircam_medmadcut

Revision 1.13  2006/12/19 13:29:06  jim
Fixed medsig and medsigcut to take care of low number stats

Revision 1.12  2006/09/08 09:22:19  jim
Fixed vircam_qmedsig so that an input BPM could be NULL

Revision 1.11  2006/06/06 13:06:48  jim
Added vircam_medsig

Revision 1.10  2006/05/26 19:33:29  jim
Fixed a doc problem

Revision 1.9  2006/05/24 13:36:48  jim
Fixed warning message

Revision 1.8  2006/05/08 14:53:16  jim
Fixed problem negative variances by making sure that floats were cast to
doubles properly

Revision 1.6  2006/03/03 14:29:46  jim
Modified definition of vir_fits and channel table

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

Revision 1.4  2006/01/23 10:30:49  jim
Mainly documentation mods

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

Revision 1.2  2005/10/14 13:19:13  jim
Added vircam_medmad

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


*/

---