No BSD License  

image thumbnail

Zernike decomposition

by

 

10 Dec 2007 (Updated )

Decomposition of a 2-D function by set of Zernike functions

deco.m
%This code was written to deal with "Zernike polynomials" code graciously
%donated by Paul Fricker via file exchange.
%Here you will find a practical example of a function decomposition by
%Zernike basis.
%The function is F below, feel free to modify
%Unlike Paul's example found in 'zernfun2.m' here the domain is the true
%unit circle, without NaN's filling it up to the unit square.

%Alex Chtchetinin, December 10, 2007


N=100;%dimension of theta-R grid
%r=linspace(1/N/2,1-1/N/2,N);t=linspace(1/N/2,2*pi*(1-1/N/2),N);
r=linspace(0,1,N);t=linspace(0,2*pi,N);
[t,r]=meshgrid(t,r);
[x,y] = pol2cart(t,r);
F=x.^3-y.^2;%Function to be decomposed
figure,surface(x,y,F);shading interp
title('Function to be decomposed');axis square
p = 0:8;%indexes of Zernike functions to be used in the decomposition
n=length(p);
z = zernfun2(p,r(:),t(:),'norm');
zz=reshape(z,N,N,length(p));
figure;
for k = 1:n
    z = zz(:,:,k);
    subplot(ceil(n^.5),ceil(n^.5),k)
    surface(x,y,z), shading interp
    set(gca,'XTick',[],'YTick',[])
    axis square
    title(['Z_' num2str(p(k))]);
end
for i=1:n
    zzz=zz(:,:,i).*F.*r;
    G = @(x,y)interp2(t,r,zzz,x,y);
    Q(i) = dblquad(G,0,2*pi,0,1);
end
disp(Q)
Z=0*zz(:,:,1);
figure;
for k = 1:n
    Z=Z+Q(k)*zz(:,:,k);
    subplot(ceil(n^.5),ceil(n^.5),k)
    surface(x,y,Z), shading interp
    set(gca,'XTick',[],'YTick',[])
    axis square
    if k==1
        title(['Z_' num2str(p(k))]);
    else
        title(['Z_0 to Z_' num2str(p(k))]);
    end
end
figure,surface(x,y,Z-F);shading interp;
title('Residual error');axis square

Contact us