function [A_nm,zmlist,cidx,V_nm] = zernike(img,n,m)if nargin>0if nargin==1n = 0;endimg=imread('A.png');d = size(img);img = double(img);xstep = 2/(d(1)-1); %每一行的步長ystep = 2/(d(2)-1);[x,y] = meshgrid(-1:xstep:1,-1:ystep:1);% 設置X為從-1到1且變換步長為xstepcircle1 = x.^2 + y.^2;inside = find(circle1<=1);% 存貯滿足該條件的所有X,Y的值mask = zeros(d);mask(inside) = ones(size(inside));[cimg,cidx] = clipimg(img,mask);z = clipimg(x+i*y,mask);p = 0.9*abs(z);theta = angle(z);c = 1;for order=1:length(n)n1 = n(order);if nargin<3m = zpossible(n1);endfor r=1:length(m)V_nmt = zpoly(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==4V_nm(:,c) = V_nmt;endc = c+1;endendelseendfunction [cimg,cindex,dim] = clipimg(img,mask)dim = size(img);cindex = find(mask~=0);cimg = img(cindex);return;function [m] = zpossible(n)if iseven(n)m = 0:2:n;elsem = 1:2:n;endreturn;function [V_nm,mag,phase] = zpoly(n,m,p,theta)R_nm = zeros(size(p));a = (n+abs(m))/2;b = (n-abs(m))/2;total = b;for s=0:totalnum = ((-1)^s)*fac(n-s)*(p.^(n-2*s));den = fac(s)*fac(a-s)*fac(b-s);R_nm = R_nm + num/den;endmag = 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:2cand = find(findex>2);candidates = findex(cand);findex(cand) = candidates-1;factorial(cand) = factorial(cand).*findex(cand);endreturn;function [verdict] = iseven(candy)verdict = zeros(size(candy));isint = find(isint(candy)==1);divided2 = candy(isint)/2;evens = (divided2==floor(divided2));verdict(isint) = evens;return;function [verdict] = isint(candy)verdict = double(round(candy))==candy;return;
1 回答
陪伴而非守候
TA貢獻1757條經驗 獲得超8個贊
% [A_nm,zmlist,cidx,V_nm] = zernike(img,n,m) 計算Zernike矩
% V_nm
為n階的Zernike多項式，定義為在極坐標系中p，theta的函數
% cidx 表示虛部值
% A_nm
為zernike矩
function [A_nm,zmlist,cidx,V_nm] =
zernike(img,n,m)
if nargin>0
if nargin==1
n=0;
end
d=size(img);
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); %
提取符合circle1<=1的數
mask=zeros(d); %
構造size（d）*size(d)的矩陣。
mask(inside)=ones(size(inside)); %
構造size(inside)*size(inside)的全為1的矩陣賦值給mask（inside）
[cimg,cidx]=clipimg(img,mask);
% 計算圖像的復數表示
int i;
z=clipimg(x+i*y,mask); % 計算Z的實部和虛部
p=0.9*abs(z); %
計算復數的模，sqrt(x,y),z=x+iy;
theta=angle(z); % 計算復數z的輻角值（tanz）
c=1;
for
order=1:length(n)
n1=n(order);
% if input arguments less than
3
if nargin<3
m=zpossible(n1);
end
for r=1:length(m)
V_nmt=zpoly(n1,m(r),p,theta); % V_nm為n階的Zernike多項式，定義為在極坐標系中p，theta的函數
zprod=cimg.*conj(V_nmt); % conj是求復數的共軛
A_nm(c)=(n1+1)*sum(sum(zprod))/pi; % (n1+1)/π*∑∑(zprod); 對于圖像而言求和代替了求積分
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;
%
判斷n是偶數還是奇數，是偶數時，m取0,2,4,6等,否則取奇數賦值m
function [m]=zpossible(n)
if
iseven(n)
m=0:2:n;
else
m=1:2:n;
end
return;
%
計算Zernike矩多項式
function
[V_nm,mag,phase]=zpoly(n,m,p,theta)
R_nm=zeros(size(p)); %
產生size(p)*size(p)的零矩陣賦給R_nm
a=(n+abs(m))/2;
b=(n-abs(m))/2;
total=b;
for
s=0:total
num=((-1)^s)*fac(n-s)*(p.^(n-2*s)); %
(-1).-1*(n-s)!r.^(n-2*s)
den=fac(s)*fac(a-s)*fac(b-s); %
s!*(a-s)!*(b-s)!
R_nm=R_nm + num/den; %
R_nm是一個實數值的徑向多項式
end
mag=R_nm; %
賦值
phase=m*theta;
V_nm=mag.*exp(i*phase); %
V_nm為n階的Zernike多項式，定義為在極坐標系中p，theta的函數
return;
% 求n的階乘
function
[factorial]=fac(n)
maxno=max(max(n));
zerosi=find(n<=0);
%取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;
function
[verdict]=iseven(candy)
verdict=zeros(size(candy));
isint=find(isint(candy)==1);
divided2=candy(isint)/2;
evens=(divided2==floor(divided2));
verdict(isint)=evens;
return;
function
[verdict]=isint(candy)
verdict =
double(round(candy))==candy;
return;
添加回答
舉報
0/150
提交
取消