histogram.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
 *              Arun Kulshreshth        (arun_2000_iitd@yahoo.com)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * 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
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef CORE_HISTOGRAM_H
#define CORE_HISTOGRAM_H

#include "multidim_array.h"
#include "multidim_array_generic.h"
#include "metadata_label.h"



class Histogram1D: public MultidimArray< double >
{
public:
    // Structure
    double hmin; // minimum value of the histogram
    double hmax; // maximum value of the histogram
    double step_size; // size of step
    double istep_size;
    int no_samples; // No. of points inside the histogram

    Histogram1D()
    {
        clear();
    }

    Histogram1D(const Histogram1D& H)
    {
        clear();
        *this = H;
    }

    void clear();

    Histogram1D& operator=(const Histogram1D& H);

    void assign(const Histogram1D& H);

    void init(double min_val, double max_val, int n_steps);

    void insert_value(double val);

#define INSERT_VALUE(histogram,value) \
{\
    if (value == histogram.hmax) { \
        size_t iii = XSIZE(histogram) - 1; \
        ++DIRECT_A1D_ELEM(histogram, iii); \
        ++histogram.no_samples; \
    } else { \
        size_t iii = (size_t) ((value - histogram.hmin) * histogram.istep_size); \
        if (iii >= 0 && iii < XSIZE(histogram)) \
        { \
            ++DIRECT_A1D_ELEM(histogram, iii); \
            ++histogram.no_samples; \
        } \
    } \
}

    double percentil(double percent_mass);

    double mass_below(double value);

    double mass_above(double value)
    {
        return no_samples - mass_below(value);
    }

    friend std::ostream& operator<<(std::ostream& o, const Histogram1D& hist);

    void write(const FileName& fn, MDLabel=MDL_X, MDLabel=MDL_COUNT);

    void val2index(double v, int& i) const
    {
        if (v == hmax)
            i = XSIZE(*this) - 1;
        else
        {
            double aux=(v - hmin) * istep_size;
            i = (int) FLOOR(aux);
        }

        if (i < 0 || i >= (int)XSIZE(*this))
            i = -1;
    }

    inline void index2val(double i, double& v) const
    {
        v = hmin + i * step_size;
    }

    double hist_min() const
    {
        return hmin;
    }

    double hist_max() const
    {
        return hmax;
    }

    double step() const
    {
        return step_size;
    }

    int stepNo() const
    {
        return XSIZE(*this);
    }

    double sampleNo() const
    {
        return no_samples;
    }

    double entropy() const;
};

 class CDF
 {
 public:
    MultidimArray<double> x;
    MultidimArray<double> probXLessThanx;
    double minVal, maxVal;
 public:
    void calculateCDF(MultidimArray<double> &V, double probStep=0.005);

    double getProbability(double x);
 };


class IrregularHistogram1D
{
public:
    Histogram1D              __hist;
    MultidimArray<double>    __binsRightLimits;
public:
    void init(const Histogram1D &oldHistogram, const MultidimArray<int> &bins);

    int val2Index(double value) const;

    void selfNormalize();

    friend std::ostream & operator << (std::ostream &_out,
                                       const IrregularHistogram1D &_h);

    inline double operator()(int i) const
    {
        return DIRECT_A1D_ELEM(__hist,i);
    }

    const Histogram1D& getHistogram() const;
};

template<typename T>
void compute_hist(const MultidimArray<T>& array, Histogram1D& hist,
                  int no_steps)
{
    double min=0, max=0;
    array.computeDoubleMinMax(min, max);
    compute_hist(array, hist, min, max, no_steps);
}

void compute_hist(const MultidimArrayGeneric& array, Histogram1D& hist,
                  int no_steps);

template<typename T>
void compute_hist(const std::vector< T > &v,
                  Histogram1D& hist,
                  int no_steps = 100)
{
    hist.clear();
    int imax=v.size();
    if (imax==0)
        return;

    // Compute minimum and maximum
    double min, max;
    min=max=v[0];
    for (int i=1; i<imax; i++)
    {
        double val=v[i];
        min=XMIPP_MIN(min,val);
        max=XMIPP_MAX(max,val);
    }
    hist.init(min, max, no_steps);

    for (int i=1; i<imax; i++)
    {
        double value=v[i];
        INSERT_VALUE(hist,value);
    }
}

template<typename T>
void compute_hist(const MultidimArray<T>& v, Histogram1D& hist,
                  double min, double max, int no_steps)
{
    hist.init(min, max, no_steps);
    T* ptr=&DIRECT_MULTIDIM_ELEM(v,0);
    size_t nmax=(MULTIDIM_SIZE(v)/4)*4;

    double value;
    for (size_t n=0; n<nmax; n+=4, ptr+=4)
    {
        value=*ptr;
        INSERT_VALUE(hist,value);
        value=*(ptr+1);
        INSERT_VALUE(hist,value);
        value=*(ptr+2);
        INSERT_VALUE(hist,value);
        value=*(ptr+3);
        INSERT_VALUE(hist,value);
    }
    for (size_t n=nmax; n<MULTIDIM_SIZE(v); ++n, ptr+=1)
    {
        value=*ptr;
        INSERT_VALUE(hist,value);
    }
}

void compute_hist(const MultidimArrayGeneric& v, Histogram1D& hist,
                  double min, double max, int no_steps);

template<typename T>
void compute_hist(const MultidimArray< T >
                  & v, Histogram1D& hist,
                  const Matrix1D< int >& corner1,
                  const Matrix1D< int >& corner2,
                  int no_steps = 100)
{
    double min, max;
    v.computeDoubleMinMax(min, max, corner1, corner2);
    hist.init(min, max, no_steps);

    Matrix1D< int > r(2);
    FOR_ALL_ELEMENTS_IN_ARRAY2D_BETWEEN(corner1, corner2)
    {
        double value=v(r);
        INSERT_VALUE(hist,value);
    }
}

double detectability_error(const Histogram1D& h1, const Histogram1D& h2);

double KLDistance(const Histogram1D& h1, const Histogram1D& h2);

template<typename T>
double effective_range(const T& v, double percentil_out = 0.25)
{
    Histogram1D hist;
    compute_hist(v, hist, 200)
    ;
    double min_val = hist.percentil(percentil_out / 2);
    double max_val = hist.percentil(100 - percentil_out / 2);
    return max_val - min_val;
}

template<typename T>
void reject_outliers(T& v, double percentil_out = 0.25)
{
    Histogram1D hist;
    compute_hist(v, hist, 400);
    double eff0 = hist.percentil(percentil_out / 2);
    double effF = hist.percentil(100 - percentil_out / 2);
    int i0, iF;
    hist.val2index(eff0, i0);
    hist.val2index(effF, iF);
    if (iF == i0) {
        hist.index2val(i0, eff0);
        hist.index2val(i0 + 1, effF);
    }

    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
    if (DIRECT_MULTIDIM_ELEM(v,n) < eff0)
        DIRECT_MULTIDIM_ELEM(v,n) = eff0;
    else if (DIRECT_MULTIDIM_ELEM(v,n) > effF)
        DIRECT_MULTIDIM_ELEM(v,n) = effF;
}

template<typename T>
void histogram_equalization(MultidimArray<T>
                            & v, int bins = 8)
{
    const int hist_steps = 200;
    Histogram1D hist;
    compute_hist(v, hist, hist_steps);

    // Compute the distribution function of the pdf
    MultidimArray<double> norm_sum(hist_steps);
    DIRECT_A1D_ELEM(norm_sum,0) = DIRECT_A1D_ELEM(hist,0);

    for (int i = 1; i < hist_steps; i++)
        DIRECT_A1D_ELEM(norm_sum,i) = DIRECT_A1D_ELEM(norm_sum,i - 1) + DIRECT_A1D_ELEM(hist,i);
    norm_sum /= MULTIDIM_SIZE(v);

    // array to store the boundary pixels of bins
    MultidimArray< double > div(bins - 1);
    int index = 0;

    for (int current_bin = 1; current_bin < bins; current_bin++)
    {
        double current_value = (double) current_bin / bins;
        while (DIRECT_A1D_ELEM(norm_sum,index) < current_value && index < FINISHINGX(norm_sum))
            index++;
        hist.index2val((double) index, DIRECT_A1D_ELEM(div,current_bin - 1));
    }

    // requantize and equalize histogram
    T* ptr=NULL;
    unsigned long int n;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY_ptr(v,n,ptr)
    {
        T vi=*ptr;
        if (vi < DIRECT_A1D_ELEM(div,0))
            *ptr = 0;
        else if (vi > DIRECT_A1D_ELEM(div,bins - 2))
            *ptr = bins - 1;
        else
        {
            index = 0;
            while (vi > DIRECT_A1D_ELEM(div,index))
                index++;
            *ptr = index;
        }
    }
}

class Histogram2D : public MultidimArray< double >
{
public:
    // Structure
    double imin; // minimum value of the i axis
    double imax; // maximum value of the i axis
    double istep_size; // size of step
    double jmin; // minimum value of the j axis
    double jmax; // maximum value of the j axis
    double jstep_size; // size of step
    int no_samples; // No. of points inside the histogram

    Histogram2D()
    {
        clear();
    }

    Histogram2D(const Histogram2D& H)
    {
        *this = H;
    }

    void clear();

    Histogram2D& operator=(const Histogram2D& H);

    void assign(const Histogram2D& H);

    void init(double imin_val, double imax_val, int in_steps,
              double jmin_val, double jmax_val, int jn_steps);

    void insert_value(double v, double u);

    friend std::ostream& operator<<(std::ostream& o, const Histogram2D& hist);

    void write(const FileName& fn);

    void val2index(double v, double u, int& i, int& j) const
    {
        if (v == imax)
            i = IstepNo() - 1;
        else
            i = (int) FLOOR((v - imin) / istep_size);

        if (u == jmax)
            j = JstepNo() - 1;

        j = (int) FLOOR((u - jmin) / jstep_size);

        if (i < 0 || i >= IstepNo())
            i = -1;

        if (j < 0 || j >= JstepNo())
            j = -1;
    }

    void index2val(double i, double j, double& v, double& u) const
    {
        v = imin + i * istep_size;
        u = jmin + j * jstep_size;
    }

    double Ihist_min() const
    {
        return imin;
    }

    double Ihist_max() const
    {
        return imax;
    }

    double Istep() const
    {
        return istep_size;
    }

    int IstepNo() const
    {
        return YSIZE(*this);
    }

    double Jhist_min() const
    {
        return jmin;
    }

    double Jhist_max() const
    {
        return jmax;
    }

    double Jstep() const
    {
        return jstep_size;
    }

    int JstepNo() const
    {
        return XSIZE(*this);
    }

    int sampleNo() const
    {
        return no_samples;
    }
};

template<typename T>
void compute_hist(const T& v1, const T& v2,
                  Histogram2D& hist, int no_steps1, int no_steps2)
{
    double min1, max1;
    v1.computeDoubleMinMax(min1, max1)
    ;

    double min2, max2;
    v2.computeDoubleMinMax(min2, max2);

    compute_hist(v1, v2, hist, min1, max1, min2, max2, no_steps1, no_steps2);
}

template<typename T>
void compute_hist(const MultidimArray<T>
                  & v1, const MultidimArray<T>& v2,
                  Histogram2D& hist,
                  double m1, double M1, double m2, double M2, int no_steps1,
                  int no_steps2)
{
    if (!v1.sameShape(v2))
        REPORT_ERROR(ERR_MULTIDIM_SIZE, "compute_hist: v1 and v2 are of different shape");

    hist.init(m1, M1, no_steps1, m2, M2, no_steps2);

    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v1)
    hist.insert_value(DIRECT_MULTIDIM_ELEM(v1, n),
                      DIRECT_MULTIDIM_ELEM(v2, n));
}

#endif