用于基因组规模模型重建、管理和分析的 RAVEN 工具箱

% Function for fitting a surface to Zernike polynomials % % Usage: % [amn, bmn, z_recon_map] = ZernikeLegendreFit(z_map, M, N, % J, K, center_x, center_y) % % Description: % This function takes as input a freeform surface and fits it to % Zernike polynomials, returning the coefficients of the fit as well % as the reconstructed map. % % Inputs: % z_map - The sag table of the freeform surface. % M - The highest azimutal fitting order. % N - The highest radial fitting order. % J - The sampling number on azimuthal direction, default: 2*M+1 % K - The sampling number on radial direction, default: 2*(N+1) % center_x, center_y - The center of the freeform surface, default: % will be calculated by FindCenter function. % % Outputs: % amn, bmn - The fitted coefficients % z_recon_map - The surface calculated by amn, bmn coefficients % % % % Acknowledgments: % This work was inspired by the research of Greg Forbes and was % funded by 1. University of Rochester, 2. Center for Freeform Optics. % % This work was partially contributed by Ilhan Kaya in 2010, who wrote % function jacobiZernike to calculate the value of Jacobi polynomials % recursively. This function was later modified by Yiwen Fan as % function jacobiZernike_table. function [amn, bmn, z_recon_map] = ZernikeLegendreFit(z_map, m_max, n_max, J, K, center_j, center_i) if nargin < 6 [center_j, center_i] = FindCenter(z_map); end if nargin < 5 K = 2 * (n_max + 1); end if nargin < 4 J = 2 * m_max + 1; end if nargin < 3 error("fit_Zernike_quadrature requires at least 3 input: z_map, M, and N."); end r = 1; intp_type = "cubic"; [x_pix, y_pix] = size(z_map); r_pix = x_pix/2 *0.99; x = linspace(-r, r, x_pix); y = linspace(-r, r, y_pix); [x_map, y_map] = meshgrid(x, y); [i_map, j_map] = meshgrid((1:x_pix), (1:y_pix)); [theta_map, r_map] = cart2pol(x_map, y_map); del = pi/J; theta = (0:(2*J))*del; theta = theta(1:end-1); % Jacobi Matrix Jm = jacobiMatrix(0, K); % Legendre matrix [V, e] = eig(Jm); u_zeropoints = diag(e)'; % zeropoints w_m = 2*(V(1,:).^2); % weights u2 = u_zeropoints; u = sqrt(u2); [x_temp, y_temp] = pol2cart(theta, u'*r_pix); x_temp = x_temp + center_j; y_temp = y_temp + center_i; Stemp_m = interp2(i_map, j_map, z_map, x_temp, y_temp, intp_type); % Azimuthal fit A_qua = zeros(m_max+1, K); B_qua = zeros(m_max+1, K); for k = 1:K Stemp = Stemp_m(k,:); FFT_S = fft(Stemp)/size(theta,2); re = real(FFT_S); im = imag(FFT_S); A_qua(:,k) = 2*re(1:m_max+1); B_qua(:,k) = -2 * im(1:m_max+1); end A_qua(1,:) = A_qua(1,:)/2; B_qua(1,:) = B_qua(1,:)/2; % Radial fit amn = zeros(m_max+1, n_max+1); bmn = zeros(m_max+1, n_max+1); for m = 0:m_max Pmns = jacobiZernike_table(n_max+1, m, u2')' .* (u.^m); Pmns_weighted = Pmns .* w_m; for n = 0:n_max Pmn = Pmns_weighted(n+1, :); amn(m+1, n+1) = dot(A_qua(m+1, :), Pmn) *(2*n+m+1)/2; bmn(m+1, n+1) = dot(B_qua(m+1, :), Pmn) *(2*n+m+1)/2; end end % Reconstructrion [i_mesh, j_mesh] = meshgrid((1:x_pix),(1:y_pix)); [theta,rho] = cart2pol(i_mesh-center_i, j_mesh-center_j); % Converting to polar coords u_map = rho(:)/r_pix; theta = theta(:); u2_map = u_map.^2; u2_map(u2_map>0.99) = NaN; z_recon_map = zeros(x_pix, y_pix); for m = 0:m_max % idx = fix((m_max-m)/m_max*n_max); idx = n_max; Pmns_map = jacobiZernike_table(idx+1, m, reshape(u2_map, 1,[])'); Pmns_map = Pmns_map' .* (reshape(u_map, 1,[]).^m); for n = 0:idx Pmn_map = reshape(Pmns_map(n+1, :), x_pix, y_pix); a = amn(m+1, n+1); b = bmn(m+1, n+1); a_map = a * cos(m * theta_map) .* Pmn_map; b_map = b * sin(m * theta_map) .* Pmn_map; z_recon_map = z_recon_map + a_map + b_map; end end z_recon_map(u_map>0.99) = NaN; end %% FindCenter function [avg_i, avg_j] = FindCenter(Z) i_index = 1:size(Z,1); j_index = 1:size(Z,2); [valid_i, valid_j] = find(~isnan(Z)); avg_i = mean(i_index(valid_i),'all'); avg_j = mean(j_index(valid_j),'all'); end %% jacobiMatrix function T_mtrx = jacobiMatrix(m, n_max) % Returns a n_max*n_max matrix for a constant m term (jacobi beta) T_mtrx = zeros(n_max, n_max); an_ = -(m+1); bn_ = m+2; for n = 1:n_max-1 s = m+2*n; an = -(s+1)*((s-n).^2+n.^2+s)./(n+1)./(s-n+1)./s; bn = (s+2)*(s+1)./(n+1)./(s-n+1); cn = (s+2)*(s-n)*n./(n+1)./(s-n+1)./s; % three terms of Jacobi matrix, them perform a diagonal similarity % transorfation T_mtrx(n, n) = -an_/bn_; T_mtrx(n, n+1) = sqrt(cn/(bn_*bn)); T_mtrx(n+1, n) = T_mtrx(n, n+1); an_ = an; bn_ = bn; end T_mtrx(n_max, n_max) = -an_/bn_; end %% jacobiZernike_table function JZ=jacobiZernike_table(k,m,xe) % k is the degree of orthogonal jacobi polynomial % k=nf=(nw-m)/2 (relation between wolf - forbes) % m is the azimuthal order % xe grid points in radial coordinate. % ilhan kaya 11/15/2010 % Yiwen modified these to export a whole table. 04/25/2023 % name change from jacobiZernike to jacobiZernike_table JZ = zeros(size(xe,1),k+1); JZ(:,1)=ones(size(xe,1),1); % n=0 JZ(:,2)=(m+2)*xe-(m+1); % n=1 --> an = -(m+1); bn = m+2; cn = 0; %recurrence Forbes Opt. Exps. Vol. 18 No. 12 June 2010 % P_(n+1)(x)= (an+bn*x)*P_n(x)-cn*P_(n-1)(x) for n = 1:k-1 s = m+2*n; an = -(s+1)*((s-n).^2+n.^2+s)./(n+1)./(s-n+1)./s; bn = (s+2)*(s+1)./(n+1)./(s-n+1); cn = (s+2)*(s-n)*n./(n+1)./(s-n+1)./s; JZ(:,n+2)=(an+bn*xe).*JZ(:,n+1)-cn*JZ(:,n); end end