www.pudn.com > P-zernike.zip > P_zernike.m, change:2006-05-23,size:1595b


function	[A_nm,zmlist,cidx,V_nm]	= P_zernike(img,n,m) 
if nargin>0 
if nargin==1 
	n	= 0; 
end 
d		= size(img);			% dimension 
img		= double(img); 
xstep		= 2/(d(1)-1); 
ystep		= 2/(d(2)-1); 
[x,y]		= meshgrid(-1:xstep:1,-1:ystep:1); 
circle1		= x.^2 + y.^2; 
inside		= find(circle1<=1); 
mask		= zeros(d); 
mask(inside)	= ones(size(inside)); 
[cimg,cidx]	= clipimg(img,mask); 
z		= clipimg(x+i*y,mask); 
p		= 0.9*abs(z);   %make sure that the p less than 0.9 according to the theory of p<=0.9; 
theta		= angle(z); 
c	= 1; 
for order=1:length(n) 
	n1	= n(order); 
	if nargin<3 
		m	= 0:n1;		 
	end 
		for r=1:length(m) 
		V_nmt		= pzpoly(n1,m(r),p,theta); 
		zprod		= cimg.*conj(V_nmt); 
		A_nm(c)		= (n1+1)*sum(sum(zprod))/pi; 
		zmlist(c,1:2)	= [n1 m(r)]; 
		  if nargout==4 
			V_nm(:,c)	= V_nmt; 
		   end 
		    c		= c+1; 
	    end 
end 
else 
end 
 
function [cimg,cindex,dim]	= clipimg(img,mask) 
 
dim	= size(img); 
cindex	= find(mask~=0); 
cimg	= img(cindex); 
return; 
 
function	[V_nm,mag,phase]	= pzpoly(n,m,p,theta) 
R_nm	= zeros(size(p)); 
a	= n+abs(m)+1; 
b	= n-abs(m); 
total	= b; 
for s=0:total 
	num	= ((-1)^s)*fac(2*n+1-s)*(p.^(n-s)); 
	den	= fac(s)*fac(a-s)*fac(b-s); 
	R_nm	= R_nm + num/den; 
end 
mag	= R_nm; 
phase	= m*theta; 
V_nm	= mag.*exp(i*phase); 
return; 
 
function [factorial]	= fac(n) 
maxno		= max(max(n)); 
zerosi		= find(n<=0); 
n(zerosi)	= ones(size(zerosi)); 
factorial	= n; 
findex		= n; 
for i=maxno:-1:2 
	cand		= find(findex>2); 
	candidates	= findex(cand); 
	findex(cand)	= candidates-1; 
	factorial(cand)	= factorial(cand).*findex(cand); 
end 
return;