zernike矩 matlab代码,3D Zernike矩计算函数 - MATLAB

% =========================================================================

% 3D Zernike 函数计算

% =========================================================================

% By Gu Jinjin 2013/09/01 - 08 程序为原创，转载请声明

%------ 参数解释 ------

%{

Input args:

n 表示3D Zernike矩的具体阶

l 表示重复度(球函数的阶), n-l 为偶数，且 0 =< l <= n

m (球函数的重复度), -l =< m <= l

X 单坐标列向量，表示单位球(US)中x的值

Y，Z 同上

Output args:

Z_nlm n阶3D Zernike 多项式的值

%}

%======================================

% 主函数

%======================================

%

function [Z_nlm] = zernikeFunction(n,l,m,X,Y,Z)

if nargin == 6

% 思想来源于文章 Precipitate shape fitting and reconstruction by means of 3D Zernike functions

v_up = (n-l)/2;

u_up = fix((l-m)/2);

v_sum = 0;

for v=0:v_up

q_nlv = calculateQ_klv(v_up,l,v);

u_sum = 0;

for u = 0:u_up

u_p1 = (-1)^(m+u)/2^(m+2*u)*combin(l,u)*combin(l-u,m+u);

a_sum = 0;

for a = 0:m

k_sum = 0;

for k=0:v

b_sum = 0;

for b = 0:k+u

b_sum = b_sum + combin(k+u,b).*X.^(2*b+a).*Y.^(2*(k+u-b)+m-a).*Z.^(2*(v-k-u)+l-m);

end

k_sum = k_sum + combin(v,k).*b_sum;

end

a_sum = a_sum + i^(m-a)*combin(m,a).*k_sum;

end

u_sum = u_sum+u_p1.*a_sum;

end

v_sum = v_sum + q_nlv.*u_sum;

end

Z_nlm = calculateC_lm(l,m).*v_sum;

% 思想来源于文章 3DZernike Descriptor for Content Based Shape Retrieval

% 6重求和计算，但重构有问题，因为论文公式表示 存在印刷错误

%{

k = (n-l)/2;

mu_up = fix((l-m)/2);

V_sum = 0;

for V=0:k

q_klV = calculateQ_klv(k,l,V);

a_sum = 0;

for a=0:V

b_sum = 0;

for b=0:V-a

u_sum = 0;

for u=0:m

mu_sum = 0;

for mu=0:mu_up

v_sum = 0;

for v=0:mu

v_sum = v_sum+combin(mu,v).*X.^(2*(v+a)+u).*Y.^(2*(mu-v+b)+m-u).*Z.^(2*(V-a-b-mu)+l-m);

end

mu_sum = mu_sum+(-1)^mu*2^(-2*mu)*combin(l,mu)*combin(l-mu,m+mu).*v_sum;

end

u_sum = u_sum+(-1)^(m-u)*combin(m,u)*i^u.*mu_sum;

end

b_sum = b_sum+combin(V-a,b).*u_sum;

end

a_sum = a_sum+combin(V,a).*b_sum;

end

V_sum = V_sum+q_klV.*a_sum;

end

Z_nlm = 2^(-m)*calculateC_lm(l,m).*V_sum;

%}

else

disp 'Wrong Input Args: No Img Info!'

end

%======================================

% 计算 C_lm 的值

%======================================

function [c_lm] = calculateC_lm(l,m)

%den = 2*sqrt(pi)*fac(l); % 分母

den = fac(l);

t = (2*l+1).*fac(l+m).*fac(l-m);

num = sqrt(t); % 分子

c_lm = num./den;

return;

%======================================

% 计算 Q_klv 的值

%======================================

function [q_klv] = calculateQ_klv(k,l,v)

p1 = (-1)^k/(2^(2*k));

p2 = sqrt((2*l+4*k+3)/3);

p3 = combin(2*k,k)*(-1)^v;

p4 = combin(k,v)*combin(2*(k+l+v)+1,2*k)/combin(k+l+v,k);

q_klv = p1*p2*p3*p4;

return;

%======================================

% 计算组合数 - Combinatorial

%======================================

% fac(n)/((fac(n-m)*fac(m))

function [combin] = combin(n,m)

num = fac(n);

den = fac(n-m).*fac(m);

combin = num./den;

return;

%======================================

% 求解n! - factorial

%======================================

% n可以为数，向量或者矩阵

function [factorial] = fac(n)

max_n = max(max(n));

zeros_n = find(n<=0); % 找到n中小于或者等于0的元素位置

n(zeros_n) = ones(size(zeros_n));

factorial = n;

n_1 = n;

% 保证n-1次乘法，1乘以任何数皆为1，所以不需要计算

for k= max_n:-1:2

cand = find(n_1>2); % 这里 1,2的阶乘皆为其本身，所以不需要计算

n_1(cand) = n_1(cand) - 1;

factorial(cand) = factorial(cand).*n_1(cand);

end

return;