operaStats.h

Go to the documentation of this file.

#ifndef LIBOPERASTATS_H
#define LIBOPERASTATS_H
/*******************************************************************
 ****                LIBRARY FOR OPERA v1.0                     ****
 *******************************************************************
 Library name: operaStats
 Version: 1.0
 Description: ThisC library implements basic statistics routines..
 Author(s): CFHT OPERA team
 Affiliation: Canada France Hawaii Telescope 
 Location: Hawaii USA
 Date: Aug/2011
 Contact: eder@cfht.hawaii.edu
 
 Copyright (C) 2011  Opera Pipeline team, Canada France Hawaii Telescope
 
 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 3 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, see:
 http://software.cfht.hawaii.edu/licenses
 -or-
 http://www.gnu.org/licenses/gpl-3.0.html
 ********************************************************************/

// $Date$
// $Id$
// $Revision$
// $Locker$
// $Log$

/*
 * Image routines.
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "operaError.h"

// for Linux...
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif

#define ELEM_SWAP(a,b) {register float t=(a);(a)=(b);(b)=t;}

float operaArrayMean(unsigned np, const float *array);
float operaArrayWeightedMean(unsigned np, const float *array, const float *weigh);
float operaArrayAvgSigClip(unsigned np, const float *array, unsigned nsig);
float operaArraySig(unsigned np, const float *array);
float operaArrayWeightedSig(unsigned np, const float *array, const float *weigh);
float operaArrayMedian(unsigned np, const float *array);
float operaArrayMedSig(unsigned np, const float *array, float median);
float operaArrayMaxValue(unsigned np, const float *array);
float operaArrayMinValue(unsigned np, const float *array);
void  operaArrayHeapSort(unsigned np, float *arr);
void operaArrayIndexSort(unsigned np, const float *array, unsigned *sindex);
float operaUniformRand(float xcen, float xmax, float xmin);
float operaGaussRand(float xcen, float sig);
unsigned operaCountPixels(unsigned np, const float *array, float minvalue, float maxvalue);


inline void operaArrayIndexedMeanQuick(unsigned np, const float *array, const float *indexmask, unsigned nb, float *meanbin) {
        
        float *sumbin = (float *)malloc(sizeof(float)*nb);
        
        memset(sumbin, 0, sizeof(float)*nb);
        memset(meanbin, 0, sizeof(float)*nb);
        
        while (np--){
                if(*indexmask) {
                        meanbin[(unsigned)*indexmask-1] += *array;
                        sumbin[(unsigned)*indexmask-1]++;
                }
                indexmask++;
                array++;
        }
        while(nb--) {
                *meanbin++ /= *sumbin++;
        }       
}

inline void operaArrayIndexedSigQuick(unsigned np, const float *array, const float *indexmask, unsigned nb, float *sigbin) {
        
        float *sumbin = (float *)malloc(sizeof(float)*nb);
        float *meanbin = (float *)malloc(sizeof(float)*nb);     
        unsigned i=np;
        unsigned j=0;
        
        memset(sumbin, 0, sizeof(float)*nb);
        memset(meanbin, 0, sizeof(float)*nb);
        memset(sigbin, 0, sizeof(float)*nb);
        
        if(np){ 
                
                while (i--){
                        if(*indexmask) {
                                j = (unsigned)*indexmask-1;
                                meanbin[j] += *array;
                                sumbin[j]++;
                        }
                        indexmask++;
                        array++;
                }
                
                j=nb;
                while(j--) {
                        meanbin[j] /= sumbin[j];
                }
                
                i=np;
                while(i--) {
                        array--;        
                        indexmask--;
                        if(*indexmask) {
                                j = (unsigned)*indexmask-1;
                                sigbin[j] += (*array - meanbin[j])*(*array - meanbin[j]);
                        }       
                }
                j=nb;
                while(j--) {
                        sigbin[j] = sqrt(sigbin[j]/sumbin[j]);
                }               
        }
}


inline float operaArrayMeanQuick(unsigned np, const float *array) {
        float xmean = 0.0;
        unsigned i=np;
        while (i--) xmean += *array++;
        if (!np)
                return 0.0;
        else
                return xmean/(float)np;
}

inline float operaArrayWeightedMeanQuick(unsigned np, const float *array, const float *weigh) {
        double xmean = 0.0;
        float npf = 0.0;        
        
        while (np--){
                xmean += (*array++)*(*weigh);
                npf += *weigh++;
        }
        if (npf == 0.0)
                return 0.0;
        else
                return xmean/npf;
}

inline float operaArraySigQuick(unsigned np, const float *array) {
        
        float xsig = 0.0;       
        float xmean = 0.0;
        unsigned i=np;
        
        if(np){
                while (i--) xmean += *array++;
                xmean/=(float)np;
                i=np;
                while (i--) {
                        array--;                        
                        xsig += (*array - xmean)*(*array - xmean);
                }       
                xsig/=(float)np;
        }
        return sqrt(xsig);      
}

inline float operaArrayWeightedSigQuick(unsigned np, const float *array, const float *weigh) {
        float xmean = 0.0;
        float xsig = 0.0;
        float npf = 0.0;        
        unsigned i=np;
        
        if (np) {
                while (i--){
                        xmean += (*array)*(*weigh);             
                        npf += *weigh;
                        array++;
                        weigh++;                        
                }       
                xmean /= npf;                           // NOTE: This may throw a floating point exception if weigh is all zeroes
                i=np;
                while (i--){
                        array--;
                        weigh--;                        
                        xsig += (*weigh)*(*array - xmean)*(*array - xmean);
                }       
                xsig /= npf;
        }
        return sqrt(xsig);      
}
inline float operaArrayAvgSigClipQuick(unsigned np, const float *array, unsigned nsig) 
{
        float xsig = 0.0;       
        float xmean = 0.0;
        float clipedmean = 0.0; 
        unsigned clipednp = 0;
        unsigned i=np;
        
        if(np){
                while (i--) xmean += *array++;
                xmean/=(float)np;
                i=np;
                while (i--) {
                        array--;                        
                        xsig += (*array - xmean)*(*array - xmean);
                }       
                xsig = sqrt(xsig/(float)np);            
                
                i=np;
                while (i--) {
                        if((*array > xmean - (float)nsig*xsig) && (*array < xmean + (float)nsig*xsig)) {
                                clipedmean += *array;
                                clipednp++;
                        }
                        array++;                        
                }                       
        }
        return(clipedmean/(float)clipednp);     
}

inline float operaArrayMedianQuick(unsigned np, /*MODIFIED!!!*/float *arr) {
        unsigned low = 0, high = np-1; 
        unsigned median = (low + high) / 2; 
        unsigned odd = (np & 1);
        unsigned middle, ll, hh; 
        
        while (1) { 
                
                if (high <= low) {/* One element only */ 
                        return arr[median] ; 
                }
                if (high == low + 1) { /* Two elements only */ 
                        if (arr[low] > arr[high]) 
                                ELEM_SWAP(arr[low], arr[high]) ; 
                        if (odd)
                                return arr[median] ; 
                        else
                                return (arr[low] + arr[high]) / 2.0;
                } 
                
                /* Find median of low, middle and high items; swap into position low */ 
                middle = (low + high) / 2; 
                if (arr[middle] > arr[high]) 
                        ELEM_SWAP(arr[middle], arr[high]) 
                        if (arr[low] > arr[high]) 
                                ELEM_SWAP(arr[low], arr[high]) 
                                if (arr[middle] > arr[low]) 
                                        ELEM_SWAP(arr[middle], arr[low]) 
                                        
                                /* Swap low item (now in position middle) into position (low+1) */ 
                                        ELEM_SWAP(arr[middle], arr[low+1]) ; 
                
                /* Nibble from each end towards middle, swapping items when stuck */ 
                ll = low + 1; 
                hh = high; 
                for (;;) { 
                        do ll++; while (arr[low] > arr[ll]) ; 
                        do hh--; while (arr[hh] > arr[low]) ; 
                        
                        if (hh < ll)
                                break;
                        ELEM_SWAP(arr[ll], arr[hh])
                } 
                
                /* Swap middle item (in position low) back into correct position */ 
                ELEM_SWAP(arr[low], arr[hh]) 
                
                /* Re-set active partition */ 
                if (hh <= median) low = ll;
                if (hh >= median) high = hh - 1; 
        }       
}

inline unsigned short operaArrayMedianQuickUSHORT(unsigned np, /*MODIFIED!!!*/unsigned short *arr) {
        unsigned low = 0, high = np-1; 
        unsigned median = (low + high) / 2; 
        unsigned odd = (np & 1);
        unsigned middle, ll, hh; 
        
        while (1) { 
                
                if (high <= low) {/* One element only */ 
                        return arr[median] ; 
                }
                if (high == low + 1) { /* Two elements only */ 
                        if (arr[low] > arr[high]) 
                                ELEM_SWAP(arr[low], arr[high]) ; 
                        if (odd)
                                return arr[median] ; 
                        else
                                return (arr[low] + arr[high]) / 2.0;
                } 
                
                /* Find median of low, middle and high items; swap into position low */ 
                middle = (low + high) / 2; 
                if (arr[middle] > arr[high]) 
                        ELEM_SWAP(arr[middle], arr[high]) 
                        if (arr[low] > arr[high]) 
                                ELEM_SWAP(arr[low], arr[high]) 
                                if (arr[middle] > arr[low]) 
                                        ELEM_SWAP(arr[middle], arr[low]) 
                                        
                                /* Swap low item (now in position middle) into position (low+1) */ 
                                        ELEM_SWAP(arr[middle], arr[low+1]) ; 
                
                /* Nibble from each end towards middle, swapping items when stuck */ 
                ll = low + 1; 
                hh = high; 
                for (;;) { 
                        do ll++; while (arr[low] > arr[ll]) ; 
                        do hh--; while (arr[hh] > arr[low]) ; 
                        
                        if (hh < ll)
                                break;
                        ELEM_SWAP(arr[ll], arr[hh])
                } 
                
                /* Swap middle item (in position low) back into correct position */ 
                ELEM_SWAP(arr[low], arr[hh]) 
                
                /* Re-set active partition */ 
                if (hh <= median) low = ll;
                if (hh >= median) high = hh - 1; 
        }       
}


#endif