vmquantc.c

有关matlab的电子书籍有一定的帮助希望有用

/* $Revision: 5.7 $ */
/* Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved. */

static char rcsid[] = "$Id: vmquantc.c,v 5.7 1997/11/24 15:57:16 eddins Exp $";

/*********************************************************************/

/*
 *
 *      vmquantc.c
 *
 *      VMQUANTC        Variance Minimization Color Quantization
 *
 *      [IM MAP] = VMQUANT( R, G, B, K, [Qr Qg Qb], DITHER, Qe ) will
 *      convert an arbitrary image comprised of RGB triples into an 
 *      indexed image IM with color map MAP.  K specifies the number 
 *      of desired entries in the target color map, and [Qr Qg Qb] 
 *      specified the number of quantization levels to assign each color 
 *      axis during color interpolation.  DITHER is a boolean that 
 *      indicates whether or not to perform error propagation
 *      dithering on output matrix.  Qe specifies the number of
 *      bits of quantization used in the error calculations.
 *
 *      1 <= Qe <= 30
 *      1 <= Qr,Qg,Qb <= 31
 *      
 *      See also: RGB2GRAY, DITHER, IMAPPROX.
 *
 *
 *      Joseph M. Winograd 7-93
 *
 *      Reference: 
 *        Xiaolin Wu, "Efficient Statistical Computation for Optimal
 *          Color Quantization," Graphics Gems II, (ed. James Arvo).
 *          Academic Press: Boston. 1991.
 *
 *      This is a MEX file for MATLAB.
 *
 */

#include "mex.h"

/*
 *      Compilation options. 
 *
 *  Define NOREG if m1box and m1base can't be in registers.
 */
#define NOREG

#define Q_CORRECTION    (1)     /* Compile with color quantization correction.
                                        (slower but more precise)       */
#define MINIMIZE_VOLUME (0)     /* Minimize volume as well as variance. */
                                /*                      (Experimental)  */


/*
 *      Macros.
 */

#define vmMAX(x,y)      (((x)>(y))?(x):(y))
#define clamp(x,l,h)    ((x)<(l)?(l):(x)>(h)?(h):(x))

#define elem3d(x,y,z)   ((x)+(y)*dim3x+(z)*dim3xy)      /* Use MATLAB/FORTRAN */
#define elem2d(x,y)     ((x)+(y)*dim2x)                 /* style indexing. */

#define v0(x)           (x[elem3d(boxll[box][R],boxll[box][G],boxll[box][B])])
#define v1(x)           (x[elem3d(boxll[box][R],boxll[box][G],boxur[box][B])])
#define v2(x)           (x[elem3d(boxll[box][R],boxur[box][G],boxll[box][B])])
#define v3(x)           (x[elem3d(boxll[box][R],boxur[box][G],boxur[box][B])])
#define v4(x)           (x[elem3d(boxur[box][R],boxll[box][G],boxll[box][B])])
#define v5(x)           (x[elem3d(boxur[box][R],boxll[box][G],boxur[box][B])])
#define v6(x)           (x[elem3d(boxur[box][R],boxur[box][G],boxll[box][B])])
#define v7(x)           (x[elem3d(boxur[box][R],boxur[box][G],boxur[box][B])])
#define region(x)                                               \
        (-x[elem3d(boxll[box][R],boxll[box][G],boxll[box][B])]  \
        + x[elem3d(boxll[box][R],boxll[box][G],boxur[box][B])]  \
        + x[elem3d(boxll[box][R],boxur[box][G],boxll[box][B])]  \
        - x[elem3d(boxll[box][R],boxur[box][G],boxur[box][B])]  \
        + x[elem3d(boxur[box][R],boxll[box][G],boxll[box][B])]  \
        - x[elem3d(boxur[box][R],boxll[box][G],boxur[box][B])]  \
        - x[elem3d(boxur[box][R],boxur[box][G],boxll[box][B])]  \
        + x[elem3d(boxur[box][R],boxur[box][G],boxur[box][B])])

/*
 *      Constants.
 */

/* Numerical constants. */

#define R               (0)     /* Parameter list/array indexing constants. */
#define G               (1)
#define B               (2)
#define K               (3)
#define Q               (4)
#define DITHER          (5)
#define QE              (6)

static int flops;
    

/*
 *      Wu's Variance Minimization Color Quantization Algorithm.
 *
 *    Return value is the number of color map entries created.
 *
 */

int32_T vmquant( uint8_T *r, uint8_T *g, uint8_T *b, uint32_T k, 
                 uint32_T qr, uint32_T qg, uint32_T qb, uint32_T qe, 
                 uint16_T *image, real64_T *map, uint32_T m, uint32_T n, 
                 int32_T **rqp, int32_T **gqp, int32_T **bqp, 
                 uint32_T *tmp_boxll, uint32_T *tmp_boxur, 
                 uint32_T *p, int32_T *m0, int32_T *m1r, int32_T *m1g, 
                 int32_T *m1b, int32_T *m2, real64_T *E, int32_T *squares ) 
{
    uint32_T **boxll;
    uint32_T **boxur;
    register int32_T e, x, y, z, box;
#ifdef NOREG /* Don't put m1box and m1base in a register. */
    register int32_T wbox, m0base;
    int32_T m1box[3], m1base[3];
#else
    register int32_T wbox, m1box[3], m0base, m1base[3];
#endif
    int32_T mindim, minind, boxes, lastbox;
    int32_T wx, m1x[3], wdiff, m1diff[3];
    int32_T *rq = *rqp, *gq = *gqp, *bq = *bqp;
    real64_T minvar, var;
    int32_T rlevels = (1 << qr) + 1, glevels = (1 << qg) + 1, 
        blevels = (1 << qb) + 1, elevels = 1 << qe;
    int32_T rshift, gshift, bshift, gmask, bmask;
    uint32_T rglevels = rlevels*glevels;
    int32_T x0, x1, x2, x3, xstep;
    uint32_T dim2x = m, dim3x = rlevels, dim3xy = rglevels;
    real64_T scale;
    /*
     * Make z volatile to force it out of a register.  Makes sure
     * that multiplication by 255 is rounded correctly on extended
     * precision machines (MAC, PC, HP300).
     */
    volatile double zr,zg,zb;

    /* initialize boxll, boxur pointers */
    boxll = (uint32_T **) mxCalloc(k, sizeof(*boxll));
    boxur = (uint32_T **) mxCalloc(k, sizeof(*boxur));
    for (x = 0; x < k; x++)
    {
        boxll[x] = tmp_boxll + x*3;
        boxur[x] = tmp_boxur + x*3;
    }

    /* Quantize r, g, b into rq, gq, bq (if not saving memory)
            and
        Calculate p (Probabilty Density Function). */

    rshift = qe - qr;           /* Set up shift variables to convert */
    gshift = qe - qg;           /* elevels to r,g,blevels. */
    bshift = qe - qb;

    if (elevels == 256)
    {
        /* we can save a lot of work */
        for (x = 0; x < m*n; x++) {
            /* Map [0,255] onto [1,xlevels-1]. */

            zr = (real64_T) r[x] + 0.5;
            zg = (real64_T) g[x] + 0.5;
            zb = (real64_T) b[x] + 0.5;
            rq[x] = (int32_T) zr;
            gq[x] = (int32_T) zg;
            bq[x] = (int32_T) zb;
            p[ elem3d( (rq[x] >> rshift) + 1, (gq[x] >> gshift) + 1, 
                       (bq[x] >> bshift) + 1 ) ]++;
        }
        flops += m * n;
    }
    else
    {
        scale = ((real64_T) elevels - 1.0) / 255.0;
        
        for (x = 0; x < m*n; x++) {
            /* Map [0,255] onto [1,xlevels-1]. */
            
            zr = r[x] * scale + 0.5;
            zg = g[x] * scale + 0.5;
            zb = b[x] * scale + 0.5;
            rq[x] = (int32_T) zr;
            gq[x] = (int32_T) zg;
            bq[x] = (int32_T) zb;
            p[ elem3d( (rq[x] >> rshift) + 1, (gq[x] >> gshift) + 1, 
                       (bq[x] >> bshift) + 1 ) ]++;
        }
        flops += m * n * 2 + 2;
    }

    /* Precalculate squares. */

    for (x = 0; x <= vmMAX(rlevels, vmMAX(glevels, blevels)); x++) {
        squares[x] = x*x;
    }

    /* Calculate m0, m1, m2.  (Zeroth, first, second moments of p) */

    /* Cumsum in R. */

    e = elem3d(1,1,1);
    for (z = 2; z <= blevels; z++) {
        for (y = 2; y <= glevels; y++) {
            for (x = 2; x <= rlevels; x++) {
                m0[e] = p[e] + m0[e-1];
                m1r[e] = x * p[e] + m1r[e-1];
                m1g[e] = y * p[e] + m1g[e-1];
                m1b[e] = z * p[e] + m1b[e-1];
                m2[e] = (squares[x] + squares[y] + squares[z]) * p[e] + 
                    m2[e-1];
                e++;
            }
            e++;
        }
        e += rlevels;
    }
            
    /* Cumsum in G. */

    for (x = 1; x <= rlevels; x++) {
        e = elem3d(x,1,1);
        for (z = 1; z < blevels; z++) {
            for (y = 1; y < glevels; y++) {
                m0[e] = m0[e] + m0[e-rlevels];
                m1r[e] = m1r[e] + m1r[e-rlevels];
                m1g[e] = m1g[e] + m1g[e-rlevels];
                m1b[e] = m1b[e] + m1b[e-rlevels];
                m2[e] = m2[e] + m2[e-rlevels];
                e += rlevels;
            }
            e += rlevels;
        }
    }
            
    /* Cumsum in B. */

    e = elem3d(1,1,1);
    for (z = 1; z < blevels; z++) {
        for (y = 1; y < glevels; y++) {
            for (x = 1; x < rlevels; x++) {
                m0[e] = m0[e] + m0[e-rglevels];
                m1r[e] = m1r[e] + m1r[e-rglevels];
                m1g[e] = m1g[e] + m1g[e-rglevels];
                m1b[e] = m1b[e] + m1b[e-rglevels];
                m2[e] = m2[e] + m2[e-rglevels];
                e++;
            }
            e++;
        }
        e += rlevels;
    }

    /* Initialize variables for first partition. */

    box = 0;
    lastbox = 0;
    boxll[0][R] = 0;  boxll[0][G] = 0;  boxll[0][B] = 0;
    boxur[0][R] = rlevels-1; boxur[0][G] = glevels-1; boxur[0][B] = blevels-1;
    wbox = region(m0);
    m1box[R] = region(m1r);
    m1box[G] = region(m1g);
    m1box[B] = region(m1b);
    E[0] = region(m2) - ((real64_T) m1box[R]*m1box[R] + 
        (real64_T) m1box[G]*m1box[G] + (real64_T) m1box[B]*m1box[B]) / 
        (real64_T) wbox;

    flops += 7;

#if MINIMIZE_VOLUME
    E[0] *= ((int32_T) boxur[0][R]-boxll[0][R]) * (boxur[0][G]-boxll[0][G]) *
                (boxur[0][B]-boxll[0][B]);
    E[0] *= ((int32_T) boxur[0][R]-boxll[0][R]) * (boxur[0][G]-boxll[0][G]) *
                (boxur[0][B]-boxll[0][B]);
    E[0] *= ((int32_T) boxur[0][R]-boxll[0][R]) * (boxur[0][G]-boxll[0][G]) *
                (boxur[0][B]-boxll[0][B]);

    flops += 3;

#endif

    /* Find k-1 partitions. */

    for (boxes = 0; boxes < k-1; boxes++) {

        /* Find the box with the greatest variance. */

        for (x = box = 0; x <= boxes; x++) {
            if (E[x] > E[box]) {
                box = x;
            }
        }

        /* If no boxes have any variance than all colors have been found.
                (Thanks, Clay!) */
        
        if (!E[box]) {
            break;
        }

        /* Find point of minimum variance in r dimension of box. */

        minvar = 0;
        mindim = R;
        minind = 0;
        m0base = -v0(m0)+v1(m0)+v2(m0)-v3(m0);
        m1base[R] = -v0(m1r)+v1(m1r)+v2(m1r)-v3(m1r);
        m1base[G] = -v0(m1g)+v1(m1g)+v2(m1g)-v3(m1g);
        m1base[B] = -v0(m1b)+v1(m1b)+v2(m1b)-v3(m1b);
        if (box != lastbox) {
            wbox = m0base+v4(m0)-v5(m0)-v6(m0)+v7(m0);
            m1box[R] = m1base[R]+v4(m1r)-v5(m1r)-v6(m1r)+v7(m1r);
            m1box[G] = m1base[G]+v4(m1g)-v5(m1g)-v6(m1g)+v7(m1g);
            m1box[B] = m1base[B]+v4(m1b)-v5(m1b)-v6(m1b)+v7(m1b);
        }

        for (y = R; y <= B; y++) {

            switch (y) {
                case R: {
                    x0 = elem3d(boxll[box][R]+1, boxll[box][G], boxll[box][B]);
                    x1 = elem3d(boxll[box][R]+1, boxll[box][G], boxur[box][B]);
                    x2 = elem3d(boxll[box][R]+1, boxur[box][G], boxll[box][B]);
                    x3 = elem3d(boxll[box][R]+1, boxur[box][G], boxur[box][B]);
                    xstep = 1;
                    break;
                }
                case G: {
                    m0base = -v0(m0)+v1(m0)+v4(m0)-v5(m0);
                    m1base[R] = -v0(m1r)+v1(m1r)+v4(m1r)-v5(m1r);
                    m1base[G] = -v0(m1g)+v1(m1g)+v4(m1g)-v5(m1g);
                    m1base[B] = -v0(m1b)+v1(m1b)+v4(m1b)-v5(m1b);
                    x0 = elem3d(boxll[box][R], boxll[box][G]+1, boxll[box][B]);
                    x1 = elem3d(boxll[box][R], boxll[box][G]+1, boxur[box][B]);
                    x2 = elem3d(boxur[box][R], boxll[box][G]+1, boxll[box][B]);
                    x3 = elem3d(boxur[box][R], boxll[box][G]+1, boxur[box][B]);
                    xstep = rlevels;
                    break;
                }
                case B: {
                    m0base = -v0(m0)+v2(m0)+v4(m0)-v6(m0);
                    m1base[R] = -v0(m1r)+v2(m1r)+v4(m1r)-v6(m1r);
                    m1base[G] = -v0(m1g)+v2(m1g)+v4(m1g)-v6(m1g);
                    m1base[B] = -v0(m1b)+v2(m1b)+v4(m1b)-v6(m1b);
                    x0 = elem3d(boxll[box][R], boxll[box][G], boxll[box][B]+1);
                    x1 = elem3d(boxll[box][R], boxur[box][G], boxll[box][B]+1);
                    x2 = elem3d(boxur[box][R], boxll[box][G], boxll[box][B]+1);
                    x3 = elem3d(boxur[box][R], boxur[box][G], boxll[box][B]+1);
                    xstep = rglevels;
                    break;
                }
            }

            for (x = boxll[box][y]+1; x < boxur[box][y]; x++) {
                wx = m0base + m0[x0] - m0[x1] - m0[x2] + m0[x3];
                m1x[R] = m1base[R] + m1r[x0] - m1r[x1] - m1r[x2] + m1r[x3];
                m1x[G] = m1base[G] + m1g[x0] - m1g[x1] - m1g[x2] + m1g[x3];
                m1x[B] = m1base[B] + m1b[x0] - m1b[x1] - m1b[x2] + m1b[x3];
                if (wx) {
                    var = ((real64_T) m1x[R]*m1x[R] + 
                                (real64_T) m1x[G]*m1x[G] + 
                                (real64_T) m1x[B]*m1x[B]) / (real64_T) wx;
                    flops += 6;
                }
                else {
                    var = 0;
                }
                if (wbox != wx) {
                    m1diff[R] = m1box[R] - m1x[R];
                    m1diff[G] = m1box[G] - m1x[G];
                    m1diff[B] = m1box[B] - m1x[B];
                    var += ((real64_T) m1diff[R]*m1diff[R] + 
                        (real64_T) m1diff[G]*m1diff[G] + 
                        (real64_T) m1diff[B]*m1diff[B]) / 
                        (real64_T) (wbox - wx);
                    flops += 7;
                }
                if (var > minvar) {     /* Maximize var to minimize variance! */


                    minvar = var;
                    minind = x;
                    mindim = y;
                }
                x0 += xstep;
                x1 += xstep;
                x2 += xstep;
                x3 += xstep;
            } /* end for each slice in y dimension */
        } /* end for each dimension */

        /* Cut the box at the point of minimum variance. */


        for (x = boxes; x >= box; x--) {
            boxll[x+1][R] = boxll[x][R];        
            boxll[x+1][G] = boxll[x][G];        
            boxll[x+1][B] = boxll[x][B];        
            boxur[x+1][R] = boxur[x][R];
            boxur[x+1][G] = boxur[x][G];
            boxur[x+1][B] = boxur[x][B];
            E[x+1] = E[x];      /* Also make room in E for new box. */
        }
        boxll[box+1][mindim] = minind;
        boxur[box][mindim] = minind;

        /* Compute the variance for two new boxes. */

        x = box;
        lastbox = box + 1;
        for (box = x; box < x+2; box++) {

            if (boxur[box][R]-boxll[box][R]+boxur[box][G]-boxll[box][G]+
                boxur[box][B]-boxll[box][B] > 3) {
                wbox = region(m0);
                m1box[R] = region(m1r);
                m1box[G] = region(m1g);
                m1box[B] = region(m1b);
                E[box] = region(m2) - ((real64_T) m1box[R]*m1box[R] + 
                    (real64_T) m1box[G]*m1box[G] +
                    (real64_T) m1box[B]*m1box[B]) / (real64_T) wbox;
                flops += 7;
#if MINIMIZE_VOLUME
                E[box] *= ((int32_T) boxur[box][R]-boxll[box][R]) * 
                        (boxur[box][G]-boxll[box][G]) *
                        (boxur[box][B]-boxll[box][B]);
                E[box] *= ((int32_T) boxur[box][R]-boxll[box][R]) *
                        (boxur[box][G]-boxll[box][G]) *
                        (boxur[box][B]-boxll[box][B]);
                E[box] *= ((int32_T) boxur[box][R]-boxll[box][R]) *
                        (boxur[box][G]-boxll[box][G]) *
                        (boxur[box][B]-boxll[box][B]);
                flops += 3;
#endif
            } /* end if box can be subdivided */
            else {
                E[box] = 0.0;
            } /* end else box cannot be subdivided */
        } /* end for each new box compute variance */

    } /* end for each color split a box */

    k = boxes + 1;

    /*  Build inverse color map on top of m2.  
                and
        If Q_CORRECTION is not being used then also calculate color map. */

    /* Build color map using power of 2 indexing scheme where
        [0,1] maps to [0,xlevels-2] */

    gmask = ~((1 << gshift) - 1);
    bmask = ~((1 << bshift) - 1);
    gshift = qr - gshift;