Histograms

Collaboration diagram for Histograms:

Classes
class	Histogram1D
class	CDF
class	IrregularHistogram1D
class	Histogram2D

Functions related to histograms 1D
template<typename T >
void	compute_hist (const MultidimArray< T > &array, Histogram1D &hist, int 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)
template<typename T >
void	compute_hist (const MultidimArray< T > &v, Histogram1D &hist, double min, double max, int no_steps)
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	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)
template<typename T >
void	reject_outliers (T &v, double percentil_out=0.25)
template<typename T >
void	histogram_equalization (MultidimArray< T > &v, int bins=8)

Functions related to histograms 2D
The vectors can be of any numerical type (short int, int, double, ...), but both of the same type. Vectors must be of the same shape, the first element of v1 and the first of v2 define the position were the first point will be inserted in the histogram, then the second of v1 and of v2, ... That is, the elements of v1 and v2 serve as coordinates within the histogram. The number of steps must always be given.
template<typename T >
void	compute_hist (const T &v1, const T &v2, Histogram2D &hist, int no_steps1, int 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)

Detailed Description

Function Documentation

compute_hist() [1/8]

template<typename T >

void compute_hist	(	const MultidimArray< T > &	array,
		Histogram1D &	hist,
		int	no_steps
	)

Compute histogram of a vector within its minimum and maximum value

Given an array as input, this function returns its histogram within the minimum and maximum of the array, in this way all the values in the array are counted. The array can be of any numerical type (short int, int, double, ...) and dimension. The number of steps must always be given.

Histogram1D hist;
compute_hist(v, hist, 100);

Definition at line 430 of file histogram.h.

{
    double min=0, max=0;
    array.computeDoubleMinMax(min, max);
    compute_hist(array, hist, min, max, no_steps);
}

compute_hist() [2/8]

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

Compute histogram of a MultidimArrayGeneric within its minimum and maximum value

Definition at line 572 of file histogram.cpp.

{
    double min=0., max=0.;
    array.computeDoubleMinMax(min, max);
    compute_hist(array, hist, min, max, no_steps);
}

compute_hist() [3/8]

template<typename T >

void compute_hist	(	const std::vector< T > &	v,
		Histogram1D &	hist,
		int	no_steps = 100
	)

Compute histogram of a vector

Definition at line 445 of file histogram.h.

{
    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);
    }
}

compute_hist() [4/8]

template<typename T >

void compute_hist	(	const MultidimArray< T > &	v,
		Histogram1D &	hist,
		double	min,
		double	max,
		int	no_steps
	)

Compute histogram of the array within two values

Given a array as input, this function returns its histogram within two values, the array values outside this range are not counted. This can be used to avoid the effect of outliers which causes a "compression" in the histogram. The array can be of any numerical type (short int, int, double, ...). The number of steps must always be given.

Histogram1D hist;
compute_hist(v, hist, 0, 1, 100);

Definition at line 486 of file histogram.h.

{
    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);
    }
}

compute_hist() [5/8]

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

Compute histogram of the MultidimArrayGeneric within two values

Definition at line 580 of file histogram.cpp.

{
    hist.init(min, max, no_steps);
#define COMPUTEHIST(type) compute_hist(MULTIDIM_ARRAY_TYPE(v,type),hist,min,max,no_steps);

    SWITCHDATATYPE(v.datatype,COMPUTEHIST)
#undef COMPUTEHIST
}

compute_hist() [6/8]

template<typename T >

void compute_hist	(	const MultidimArray< T > &	v,
		Histogram1D &	hist,
		const Matrix1D< int > &	corner1,
		const Matrix1D< int > &	corner2,
		int	no_steps = 100
	)

Compute histogram within a region (2D or 3D)

The region is defined by its corners

Definition at line 522 of file histogram.h.

{
    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);
    }
}

compute_hist() [7/8]

template<typename T >

void compute_hist	(	const T &	v1,
		const T &	v2,
		Histogram2D &	hist,
		int	no_steps1,
		int	no_steps2
	)

Compute histogram of two arrays within their minimum and maximum values

Given two arrays as input, this function returns their joint histogram within the minimum and maximum of the arrays, in this way all the values in the arrays are counted. Both arrays must have the same shape

Definition at line 966 of file histogram.h.

{
    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);
}

compute_hist() [8/8]

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
	)

Compute histogram of two arrays within given values

Given two arrays as input, this function returns their joint histogram within the specified values, all the values lying outside are not counted

Definition at line 985 of file histogram.h.

{
    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));
}

detectability_error()

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

Compute the detectability error between two pdf's

The input histograms are expressed as probability density functions representing two different objects with the same parameter of variation, for instance, the protein and background both defined by a grey-level. The first histogram would be the histogram of the grey-levels associated to voxels which we know to belong to the protein area, while the second histogram is for the grey-level of voxels which we know to belong to the background. Then the intersection area between the both associated pdf's represents the probability of committing an error when classifying a voxel, ie, the probability of saying that a voxel is protein when it really is background and viceversa. This function allows you to compute this probability error when the two histograms are provided. Be careful that this probability usually has got a very low value.

detect_error = detectability_error(hist1, hist2);

Definition at line 324 of file histogram.cpp.

{
    double hmin, hmax;
    double step;
    double v;
    double error = 0;
    int ih1, ih2; // Indexes within the histograms
    double p1, p2; // Probability associated

    // Find global range
    hmin = XMIPP_MAX(h1.hmin, h2.hmin);
    hmax = XMIPP_MIN(h1.hmax, h2.hmax);
    step = XMIPP_MIN(h1.step_size, h2.step_size) / 2;

    // Go over the range computing the errors
    v = hmin;
    int N = 0;
    while (v <= hmax)
    {
        h1.val2index(v, ih1);
        p1 = A1D_ELEM(h1, ih1) / h1.no_samples;
        h2.val2index(v, ih2);
        p2 = A1D_ELEM(h2, ih2) / h2.no_samples;
        //#define DEBUG
#ifdef DEBUG

        std::cout << "Comparing at " << v << " (" << ih1 << ") p1=" << p1 << " p2= " << p2 << std::endl;
        std::cout << "   hmin " << hmin << " hmax " << hmax << " stepsize " << h1.step_size << std::endl;
#endif//;

        if (p1 != 0 && p2 != 0)
        {
            if (p1 > p2)
                error += p2;
            else
                error += p1;
        }
        v += step;
        N++;
    }

    // Normalise such that the result is the area of a probability function
    error *= step / (hmax - hmin);
    error /= N;
#ifdef DEBUG

    std::cout << "Total error = " << error << std::endl;
#endif

    return error;
}

effective_range()

template<typename T >

double effective_range	(	const T &	v,
		double	percentil_out = 0.25
	)

Returns the effective range of a multidimensional array

The effective range is defined as the difference of those two values comprising a given central percentage of the array histogram. This function is used to compute the range removing outliers. The default central percentage is 99.75%, although this value should be increased as the number of values in the array decreases. For the default, for instance, the 0.125% of the smaller values are left out as well as the 0.125% of the higher values. The range is given always as a double number.

double range = v.effective_range();
// range for the 99.75% of the mass

double range = v.effective_range(1);
// range for the 99% of the mass

Definition at line 590 of file histogram.h.

{
    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;
}

histogram_equalization()

template<typename T >

void histogram_equalization	(	MultidimArray< T > &	v,
		int	bins = 8
	)

Histogram equalization and re-quantization

This function equalizes the histogram of the input multidimensional array, and re-quantize the input array to a specified number of bins. The output array is defined between 0 and bins-1.

Definition at line 633 of file histogram.h.

{
    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;
        }
    }
}

KLDistance()

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

Compute the Kullback-Leibler distance between two pdf's

The input histograms are expressed as probability density functions.

distance = KLDistance(hist1, hist2);

Definition at line 377 of file histogram.cpp.

{
    if (XSIZE(h1) != XSIZE(h2))
        REPORT_ERROR(ERR_MULTIDIM_SIZE,"KLDistance: Histograms of different sizes");

    double retval = 0;
    FOR_ALL_ELEMENTS_IN_ARRAY1D(h1)
    if (A1D_ELEM(h2,i) >1e-180 && A1D_ELEM(h1,i) >1e-180)
        retval += A1D_ELEM(h1,i) * log10(A1D_ELEM(h1,i) / A1D_ELEM(h2,i));
    return retval;
}

reject_outliers()

template<typename T >

void reject_outliers	(	T &	v,
		double	percentil_out = 0.25
	)

Clips the array values within the effective range

Look at the documentation of effective_rage

Definition at line 605 of file histogram.h.

{
    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;
}