rt_band_get_histogram()

rt_histogram rt_band_get_histogram	(	rt_bandstats	stats,
		uint32_t	bin_count,
		double *	bin_width,
		uint32_t	bin_width_count,
		int	right,
		double	min,
		double	max,
		uint32_t *	rtn_count
	)

Count the distribution of data.

Parameters

stats	: a populated stats struct for processing
bin_count	: the number of bins to group the data by
bin_width	: the width of each bin as an array
bin_width_count	: number of values in bin_width
right	: evaluate bins by (a,b] rather than default [a,b)
min	: user-defined minimum value of the histogram a value less than the minimum value is not counted in any bins if min = max, min and max are not used
max	: user-defined maximum value of the histogram a value greater than the max value is not counted in any bins if min = max, min and max are not used
rtn_count	: set to the number of bins being returned

Returns

the histogram of the data or NULL

Definition at line 412 of file rt_statistics.c.

  {
        rt_histogram bins = NULL;
        int init_width = 0;
        uint32_t i;
        uint32_t j;
        double tmp;
        double value;
        int sum = 0;
        double qmin;
        double qmax;
 
#if POSTGIS_DEBUG_LEVEL > 0
        clock_t start, stop;
        double elapsed = 0;
#endif
 
        RASTER_DEBUG(3, "starting");
#if POSTGIS_DEBUG_LEVEL > 0
        start = clock();
#endif
 
        assert(NULL != stats);
        assert(NULL != rtn_count);
 
        if (stats->count < 1 || NULL == stats->values) {
                rterror("rt_util_get_histogram: rt_bandstats object has no value");
                return NULL;
        }
 
        /* bin width must be positive numbers and not zero */
        if (NULL != bin_width && bin_width_count > 0) {
                for (i = 0; i < bin_width_count; i++) {
                        if (bin_width[i] < 0 || FLT_EQ(bin_width[i], 0.0)) {
                                rterror("rt_util_get_histogram: bin_width element is less than or equal to zero");
                                return NULL;
                        }
                }
        }
 
        /* ignore min and max parameters */
        if (FLT_EQ(max, min)) {
                qmin = stats->min;
                qmax = stats->max;
        }
        else {
                qmin = min;
                qmax = max;
                if (qmin > qmax) {
                        qmin = max;
                        qmax = min;
                }
        }
 
        /* # of bins not provided */
        if (bin_count <= 0) {
                /*
                        determine # of bins
                        http://en.wikipedia.org/wiki/Histogram
 
                        all computed bins are assumed to have equal width
                */
                /* Square-root choice for stats->count < 30 */
                if (stats->count < 30)
                        bin_count = ceil(sqrt(stats->count));
                /* Sturges' formula for stats->count >= 30 */
                else
                        bin_count = ceil(log2((double) stats->count) + 1.);
 
                /* bin_width_count provided and bin_width has value */
                if (bin_width_count > 0 && NULL != bin_width) {
                        /* user has defined something specific */
                        if (bin_width_count > bin_count)
                                bin_count = bin_width_count;
                        else if (bin_width_count > 1) {
                                tmp = 0;
                                for (i = 0; i < bin_width_count; i++) tmp += bin_width[i];
                                bin_count = ceil((qmax - qmin) / tmp) * bin_width_count;
                        }
                        else
                                bin_count = ceil((qmax - qmin) / bin_width[0]);
                }
                /* set bin width count to zero so that one can be calculated */
                else {
                        bin_width_count = 0;
                }
        }
 
        /* min and max the same */
        if (FLT_EQ(qmax, qmin)) bin_count = 1;
 
        RASTER_DEBUGF(3, "bin_count = %d", bin_count);
 
        /* bin count = 1, all values are in one bin */
        if (bin_count < 2) {
                bins = rtalloc(sizeof(struct rt_histogram_t));
                if (NULL == bins) {
                        rterror("rt_util_get_histogram: Could not allocate memory for histogram");
                        return NULL;
                }
 
                bins->count = stats->count;
                bins->percent = -1;
                bins->min = qmin;
                bins->max = qmax;
                bins->inc_min = bins->inc_max = 1;
 
                *rtn_count = bin_count;
                return bins;
        }
 
        /* establish bin width */
        if (bin_width_count == 0) {
                bin_width_count = 1;
 
                /* bin_width unallocated */
                if (NULL == bin_width) {
                        bin_width = rtalloc(sizeof(double));
                        if (NULL == bin_width) {
                                rterror("rt_util_get_histogram: Could not allocate memory for bin widths");
                                return NULL;
                        }
                        init_width = 1;
                }
 
                bin_width[0] = (qmax - qmin) / bin_count;
        }
 
        /* initialize bins */
        bins = rtalloc(bin_count * sizeof(struct rt_histogram_t));
        if (NULL == bins) {
                rterror("rt_util_get_histogram: Could not allocate memory for histogram");
                if (init_width) rtdealloc(bin_width);
                return NULL;
        }
        if (!right)
                tmp = qmin;
        else
                tmp = qmax;
        for (i = 0; i < bin_count;) {
                for (j = 0; j < bin_width_count; j++) {
                        bins[i].count = 0;
                        bins->percent = -1;
 
                        if (!right) {
                                bins[i].min = tmp;
                                tmp += bin_width[j];
                                bins[i].max = tmp;
 
                                bins[i].inc_min = 1;
                                bins[i].inc_max = 0;
                        }
                        else {
                                bins[i].max = tmp;
                                tmp -= bin_width[j];
                                bins[i].min = tmp;
 
                                bins[i].inc_min = 0;
                                bins[i].inc_max = 1;
                        }
 
                        i++;
                }
        }
        if (!right) {
                bins[bin_count - 1].inc_max = 1;
 
                /* align last bin to the max value */
                if (bins[bin_count - 1].max < qmax)
                        bins[bin_count - 1].max = qmax;
        }
        else {
                bins[bin_count - 1].inc_min = 1;
 
                /* align first bin to the min value */
                if (bins[bin_count - 1].min > qmin)
                        bins[bin_count - 1].min = qmin;
        }
 
        /* process the values */
        for (i = 0; i < stats->count; i++) {
                value = stats->values[i];
 
                /* default, [a, b) */
                if (!right) {
                        for (j = 0; j < bin_count; j++) {
                                if (
                                        (!bins[j].inc_max && value < bins[j].max) || (
                                                bins[j].inc_max && (
                                                        (value < bins[j].max) ||
                                                        FLT_EQ(value, bins[j].max)
                                                )
                                        )
                                ) {
                                        bins[j].count++;
                                        sum++;
                                        break;
                                }
                        }
                }
                /* (a, b] */
                else {
                        for (j = 0; j < bin_count; j++) {
                                if (
                                        (!bins[j].inc_min && value > bins[j].min) || (
                                                bins[j].inc_min && (
                                                        (value > bins[j].min) ||
                                                        FLT_EQ(value, bins[j].min)
                                                )
                                        )
                                ) {
                                        bins[j].count++;
                                        sum++;
                                        break;
                                }
                        }
                }
        }
 
        for (i = 0; i < bin_count; i++) {
                bins[i].percent = ((double) bins[i].count) / sum;
        }
 
#if POSTGIS_DEBUG_LEVEL > 0
        stop = clock();
        elapsed = ((double) (stop - start)) / CLOCKS_PER_SEC;
        RASTER_DEBUGF(3, "elapsed time = %0.4f", elapsed);
 
        for (j = 0; j < bin_count; j++) {
                RASTER_DEBUGF(5, "(min, max, inc_min, inc_max, count, sum, percent) = (%f, %f, %d, %d, %d, %d, %f)",
                        bins[j].min, bins[j].max, bins[j].inc_min, bins[j].inc_max, bins[j].count, sum, bins[j].percent);
        }
#endif
 
        if (init_width) rtdealloc(bin_width);
        *rtn_count = bin_count;
        RASTER_DEBUG(3, "done");
        return bins;
}

References rt_bandstats_t::count, rt_histogram_t::count, genraster::count, FLT_EQ, rt_histogram_t::inc_max, rt_histogram_t::inc_min, rt_bandstats_t::max, rt_histogram_t::max, rt_bandstats_t::min, rt_histogram_t::min, rt_histogram_t::percent, RASTER_DEBUG, RASTER_DEBUGF, rtalloc(), rtdealloc(), rterror(), genraster::value, and rt_bandstats_t::values.

Referenced by RASTER_histogram(), and test_band_stats().

Here is the call graph for this function:

Here is the caller graph for this function: