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);endreturn;