% 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
没有合适的资源?快使用搜索试试~ 我知道了~
用于基因组规模模型重建、管理和分析的 RAVEN 工具箱
共1个文件
m:1个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 171 浏览量
2023-08-21
09:44:07
上传
评论
收藏 2KB ZIP 举报
温馨提示
用于基因组规模模型重建、管理和分析的 RAVEN 工具箱。.zip
资源推荐
资源详情
资源评论
收起资源包目录
ZernikeLegendreFit.zip (1个子文件)
ZernikeLegendreFit.m 6KB
共 1 条
- 1
资源评论
小风飞子
- 粉丝: 369
- 资源: 1962
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
最新资源
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功