Zernike Moments Computation and Reconstruction

[Moon2]

Hi all I have the following c++ code, using opencv, that I wrote to
calculate the Zernike moments I have defined the radial polynomial as:

//Radial polynomial
//p-|q| must be even and |m|<=n
double R(int p,int q, double rho)
{
double val=0.0;
for(int s=0; s<=((p-q)/2);s++)
{
double num = pow(-1,s)*fact(p-s);
double denom = fact(s) * fact( ((p+q)/2) -s ) *
fact( ((p-q)/2)-s );
//double denom = fact(s) * fact( (p-2*s+q)/2) *
fact( (p-2*s-q)/2);
val += (num/denom)*pow(rho,p-2*s);
}
return val;
}

Then the Zernike function to calculate the zernike moments for an NxN
image are defined as....

double Zernike::CalculateMoments(IplImage* image, int n, int m)
{
ImageSize = image->width;
uchar*imagedata = (uchar*)image->imageData;
int step=image->widthStep;

double rho, theta;
cmplx zernike(0,0);

for(int i=0;i<ImageSize;i++)
{
//map y to unit circle
for(int j=0;j<ImageSize;j++)
{
//map x, y to unit ci rcle
double y = -1 + ((2*(j+0.5))/ImageSize;
double x = -1 + ((2*(j+0.5))/ImageSize;
theta=0.0;
rho = (double)sqrt(x*x+y*y);

if(rho<=1.0)
{
if(x!=0)
theta = atan(y/x);

cmplx pixel = (int)imagedata[i*step+j]; //pixel + i0
zernike+= conj( R(n,m,rho)* polar(1.0, m*theta))*pixel //Znm
= [Rnm(r)*exp(-jmO]]*f(x,y)
}
}
}
cmplx val = (n+1)/M_PI;
zernike= zernike*val;

double result = abs(zernike/ImageSize); //return normalised
Zernike moment
}

The Zernike function doesnt return Z11 to =0, and Z00 to be 1/pi... I
dont know where the problem is there!

Now the reconstruction function is defined as:

IplImage* Zernike::Recon(int moments)
{
//create a dummy NxN image
IplImage* reconimg =
cvCreateImage(cvSize(ImageSize,ImageSize),IPL_DEPTH_8U,1);
uchar* imagedata = (uchar*)reconimg->imageData;
int step=reconimg->widthStep;

double rho, theta;
cmplx zernVal(0,0);

for(int i=0;i<ImageSize;i++)
{
for(int j=0;j<ImageSize;j++)
{
//map x, y to unit circle
double y = -1 + ((2*(j+0.5))/ImageSize;
double x = -1 + ((2*(j+0.5))/ImageSize;
theta=0.0;

rho = (double)sqrt(x*x+y*y);

//get each zernike and multiply by Vpq(x,y)
if(rho<=1.0)
{
if(x!=0)
theta = atan(y/x);

for(int n=0;n<=moments;n++)
{
if(n%2==0)
{
for(int m=0;m<=n;m+=2)
{
//ZernikeMom
array that holds Znm
zernVal+= ZernikeMom[n][m]*R(n,m,rho)* polar(1.0, n*theta);
}
}
else
{
for(int m=1;m<=n;m+=2)
{
zernVal+= ZernikeMom[n][m]*R(n,m,rho)* polar(1.0,
n*theta);
}
}
}
}
imagedata[i*step+j] = abs(zernVal);//zernike;
}
}
return reconimg;
}

Well my question boils down to I have the zernike moment function and
im trying to reconstruct the original image with no success if any one
could please take the time to read this and help I would be greatly
appreciated.

Mike

[Moon2]

sorry i have the following typedef complex<double> cmplx