histo_eq.c

Go to the documentation of this file.

/**************************************************************
* Rast_histogram_eq (histo, map, min, max)
*
*   struct Histogram *histo;    histogram as returned by Rast_read_histogram()
*   unsigned char **map;        equalized category mapping
*   CELL *min, *max;            min,max category for map
*
* perform histogram equalization
* inputs are histo, output is map,min,max
****************************************************************/
#include <grass/gis.h>
#include <grass/raster.h>

void Rast_histogram_eq(const struct Histogram *histo,
                       unsigned char **map, CELL * min, CELL * max)
{
    int i;
    int x;
    CELL cat, prev;
    double total;
    double sum;
    double span;
    int ncats;
    long count;
    unsigned char *xmap;
    int len;
    int first, last;

    ncats = Rast_get_histogram_num(histo);
    if (ncats == 1) {
        *min = *max = Rast_get_histogram_cat(0, histo);
        *map = xmap = (unsigned char *)G_malloc(1);
        *xmap = 0;
        return;
    }
    if ((*min = Rast_get_histogram_cat(first = 0, histo)) == 0)
        *min = Rast_get_histogram_cat(++first, histo);
    if ((*max = Rast_get_histogram_cat(last = ncats - 1, histo)) == 0)
        *max = Rast_get_histogram_cat(--last, histo);
    len = *max - *min + 1;
    *map = xmap = (unsigned char *)G_malloc(len);

    total = 0;
    for (i = first; i <= last; i++) {
        if (Rast_get_histogram_cat(i, histo) == 0)
            continue;
        count = Rast_get_histogram_count(i, histo);
        if (count > 0)
            total += count;
    }
    if (total <= 0) {
        for (i = 0; i < len; i++)
            xmap[i] = 0;
        return;
    }

    span = total / 256;

    sum = 0.0;
    cat = *min - 1;
    for (i = first; i <= last; i++) {
        prev = cat + 1;
        cat = Rast_get_histogram_cat(i, histo);
        count = Rast_get_histogram_count(i, histo);
        if (count < 0 || cat == 0)
            count = 0;
        x = (sum + (count / 2.0)) / span;
        if (x < 0)
            x = 0;
        else if (x > 255)
            x = 255;
        sum += count;

        while (prev++ <= cat)
            *xmap++ = x;
    }
}