Next: Feature Extraction - Haralick Up: HeLa Data Previous: mb_zernike.m

### mb_Znl.cpp

``` /*/////////////////////////////////////////////////////////////////////////
//
//
//                            mb_Znl.cpp
//
//
//                           Michael Boland
//                            09 Dec 1998
//
//  Revisions:
//
/////////////////////////////////////////////////////////////////////////*/

#include "mex.h"
#include "matrix.h"
#include
#include
#include

#define row 0
#define col 1

//
// Calculates n! (uses double arithmetic to avoid overflow)
//
double factorial(double n)
{
if(n < 0)
return(0.0) ;
if(n == 0.0)
return(1.0) ;
else
return(n * factorial(n-1.0)) ;
}

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

int n ;                        /* degree (n) of the Zernike moment */
int l ;                        /* angular dependence of the moment */
double* X ;                    /* list of X coordinates of pixels */
double* Y ;                    /* list of Y coordinates of pixels */
double* P ;                    /* list of values of pixels */
int outputsize[2] = {1,1} ;    /* Dimensions of return variable */
register double x, y, p ;      /* individual values of X, Y, P */
int i ;
complex sum = 0.0 ;    /* Accumulator for complex moments */
complex Vnl = 0.0 ;    /* Inner sum in Zernike calculations */
double* preal ;                /* Real part of return value */
double* pimag ;                /* Imag part of return value */
int m ;

if (nrhs != 5) {
mexErrMsgTxt("
MB_Znl(N, L, X, Y, P),
Zernike moment generating function.  The moment of degree n and
angular dependence l for the pixels defined by coordinate vectors
X and Y and intensity vector P.  X, Y, and P must have the same
length.") ;

} else if (nlhs != 1) {
mexErrMsgTxt("mb_Znl returns a single output.\n") ;
}

if ( !mxIsNumeric(prhs[0]) || (mxGetM(prhs[0]) != 1) ||
(mxGetN(prhs[0]) != 1) ) {
mexErrMsgTxt("The first argument (n) should be a scalar\n") ;
}

if ( !mxIsNumeric(prhs[1]) || (mxGetM(prhs[1]) != 1) ||
(mxGetN(prhs[1]) != 1) ) {
mexErrMsgTxt("The second argument (l) should be a scalar\n") ;
}

if ( !mxIsNumeric(prhs[2]) || (mxIsComplex(prhs[2])) ) {
mexErrMsgTxt("The 3d argument (X) should be numeric and not complex.");
}

if ( !mxIsNumeric(prhs[3]) || (mxIsComplex(prhs[3])) ) {
mexErrMsgTxt("The 3d argument (Y) should be numeric and not complex.");
}

if ( !mxIsNumeric(prhs[4]) || (mxIsComplex(prhs[4])) ) {
mexErrMsgTxt("The 3d argument (P) should be numeric and not complex.");
}

if (mxGetM(prhs[2])!=mxGetM(prhs[3]) ||
(mxGetM(prhs[3])!=mxGetM(prhs[4]))){
mexErrMsgTxt("X, Y, and P must have the same number of rows.") ;
}

if (mxGetM(prhs[2]) < mxGetM(prhs[2])) {
mexErrMsgTxt("X, Y, and P should be column vectors.") ;
}

n = (int)mxGetScalar(prhs[0]) ;
l = (int)mxGetScalar(prhs[1]) ;

X = mxGetPr(prhs[2]) ;
Y = mxGetPr(prhs[3]) ;
P = mxGetPr(prhs[4]) ;

for(sum = 0.0, i = 0 ; i < mxGetM(prhs[2]) ; i++) {
x = X[i] ;
y = Y[i] ;
p = P[i] ;

for(Vnl = 0.0, m = 0; m <= (n-l)/2; m++) {
Vnl += (pow((double)-1.0,(double)m)) * ( factorial(n-m) ) /
( factorial(m) * (factorial((n - 2.0*m + l) / 2.0)) *
(factorial((n - 2.0*m - l) / 2.0)) ) *
( pow( sqrt(x*x + y*y), (double)(n - 2*m)) ) *
polar(1.0, l*atan2(y,x)) ;
//
// NOTE: This function did not work with the following:
//  ...pow((x*x + y*y), (double)(n/2 -m))...
//  perhaps pow does not work properly with a non-integer
//  second argument.
// 'not work' means that the output did not match the 'old'
//  Zernike calculation routines.
//
}

sum += p * conj(Vnl) ;
}

sum *= (n+1)/PI ;

/* Assign the returned value */

plhs[0] = mxCreateNumericArray(2, outputsize, mxDOUBLE_CLASS, mxCOMPLEX) ;
preal = mxGetPr(plhs[0]) ;
pimag = mxGetPi(plhs[0]) ;

preal[0] = real(sum) ;
pimag[0] = imag(sum) ;

}
```