/*
    This program is copyright (c) 2012 Volker Schatz.

    It 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 can obtain a copy of the GNU General Public License from the following
    URL: http://www.gnu.org/licenses/gpl.html .

    
    To compile, run the following command: (or just: grep ^gcc pnmclahe.c | sh)

gcc -o pnmclahe pnmclahe.c -lm
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <netinet/in.h>


void usageexit(void);
char *addbeforeext( const char *fname, const char *suffix );
void skipcomment( char **read );
unsigned char *load_pnm(const char *filename, int *width, int *height, int *maxval);
void *loadfile(const char *filename, unsigned *filesize);
void write_pnm( const char *fname, unsigned char *buf, int width, int height, int maxval );
void clahescale(unsigned *data, int width, int height, int max, int ntilesx, int ntilesy, double normlimit, unsigned *scaleunit);
unsigned *touns(const unsigned char *orig, int size, int origmax, int *monomax);
void scaleuns(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit);
unsigned ushortmax3(const unsigned short *data);
unsigned ucharmax3(const unsigned char *data);
unsigned *getbrightness(const unsigned char *orig, int size, int origmax, int *monomax);
unsigned satscale(unsigned orig, unsigned scale, unsigned scaleunit, unsigned max);
void scalebrightness(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit);
unsigned *getsaturation(const unsigned char *orig, int size, int origmax, int *monomax);
void scalesaturation(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit);


int main(int argc, char **argv)
{
    unsigned *(*getmono)(const unsigned char *orig, int size, int origmax, int *monomax);
    void (*scalemono)(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit);
    unsigned *mono;
    unsigned char *pic;
    char *parseend, *destname;
    double limit= 3.0;
    unsigned unit;
    int argind, width, height, size, origmax, monomax;
    int ntilesx= 1, ntilesy= 1;

    if( argc < 2 )
        usageexit();
    argind= 1;
    getmono= getbrightness;
    scalemono= scalebrightness;
    while( argv[argind][0] == '-' ) {
        if( !strcmp(argv[argind], "-l") ) {
            if( argc <= argind+1 )
                usageexit();
            limit= strtod(argv[argind+1], NULL);
            argind += 2;
        }
        else if( !strcmp(argv[argind], "-t") ) {
            if( argc <= argind+1 )
                usageexit();
            ntilesx= strtol(argv[argind+1], &parseend, 10);
            if( ntilesx <= 0 || parseend == argv[argind+1] || *parseend != 'x' )
                usageexit();
            argv[argind+1]= parseend+1;
            ntilesy= strtol(argv[argind+1], &parseend, 10);
            if( ntilesy <= 0 || parseend == argv[argind+1] || *parseend )
                usageexit();
            argind += 2;
        }
        else if( !strcmp(argv[argind], "-s") ) {
            getmono= getsaturation;
            scalemono= scalesaturation;
            ++argind;
        }
        else if( !strcmp(argv[argind], "-h") )
            usageexit();
        else {
            fprintf(stderr, "Unknown option %s.\n", argv[argind]);
            usageexit();
        }
    }
    if( argc != argind+1 && argc != argind+2 )
        usageexit();
    pic= load_pnm(argv[argind], &width, &height, &origmax);
    size= width*height;
    if( origmax > 0 ) {
        getmono= touns;
        scalemono= scaleuns;
    }
    mono= (*getmono)(pic, size, origmax, &monomax);
    clahescale(mono, width, height, monomax, ntilesx, ntilesy, limit, &unit);
    (*scalemono)(pic, size, origmax, mono, unit);
    if( argc == argind + 2 )
        write_pnm(argv[argind+1], pic, width, height, origmax);
    else {
        destname= addbeforeext(argv[argind], "_clahe");
        if( !destname ) {
            fprintf(stderr, "Could not allocate memory for destination file name.  Aborting.\n");
            exit(1);
        }
        write_pnm(destname, pic, width, height, origmax);
        free(destname);
    }
    free(mono);
    free(pic);
    return 0;
}


void usageexit(void)
{
    printf("usage: pnmclahe [-s] [-t <xtiles>x<ytiles>] [-l <limit>] <source.pnm> [ <dest.pnm> ]\n"
           "Apply contrast-limited adaptive histogram equalisation to a greyscale or\n"
           "colour picture.  Defaults are 1x1 tiles (i.e., non-adaptive CLHE) and a\n"
           "normalised clip limit of 3.  -s operates on saturation, not brightness.\n");
    exit(2);
}


/*
  addbeforeext  Add a string to a file name just before the extension (if any)

  ->    fname   original file name
        suffix  string to append to trunk of file name
  <-  New file name, or NULL on failure.  Has to be freed by the caller.
*/
char *addbeforeext( const char *fname, const char *suffix )
{
    char *dest;
    const char *read;

    dest= malloc(strlen(fname) + strlen(suffix) + 1);
    if( !dest )
        return NULL;
    read= strrchr(fname, '.');
    if( read ) {
        strncpy(dest, fname, read-fname);
        strcpy(dest + (read-fname), suffix);
        strcpy(dest + (read-fname) + strlen(suffix), read);
    }
    else {
        strcpy(dest, fname);
        strcat(dest, suffix);
    }
    return dest;
}


/*
  skipcomment   Skip comment in PNM header - skips optional white space, then
                '#' up to the next newline

  <->   read    pointer to read position
*/
void skipcomment( char **read )
{
    while( isspace(**read) )      ++*read;
    if( **read == '#' )
        while( **read != '\r' && **read != '\n' )   ++*read;
}


/*
    Load a .pnm image.  maxval < 0 means colour (3 components); |maxval| > 255
    means the return pointer is really an unsigned short *.
*/
unsigned char *load_pnm(const char *filename, int *width, int *height, int *maxval)
{
    unsigned char *data;
    char *scan;
    unsigned linesize;
    int filesize, w, h, max, count;

    data= loadfile(filename, &filesize);
    if( filesize < 10 || data[0] != 'P' || (data[1] != '5' && data[1] != '6') ) {
        fprintf(stderr, "%s is not a binary PNM file\n", filename);
        exit(1);
    }
    scan= data + 2;
    skipcomment(&scan);
    w= strtol(scan, &scan, 10);
    skipcomment(&scan);
    h= strtol(scan, &scan, 10);
    skipcomment(&scan);
    max= strtol(scan, &scan, 10);
    if( !w || !h || !max ) {
        fprintf(stderr, "Illegal header values in %s\n", filename);
        exit(1);
    }
    if( data[1] == '6' )
        max= -max;          /* indicate colour */
    filesize -= scan + 1 - (char*)data;
    memmove(data, scan + 1, filesize);
    linesize= max > 0 ? w : 3*w;
    if( abs(max) > 255 ) {
        linesize = 2*linesize;
        for( count= 0; count< filesize/2; ++count )
            ((short*)data)[count]= ntohs(((short*)data)[count]);
    }
    if( (unsigned)h * linesize > filesize )
        h= filesize / linesize;
    if( width )
        *width= w;
    if( height )
        *height= h;
    if( maxval )
        *maxval= max;
    return data;
}


void *loadfile(const char *filename, unsigned *filesize)
{
    FILE *in;
    void *buf;
    unsigned size;

    in= fopen(filename, "rb");
    if( !in ) {
        fprintf(stderr, "Could not open `%s' for reading.\n", filename);
        exit(1);
    }
    fseek(in, 0, SEEK_END);
    size= ftell(in);
    fseek(in, 0, SEEK_SET);
    buf= malloc(size);
    if( !buf ) {
        fprintf(stderr, "Could not allocate memory for loading `%s'.\n", filename);
        fclose(in);
        exit(1);
    }
    if( ferror(in) || fread(buf, 1, size, in) != size ) {
        fprintf(stderr, "Error reading `%s'.\n", filename);
        free(buf);
        fclose(in);
        exit(1);
    }
    fclose(in);
    if( filesize )
        *filesize= size;
    return buf;
}


/*
    Write a .pnm image.  16-bit pixel values will have their bytes swapped.
*/
void write_pnm( const char *fname, unsigned char *buf, int width, int height, int maxval )
{
    FILE *out;
    int count, size, samplebytes, colour;

    size= width*height;
    if( colour= maxval < 0 ) {
        size *= 3;
        maxval= -maxval;
    }
    samplebytes= maxval > 255 ? 2 : 1;
    if( maxval > 255 )
        for( count= 0; count< size; ++count )
            ((short*)buf)[count]= htons(((short*)buf)[count]);
    out= fopen(fname, "wb");
    if( !out ) {
        fprintf(stderr, "Could not create output file `%s'.  Aborting.\n", fname);
        exit(1);
    }
    fprintf(out, "P%c\n%u %u\n%u\n", colour ? '6' : '5', width, height, maxval);
    fwrite(buf, samplebytes, size, out);
    if( ferror(out) ) {
        fprintf(stderr, "Error writing `%s'.\n", fname);
        fclose(out);
        exit(1);
    }
    fclose(out);
}



#define MINAVGBIN   2
#define REDISTPREC  8
#define REDISTLOOPS 5

void clahescale(unsigned *data, int width, int height, int max, int ntilesx, int ntilesy, double normlimit, unsigned *scaleunit)
{
    unsigned **tilehists, **tilexforms, *xformmax, *hist, *row, **xforms, *xmax;
    size_t allocsize;
    unsigned long long unit;
    unsigned xformval, limit, histval, excess, nbelow, excperbin, cumredist, prevredist;
    int tilew, tileh, halftilew, halftileh, tilesize;
    int xtile, ytile, intilex, intiley, x, y;
    int count, histbits, histshift, histmax, histind, origval, redistcount;

    tilew= (width + ntilesx/2) / ntilesx;
    tileh= (height + ntilesy/2) / ntilesy;
    if( tilew * ntilesx != width || tileh * ntilesy != height )
        fprintf(stderr, "Warning: Requested number of tiles does not divide width or height - some tiles will have diverging sizes.\n");
    halftilew= tilew/2;
    halftileh= tileh/2;
    tilesize= tilew*tileh;

    /* choose histogram bin size so that average bin fill >= MINAVGBIN */
    /* histbits= (int)floor(log2(tilesize / MINAVGBIN));  */
    histshift= (int)floor(log2((max+1) * MINAVGBIN / tilesize));
    /* histdiv= (double)(max+1) * MINAVGBIN / tilesize; */
    if( histshift < 0 )
        histshift= 0;
    histmax= max >> histshift;

    /* allocate and initialise tile histograms and transformation functions */
    allocsize= 2 * ntilesx*ntilesy * (sizeof(unsigned*) + sizeof(unsigned)*(histmax+1) + sizeof(unsigned));
    tilehists= malloc(allocsize);
    if( !tilehists ) {
        fprintf(stderr, "Could not allocate memory for tile histograms.\n");
        exit(1);
    }
    memset(tilehists, 0, allocsize);
    tilexforms= (unsigned **)((char*)tilehists + ntilesx*ntilesy * (sizeof(unsigned*) + sizeof(unsigned)*(histmax+1)));
    tilehists[0]= (unsigned*)(tilehists + ntilesx*ntilesy);
    tilexforms[0]= (unsigned*)(tilexforms + ntilesx*ntilesy);
    for( count= 1; count< ntilesx*ntilesy; ++count ) {
        tilehists[count]= tilehists[count-1] + histmax+1;
        tilexforms[count]= tilexforms[count-1] + histmax+1;
    }
    xformmax= tilexforms[ntilesx*ntilesy-1] + histmax+1;

    /* compute tile histograms */
    ytile= 0;
    intiley= tileh;
    row= data;
    for( y= 0; y < height; ++y, row += width ) {
        hist= tilehists[ytile*ntilesx];
        intilex= tilew;
        for( x= 0; x < width; ++x ) {
            histind= row[x] >> histshift;
            ++hist[histind];
            if( --intilex == 0 && x+halftilew < width ) {
                intilex= tilew;
                hist += histmax+1;
            }
        }
        if( --intiley == 0 && y+halftileh < height ) {
            intiley= tileh;
            ++ytile;
            hist += histmax + 1;
        }
        else
            hist -= (ntilesx-1)*(histmax+1);
    }

    /* clip, redistribute and integrate histogram to obtain transformation function */
    if( normlimit*tilesize/(histmax+1) >= UINT_MAX )
        limit= UINT_MAX;
    else
        limit= (unsigned)floor(normlimit*tilesize/(histmax+1) + 0.5);
    for( count= 0; count < ntilesx*ntilesy; ++count ) {
        nbelow= excess= 0;
        for( histind= 0; histind <= histmax; ++histind ) {
            histval= tilehists[count][histind];
            if( histval > limit ) {
                excess += histval - limit;
                tilehists[count][histind]= limit;
            }
            else if( histval < limit )
                ++nbelow;
        }
        for( redistcount= 0; excess && nbelow && redistcount < REDISTLOOPS; ++redistcount ) {
            /* we add REDISTPREC fractional bits to the excess calculation to
               redistribute it more evenly.  */
            excperbin= (excess << REDISTPREC) / nbelow;
            cumredist= prevredist= 0;
            nbelow= excess= 0;
            for( histind= 0; histind <= histmax; ++histind ) {
                histval= tilehists[count][histind];
                if( histval == limit )
                    continue;
                cumredist += excperbin;
                histval += (cumredist >> REDISTPREC) - prevredist;
                prevredist= cumredist >> REDISTPREC;
                if( histval > limit ) {
                    excess += histval - limit;
                    tilehists[count][histind]= limit;
                }
                else {
                    tilehists[count][histind]= histval;
                    if( histval < limit )
                        ++nbelow;
                }
            }
        }
        xformval= 0;
        for( histind= 0; histind <= histmax; ++histind ) {
            tilexforms[count][histind]= xformval;
            xformval += tilehists[count][histind];
        }
        xformmax[count]= xformval - histval + 1;
    }

    /* apply transformation, but save scaling factor, not result */
    /* choose unit so that intermediate quantities fit into 64 bits, see below */
    if( ntilesx > 1 && ntilesy > 1 )
        unit= ULLONG_MAX / tilesize / max / tilesize;
    else
        unit= ULLONG_MAX / (ntilesx > 1 ? tilew : tileh) / max / tilesize;
    if( unit > 0x10000 )
        unit= 0x10000;
    else if( !unit ) {
        fprintf(stderr, "Warning: Insufficient headroom for bilinear interpolation; contrast amplification will be inaccurate and intermediate results may overflow.\n");
        unit= 1;
    }
    else if( unit < max )
        fprintf(stderr, "Warning: Constrained headroom for bilinear interpolation; contrast amplification will be inaccurate.\n");
    *scaleunit= (unsigned)unit;
    unit *= max;
    /* first, the top half-tile row */
    for( y= 0; y< halftileh; ++y ) {
        for( x= 0; x< halftilew; ++x, ++data ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= unit * tilexforms[0][origval] / origval / xformmax[0];
        }
        xtile= 0;
        intilex= 0;
        if( ntilesx > 1 ) while( 13 ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= (unit*intilex*tilexforms[xtile+1][origval] / xformmax[xtile+1]
                + unit*(tilew-1-intilex)*tilexforms[xtile][origval] / xformmax[xtile])
                / tilew / origval;
            ++x;
            ++data;
            if( ++intilex >= tilew ) {
                if( ++xtile == ntilesx-1 )
                    break;
                intilex= 0;
            }
        }
        for( ; x< width; ++x, ++data ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= unit * tilexforms[xtile][origval] / origval / xformmax[xtile];
        }
    }
    /* full tile rows in the bulk of the image */
    ytile= 0;
    intiley= 0;
    xforms= tilexforms;
    xmax= xformmax;
    if( ntilesy > 1 ) while( 13 ) {
        for( x= 0; x< halftilew; ++x, ++data ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= (unit*intiley*xforms[ntilesx][origval] / xmax[ntilesx]
                    + unit*(tileh-1-intiley)*xforms[0][origval] / xmax[0])
                    / tileh / origval;
        }
        xtile= 0;
        intilex= 0;
        if( ntilesx > 1 ) while( 13 ) {
            origval= *data;
            if( !origval ) *data= unit; else
            /* unit was chosen so that these products still fits into 64 bits */
            *data= (unit*intilex*intiley*xforms[ntilesx+1][origval] / xmax[ntilesx+1]
                + unit*(tilew-1-intilex)*intiley*xforms[ntilesx][origval] / xmax[ntilesx]
                + unit*intilex*(tileh-1-intiley)*xforms[1][origval] / xmax[1]
                + unit*(tilew-1-intilex)*(tileh-1-intiley)*xforms[0][origval] / xmax[0])
                / tilesize / origval;
            ++x;
            ++data;
            if( ++intilex >= tilew ) {
                ++xforms;
                ++xmax;
                if( ++xtile == ntilesx-1 )
                    break;
                intilex= 0;
            }
        }
        for( ; x< width; ++x, ++data ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= (unit*intiley*xforms[ntilesx][origval] / xmax[ntilesx]
                    + unit*(tileh-1-intiley)*xforms[0][origval] / xmax[0])
                    / tileh / origval;
        }
        ++y;
        if( ++intiley >= tileh ) {
            ++xforms;
            ++xmax;
            if( ++ytile == ntilesy-1 )
                break;
            intiley= 0;
        }
        else {
            xforms -= ntilesx - 1;
            xmax -= ntilesx - 1;
        }
    }
    /* last, the bottom half-tile row */
    for( ; y< height; ++y ) {
        for( x= 0; x< halftilew; ++x, ++data ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= unit * xforms[0][origval] / origval / xmax[0];
        }
        xtile= 0;
        intilex= 0;
        if( ntilesx > 1 ) while( 13 ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= (unit*intilex*xforms[xtile+1][origval] / xmax[xtile+1]
                    + unit*(tilew-1-intilex)*xforms[xtile][origval] / xmax[xtile])
                    / tilew / origval;
            ++x;
            ++data;
            if( ++intilex >= tilew ) {
                if( ++xtile == ntilesx-1 )
                    break;
                intilex= 0;
            }
        }
        for( ; x< width; ++x, ++data ) {
            origval= *data;
            if( !origval ) *data= unit; else
            *data= unit * xforms[xtile][origval] / origval / xmax[xtile];
        }
    }
}


/*************************************************************************
    Different functions for extracting and scaling per-pixel quantities
    that CLAHE can be applied to.
*************************************************************************/

/************************************************************************
    For greyscale images just copy data to unsigned array, and scale
    original value.
*/

unsigned *touns(const unsigned char *orig, int size, int origmax, int *monomax)
{
    unsigned *mono;
    const unsigned short *orig16= (const unsigned short *)orig;
    int count;

    mono= malloc(size * sizeof(unsigned));
    if( !mono ) {
        fprintf(stderr, "Could not allocate memory for intermediate image.\n");
        exit(1);
    }
    if( origmax > 255 )
        for( count= 0; count< size; ++count )
            mono[count]= orig16[count];
    else
        for( count= 0; count< size; ++count )
            mono[count]= orig[count];
    *monomax= origmax;
    return mono;
}


void scaleuns(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit)
{
    unsigned short *dest16= (unsigned short *)dest;
    int count;

    origmax= abs(origmax);
    if( origmax > 255 )
        for( count= 0; count< size; ++count )
            dest16[count]= satscale(dest16[count], scale[count], scaleunit, origmax);
    else
        for( count= 0; count< size; ++count )
            dest[count]= satscale(dest[count], scale[count], scaleunit, origmax);
}


/************************************************************************
    For colour images there is a choice of quantities.

    First, brightness = the maximum RGB colour component.  We scale the other
    components by the same factor.
*/

unsigned ushortmax3(const unsigned short *data)
{
    if( data[0] > data[1] )
        if( data[0] > data[2] )     return (unsigned)data[0];
        else                        return (unsigned)data[2];
    else
        if( data[1] > data[2] )     return (unsigned)data[1];
        else                        return (unsigned)data[2];
}


unsigned ucharmax3(const unsigned char *data)
{
    if( data[0] > data[1] )
        if( data[0] > data[2] )     return (unsigned)data[0];
        else                        return (unsigned)data[2];
    else
        if( data[1] > data[2] )     return (unsigned)data[1];
        else                        return (unsigned)data[2];
}


unsigned *getbrightness(const unsigned char *orig, int size, int origmax, int *monomax)
{
    unsigned *mono;
    const unsigned short *orig16= (const unsigned short *)orig;
    int count, min, max;

    mono= malloc(size * sizeof(unsigned));
    if( !mono ) {
        fprintf(stderr, "Could not allocate memory for intermediate image.\n");
        exit(1);
    }
    origmax= abs(origmax);
    if( origmax > 255 )
        for( count= 0; count< size; ++count )
            mono[count]= ushortmax3(orig16 + 3*count);
    else
        for( count= 0; count< size; ++count )
            mono[count]= ucharmax3(orig + 3*count);
    *monomax= origmax;
    return mono;
}


/*
    Scale a value while observing a saturation threshold.
*/
unsigned satscale(unsigned orig, unsigned scale, unsigned scaleunit, unsigned max)
{
    unsigned scaled= orig*scale;

    if( scaled > max*scaleunit )    return max;
    else                            return scaled / scaleunit;
}


void scalebrightness(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit)
{
    unsigned short *dest16= (unsigned short *)dest;
    int count;

    origmax= abs(origmax);
    if( origmax > 255 )
        for( count= 0; count< size; ++count ) {
            dest16[3*count]= (unsigned short)satscale(dest16[3*count], scale[count], scaleunit, origmax);
            dest16[3*count+1]= (unsigned short)satscale(dest16[3*count+1], scale[count], scaleunit, origmax);
            dest16[3*count+2]= (unsigned short)satscale(dest16[3*count+2], scale[count], scaleunit, origmax);
        }
    else
        for( count= 0; count< size; ++count ) {
            dest[3*count]= (unsigned char)satscale(dest[3*count], scale[count], scaleunit, origmax);
            dest[3*count+1]= (unsigned char)satscale(dest[3*count+1], scale[count], scaleunit, origmax);
            dest[3*count+2]= (unsigned char)satscale(dest[3*count+2], scale[count], scaleunit, origmax);
        }
}


/************************************************************************
    Saturation = the difference between maximum and minimum RGB colour
    component.

    Scaling the RGB pixels is done by transforming the three colour components
    so that the ratio between the minimum component and the headroom above the
    maximum component remains the same.
*/

unsigned *getsaturation(const unsigned char *orig, int size, int origmax, int *monomax)
{
    unsigned *mono;
    const unsigned short *orig16= (const unsigned short *)orig;
    int count, min, max;

    mono= malloc(size * sizeof(unsigned));
    if( !mono ) {
        fprintf(stderr, "Could not allocate memory for intermediate image.\n");
        exit(1);
    }
    origmax= abs(origmax);
    if( origmax > 255 )
        for( count= 0; count< size; ++count, orig16 += 3 )
            mono[count]= (abs((int)orig16[0]-(int)orig16[1])
                        + abs((int)orig16[1]-(int)orig16[2])
                        + abs((int)orig16[2]-(int)orig16[0]))/2;
    else
        for( count= 0; count< size; ++count, orig += 3 )
            mono[count]= (abs((int)orig[0]-(int)orig[1])
                        + abs((int)orig[1]-(int)orig[2])
                        + abs((int)orig[2]-(int)orig[0]))/2;
    *monomax= origmax;
    return mono;
}


void scalesaturation(unsigned char *dest, int size, int origmax, unsigned *scale, unsigned scaleunit)
{
    unsigned short *dest16= (unsigned short *)dest;
    int result[3];
    int count, minind, midind, maxind, oldrange, newrange;
    int chgpix, chgsum;

    origmax= abs(origmax);
    if( origmax > 255 )
        for( count= 0; count< size; ++count, dest16 += 3, ++scale ) {
            if( dest16[0] > dest16[1] )
                if( dest16[1] > dest16[2] ) {
                    maxind= 0;  midind= 1;  minind= 2;
                }
                else {
                    minind= 1;
                    if( dest16[0] > dest16[2] ) { maxind= 0;  midind= 2; }
                    else                        { maxind= 2;  midind= 0; }
                }
            else
                if( dest16[2] > dest16[1] ) {
                    maxind= 2;  midind= 1;  minind= 0;
                }
                else {
                    maxind= 1;
                    if( dest16[2] > dest16[0] ) { minind= 0;  midind= 2; }
                    else                        { minind= 2;  midind= 0; }
                }
            oldrange= dest16[maxind]-dest16[minind];
            if( oldrange == origmax )
                continue;
            newrange= satscale(oldrange, *scale, scaleunit, origmax);
            result[minind]= dest16[minind] * (origmax - newrange) / (origmax - oldrange);
            result[midind]= dest16[minind] + satscale(dest16[midind] - dest16[minind], *scale, scaleunit, origmax - dest16[minind]);
            result[maxind]= dest16[minind] + newrange;
            if( result[maxind] > origmax )
                result[maxind]= origmax;
            dest16[0]= result[0];
            dest16[1]= result[1];
            dest16[2]= result[2];
        }
    else
        for( count= 0; count< size; ++count, dest += 3, ++scale ) {
            if( dest[0] > dest[1] )
                if( dest[1] > dest[2] ) {
                    maxind= 0;  midind= 1;  minind= 2;
                }
                else {
                    minind= 1;
                    if( dest[0] > dest[2] ) { maxind= 0;  midind= 2; }
                    else                    { maxind= 2;  midind= 0; }
                }
            else
                if( dest[2] > dest[1] ) {
                    maxind= 2;  midind= 1;  minind= 0;
                }
                else {
                    maxind= 1;
                    if( dest[2] > dest[0] ) { minind= 0;  midind= 2; }
                    else                    { minind= 2;  midind= 0; }
                }
            oldrange= dest[maxind]-dest[minind];
            if( oldrange == origmax )
                continue;
            newrange= satscale(oldrange, *scale, scaleunit, origmax);
            chgsum += newrange - oldrange;
            ++chgpix;
            result[minind]= dest[minind] * (origmax - newrange) / (origmax - oldrange);
            result[midind]= dest[minind] + satscale(dest[midind] - dest[minind], *scale, scaleunit, origmax - dest[minind]);
            result[maxind]= dest[minind] + newrange;
            if( result[maxind] > origmax )
                result[maxind]= origmax;
            dest[0]= result[0];
            dest[1]= result[1];
            dest[2]= result[2];
        }
}