ordf.c

A toolbox for the non-local means algorithm

#include "ordf_ordfilt2.h"

#define TYPE int8_T
#define SELECT select_int8
void ordfilt2_int8
#include "ordf_ordfilt2.h"

#define TYPE int16_T
#define SELECT select_int16
void ordfilt2_int16
#include "ordf_ordfilt2.h"

#define TYPE int32_T
#define SELECT select_int32
void ordfilt2_int32
#include "ordf_ordfilt2.h"

#define TYPE mxLogical
#define SELECT select_logical
void ordfilt2_logical
#include "ordf_ordfilt2.h"

#define TYPE double
#define SELECT select_double
void ordfilt2_double
#include "ordf_ordfilt2.h"

#define TYPE float
#define SELECT select_single
void ordfilt2_single
#include "ordf_ordfilt2.h"

/* method which handles additive offsets */
#define TYPE double
#define SELECT select_double
#define ADD_OFFSET
void ordfilt2_add_offset
#include "ordf_ordfilt2.h"

/*
 * ordfilt_hist_uint8
 *
 * Implements histogram approach to computing order statistic.
 */
#define TYPE uint8_T
#define MAX_VAL_TYPE MAX_uint8_T
#define MIN_VAL_TYPE MIN_uint8_T
void ordfilt_hist_uint8
#include "ordf_hist.h"

#define TYPE int8_T
#define MAX_VAL_TYPE MAX_int8_T
#define MIN_VAL_TYPE MIN_int8_T
#define SIGNED
void ordfilt_hist_int8
#include "ordf_hist.h"

#define TYPE uint16_T
#define MAX_VAL_TYPE MAX_uint16_T
#define MIN_VAL_TYPE MIN_uint16_T
void ordfilt_hist_uint16
#include "ordf_hist.h"

#define TYPE int16_T
#define MAX_VAL_TYPE MAX_int16_T
#define MIN_VAL_TYPE MIN_int16_T
#define SIGNED
void ordfilt_hist_int16
#include "ordf_hist.h"


/* These defines specify the size of the domain matrix beyond which
   the histogram approach (faster for a larger domain) of calculating
   the order-statistic is used instead of quicksort style method.
   These values were determined experimentally.
*/
#define HIST8_REQ_ROWS 7
#define HIST8_REQ_COLS 1

#define HIST16_REQ_ROWS     3
#define HIST16_REQ_ELEMENTS 520

/* This method determines if it's more efficient to use the 
   histogram based filtering in case of INT8 and UINT8
   data type.
*/
bool useHistMethodForINT8(int numOffsets, int Md, int Nd)
{
    /* Invoke histogram based method only if a full domain 
       matrix was specified and it contained more than a specified
       number of rows and columns.
    */
    if( (numOffsets == Md*Nd) && 
        (Md >= HIST8_REQ_ROWS) && 
        (Nd >= HIST8_REQ_COLS) )
    {    
        return true;
    }
    else
    {
        return false;
    }
} 

/* This method determines if it's more efficient to use the 
   histogram based filtering in case of INT16 and UINT16
   data type.
*/
bool useHistMethodForINT16(int numOffsets, int Md, int Nd)
{
    /* Invoke histogram based method only if a full domain 
       matrix was specified and it contained more than a specified
       number of rows and more than a specified number of elements.
    */
    if( (numOffsets == Md*Nd) && 
        (Md >= HIST16_REQ_ROWS) && 
        (numOffsets >= HIST16_REQ_ELEMENTS) )
    {
        return true;
    }
    else
    {
        return false;
    }
} 

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    int startRow;  /* zero-based */
    int startCol;  /* zero-based */
    int order;     /* zero-based */
    int Mb;
    int Nb;
    int Md;
    int Nd;
    int Ma;
    mxClassID classA;
    int sizeB[2];
    int *offsets;
    int numOffsets;
    const mxArray *A;
    mxArray *B;
    double *add;

    ValidateInputs(nlhs, plhs, nrhs, prhs, &A, &startRow, &startCol,
                   &order, &offsets, &numOffsets, &add, &Mb, &Nb,
                   &Md, &Nd);

    /* Create output */
    classA = mxGetClassID(A);
    sizeB[0] = Mb;
    sizeB[1] = Nb;

    if (mxIsLogical(A))
    {
        B = mxCreateLogicalArray(2, sizeB);
    }
    else
    {
        B = mxCreateNumericArray(2, sizeB, classA, mxREAL);
    }
    
    Ma = mxGetM(A);

    switch (classA)
    {
      case mxUINT8_CLASS:
        if( useHistMethodForINT8(numOffsets, Md, Nd) )
        {
            ordfilt_hist_uint8((uint8_T *)mxGetData(A), 
                               (uint8_T *)mxGetData(B),
                               startRow, startCol, Mb, Nb,
                               Ma, order, offsets,
                               numOffsets, Md, Nd);
        }
        else
        {
            ordfilt2_uint8((uint8_T *) mxGetData(A), (uint8_T *) mxGetData(B),
                           startRow, startCol, Mb, Nb, Ma, order, offsets,
                           numOffsets);
        }
        break;
        
      case mxUINT16_CLASS:
        if( useHistMethodForINT16(numOffsets, Md, Nd) )
        {
            ordfilt_hist_uint16((uint16_T *)mxGetData(A), 
                                (uint16_T *)mxGetData(B),
                                startRow, startCol, Mb, Nb,
                                Ma, order, offsets,
                                numOffsets, Md, Nd);
        }
        else
        {
            ordfilt2_uint16((uint16_T *) mxGetData(A), 
                            (uint16_T *) mxGetData(B),
                            startRow, startCol, Mb, Nb, Ma, 
                            order, offsets, numOffsets);
        }
        
        break;

      case mxUINT32_CLASS:
        ordfilt2_uint32((uint32_T *) mxGetData(A), (uint32_T *) mxGetData(B),
                        startRow, startCol, Mb, Nb, Ma, order, offsets,
                        numOffsets);
        break;

      case mxINT8_CLASS:
        if( useHistMethodForINT8(numOffsets, Md, Nd) )
        {
            ordfilt_hist_int8((int8_T *)mxGetData(A), 
                              (int8_T *)mxGetData(B),
                              startRow, startCol, Mb, Nb,
                              Ma, order, offsets,
                              numOffsets, Md, Nd);
        }
        else
        {
            ordfilt2_int8((int8_T *) mxGetData(A), (int8_T *) mxGetData(B),
                          startRow, startCol, Mb, Nb, Ma, order, offsets,
                          numOffsets);
        }
        break;
        
      case mxINT16_CLASS:
        if( useHistMethodForINT16(numOffsets, Md, Nd) )
        {
            ordfilt_hist_int16((int16_T *)mxGetData(A), 
                               (int16_T *)mxGetData(B),
                               startRow, startCol, Mb, Nb,
                               Ma, order, offsets,
                               numOffsets, Md, Nd);
        }
        else
        {
            ordfilt2_int16((int16_T *) mxGetData(A), 
                           (int16_T *) mxGetData(B),
                           startRow, startCol, Mb, Nb, Ma, 
                           order, offsets, numOffsets);
        }

        break;

      case mxINT32_CLASS:
        ordfilt2_int32((int32_T *) mxGetData(A), (int32_T *) mxGetData(B),
                       startRow, startCol, Mb, Nb, Ma, order, offsets,
                       numOffsets);
        break;

      case mxLOGICAL_CLASS:
        ordfilt2_logical(mxGetLogicals(A), mxGetLogicals(B),
                         startRow, startCol, Mb, Nb, Ma, order, offsets,
                         numOffsets);
        break;
      
      case mxDOUBLE_CLASS:
        if(add==NULL)
	{
	    ordfilt2_double((double *) mxGetData(A), (double *) mxGetData(B),
			    startRow, startCol, Mb, Nb, Ma, order, offsets,
			    numOffsets);
	}
	else /* additive offsets were specified */
	{   
	    ordfilt2_add_offset((double *) mxGetData(A), 
                                (double *) mxGetData(B),
				startRow, startCol, Mb, Nb, 
                                Ma, order, offsets,
				numOffsets, add);
	}
	break;
        
      case mxSINGLE_CLASS:
        ordfilt2_single((float *) mxGetData(A), (float *) mxGetData(B),
                        startRow, startCol, Mb, Nb, Ma, order, offsets,
                        numOffsets);
        break;
  
      default:
        /* Should have errored in ValidateInputs() */
        mexErrMsgIdAndTxt("Images:ordf:invalidInputClass",
                          "%s",
                          "Invalid input class.");
    }

    mxFree((void *) offsets);

    plhs[0] = B;
}