balance_lib.c

Go to the documentation of this file.

/*
 * Copyright 2009-2011 IPOL Image Processing On Line http://www.ipol.im/
 *
 * 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://www.gnu.org/licenses/>.
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <string.h>

/* ensure consistency */
#include "balance_lib.h"

static void minmax_u8(const unsigned char *data, size_t size,
                      unsigned char *ptr_min, unsigned char *ptr_max)
{
    unsigned char min, max;
    size_t i;

    /* compute min and max */
    min = data[0];
    max = data[0];
    for (i = 1; i < size; i++) {
        if (data[i] < min)
            min = data[i];
        if (data[i] > max)
            max = data[i];
    }

    /* save min and max to the returned pointers if available */
    if (NULL != ptr_min)
        *ptr_min = min;
    if (NULL != ptr_max)
        *ptr_max = max;
    return;
}

static void minmax_f32(const float *data, size_t size,
                       float *ptr_min, float *ptr_max)
{
    float min, max;
    size_t i;

    /* compute min and max */
    min = data[0];
    max = data[0];
    for (i = 1; i < size; i++) {
        if (data[i] < min)
            min = data[i];
        if (data[i] > max)
            max = data[i];
    }

    /* save min and max to the returned pointers if available */
    if (NULL != ptr_min)
        *ptr_min = min;
    if (NULL != ptr_max)
        *ptr_max = max;
    return;
}

static void quantiles_u8(const unsigned char *data, size_t size,
                         size_t nb_min, size_t nb_max,
                         unsigned char *ptr_min, unsigned char *ptr_max)
{
    /*
     * the histogram must hold all possible "unsigned char" values,
     * including 0
     */
    size_t h_size = UCHAR_MAX + 1;
    size_t histo[UCHAR_MAX + 1];
    size_t i;

    /* make a cumulative histogram */
    memset(histo, 0x00, h_size * sizeof(size_t));
    for (i = 0; i < size; i++)
        histo[(size_t) data[i]] += 1;
    for (i = 1; i < h_size; i++)
        histo[i] += histo[i - 1];

    /* get the new min/max */

    if (NULL != ptr_min) {
        /* simple forward traversal of the cumulative histogram */
        /* search the first value > nb_min */
        i = 0;
        while (i < h_size && histo[i] <= nb_min)
            i++;
        /* the corresponding histogram value is the current cell position */
        *ptr_min = (unsigned char) i;
    }

    if (NULL != ptr_max) {
        /* simple backward traversal of the cumulative histogram */
        /* search the first value <= size - nb_max */
        i = h_size - 1;
        /* i is unsigned, we check i<h_size instead of i>=0 */
        while (i < h_size && histo[i] > (size - nb_max))
            i--;
        /*
         * if we are not at the end of the histogram,
         * get to the next cell,
         * the last (backward) value > size - nb_max
         */
        if (i < h_size - 1)
            i++;
        *ptr_max = (unsigned char) i;
    }
    return;
}

static int cmp_f32(const void *a, const void *b)
{
    if (*(const int *) a > *(const int *) b)
        return 1;
    if (*(const int *) a < *(const int *) b)
        return -1;
    return 0;
}

static void quantiles_f32(const float *data, size_t size,
                          size_t nb_min, size_t nb_max,
                          float *ptr_min, float *ptr_max)
{
    float *data_tmp;

    data_tmp = (float *) malloc(size * sizeof(float));

    /* copy the data and sort */
    memcpy(data_tmp, data, size * sizeof(float));
    qsort(data_tmp, size, sizeof(float), &cmp_f32);

    /* get the min/max */
    if (NULL != ptr_min)
        *ptr_min = data_tmp[nb_min];
    if (NULL != ptr_max)
        *ptr_max = data_tmp[size - 1 - nb_max];

    free(data_tmp);
    return;
}

static unsigned char *rescale_u8(unsigned char *data, size_t size,
                                 unsigned char min, unsigned char max)
{
    size_t i;

    if (max <= min)
        for (i = 0; i < size; i++)
            data[i] = UCHAR_MAX / 2;
    else {
        /* build a normalization table */
        unsigned char norm[UCHAR_MAX + 1];
        for (i = 0; i < min; i++)
            norm[i] = 0;
        for (i = min; i < max; i++)
            /*
             * we can't store and reuse UCHAR_MAX / (max - min) because
             *     105 * 255 / 126.            -> 212.5, rounded to 213
             *     105 * (double) (255 / 126.) -> 212.4999, rounded to 212
             */
            norm[i] = (unsigned char) ((i - min) * UCHAR_MAX
                                       / (double) (max - min) + .5);
        for (i = max; i < UCHAR_MAX + 1; i++)
            norm[i] = UCHAR_MAX;
        /* use the normalization table to transform the data */
        for (i = 0; i < size; i++)
            data[i] = norm[(size_t) data[i]];
    }
    return data;
}

static float *rescale_f32(float *data, size_t size, float min, float max)
{
    size_t i;

    if (max <= min)
        for (i = 0; i < size; i++)
            data[i] = .5;
    else
        for (i = 0; i < size; i++)
            data[i] = (min > data[i] ? 0. :
                       (max < data[i] ? 1. : (data[i] - min) / (max - min)));

    return data;
}

unsigned char *balance_u8(unsigned char *data, size_t size,
                          size_t nb_min, size_t nb_max)
{
    unsigned char min, max;

    /* sanity checks */
    if (NULL == data) {
        fprintf(stderr, "a pointer is NULL and should not be so\n");
        abort();
    }
    if (nb_min + nb_max >= size) {
        nb_min = (size - 1) / 2;
        nb_max = (size - 1) / 2;
        fprintf(stderr, "the number of pixels to flatten is too large\n");
        fprintf(stderr, "using (size - 1) / 2\n");
    }

    /* get the min/max */
    if (0 != nb_min || 0 != nb_max)
        quantiles_u8(data, size, nb_min, nb_max, &min, &max);
    else
        minmax_u8(data, size, &min, &max);

    /* rescale */
    (void) rescale_u8(data, size, min, max);

    return data;
}

float *balance_f32(float *data, size_t size, size_t nb_min, size_t nb_max)
{
    float min, max;

    /* sanity checks */
    if (NULL == data) {
        fprintf(stderr, "a pointer is NULL and should not be so\n");
        abort();
    }
    if (nb_min + nb_max >= size) {
        nb_min = (size - 1) / 2;
        nb_max = (size - 1) / 2;
        fprintf(stderr, "the number of pixels to flatten is too large\n");
        fprintf(stderr, "using (size - 1) / 2\n");
    }

    /* get the min/max */
    if (0 != nb_min || 0 != nb_max)
        quantiles_f32(data, size, nb_min, nb_max, &min, &max);
    else
        minmax_f32(data, size, &min, &max);

    /* rescale */
    (void) rescale_f32(data, size, min, max);

    return data;
}