ditherc.c

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

/* Copyright 1993-1998 The MathWorks, Inc. All Rights Reserved. */

/* $Revision: 5.5 $  $Date: 1997/11/24 15:56:55 $  $Author: eddins $ */

static char rcsid[] = "$Id: ditherc.c,v 5.5 1997/11/24 15:56:55 eddins Exp $";

/*
 *
 *	ditherc.c
 *
 *	DITHERC.MEX		Error Propagation Image Dithering 
 *
 *	IM = DITHER( R, G, B, MAP, Qm, Qe ) will dither the R, G, B encoded
 *	image against the colormap MAP to create the dithered indexed image
 *	IM.  Qm specifies the number of quantization bits to use along each 
 *	color axis of the inverse color map.  Qe specifies the number of
 *	quantization bits to use for color space error calculations.
 *
 *      Qm must be at least 1 and no greater than 10.
 *      Qe must be at least 1 and no greater than 31.
 *
 *      This function requires temporary memory storage proportional to
 *      c1*2^(3*Qm) + c2*2^(Qe+1) + c3*prod(size(R)).  c1, c2 and c3 are
 *      architecture dependant.  Memory usage is typically around 256Kb for
 *      a 100x100 image using the default parameters.
 *
 *	Joseph M. Winograd 7-93
 *
 *	References: 
 *	  R. W. Floyd and L. Steinberg, "An Adaptive Algorithm for
 *	    Spatial Gray Scale," International Symposium Digest of Technical
 *	    Papers, Society for Information Displays, 36.  1975.
 *        Spencer W. Thomas, "Efficient Inverse Color Map Computation",
 *          Graphics Gems II, (ed. James Arvo), Academic Press: Boston.
 *          1991. (includes source code)
 *
 */
 
#include "mex.h"
#include "invcmap.h"
 
#define clamp(x,l,h)    ((x)<(l)?(l):(x)>(h)?(h):(x))
#define maxx(x,y)	((x)>(y)?(x):(y))
#define elem2d(x,y)     ((x)+(y)*dim2x)                 

#define R		(0)	/* Parameter list/array indexing constants. */
#define G		(1)
#define B		(2)
#define MAP		(3)
#define QM		(4)	/* Quantization levels in map. */
#define QE		(5)	/* Quantization levels in error calculation. */

static uint32_T flops;

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

/*
 *	Floyd-Steinberg Error Propogation Dithering
 */

void fsdither( uint8_T *r, uint8_T *g, uint8_T *b, 
               int32_T *mapq, int32_T k, int32_T *inverse_colormap, 
               uint16_T *dithered_image, int32_T m, int32_T n, 
               int32_T qm, int32_T qe, int32_T *rq, int32_T *gq, int32_T *bq, 
               int32_T *tab1, int32_T *tab3, int32_T *tab5, int32_T *tab7)
{
    real64_T scale;
    register int32_T x, y, z, Rval, Gval, Bval;
    int32_T pixel_value;
    int32_T dim2x = k, elevels, mask, bitshift0,
	bitshift1, bitshift2, Roffset, Goffset,
	Boffset;
    /*
     * 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 real64_T zr,zg,zb;

    /* Handle special case of no dithering possible. */

    if (qe < qm) {
	elevels = (1 << qm) - 1;	/* Really mlevels-1 */
        scale = (real64_T) elevels / 255.0;
	for (z = 0; z < m*n; z++) {
            zr = (real64_T) r[z] * scale + 0.5;
            zg = (real64_T) g[z] * scale + 0.5;
            zb = (real64_T) b[z] * scale + 0.5;
	    dithered_image[z] = (uint16_T) (inverse_colormap[
		(((int32_T) zr) << (qm+qm)) +
		(((int32_T) zg) << qm) +
                ((int32_T) zb) ]); 
	}
        flops += m * n * 2 + 1;
	return;
    }

    /* Set up for dither. */

    elevels = (1 << qe);
    bitshift0 = qe - qm;
    bitshift1 = qm - bitshift0;
    bitshift2 = qm + bitshift1;
    mask = ~((1 << bitshift0) - 1);
    Roffset = elem2d(0,R);
    Goffset = elem2d(0,G);
    Boffset = elem2d(0,B);

    /* Calculate fast math table values. */
 
    y = elevels / 32;
    for (x = -elevels; x < 0; x++) {
        tab1[x+elevels] = -((y-x) / 16);
        tab3[x+elevels] = -((y-3*x) / 16);
        tab5[x+elevels] = -((y-5*x) / 16);
        tab7[x+elevels] = -((y-7*x) / 16);
    }
    for (x = 0; x < elevels; x++) {
        tab1[x+elevels] = ((y+x) / 16);
        tab3[x+elevels] = ((y+3*x) / 16);
        tab5[x+elevels] = ((y+5*x) / 16);
        tab7[x+elevels] = ((y+7*x) / 16);
    }
 
 
    for (x = 0; x < 3*k; x++) {
        mapq[x] -= elevels;
    }
 
    /* Quantize r, g, b, color map entries. */

    /* When elevels = 256, we can save a lot of work here */
    if (elevels == 256)
    {
        for (x = 0; x < m*n; x++)
        {
            rq[x] = (int32_T) r[x];
            gq[x] = (int32_T) g[x];
            bq[x] = (int32_T) b[x];
        }
    }
    else
    {
        scale = (elevels - 1.0) / 255.0;
        for (x = 0; x < m*n; x++) {
            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;
        }
        flops += m * n * 2 + 2;
    }

 
    for (x = z = 0; x < n; x++) {
 
        for (y = 0; y < m; y++, z++) {
 
            /* Place this pixel. */
 
            pixel_value = inverse_colormap[
		 (((rq[z]=clamp(rq[z],0,elevels-1)) & mask)<<bitshift2) +
                 (((gq[z]=clamp(gq[z],0,elevels-1)) & mask)<<bitshift1) +
		 ((bq[z]=clamp(bq[z],0,elevels-1))>>bitshift0) ];
	    dithered_image[z] = (uint16_T) pixel_value;
 
            /* Calculate errors. */
 
            Rval = rq[z] -= mapq[pixel_value + Roffset];
            Gval = gq[z] -= mapq[pixel_value + Goffset];
            Bval = bq[z] -= mapq[pixel_value + Boffset];
 
 
            /* Apply error to next pixel. */
 
            if (y != m-1) {
            	rq[z+1] += tab7[Rval];
            	gq[z+1] += tab7[Gval];
            	bq[z+1] += tab7[Bval];
            }
 
            /* Apply error to right pixel. */
 
            if (x != n-1) {
            	rq[z+m] += tab5[Rval];
            	gq[z+m] += tab5[Gval];
            	bq[z+m] += tab5[Bval];
 
            	/* Apply error to upper right pixel. */
 
            	if (y) {
        	    rq[z+m-1] += tab3[Rval];
        	    gq[z+m-1] += tab3[Gval];
        	    bq[z+m-1] += tab3[Bval];
            	}
 
            	/* Apply error to lower right pixel. */
 
            	if (y != m-1) {
            	    rq[z+m+1] += tab1[Rval];
            	    gq[z+m+1] += tab1[Gval];
            	    bq[z+m+1] += tab1[Bval];
                }
	    }
        }
    }

}
 
 
/*
 *      Gateway routine.
 */

void ValidateInputs(int nlhs,
                    mxArray *plhs[],
                    int nrhs,
                    const mxArray *prhs[],
                    int32_T *qm,
                    int32_T *qe)
{
    int32_T m;
    int32_T n;
    int32_T p;
    int32_T q;
    int32_T i;
    real64_T *map;
                
    if (nrhs != 6) 
    {
        mexErrMsgTxt("ditherc expects six input arguments");
    }
    
    m = mxGetM(prhs[R]);
    n = mxGetN(prhs[R]);
    
    if ((m != mxGetM(prhs[G])) || (n != mxGetN(prhs[G])) ||
        (m != mxGetM(prhs[B])) || (n != mxGetN(prhs[B])))
    {
        mexErrMsgTxt("R, G, and B matrices must be the same size");
    }
    
    if (!mxIsUint8(prhs[R]) || mxIsComplex(prhs[R]) ||
        !mxIsUint8(prhs[G]) || mxIsComplex(prhs[G]) ||
        !mxIsUint8(prhs[B]) || mxIsComplex(prhs[B]))
    {
        mexErrMsgTxt("R, G, and B must be real uint8 matrices");
    }
    
    p = mxGetM(prhs[MAP]);
    q = mxGetN(prhs[MAP]);
    if (p < 2)
    {
        mexErrMsgTxt("Colormap must have at least two entries");
    }
    if (p > 65536)
    {
        mexErrMsgTxt("Colormap must not have more than 65536 entries");
    }
    if (q != 3)
    {
        mexErrMsgTxt("Colormap must be an M-by-3 matrix");
    }
    if (!mxIsDouble(prhs[MAP]) || mxIsComplex(prhs[MAP]))
    {
        mexErrMsgTxt("Colormap must be real");
    }
    
    if (!mxIsDouble(prhs[QM]) || (mxGetNumberOfElements(prhs[QM]) != 1) ||
        mxIsComplex(prhs[QM]))
    {
        mexErrMsgTxt("Qm must be a real double scalar");
    }
    
    if (!mxIsDouble(prhs[QE]) || (mxGetNumberOfElements(prhs[QE]) != 1) ||
        mxIsComplex(prhs[QE]))
    {
        mexErrMsgTxt("Qe must be a real double scalar");
    }
    
    *qm = (int32_T) mxGetScalar(prhs[QM]);
    *qe = (int32_T) mxGetScalar(prhs[QE]);

    if (*qm < 1)
    {
        mexErrMsgTxt("Qm must be at least 1");
    }
    if (*qm > ((sizeof(int32_T) * 8) / 3))
    {
        mexErrMsgTxt("Qm must be no greater than 10");
    }
    
    if (*qe < 1)
    {
        mexErrMsgTxt("Qe must be at least 1");
    }
    if (*qe > (sizeof(int32_T) * 8 - 1))
    {
        mexErrMsgTxt("Qe must be no greater than 31");
    }

    map = mxGetPr(prhs[MAP]);
    for (i = 0; i < p*q; i++)
    {
        if ((map[i] < 0.0) || (map[i] > 1.0))
        {
            mexErrMsgTxt("Colormap contains out-of-range values");
        }
    }
    
    /* hey, if we made it all the way here, everything's ok */
}
    
 
void mexFunction(
	int nlhs,
	mxArray *plhs[],
	int nrhs,
	const mxArray *prhs[])
{
    /* Declare input variable data. */

    uint8_T     *r, *g, *b;
    real64_T	*map;
    int32_T	qm, qe;
    int32_T     out_size[2];
        
    /* Declare output variable data. */

    mxArray     *result;
    uint16_T   *dithered_image;

    /* Declare working storage data. */

    uint32_T *scratchpad;
    int32_T *inverse_colormap;
    int32_T m, n, k, x, elevels, mlevels;
    int32_T *rq, *gq, *bq, *mapq, *tab1, *tab3, *tab5, *tab7;
    int32_T i;
    uint8_T *pu8;
    real64_T *pr;
    /*
     * 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 real64_T z;

    /* initialize flops counter */
    flops = 0;
    
    /* Check input arguments. */
    ValidateInputs(nlhs, plhs, nrhs, prhs, &qm, &qe);

    elevels = 1 << maxx(qm, qe);
    mlevels = 1 << qm;

    m = mxGetM(prhs[R]);
    n = mxGetN(prhs[R]);

    map = mxGetPr(prhs[MAP]);
    k = mxGetM(prhs[MAP]);

    /* create output array */
    out_size[0] = m;
    out_size[1] = n;
    if (k <= 256)
    {
        result = mxCreateNumericArray(2, out_size, mxUINT8_CLASS, mxREAL);
    }
    else
    {
        result = mxCreateNumericArray(2, out_size, mxDOUBLE_CLASS, mxREAL);
    }
    
    if (mxIsEmpty(result))
    {
        /* nothing to do */
        plhs[0] = result;
        return;
    }

    r = (uint8_T *) mxGetData(prhs[R]);
    g = (uint8_T *) mxGetData(prhs[G]);
    b = (uint8_T *) mxGetData(prhs[B]);

    /* Allocate working memory. */
    dithered_image = mxCalloc(mxGetNumberOfElements(result), 
                              sizeof(*dithered_image));
    inverse_colormap = mxCalloc( mlevels*mlevels*mlevels, 
                                 sizeof(*inverse_colormap) );
    scratchpad = mxCalloc( mlevels*mlevels*mlevels, 
                           sizeof(*scratchpad) );
    rq = mxCalloc( m*n, sizeof( *rq ) );
    gq = mxCalloc( m*n, sizeof( *gq ) );
    bq = mxCalloc( m*n, sizeof( *bq ) );
    mapq = mxCalloc( k*3, sizeof( *mapq ) );
    tab1 = mxCalloc( elevels*2, sizeof( *tab1 ) );
    tab3 = mxCalloc( elevels*2, sizeof( *tab3 ) );
    tab5 = mxCalloc( elevels*2, sizeof( *tab5 ) );
    tab7 = mxCalloc( elevels*2, sizeof( *tab7 ) );

    for (x = 0; x < 3*k; x++) {
        z = map[x] * (elevels-1) + 0.5;
        mapq[x] = (int32_T) z;

    }
    flops += 9 * k;
 
    /* Call algorithms. */

    inv_cmap_2( k, mapq, qm, maxx(qm,qe), scratchpad, inverse_colormap, 
		&flops );


    fsdither( r, g, b, mapq, k, inverse_colormap, dithered_image, 
		m, n, qm, qe, rq, gq, bq, 
		tab1, tab3, tab5, tab7 );

    if (k <= 256)
    {
        /* uint8 zero-based output */
        pu8 = (uint8_T *) mxGetData(result);
        for (i = 0; i < m*n; i++)
        {
            *pu8++ = (uint8_T) dithered_image[i];
        }
    }
    else
    {
        /* double one-based output */
        pr = (real64_T *) mxGetData(result);
        for (i = 0; i < m*n; i++)
        {
            *pr++ = (real64_T) dithered_image[i] + 1;
        }
    }
    
    plhs[0] = result;
    
    /* Report flop count. */
    mexAddFlops(flops);

    /* Free working memory */
    mxFree(tab7);
    mxFree(tab5);
    mxFree(tab3);
    mxFree(tab1);
    mxFree(bq);
    mxFree(gq);
    mxFree(rq);
    mxFree(scratchpad);
    mxFree(inverse_colormap);
    mxFree(dithered_image);
    mxFree(mapq);

} /* mexFunction() */