applylut.c

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

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

/*

   applylut.c  .MEX file
               Applies a 16-by-1 lookup table to 2-by-2 sliding
               neighborhoods in a binary image; or applies a
               512-by-1 lookup table to 3-by-3 sliding neighborhoods
               in a binary image.

               For 2-by-2 neighborhoods, the neighborhood examined
               for each pixel is the lower right one, and lookup
               indices are determined by applying this mask:

                        1 4
                        2 8

               For 3-by-3 neighborhoods, the lookup indices are
               determined by applying this mask:

                        1  8  64
                        2 16 128
                        4 32 256

               A = APPLYLUT(BW, LUT)

               BW and LUT and be either numeric double or numeric
               uint8.  BW and LUT must be real.

               Case 1: LUT is uint8
                 Output will be uint8

               Case 2: LUT is double
                 If all LUT values are integers between 0 and 255,
                 output will be uint8; otherwise output will be
                 double.           

               If output contains only 0's and 1's, it will be logical.

               Steven L. Eddins
               March 1996

*/

static char rcsid[] = "$Revision: 1.8 $";

#include "mex.h"

#define MATRIX_REF(PR, NUMROWS, R, C) \
   (*((PR) + (NUMROWS)*(C) + (R)))

mxArray *ConvertToDouble(const mxArray *A)
{
    int nrhs = 1;
    int nlhs = 1;
    mxArray *prhs[1];
    mxArray *plhs[1];
    mxArray *result;
    int errorFlag;

    if (A == NULL)
    {
        mexErrMsgTxt("applylut internal error 3");
    }

    prhs[0] = (mxArray *) A;
    errorFlag = mexCallMATLAB(nlhs, plhs, nrhs, prhs, "double");
    if (errorFlag == 0) {
        /* Successful operation */
        result = plhs[0];
    } else {
        /* double conversion failed */
        /* We should never get here; error should have long-jumped */
        /* back to the command prompt */
        mexErrMsgTxt("applylut internal error 4");
    }

    return(result);

}

mxArray *ConvertToUint8(const mxArray *A)
{
    int nrhs = 1;
    int nlhs = 1;
    mxArray *prhs[1];
    mxArray *plhs[1];
    mxArray *result;
    int errorFlag;

    if (A == NULL)
    {
        mexErrMsgTxt("applylut internal error 1");
    }

    prhs[0] = (mxArray *) A;
    errorFlag = mexCallMATLAB(nlhs, plhs, nrhs, prhs, "uint8");
    if (errorFlag == 0) {
        /* Successful operation */
        result = plhs[0];
    } else {
        /* double conversion failed */
        /* We should never get here; error should have long-jumped */
        /* back to the command prompt */
        mexErrMsgTxt("applylut internal error 2");
    }

    return(result);

}

static int16_T weights2[2][2] = {1, 4, 2, 8};

int16_T Nhood2Offset(uint8_T *pBWin, int numRows, int numCols, 
                     int r, int c) {

    int maxR, maxC;
    int rr, cc;
    int16_T result = 0;

    /* Figure out the neighborhood extent that does not go past */
    /* image boundaries */
    if (r == (numRows-1)) {
        maxR = 0;
    } else {
        maxR = 1;
    }
    if (c == (numCols-1)) {
        maxC = 0;
    } else {
        maxC = 1;
    }

    for (rr = 0; rr <= maxR; rr++) {
        for (cc = 0; cc <= maxC; cc++) {
            result += weights2[rr][cc] *
              (MATRIX_REF(pBWin, numRows, r+rr, c+cc) != 0);
        }
    }

    return(result);
}

static int16_T weights3[3][3] = {1, 8, 64, 2, 16, 128, 4, 32, 256};

int16_T Nhood3Offset(uint8_T *pBWin, int numRows, int numCols,
                     int r, int c) {

    int minR, maxR, minC, maxC;
    int rr, cc;
    int16_T result = 0;
    
    /* Figure out the neighborhood extent that does not go past */
    /* image boundaries */
    if (r == 0) {
        minR = 1;
    } else {
        minR = 0;
    }
    if (r == (numRows-1)) {
        maxR = 1;
    } else {
        maxR = 2;
    }
    if (c == 0) {
        minC = 1;
    } else {
        minC = 0;
    }
    if (c == (numCols-1)) {
        maxC = 1;
    } else {
        maxC = 2;
    }

    for (rr = minR; rr <= maxR; rr++) {
        for (cc = minC; cc <= maxC; cc++) {
            result += weights3[rr][cc] * 
              (MATRIX_REF(pBWin, numRows, r + rr - 1, c + cc - 1) != 0);
        }
    }

    return(result);
}

void Compute2by2DoubleResult(mxArray *BWout, 
                             mxArray *BWin, 
                             mxArray *lut) {
    
    int numRows, numCols;
    int r, c;
    uint8_T *pBWin;
    double *plut;
    double *pBWout;

    pBWin = (uint8_T *) mxGetData(BWin);
    plut = (double *) mxGetData(lut);
    pBWout = (double *) mxGetData(BWout);

    numRows = mxGetM(BWin);
    numCols = mxGetN(BWin);

    for (c = 0; c < numCols; c++) {
        for (r = 0; r < numRows; r++) {
            MATRIX_REF(pBWout, numRows, r, c) =
              *(plut + Nhood2Offset(pBWin, numRows, numCols, r, c));
        }
    }
}


void Compute3by3DoubleResult(mxArray *BWout, 
                             mxArray *BWin, 
                             mxArray *lut) {
    
    int numRows, numCols;
    int r, c;
    uint8_T *pBWin;
    double *plut;
    double *pBWout;

    pBWin = (uint8_T *) mxGetData(BWin);
    plut = (double *) mxGetData(lut);
    pBWout = (double *) mxGetData(BWout);

    numRows = mxGetM(BWin);
    numCols = mxGetN(BWin);

    for (c = 0; c < numCols; c++) {
        for (r = 0; r < numRows; r++) {
            MATRIX_REF(pBWout, numRows, r, c) =
              *(plut + Nhood3Offset(pBWin, numRows, numCols, r, c));
        }
    }
}

void Compute2by2Uint8Result(mxArray *BWout, 
                            mxArray *BWin, 
                            mxArray *lut) {
    
    int numRows, numCols;
    int r, c;
    uint8_T *pBWin;
    double *plut;
    uint8_T *pBWout;

    pBWin = (uint8_T *) mxGetData(BWin);
    plut = (double *) mxGetData(lut);
    pBWout = (uint8_T *) mxGetData(BWout);

    numRows = mxGetM(BWin);
    numCols = mxGetN(BWin);

    for (c = 0; c < numCols; c++) {
        for (r = 0; r < numRows; r++) {
            MATRIX_REF(pBWout, numRows, r, c) = (uint8_T)
              *(plut + Nhood2Offset(pBWin, numRows, numCols, r, c));
        }
    }
}

void Compute3by3Uint8Result(mxArray *BWout, 
                            mxArray *BWin, 
                            mxArray *lut) {
    
    int numRows, numCols;
    int r, c;
    uint8_T *pBWin;
    double *plut;
    uint8_T *pBWout;

    pBWin = (uint8_T *) mxGetData(BWin);
    plut = (double *) mxGetData(lut);
    pBWout = (uint8_T *) mxGetData(BWout);

    numRows = mxGetM(BWin);
    numCols = mxGetN(BWin);

    for (c = 0; c < numCols; c++) {
        for (r = 0; r < numRows; r++) {
            MATRIX_REF(pBWout, numRows, r, c) = (uint8_T)
              *(plut + Nhood3Offset(pBWin, numRows, numCols, r, c));
        }
    }
}

/* In addition to validating the inputs, this function */
/* converts the input BW array to uint8 and the input LUT */
/* array to double if necessary. */
void ValidateInputs(int nlhs, int nrhs, const mxArray *prhs[],
                    mxArray **BW, bool *freeBW, mxArray **lut, 
                    bool *freelut, mxClassID *outputClass, 
                    bool *outputIsLogical, int *nhoodSize) {

    int i;
    int lutLength;
    double *plut;   /* pointer to lut array */

    if (nrhs != 2) {
        mexErrMsgTxt("Two inputs required");
    }

    for (i = 0; i < nrhs; i++) {
        if (mxIsComplex(prhs[i])) {
            mexErrMsgTxt("Inputs must be real");
        }
        if (! mxIsNumeric(prhs[i])) {
            mexErrMsgTxt("Inputs must be numeric");
        }
        if ((! mxIsDouble(prhs[i])) && (! mxIsUint8(prhs[i]))) {
            mexErrMsgTxt("Inputs must be double or uint8");
        }
        if (mxGetNumberOfDimensions(prhs[i]) > 2) {
            mexErrMsgTxt("Inputs must be 2-D");
        }
    }

    if (nlhs > 1) {
        mexErrMsgTxt("Two many output arguments");
    }

    /* Determine output type */
    *outputClass = mxUINT8_CLASS;
    if (mxIsDouble(prhs[1])) {
        plut = (double *) mxGetData(prhs[1]);
        for (i = 0; i < mxGetNumberOfElements(prhs[1]); i++) {
            if (!mxIsFinite(*plut) || (*plut != (uint8_T) *plut)) {
                *outputClass = mxDOUBLE_CLASS;
                 break;
            }
            plut++;
        }
    }

    if (! mxIsUint8(prhs[0])) {
        *BW = ConvertToUint8(prhs[0]);
        *freeBW = true;
    } else {
        *BW = (mxArray *) prhs[0];
        *freeBW = false;
    }

    if (! mxIsDouble(prhs[1])) {
        *lut = ConvertToDouble(prhs[1]);
        *freelut = true;
    } else {
        *lut = (mxArray *) prhs[1];
        *freelut = false;
    }

    lutLength = mxGetNumberOfElements(*lut);
    if (lutLength == 16) {
        *nhoodSize = 2;
    } else if (lutLength == 512) {
        *nhoodSize = 3;
    } else {
        mexErrMsgTxt("Invalid lookup table");
    }
    
    /* Should we set the logical flag on the output? */
    /* if output contains no values other than 0 or 1, then */
    /* set the logical flag */
    *outputIsLogical = true;
    plut = (double *) mxGetData(*lut);
    for (i = 0; i < lutLength; i++) {
        if ((*plut != 0.0) && (*plut != 1.0)) {
            *outputIsLogical = false;
            break;
        }
        plut++;
    }
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    uint16_T **indices;         /* neighborhood indices */
    mxClassID outputClass;      /* output array class */
    bool outputIsLogical;       /* Set output logical flag? */
    bool freeBW;                /* Free BW at the end? */
    bool freelut;               /* Free lut at the end? */
    mxArray *BW;                /* Input image */
    mxArray *lut;               /* Input lookup table */
    mxArray *BW2;               /* Output image */
    int nhoodSize;              /* Neighborhood size; 2 or 3 */

    /* ValidateInputs always returns a uint8 BW and a double lut */
    ValidateInputs(nlhs, nrhs, prhs, &BW, &freeBW, &lut, &freelut, 
                   &outputClass, &outputIsLogical, &nhoodSize);

    BW2 = mxCreateNumericArray(mxGetNumberOfDimensions(BW),
                               mxGetDimensions(BW), outputClass, mxREAL);
    if (outputIsLogical) {
        mxSetLogical(BW2);
    }
    
    if (! mxIsEmpty(BW2)) {
        /* Output is not empty, so we actually have to do some work. */

        if (nhoodSize == 2) {
            if (outputClass == mxDOUBLE_CLASS) {
                Compute2by2DoubleResult(BW2, BW, lut);
            } else {
                Compute2by2Uint8Result(BW2, BW, lut);
            }
        } else {
            if (outputClass == mxDOUBLE_CLASS) {
                Compute3by3DoubleResult(BW2, BW, lut);
            } else {
                Compute3by3Uint8Result(BW2, BW, lut);
            }
        }

    }

    /* If we made it this far, all is OK.  Might as well let the */
    /* left-hand side have it. */

    plhs[0] = BW2;

    if (freeBW) {
        mxDestroyArray(BW);
    }
    if (freelut) {
        mxDestroyArray(lut);
    }
}