/*
 *            FAST IMAGE HISTAGRAM EQUILIZATION ALGORITHM
 *
 *                         Copywrite 1998
 *                         by Melinda Green
 *                         melinda@superliminal.com
 *                         all rights reserved
 *
 * Permission is granted for use in all non-commercial applications.
 * All commercial use of this code or algorithm requires written 
 * permission.
 *
 * DESCRIPION:
 *
 *    When given an array of integer values need to be converted
 *    into an image, colors must be assigned for each value. This
 *    module assumes that a color ramp will be selected, and
 *    that the lowest integer values should be assigned color map
 *    indices at one end of that ramp, and the higest values indices
 *    at the other end. The problem is that there are any number of
 *    ways that those indices can be assigned. The most obvious is to
 *    partition the ramp linearly from the lowest value to the highest,
 *    but that only gives reasonable results when the input values are
 *    equally distributed. In the worst case, almost all pixels in
 *    an image would have very similar integer values, but a few would
 *    be very different. Using linear interpolation, all the similar
 *    values would be assigned only one or two unique color indices and
 *    would be indistinguishable in the final image. All the fine
 *    detail would be lost. Histagram equalization attempts to allocate
 *    ranges of the color ramp proportional to the image area ocupied
 *    by various pixel value ranges so that more colors are allocated
 *    to the pixel ranges that most need it, and only few colors are
 *    allocated for those rare pixel values.
 *
 *    A proper histagram equalization algorithm would require a careful
 *    analysis of the input image in order to allocate colors in the
 *    most ideal way. The technique used here does not perform an ideal
 *    allocation of color ranges, but it has the advantages of being
 *    very simple and extremely fast.
 *
 *    This module contains the following public methods:
 *         select_samples()
 *         lookup_color()
 *         dump_samples()
 *
 *    select_samples() takes an array of integer values representing an
 *    image to be colored and sets up the data structures needed to
 *    perform subsequent coloring of pixel values. It must be called
 *    before lookup_color(). It returns the number of colors lookup_colors
 *    requires up to the maximum (currently hardcoded at 256 but easily
 *    changed).
 *
 *    lookup_color() takes a single raw pixel integer to be colored,
 *    and returns the colormap index that should be used to color that
 *    pixel.
 *
 *    dump_samples() is simply a debugging aid which can be used to
 *    print out the internal sample array after select_samples has run.
 */

#include <stdio.h>
#include <math.h> /* for qsort */

/*
 * Largest number of colors allowed.
 */
#define MAX_MAP 256


/*
 * When choosing color samples, this is the maximum number of pixels
 * to randomly sample.
 */
#define MAX_COLOR_TRYS (MAX_MAP * 10)



const unsigned long MAX_VAL = 65931;



/*
 * A sampled table of exit values used to partition the color map
 * into color ranges.
 */
int Samples[MAX_MAP];
int NSamples; /* The number of samples in the 'Samples' array.  */



int lookup_color(val)
int val;
{
    int i;

    for(i=0; i<NSamples; i++)
		if(val <= Samples[i])
			return(i);
    return(0);
}




static int comp_func(a, b)
int *a, *b;
{
    if(*a < *b)
		return(-1);
    else if(*b < *a)
		return(1);
    return(0);
}


void dump_sample(char *str)
{
    int i;
    fprintf(stderr, "%s\n", str);
    for(i=0; i<NSamples; i++) {
		fprintf(stderr, "%d ", Samples[i]);
		if(9 == (i % 10))
			fprintf(stderr, "\n");
    }
    fprintf(stderr, "\n");
}


int in_set(int val, int set_size, int set[])
{
    int i;
    for(i=0; i<set_size; i++)
		if(set[i] == val)
			return 1;
    return 0;
}


/*
 * Builds the `Samples' array.  If all goes well, this array will contain
 * 'goal' unique exit values of random values from the 'vals' array.
 * 'nvals' is the length of that array.
 * This will not be possible when enough unique values can't be found,
 * in which case one color will be assigned for each unique exit value.
 *
 * The global 'NSamples' value is set to be the number of values
 * loaded into the global 'Samples' array.
 */

int select_samples(int goal, int nvals, int vals[])
{
    int i, j, attempts;
    int test_samples[MAX_MAP], nunique;

    printf("trying to select %d from %d vals\n", goal, nvals);
    nunique = 0;
    for(i=0; i<nvals; i++) {
		int is_unique=1;
		for(j=0; j<nunique; j++) /* see if we have this one yet */
			if(test_samples[j] == vals[i])
				is_unique = 0;
		if(is_unique)
			test_samples[nunique++] = vals[i];
		if(nunique >= goal)
			goto enough_colors;
    }

    /*
     * if here, we must not be enough colors to fill the whole map
     * so we'll take what we can get.
     */
    for(i=0; i<nunique; i++)
		Samples[i] = test_samples[i];
    NSamples = nunique;
    fprintf(stderr, "only %d values\n", NSamples);
    goto done_sampling;

    enough_colors:

    fprintf(stderr, "found enough samples to fill color map\n");

    /*
     * there exist enough colors, so fill the map with random samples
     */
    NSamples = 0;
    for(attempts=0; NSamples<goal-1; attempts++)
    {
		/*
		 * After examining MAX_COLOR_TRYS number of pixels, give up.
		 */
		if(attempts > MAX_COLOR_TRYS) {
			fprintf(stderr, "select_samples: quitting random selection\n");
			break;
		}

		/* choose a random pixel */
		i = random() % nvals;
		if( ! in_set(vals[i], NSamples, Samples))
			Samples[NSamples++] = vals[i];
    }

    /*
     * fill any remaining slots which occure when random sampling quit.
     */
    for(i=0; i<nunique && NSamples<goal-1; i++)
		if(!in_set(test_samples[i], NSamples, Samples))
			Samples[NSamples++] = test_samples[i];

    done_sampling:

    Samples[NSamples++] = MAX_VAL;
    qsort(Samples, NSamples, sizeof(int), comp_func);
    /* dump_sample("samples"); */

	return NSamples;
}