image procssing toolbox 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++) {