vmquantc.c

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

    bshift = qr + qg - bshift;

    for (box = 0; box < k; box++) {
        for (z = boxll[box][B]; z <= boxur[box][B]-1; z++) {
            for (y = boxll[box][G]; y <= boxur[box][G]-1; y++) {
                for (x = boxll[box][R]; x <= boxur[box][R]-1; x++) {
                    m2[x + (y << qr) + (z << (qr+qg))] = box;
                } /* for all R in box */
            } /* for all G in box */
        } /* for all B in box */
#if !Q_CORRECTION
        /* Assign color map entry to mean of quantized values. */
        /* (Inefficient?) */

        wbox = region(m0);
        if (!(m1x[R] = rlevels - 2)) {
            m1x[R] = 1;
        }
        if (!(m1x[G] = glevels - 2)) {
            m1x[G] = 1;
        }
        if (!(m1x[B] = blevels - 2)) {
            m1x[B] = 1;
        }
        map[elem2d(box,R)] = ((region(m1r) / wbox) - 2) / 
                                (real64_T) m1x[R];
        map[elem2d(box,G)] = ((region(m1g) / wbox) - 2) / 
                                (real64_T) m1x[G];
        map[elem2d(box,B)] = ((region(m1b) / wbox) - 2) / 
                                (real64_T) m1x[B];
        flops += 3;
#endif /* !Q_CORRECTION */
    } /* for each box build inverse color map */

    /* Build image from inverse color map. 
       (This step is not necessary if !Q_CORRECTION and dithering, but 
       is done anyway) */

    for (x = 0; x < m*n; x++) {
#if Q_CORRECTION
        image[x] = (uint16_T) m2[ (rq[x] >> rshift) + 
                                ((gq[x] & gmask) << gshift) +
                                ((bq[x] & bmask) << bshift) ];
#else
        image[x] = (uint16_T) m2[ (rq[x] >> rshift) + 
                                 ((gq[x] & gmask) << gshift) +
                                 ((bq[x] & bmask) << bshift) + 1];
#endif /* Q_CORRECTION */
    }
    
    /* If Q_CORRECTION, calculate corrected color map. */
    
#if Q_CORRECTION
    
    dim2x = k;          /* Change array size for manual indexing macros. */
    
    for (x = 0; x < k; x++) {
        E[x] = 0;
    }
    for (x = 0; x < m*n; x++) {
        map[elem2d(y = (int32_T) image[x]++, R)] += (real64_T) r[x] / 255.0;
        map[elem2d(y, G)] += (real64_T) g[x] / 255.0;
        map[elem2d(y, B)] += (real64_T) b[x] / 255.0;
        E[y] += 1.0;
    }

    flops += m * n * 7;

    for (x = 0; x < k; x++) {
        map[elem2d(x,R)] /= E[x];
        map[elem2d(x,G)] /= E[x];
        map[elem2d(x,B)] /= E[x];
    }

    mxFree(boxll);
    mxFree(boxur);

    flops += k * 3;
    
#endif /* Q_CORRECTION */

    return k;

} /* vmquant() */


/*
 *      Floyd-Steinberg Dithering Algorithm.
 *      (This one is not exactly the same as in ditherc.c
 */

void fsdither( int32_T *rq, int32_T *gq, int32_T *bq, int32_T k, 
               int32_T qr, int32_T qg, int32_T qb, int32_T qe, 
               uint16_T *image, real64_T *map, 
               int32_T *mapq, int32_T m, int32_T n, int32_T *inverse_colormap, 
               int32_T *tab1, int32_T *tab3, int32_T *tab5, int32_T *tab7 )
{
    int32_T x, y, z, pixel_value, elevels = (1 << qe);
    int32_T dim2x = k, dim3x, dim3xy;
    int32_T rshift = qe - qr, gshift = qr - (qe - qg), 
        bshift = qr + qg - (qe - qb);
    int32_T gmask = ~((1 << (qe-qg))-1), bmask = ~((1 << (qe-qb))-1);
    int32_T rval, gval, bval;

    dim3x = (1 << qr) + 1;
    dim3xy = dim3x * ((1 << qg) + 1);

    /* 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);
    }

    /* Quantize color map entries. */
    
    for (x = 0; x < k; x++) {
        mapq[elem2d(x,R)] = ((int32_T) (map[elem2d(x,R)] * 
                                        (elevels-1) + 0.5)) - elevels;
        mapq[elem2d(x,G)] = ((int32_T) (map[elem2d(x,G)] * 
                                        (elevels-1) + 0.5)) - elevels;
        mapq[elem2d(x,B)] = ((int32_T) (map[elem2d(x,B)] * 
                                        (elevels-1) + 0.5)) - elevels;
    }
    flops += 6*k;

    for (x = z = 0; x < n; x++) {

        for (y = 0; y < m; y++, z++) {

            /* Place this pixel. */

            rq[z] = clamp( rq[z], 0, elevels-1 );
            gq[z] = clamp( gq[z], 0, elevels-1 );
            bq[z] = clamp( bq[z], 0, elevels-1 );

            pixel_value = inverse_colormap[ (rq[z] >> rshift) +
                        ((gq[z] & gmask) << gshift) + 
                        ((bq[z] & bmask) << bshift) ];
            image[z] = (uint16_T) pixel_value;

            rval = rq[z] -= mapq[elem2d(pixel_value,R)];
            gval = gq[z] -= mapq[elem2d(pixel_value,G)];
            bval = bq[z] -= mapq[elem2d(pixel_value,B)];

            /* 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];
                }
            }
        }
    }

}

void ValidateInputs(int32_T nlhs,
                    mxArray *plhs[],
                    int32_T nrhs,
                    const mxArray *prhs[],
                    int32_T *dither,
                    uint32_T *qe, uint32_T *qr, uint32_T *qg, uint32_T *qb,
                    uint32_T *k)
{
    int32_T m;
    int32_T n;
    real64_T *pr;
        
    if (nrhs != 7)
    {
        mexErrMsgTxt("vmquantc expects six input arguments");
    }
    
    if (nlhs != 2)
    {
        mexErrMsgTxt("vmquantc expects two output 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");
    }

    if (!mxIsDouble(prhs[K]))
    {
        mexErrMsgTxt("K must be a scalar double");
    }
    *k = (uint32_T) mxGetScalar(prhs[K]);
    if (*k < 2)
    {
        mexErrMsgTxt("K must be greater than 1");
    }
    if (*k > 65536)
    {
        mexErrMsgTxt("K must be no greater than 65536");
    }

    if ((mxGetNumberOfElements(prhs[Q]) != 3) || !mxIsDouble(prhs[Q]))
    {
        mexErrMsgTxt("Q must be a 3-element double vector");
    }
    pr = mxGetPr(prhs[Q]);
    *qr = (uint32_T) *pr;
    *qg = (uint32_T) *(pr+1);
    *qb = (uint32_T) *(pr+2);
    if ((*qr < 1) || (*qr > (8 * sizeof(int32_T) - 1)))
    {
        mexErrMsgTxt("Qr out of range");
    }
    if ((*qg < 1) || (*qg > (8 * sizeof(int32_T) - 1)))
    {
        mexErrMsgTxt("Qg out of range");
    }
    if ((*qb < 1) || (*qb > (8 * sizeof(int32_T) - 1)))
    {
        mexErrMsgTxt("Qb out of range");
    }
    
    if (!mxIsDouble(prhs[DITHER]))
    {
        mexErrMsgTxt("DITHER must be a double scalar");
    }
    *dither = (int32_T) (mxGetScalar(prhs[DITHER]) != 0.0);
    
    if (!mxIsDouble(prhs[QE]))
    {
        mexErrMsgTxt("Qe must be a double scalar");
    }
    *qe = (uint32_T) mxGetScalar(prhs[QE]);
    if ((*qe < 1) || (*qe > (8 * sizeof(int32_T) - 2)))
    {
        mexErrMsgTxt("Qe out of range");
    }
}


/*
 *      Gateway routine.
 */

void mexFunction( int32_T nlhs, 
                  mxArray *plhs[], 
                  int32_T nrhs, 
                  const mxArray *prhs[] )
{

    int32_T x, y;
    real64_T *dptr;

    /* Declare input variable data. */

    uint8_T *r, *g, *b;
    real64_T *q;
    uint32_T k, qr, qg, qb, qe;
    int32_T dither;
    uint32_T m, n;

    /* Declare output variable data. */

    uint16_T *image;
    real64_T *map;
    mxArray *X;
    mxArray *mapArray;
    int32_T out_size[2];
    uint8_T *pu8;
    real64_T *pr;
            
    /* Declare working storage data. */

    int32_T *rq, *gq, *bq, i;
    uint32_T *p;
    int32_T *m0, *m1r, *m1g, *m1b, *m2, *squares, flops = 0;
    uint32_T *boxll, *boxur;
    int32_T *mapq, *tab1, *tab3, *tab5, *tab7;
    int32_T rlevels, glevels, blevels, elevels;
    real64_T *E;
    
    /* Check input arguments. */
    ValidateInputs(nlhs, plhs, nrhs, prhs, &dither, &qe, &qr, &qg, &qb, &k);

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

    qe = vmMAX( qe, vmMAX( qr, vmMAX( qg, qb ) ) );
    elevels = 1 << qe;
    rlevels = (1 << qr) + 1;
    glevels = (1 << qg) + 1;
    blevels = (1 << qb) + 1;

    m = mxGetM(prhs[R]);
    n = mxGetN(prhs[R]);
    
    out_size[0] = m;
    out_size[1] = n;
    if (k <= 256)
    {
        X = mxCreateNumericArray(2, out_size, mxUINT8_CLASS, mxREAL);
    }
    else
    {
        X = mxCreateNumericArray(2, out_size, mxDOUBLE_CLASS, mxREAL);
    }
    
    if (mxIsEmpty(X))
    {
        /* we're done! */
        plhs[0] = X;
        plhs[1] = mxCreateDoubleMatrix(0, 0, mxREAL);
        return;
    }


    /* Allocate working memory. */
    /* NOTE: The order in which these arrays are allocated is important.
        Adjust the order of the mxFree calls below after changing. */

    image = mxCalloc( m*n, sizeof( *image ) );
    map = mxCalloc( 3*k, sizeof( *map ) );
    rq = mxCalloc( m*n, sizeof( *rq ) );
    gq = mxCalloc( m*n, sizeof( *gq ) );
    bq = mxCalloc( m*n, sizeof( *bq ) );
    m2 = mxCalloc( rlevels*glevels*blevels, sizeof( *m2 ) );
    E = mxCalloc( k, sizeof( *E ) );
    boxll = mxCalloc( 3*k, sizeof( *boxll ) );          
    boxur = mxCalloc( 3*k, sizeof( *boxur ) );          
    p = mxCalloc( rlevels*glevels*blevels, sizeof( *p ) );      
    m0 = mxCalloc( rlevels*glevels*blevels, sizeof( *m0 ) );
    m1r = mxCalloc( rlevels*glevels*blevels, sizeof( *m1r ) );
    m1g = mxCalloc( rlevels*glevels*blevels, sizeof( *m1g ) );
    m1b = mxCalloc( rlevels*glevels*blevels, sizeof( *m1b ) );
    squares = mxCalloc( vmMAX( rlevels, vmMAX( glevels, blevels ) ) + 1, 
                        sizeof( *squares ) ); 

    /* Call algorithm. */

    k = vmquant( r, g, b, k, qr, qg, qb, qe, image, map, m, n, 
        &rq, &gq, &bq, boxll, boxur, p, m0, m1r, 
                 m1g, m1b, m2, E, squares );

    /* Prepare return matrix now so dither can quantize. */

    mapArray = mxCreateDoubleMatrix(k, 3, mxREAL);
    dptr = mxGetPr( mapArray );
    for ( i = 3*k; i--; ) {
        *dptr++ = *map++;
    }
    map -= 3*k;
    
    if (dither) {

        /* Make allocations for dither routine. */

        mapq = mxCalloc( k*3, sizeof( int32_T ) );
        tab1 = mxCalloc( elevels*2, sizeof( int32_T ) );
        tab3 = mxCalloc( elevels*2, sizeof( int32_T ) );
        tab5 = mxCalloc( elevels*2, sizeof( int32_T ) );
        tab7 = mxCalloc( elevels*2, sizeof( int32_T ) );
        
        fsdither( rq, gq, bq, k, qr, qg, qb, qe, image, map, mapq, 
                m, n, m2, tab1, tab3, tab5, tab7 );

        mxFree( tab7 );
        mxFree( tab5 );
        mxFree( tab3 );
        mxFree( tab1 );
        mxFree( mapq );
    }

    if (k <= 256)
    {
        pu8 = (uint8_T *) mxGetData(X);
        for (i = 0; i < m*n; i++)
        {
            *pu8++ = (uint8_T) (image[i] - 1);
        }
    }
    else
    {
        pr = mxGetPr(X);
        for (i = 0; i < m*n; i++)
        {
            *pr++ = image[i];
        }
        flops += m*n;
    }

    plhs[0] = X;
    plhs[1] = mapArray;
        
    /* Report flop count. */
    mexAddFlops(flops);
    
    mxFree( squares );
    mxFree( m1b );
    mxFree( m1g );
    mxFree( m1r );
    mxFree( m0 );
    mxFree( p );
    mxFree( boxur );
    mxFree( boxll );
    mxFree( E );
    mxFree( m2 );
    mxFree( bq );
    mxFree( gq );
    mxFree( rq );
    mxFree( map );
    mxFree( image );
        
} /* mexFunction() */