利用前36阶zernike函数拟合曲面：
脚本程序
clc;clear;
load unwrap_ph.mat
unwrap_ph=max(max(unwrap_ph))-unwrap_ph;
unwrap_ph=unwrap_ph(:,81:560);
x=linspace(-1,1,size(unwrap_ph,1));
y=linspace(-1,1,size(unwrap_ph,2));
xy=[x;y];
a0=ones(1,36);
unwrap_ph(isnan(unwrap_ph)) = 0;
a=lsqcurvefit('SH',a0,xy,unwrap_ph)
FIT=SH(a,xy);
FIT(FIT==0)=nan;
figure(1)
imagesc(FIT)
colormap('jet')
grid off
shading interp
figure(2)
imagesc(unwrap_ph)
colormap('jet')
grid off
shading interpzernike函数
function z = zernfun(n,m,r,theta,nflag)
m=m*-1;
if (~any(size(n)==1))||( ~any(size(m)==1))
error('zernfun:NMvectors','N and M must be vectors.')
end
if length(n)~=length(m)
error('zernfun:NMlength','N and M must be the same length.')
end
n=n(:);
m=m(:);
if any(mod(n-m,2))
error('zernfun:NMmultiplesof2', ...
'All N and M must differ by multiples of 2 (including 0).')
end
if any(m>n)
error('zernfun:MlessthanN', ...
'Each M must be less than or equal to its corresponding N.')
end
if any(r>1|r<0)
error('zernfun:Rlessthan1','All R must be between 0 and 1.')
end
if (~any(size(r)==1))||(~any(size(theta)==1))
error('zernfun:RTHvector','R and THETA must be vectors.')
end
r=r(:);
theta=theta(:);
length_r=length(r);
if length_r~=length(theta)
error('zernfun:RTHlength', ...
'The number of R- and THETA-values must be equal.')
end
% Check normalization:
% --------------------
if nargin==5&&ischar(nflag)
isnorm=strcmpi(nflag,'norm');
if ~isnorm
error('zernfun:normalization','Unrecognized normalization flag.')
end
else
isnorm=false;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Zernike Polynomials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the required powers of r:
% -----------------------------------
m_abs=abs(m);
rpowers=[];
for j=1:length(n)
rpowers=[rpowers m_abs(j):2:n(j)];
end
rpowers=unique(rpowers);
% Pre-compute the values of r raised to the required powers,
% and compile them in a matrix:
% -----------------------------
if rpowers(1)==0
rpowern=arrayfun(@(p)r.^p,rpowers(2:end),'UniformOutput',false);
rpowern=cat(2,rpowern{:});
rpowern=[ones(length_r,1) rpowern];
else
rpowern=arrayfun(@(p)r.^p,rpowers,'UniformOutput',false);
rpowern=cat(2,rpowern{:});
end
% Compute the values of the polynomials:
% --------------------------------------
y=zeros(length_r,length(n));
for j=1:length(n)
s=0:(n(j)-m_abs(j))/2;
pows=n(j):-2:m_abs(j);
for k=length(s):-1:1
p=(1-2*mod(s(k),2))* ...
prod(2:(n(j)-s(k)))/ ...
prod(2:s(k))/ ...
prod(2:((n(j)-m_abs(j))/2-s(k)))/ ...
prod(2:((n(j)+m_abs(j))/2-s(k)));
idx=(pows(k)==rpowers);
y(:,j)=y(:,j) + p*rpowern(:,idx);
end
if isnorm
y(:,j)=y(:,j)*sqrt((1+(m(j)~=0))*(n(j)+1)/pi);
end
end
% END: Compute the Zernike Polynomials
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Zernike functions:
% ------------------------------
idx_pos=m>0;
idx_neg=m<0;
z = y;
if any(idx_pos)
if (m>=0&&rem(n,2)==1)
z(:,idx_pos)=1-y(:,idx_pos).*sin(theta*m(idx_pos)');
else
z(:,idx_pos)=y(:,idx_pos).*sin(theta*m(idx_pos)');
end
end
if any(idx_neg)
if (m>=0&&rem(n,2)==1)
z(:,idx_neg)=1-y(:,idx_neg).*cos(theta*m(idx_neg)');
else
z(:,idx_neg)=y(:,idx_neg).*cos(theta*m(idx_neg)');
end
end
% EOF zernfun拟合函数
function z=SH(c,data)
x=data(1,:);
y=data(2,:);
[X,Y]=meshgrid(x,y);
[theta,r]=cart2pol(X,Y);
idx=r<=1;
zz=zeros(size(X));
b=c(1:4);
a=c(5:end);
zz(idx)=b(1)*zernfun(0,0,r(idx),theta(idx))+b(2)*zernfun(1,1,r(idx),theta(idx))+b(3)*zernfun(1,-1,r(idx),theta(idx))+b(4)*zernfun(2,0,r(idx),theta(idx))+...
a(1)*zernfun(2,2,r(idx),theta(idx))+a(2)*zernfun(2,-2,r(idx),theta(idx))+a(3)*zernfun(3,1,r(idx),theta(idx))+...
a(4)*zernfun(3,-1,r(idx),theta(idx))+a(5)*zernfun(3,3,r(idx),theta(idx))+a(6)*zernfun(3,-3,r(idx),theta(idx))+...
a(7)*zernfun(4,0,r(idx),theta(idx))+a(8)*zernfun(4,2,r(idx),theta(idx))+a(9)*zernfun(4,-2,r(idx),theta(idx))+...
a(10)*zernfun(4,4,r(idx),theta(idx))+a(11)*zernfun(4,-4,r(idx),theta(idx))+a(12)*zernfun(5,1,r(idx),theta(idx))+...
a(13)*zernfun(5,-1,r(idx),theta(idx))+a(14)*zernfun(5,3,r(idx),theta(idx))+a(15)*zernfun(5,-3,r(idx),theta(idx))+...
a(16)*zernfun(5,5,r(idx),theta(idx))+a(17)*zernfun(5,-5,r(idx),theta(idx))+a(18)*zernfun(6,0,r(idx),theta(idx))+...
a(19)*zernfun(6,2,r(idx),theta(idx))+a(20)*zernfun(6,-2,r(idx),theta(idx))+a(21)*zernfun(6,4,r(idx),theta(idx))+...
a(22)*zernfun(6,-4,r(idx),theta(idx))+a(23)*zernfun(6,6,r(idx),theta(idx))+a(24)*zernfun(6,-6,r(idx),theta(idx))...
+a(25)*zernfun(7,1,r(idx),theta(idx))+a(26)*zernfun(7,-1,r(idx),theta(idx))+a(27)*zernfun(7,3,r(idx),theta(idx))...
+a(28)*zernfun(7,-3,r(idx),theta(idx))+a(29)*zernfun(7,5,r(idx),theta(idx))+a(30)*zernfun(7,-5,r(idx),theta(idx))...
+a(31)*zernfun(7,7,r(idx),theta(idx))+a(32)*zernfun(7,-7,r(idx),theta(idx));
z=zz;
原图
拟合结果
应当注意的是zernike是圆对称函数，也就意味着其拟合的是一个圆，所以原图圆外的数值并未拟合。
最近发现了一个zernike拟合的实现，这里给出分享吧：
function [output1 output2] = ZernikeCalc( ZernikeList, Zdata, mask, ...
ZernikeDef, ShapeParameter, ...
unitCircle, MMESorC)
%ZERNIKECALC Uses 'mask' region to fit circular (or other shape) Zernikes to surface data.
%
% VERSION: 2013-02-07 (YYYY-MM-DD)
%
% Fits circular, hexagonal, rectangular, square, elliptical, or annulus
% orthogonal Zernike polynomials to the surface data provided. If no surface
% data is provided then plots of these polynomial functions over the
% region specified are produced. Several different convesions (sign,
% normalization) can be explicitly specified.
%
% function [output1 output2] = ZernikeCalc( ZernikeList, Zdata, mask, ...
% ZernikeDef, ShapeParameter, ...
% unitCircle, MMESorC)
%
% INPUT:
% ZernikeList = if size() = (1,k) then this is a Zernike j-ordering
% list of numbers (j).
% if size() = (2,k) then this is Zernike (n, m) list.
% row 1 is the list of n values, row 2 is the list of
% m values. m is positive or negative values to indicate
% whether sin() or cos() factor for the term. See MMESorC
% parameter for the meaning of a minus m value.
% DEFAULT: [1:15]. First 15 j-ordered Zernikes.
%
% Zdata = If this is a single column then it contains the
% Zernike coefficients and we are to calculate the
% surface data.
% If this is more than a single column then it is a
% surface data and we are to calculate the
% coefficients for the Zernike polynomials specified.
% The surface data must be a phase unwrapped surface data.
% DEFAULT: ones(15,1).
%
% mask = Matrix (same size as surface data) indicating
% the individual pixels in the surface data that are
% to be used for fitting of the Zernike polynomials.
% '0' = don't use this pixel, '1' = use this pixel.
% If mask is a scalar number n, then an nxn mask matrix
% is created with an n diameter circle as the mask.
% DEFAULT: 100x100 with circluar mask.
%
% ZernikeDef = One of 'FRINGE', 'ISO', 'WYANT', 'MAHAJAN', 'NOLL',
% 'B&W', 'STANDARD'
% 'HEXAGON', 'HEXAGONR30', 'ELLIPSE',
% 'RECTANGLE', 'SQUARE', 'ANNULUS'
% See table below for possible j-value ranges.
% NOTE: 'MAHAJAN' and 'NOLL' are programmed to be the same.
% NOTE: 'B&W' and 'STANDARD' are programmed to be the same.
% DEFAULT: 'FRINGE'.
%
% ShapeParameter = For 'ELLIPSE' this is the aspect ratio (b).
% For 'RECTANGLE' this is half-width along the x-axis (a).
% For 'ANNULUS' this is the obscuration ratio (e).
% For all other shapes this is ignored.
% DEFAULT: 0, which is not valid for 'ELLIPSE' and 'RECTANGLE'.
%
% unitCircle = Gives the unit circle around the mask which the circular
% Zernikes are defined on. It is at most a 1x3 vector
% specifying centerRow, centerCol, radiusInPixels. The
% permutations and ording allowed are as follows:
% [] empty so all 3 defaults are calculated.
% [radiusInPixels] the default centerRow and centerCol
% are calculated.
% [centerRow, centerCol] the default radiusInPixels is
% calculated.
% [centerRow, centerCol, radiusInPixels] no defaults.
% DEFAULT: The defaults for centerRow and centerCol are
% calculated by calculating the weighted
% average of row and column for the mask's 1's.
% The default radiusInPixels is calculated to
% be the largest distance from (centerRow,centerCol)
% to all mask's pixels with a '1' value.
%
% MMESorC = Either '-m=sin' or '-m=cos'. Indicates, for (n, m) ordering
% what a minus m value represents, a sine factor or a
% cosine factor.
% DEFAULT: '-m=sin'. Ignored when j ordering is used.
%
% OUTPUT:
% - When no output result is requested (ZernikeCalc(...)) this function
% produces plots of the Zernike polynomials requested.
% If Zdata is a single column specifying the coefficients then
% a plot of the sum of each of the Zernikes specified is produced.
% If Zdata is an n x m matrix of surface data then 3 plots are
% produced: 1) of the surface data, 2) of the fitted surface, 3) the
% difference between the surface data and the fitted surface.
% - When one output result is requested (results = ZernikeCalc(...)) then 1
% of 2 possible results are returned:
% if the input Zdata is a single column specifying the coefficients
% to multiply the Zernike polynomials by then the result retruned
% is the Zernike polynomial value matrices (across the mask area).
% if the input Zdata is a matrix for which the Zernikes are to be
% fit then the results returned is a column vector of the Zernike
% coefficients corrseponding to (and in the same order as) the
% polynomials identified by ZernikeList.
% - If 2 output results are requested ([out1 out2] = ZernikeCalc(...)) then
% if the input for Zdata is a single column, giving the coefficients
% to multiply the Zernike polynomials by, then out1 is the sum
% of the Zernike polynomials requested, and out2 is a 3 dimensional
% matrix of all the n Zernike polynomials requested [:,:,n].
% if the input for Zdata is not a column vector then out1 is the
% 3 dimensional data cube of the n Zernike polynomials requested [:,:,n]
% and out2 are the coefficients used.
%
%
% Examples:
%
% ZernikeCalc
% - Displays the first 15 Fringe Zernikes in 15 color plots.
%
% ZernikeCalc([4 5 7 8], [0.4; -0.6; 1.2; 0.25])
% - Displays Fringe Zernikes Zj=4, Zj=5, Zj=7, Zj=8 multiplied by the
% scaling factors 0.4, -0.6, 1.2 and 0.25, respectively in 4 separate
% color plots.
%
% ZernikeCalc([4 5 7 8], [0.4; -0.6; 1.2; 0.25], [], 'standard')
% - Same as the last case only using standard Zernikes rather than Fringe
% Zernikes.
%
% ZernikeCalc([2 2; 2 0; 3 3; 3 1]', [0.4; -0.6; 1.2; 0.25], [], 'standard')
% - Same as last case now using Z(n,m) notation to specify which Zernikes
% to use.
%
% Let SD be an n x m matrix of surface data to which the specified
% Zernikes are to be fit. Then
%
% coeff = ZernikeCalc([2 2; 2 0; 3 3; 3 1]', SD, [], 'standard')
%
% returns a column consisting of the calculated fitting coefficients
% in coeff. No plots are produced.
%
% [DC, coeff] = ZernikeCalc([2 2; 2 0; 3 3; 3 1]', SD, [], 'standard')
%
% returns a column consisting of the calculated fitting coefficients
% in coeff and a n x m x 4 data cube. ('4' because 4 Zernikes were
% specified) Each DC(:, :, i) (i=1,2,3,4) is the ith specified Zernike
% fitted to the surface data SD across the (default) mask area.
%
% [DC, coeff] = ZernikeCalc([4 5 7 8], SD, [], 'annulus', 0.25)
%
% This uses the annular Zernikes, with a central obscuration radius ratio of 0.25
% to fit the surface data. See Ref. 1 for details on noncircular Zernikes.
%
%
% Circular Zernike polynomials are available in several optical design
% software packages, including Code V(R), OSLO(R), Zemax(R), etc.
%
% Table 1 Software Conventions
%--------------------------------------------
%|INPUT PARAM |APERTURE |SOFTWARE|ORDER & |
%|ZernikeDef | SHAPE | | RANGE |
%|------------|---------|--------|----------|
%|'B&W', |Circle |CODE V | (n, m), |
%|'STANDARD' | | | j=1... |
%|------------|---------|--------|----------|
%|'MAHAJAN', |Circle |ZEMAX | (n, m), |
%|'NOLL' | | | j=1... |
%|------------|---------|--------|----------|
%|'FRINGE' |Circle |CODE V, | j=1...37 |
%| | |ZEMAX | |
%|------------|---------|--------|----------|
%|'ISO' |Circle | | j=0..35 |
%|------------|---------|--------|----------|
%|'WYANT' |Circle | OSLO | j=0...36 |
%|------------|---------|--------|----------|
%|'HEXAGON' |Hexagon | | j=1...45 |
%|------------|---------|--------|----------|
%|'HEXAGONR30'|Heaxgon | | j=1...28 |
%| |rotated | | |
%| |30 deg. | | |
%|------------|---------|--------|----------|
%|'ELLIPSE' |Ellipse | | j=1..15 |
%|------------|---------|--------|----------|
%|'RECTANGLE' |Rectangle| | j=1...15 |
%|------------|---------|--------|----------|
%|'SQUARE' |Square | | j=1...45 |
%|------------|---------|--------|----------|
%|'ANNULUS' |Annulus | | j=1...35,|
%| | | | j<>29,30,|
%| | | | 31,32 |
%--------------------------------------------
%
% Ref. 1: Mahajan, V.N., G.-m. Dai, "Orthonormal polynomials in wavefront
% analysis: analytical solution," J. Opt. Soc. Am. A, Vol. 24, No. 9
% Sept. 2007.
%
%
% Updates: 2012-01-08 (YYYY-MM-DD)
% RWG - Added default mask shapes for the different ZernikeDef
% input parameter values.
%
% 2012-01-08 (YYYY-MM-DD)
% RWG - When no output requested ZernikeCalc will print all
% Zernike polynomials specified.
%
%
% Code Copyrighted, 2011-2013 by Robert W. Gray and Joseph M. Howard. All
% rights reserved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For validation of the equations, uncomment the next 2 lines.
% Then, it doesn't matter what input parameters are specified.
%
% validateZ();
% return;
%
% end validation of equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assign default values to input parameters that are empty.
% Empty is '' for strings and [] for vectors and matrices.
%
if nargin == 0 || isempty(ZernikeList)
ZernikeList = 1:15; % default is first 15 fringe Zernikes.
end % if statement
if nargin <= 1 || isempty(Zdata)
theSize = size(ZernikeList);
Zdata = ones(theSize(1,2),1); % all coefficients are 1.
end % if statement
if nargin <= 3 || isempty(ZernikeDef)
ZernikeDef = 'FRINGE';
end % if statement
% Convert to upper case
ZernikeDef = upper(ZernikeDef);
if nargin <= 4 || isempty(ShapeParameter)
ShapeParameter = 0;
end % if statement
if nargin <= 2 || isempty(mask)
% make a default mask
defaultRows = 100;
defaultCols = 100;
theZdataSize = size(Zdata);
if (theZdataSize(1,1) > 1) && (theZdataSize(1,2) > 1)
defaultRows = theZdataSize(1,1);
defaultCols = theZdataSize(1,2);
end % if statement
mask = makeDefaultMask(ZernikeDef, defaultRows, defaultCols, ShapeParameter);
end % if no mask statement
sm = size(mask);
if (sm(1,1) == 1) && (sm(1,2) == 1)
% if mask is a scalar n then create nxn matrix
% make a circular mask
defaultRows = mask(1,1);
defaultCols = mask(1,1);
mask = makeDefaultMask(ZernikeDef, defaultRows, defaultCols, ShapeParameter);
end % end if statement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate default centerRow and centerCol values.
% These might have been specified, but if they haven't
% use these values.
radiusInPixels = 0;
theSize = size(mask);
numrows = theSize(1,1);
numcols = theSize(1,2);
% calculate center by weighted averages.
sumMaskRow = sum(mask,2);
sumMaskCol = sum(mask,1);
sumMaskAll = sum(sum(mask));
centerRow = sum(sumMaskRow .* ((1:numrows)')) / sumMaskAll;
centerCol = sum(sumMaskCol .* (1:numcols)) / sumMaskAll;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin <= 5 || isempty(unitCircle)
% The unit circle associated with the mask has not been specified.
% unitCircle is empty.
unitCircle = [centerRow, centerCol];
end % if statement
% At this point, unitCircle is not empty.
theSize = size(unitCircle);
nc = theSize(1,2);
switch nc
case 1
% only the radiusInPixels has been specified
radiusInPixels = unitCircle(1,1);
case 2
% the centerRow and centerCol have been specified
% so can now calculate the radius in Pixels.
centerRow = unitCircle(1,1);
centerCol = unitCircle(1,2);
% a matrix such that each element (r,c) has the value r-centerRow
rm = (((1:numrows)-centerRow)'*ones(1,numcols)).*mask;
% a matrix such that each element (r,c) has the value c-centerCol
cm = (ones(numrows,1)*((1:numcols)-centerCol)).*mask;
% sqrt(rm.^2 + cm.^2) is a matrix such that (r,c) contains the distance
% of (r,c) to the center (centerRow, centerCol).
radiusInPixels = max(max(sqrt(cm.^2 + rm.^2)));
case 3
% the centerRow, centerCol, radiusInPixels have been specified
centerRow = unitCircle(1,1);
centerCol = unitCircle(1,2);
radiusInPixels = unitCircle(1,3);
otherwise
% error.
end % switch statement
if nargin <= 6 || isempty(MMESorC)
MMESorC = '-m=sin';
end % if stateemnt
%
% end of section on default input assignments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input validation
% Too much code is distracting, so put validation in separate function.
validateInput(ZernikeList, Zdata, mask, ...
ZernikeDef, ShapeParameter, ...
centerRow, centerCol, radiusInPixels, MMESorC);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Zernike polynomials are caluclated using (n, m) ordering
% so need to calculate (n,m,sc) even when j ordering is specified.
% sc = 's' means sin() factor
% sc = 'c' means cos() factor
% sc = ' ' means m = 0 so has no sin() or cos() factor
%
theSize = size(ZernikeList);
maxNumOfZernikeTerms = theSize(1,2);
n = zeros(1, theSize(1,2));
m = zeros(1, theSize(1,2));
sc = ' ';
if theSize(1,1) == 1
% the ZernikeList is list of j order values.
for k=1:theSize(1,2)
% convert the j ordering of type ZernikeDef to (n,m,sc) of
% same ZernikeDef.
[n(k) m(k) sc(k)] = convertj(ZernikeList(1, k), ZernikeDef);
end % for statement
else
% the ZernikeList is list of (n,m) pairs.
for k=1:theSize(1,2)
% convert (n,m) to (n,m,sc) using MMESorC
n(k) = ZernikeList(1, k);
m(k) = abs(ZernikeList(2, k));
sc(k) = ' ';
switch MMESorC
case '-m=sin'
if ZernikeList(2, k) < 0
sc(k) = 's';
end % if statement
if ZernikeList(2, k) > 0
sc(k) = 'c';
end % if statement
case '-m=cos'
if ZernikeList(2, k) < 0
sc(k) = 'c';
end % if statement
if ZernikeList(2, k) > 0
sc(k) = 's';
end % if statement
end % switch statement
end % for k statement
end % if j or (n,m) ordering
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Preallocate some of the matrices.
%
% We'll need the size of the mask matrix (same as surface data).
theSize = size(mask);
rows = theSize(1,1);
cols = theSize(1,2);
% pre-allocate vectors and matrices.
numOfMaskPixels = sum(sum(mask));
zcoeff = zeros(maxNumOfZernikeTerms, 1);
xMatrices = zeros(rows, cols, maxNumOfZernikeTerms);
xVectors = zeros(numOfMaskPixels, maxNumOfZernikeTerms);
yVector = zeros(numOfMaskPixels, 1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check if surface data is the input or if a coefficent vector is
% the input.
interferogramInput = true;
theSize = size(Zdata);
if theSize(1,2) == 1
% There is one column in Zdata so it is a column vector of coefficents.
zcoeff = Zdata;
interferogramInput = false;
end % if statement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for each pixel (matrix element) we need the distance from the
% center (centerRow, centerCol). So we make a matrix that has
% as element values the distance of that (row, col) element to
% (centerRow, centerCol).
%
% a matrix such that each element (r,c) has the value r-centerRow
rm = ((1:rows)-centerRow)'*ones(1,cols);
% a matrix such that each element (r,c) has the value c-centerCol
cm = ones(rows,1)*((1:cols)-centerCol);
% sqrt(rm.^2 + cm.^2)./radiusInPixels is a matrix such that (r,c) contains
% the normalized distance of (r,c) to the center (centerRow, centerCol).
% Then we reshape this into a column vector.
rho = reshape(sqrt(cm.^2 + rm.^2)./radiusInPixels, rows*cols, 1);
% atan2(rm, cm) is a matrix such that each element (r,c) contains
% the angle (radians) from the centerRow axis to (r,c). Then we
% reshape this into a vector.
theta = reshape(atan2(rm, cm), rows*cols, 1);
% reshape the mask into a vector.
vmask = reshape(mask, rows*cols, 1);
if interferogramInput
% we have surface data so reshape it just like rho and theta.
yVector = reshape(Zdata, rows*cols, 1);
% and remove all the elements that don't have a corresponding '1'
% in the mask.
yVector(vmask ==0) = [];
end % if statement
for i=1:maxNumOfZernikeTerms
% calculate the ith Zernike polynomial value for each pixel
% in the mask matrix (now a vector).
switch ZernikeDef
% Handle the special shapes.
case 'HEXAGON'
hldZ = ZHexagon(ZernikeList(1,i), rho, theta);
case 'HEXAGONR30'
hldZ = ZHexagonR30(ZernikeList(1,i), rho, theta);
case 'ELLIPSE'
hldZ = ZEllipse(ZernikeList(1,i), rho, theta, ShapeParameter);
case 'RECTANGLE'
hldZ = ZRectangle(ZernikeList(1,i), rho, theta, ShapeParameter);
case 'SQUARE'
hldZ = ZSquare(ZernikeList(1,i), rho, theta);
case 'ANNULUS'
hldZ = ZAnnulus(ZernikeList(i), n(i), m(i), rho, theta, ShapeParameter);
otherwise
% Otherwise, its a circle Zernike.
hldZ = calculateZernike(n(i), m(i), sc(i), rho, theta, ZernikeDef);
end % switch ZernikeShape
% reshape the column vector result into a (rows, cols) matrix.
xMatrices(:, :, i) = reshape(hldZ,rows,cols);
% remove the elements from the Zernike calculation that do not
% have a corresponding '1' in the mask.
hldZ(vmask == 0) = [];
% this is one of the Zernike polynomial results for each pixel
% for which the mask is a '1'.
xVectors(:, i) = hldZ;
end % for i statement
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the regression (least squares fit) only if the
% input Zdata is surface data.
%
if interferogramInput
% Use least squares fit to determine the coefficients.
zcoeff = xVectors\yVector;
end % if statement
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiply Zernike polynomial matrices by the calculated coefficients.
% since we already have the Zernike matrices calculated, we do this here.
%
zm = zeros(rows, cols, maxNumOfZernikeTerms);
for i=1:maxNumOfZernikeTerms
% multiply the Zernike matrix by the corresponding coefficient
% and by the mask to zero out pixels that we don't care about.
zm(:,:,i) = xMatrices(:, :, i) .* zcoeff(i) .* mask;
end % for statement
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Change the output depending on input and output specified.
%
if nargout < 1
% plot Zernike figures
sizeZD = size(Zdata);
if sizeZD(1,2) == 1
% plot each of the Zernike polynomials
figs(zm);
% plot the sum of the request Zernike Polynomials
% figs(sum(zm,3));
return;
end % if statement
% plot 1) input surface data, 2) the fit results, 3) their difference
toPlot(:,:,1) = Zdata;
theFit = sum(zm, 3);
toPlot(:,:,2) = theFit;
theDiff = Zdata - theFit;
toPlot(:,:,3) = theDiff;
figs(toPlot, mask, 3, {'Input Data', 'Zernike Fit', 'Input Data - Zernike Fit'});
return;
end % if nargout < 1
% At this point we know nargout >= 1
if nargout < 2
% This means nargout == 1.
if interferogramInput
% return calculated coefficients only
output1 = zcoeff;
return;
end % if interferogramInput
% Not an interferogram as Input. Must be coefficients as input
% so return only the interferogram. They already know the
% coefficients.
output1 = zm;
return;
end % if nargout < 2
% At this point we know the output requested is 2 (or more).
theZdataSize = size(Zdata);
if theZdataSize(1,2) == 1
output1 = sum(zm,3);
output2 = zm;
return;
end % if statement
output1 = zm;
output2 = zcoeff;
end % ZernikeCalc
function handle = figs(data,mask,numplotcols,titles,labeloption,labelprecision,handle)
%displays array data using imagesc function
% function h = figs(data,mask,numplotcols,titles,labeloption,labelprecision)
% INPUT: data, 2, 3 or 4 dimensional stack of 2d arrays
% 2-d data of column vectors (i.e. non-square) are reshaped
% into 3-d stack of square figures
% 3-d data is subplotted with max of 6 cols and rows, with multiple
% figures created.
% 4-d data is plotted in a single figure with the 3rd and 4th
% dimension determining the subplot array size.
% mask is optional, either one array for all, or a stack similar
% in size to data.
% numplotcols = set number of plot columns desired (default = 6 max)
% titles = cell array of text titles for each plot (optional)
% labeloption = 1, RMS data given
% = 2, Avg, RMS, and Peak-to-Valley given
% labelprecision = number of digits in label (default = 4)
% handle = handle of current figure to use
%
% OUTPUT: handle = handles for each figure
% To do: 1. add 5d capability
% 2. ensure mask is permuted along with data for 4d input
if nargin<1, handle = figure; return; end
if nargin<2, mask = []; end
if nargin<3, numplotcols = []; end
if nargin<4, titles = {}; end; if ~iscell(titles), titles = {}; end
if nargin<5, labeloption = []; end; if isempty(labeloption), labeloption = 1; end
if nargin<6, labelprecision = []; end; if isempty(labelprecision), labelprecision = 4; end
if nargin<7, handle = []; end
sizedata = size(data);
%extract data from structure if given as such
if isstruct(data)
datastruct = data; clear data;
for i=1:length(datastruct)
if isfield(datastruct,'opd'), data(:,:,i) = datastruct(i).opd; end
if isfield(datastruct,'mask'), mask(:,:,i) = datastruct(i).mask; end
if isfield(datastruct,'psf'), data(:,:,i) = datastruct(i).psf; end
end
end
% reshape data if 2d column or row input and square (e.g. size(col) = 100^2)
if size(data,1) ~= size(data,2) % if non-square input matrix
sqrtval = sqrt(length(data(:,1)));
if mod(length(data(:,1)),10000) == 0 %if data is integer lengths of 100x100
disp('Data appears to be integer number of 100x100 matrices in a column vector.')
disp('Reshaping data into 2-D stack.')
num100x100 = length(data(:,1))/10000;
if num100x100 == 1
data = reshape(data,100,100,size(data,2)); %3d data
else
data = reshape(data,100,100,num100x100,size(data,2)); %4d data
end
elseif mod(length(data(:,1)), sqrtval) == 0 % if data square
disp('Data appears to be square in column format, reshaping into 2d.');
data = reshape(data,sqrtval,sqrtval,size(data,2));
else
sqrtval = sqrt(length(data(:,2)));
if mod(length(data(:,1)), sqrtval) == 0 % if data square
disp('Data appears to be square in row format, reshaping into 2d.');
data = reshape(data',sqrtval,sqrtval,size(data',2));
end
end
end
% determine number of plots, and convert 5d and 4d data to 3d stack
if ndims(data)==5
[s1,s2,s3,s4,s5] = size(data);
data = reshape( permute(data,[1 2 5 4 3]),s1,s2,s3*s4*s5);
end
if ndims(data)==4
[s1,s2,s3,s4] = size(data);
numplots = s3*s4;
data = reshape( permute(data,[1 2 4 3]),s1,s2,s3*s4);
data4d = s4;
else
numplots = size(data,3);
data4d = 0;
end
%create mask data if needed
if isempty(mask), mask = ~isnan(data) & data~=0; end
% determine subplot format: rows and cols
if isempty(numplotcols)
if data4d, numplotcols = data4d; %let 4d data size determine rows and cols
elseif numplots==1, numplotcols=1;
elseif numplots<5, numplotcols = 2;
elseif numplots<7, numplotcols = 3;
elseif numplots<9, numplotcols = 4;
elseif numplots<11, numplotcols = 5;
else numplotcols = 6;
end
end
numplotrows = ceil(numplots/numplotcols);
if numplotrows>6 && data4d<1, numplotrows=6; end
if numplots>numplotrows*numplotcols
numfigs = ceil(numplots/(numplotrows*numplotcols)); %number of figures
else
numfigs = 1;
end
if ~iscell(titles) %array of text given
for i=1:size(titles,1)
t{i} = titles(i,:);
end
titles = t;
end
numtitles = length(titles(:));
disp(['Total number of figures = ' int2str(numfigs)]);
disp(['Number of plots per figure = ' int2str(numplotrows*numplotcols)]);
disp(['Number of plots rows = ' int2str(numplotrows)]);
disp(['Number of plots cols = ' int2str(numplotcols)]);
% plot data
for k=1:numfigs
scrsz = get(0,'ScreenSize');
if numfigs>1 || isempty(handle)
handle = figure('Position',[scrsz(3)/4 scrsz(4)/4 0.6*scrsz(3) 0.6*scrsz(4)],'color',[1 1 1]);
else figure(handle);
end
set(gcf,'Name',inputname(1));
for i=1:numplotrows
for j=1:numplotcols
plotnum = (i-1)*numplotcols+j+(k-1)*numplotrows*numplotcols; %data to plot
subplotnum = (i-1)*numplotcols+j; %subplot location to plot
if plotnum<numplots+1,
plotdata = data(:,:,plotnum)*1; % *1 converts logical to doubles
if size(mask,3)>1, maskdata = mask(:,:,plotnum); else maskdata = mask(:,:,1); end
subplot(numplotrows,numplotcols,subplotnum);
plotdata(maskdata==0)=nan;
imagesc(plotdata);
axis square;
axis xy;
alpha(1*maskdata);
set(gca,'XTick',[]);
set(gca,'YTick',[]);
set(gca,'ZTick',[]);
if numtitles<1
title(int2str(plotnum));
elseif plotnum>numtitles %do nothing if partial title list given
else title(titles{plotnum});
end
put_xlabel(plotdata,maskdata,labeloption,labelprecision); %subfunction included below
end
if nargin>3, if length(titles)>=plotnum, title(titles(plotnum)); end; end
end
end
end
end % function figs
% X Label for each plot
function put_xlabel(plotdata,maskdata,labeloption,labelprecision)
if nargin<2
maskdata = ~isnan(maskdata) & data~=0;
datavals = plotdata(maskdata);
else
if ~islogical(maskdata), maskdata = logical(maskdata); end
datavals = plotdata(maskdata);
end
if ~isempty(datavals) %only analyze statistics on a real data set
stdval = std(plotdata(maskdata));
% pk2val = pv(plotdata,maskdata);
maxval = max(datavals);
minval = min(datavals);
meanval = mean(datavals);
p95 = sort(datavals);
p95 = p95(ceil(.95*length(datavals)));
else
stdval = 0; pk2val = 0; maxval = 0; minval = 0; meanval = 0; p95 = 0;
end
if labeloption==1
xlabel(['RMS = ' num2str(stdval,labelprecision)]);
elseif labeloption==2
xlabel({['AVG = ' num2str(meanval,labelprecision)];['RMS = ' num2str(stdval,labelprecision)];['PV = ' num2str(pk2val,labelprecision)]});
elseif labeloption==3
xlabel(['RMS = ' num2str(stdval*1e9,labelprecision) ' nm']);
end
end % function put_xlabel
function result = calculateZernike(n, m, sc, rho, theta, ZernikeDef)
%
% Calculates the Zernike polynomial value for the given pixel location.
%
% INPUT:
% n = radial order.
% m = azimuthal order.
% sc = ' ' for m = 0, = 's' for sin() term, = 'c' for cos() term.
% rho = normalized radial distance to pixel location.
% theta = angle from the x-axis in radians of pixel location.
% ZernikeDef = One of 'MAHAJAN', 'NOLL', 'FRINGE', 'ISO', 'WYANT',
% 'B&W', 'CIRCLE', 'HEXAGON', 'HEXAGONR30',
% 'RECTANGLE', 'SQUARE', 'ANNULUS'
%
% OUTPUT
% results = The Zernike polynomial (n,m,sc) value for the given pixel (rho, theta).
%
% calculate radial part Rnm
Rnm = zeros(size(rho));
for s=0:(n-m)/2
numerator = (-1)^s * factorial(n-s);
denominator = factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s);
Rnm = Rnm + (numerator / denominator) * rho.^(n-2*s);
end % for s statement
% 3 cases. sc=' ', 's', 'c'.
theFactor = 1;
switch sc
case ' '
% means m=0
theFactor = sqrt(n+1);
result = Rnm;
case 's'
theFactor = sqrt(2*(n+1));
result = Rnm .* sin(m*theta);
case 'c'
theFactor = sqrt(2*(n+1));
result = Rnm .* cos(m*theta);
end % switch sc
switch ZernikeDef
case {'MAHAJAN', 'NOLL', ...
'HEXAGON', 'HEXAGONR30', 'RECTANGLE', 'SQUARE', 'ANNULUS'}
result = theFactor * result;
end % switch
end % functon calculateZernike
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The following functions are used to calculate (n,m,sc) from a j-order
% number.
%
%
function [n, m, sc] = convertj(j, ztype)
%CONVERTJ Convert ordering parameter j to n, m, sc parameters.
%
% function [n, m, sc] = convertj(j, ztype)
%
% INPUT:
% j = Zernike polynomial order number. j >= 0 and integer.
% For 'FRINGE' 1 <= j <= 37.
% For 'ISO' 0 <= j <= 35.
% ztype = type of odering: 'FRINGE', 'ISO', 'WYANT', 'B&W',
% 'STANDARD', 'MAHAJAN', 'NOLL',
% 'HEXAGON', 'HEXAGONR30', 'ELLIPSE',
% 'RECTANGLE', 'SQUARE', 'ANNULUS'
%
% OUTPUT:
% n = radial order.
% m = azimuthal order
% sc = indicate whether sine or cosine or neither factor
% 's' means sin() factor
% 'c' means cos() factor
% ' ' means m = 0 so has no sin() or cos() factor
%
%
% DANGER: Validation of input is NOT performed.
% This function should only be called from ZernikeCalc
% which does the input validation.
%
switch ztype
case 'FRINGE'
[n,m,sc] = fringe(j);
case {'ISO', 'WYANT'}
[n,m,sc] = fringe(j+1);
case {'B&W', 'STANDARD'}
[n,m,sc] = bw(j);
case {'MAHAJAN','NOLL', ...
'HEXAGON', 'HEXAGONR30', 'ELLIPSE', ...
'RECTANGLE', 'SQUARE', 'ANNULUS'}
[n,m,sc] = mahajan(j);
end % switch ztype
end % function convertj
function [n, m, sc] = fringe(j)
%
% Note that 'FRINGE', 'ISO', 'WYANT' use this function
% to assign (n, m) pairs to the j values.
%
% sc = 's' = sin(), sc = 'c' = cos(), sc = ' ' for neither.
%
switch j
case 1
n = 0; m = 0; sc = ' ';
case 2
n = 1; m = 1; sc = 'c';
case 3
n = 1; m = 1; sc = 's';
case 4
n = 2; m = 0; sc = ' ';
case 5
n = 2; m = 2; sc = 'c';
case 6
n = 2; m = 2; sc = 's';
case 7
n = 3; m = 1; sc = 'c';
case 8
n = 3; m = 1; sc = 's';
case 9
n = 4; m = 0; sc = ' ';
case 10
n = 3; m = 3; sc = 'c';
case 11
n = 3; m = 3; sc = 's';
case 12
n = 4; m = 2; sc = 'c';
case 13
n = 4; m = 2; sc = 's';
case 14
n = 5; m = 1; sc = 'c';
case 15
n = 5; m = 1; sc = 's';
case 16
n = 6; m = 0; sc = ' ';
case 17
n = 4; m = 4; sc = 'c';
case 18
n = 4; m = 4; sc = 's';
case 19
n = 5; m = 3; sc = 'c';
case 20
n = 5; m = 3; sc = 's';
case 21
n = 6; m = 2; sc = 'c';
case 22
n = 6; m = 2; sc = 's';
case 23
n = 7; m = 1; sc = 'c';
case 24
n = 7; m = 1; sc = 's';
case 25
n = 8; m = 0; sc = ' ';
case 26
n = 5; m = 5; sc = 'c';
case 27
n = 5; m = 5; sc = 's';
case 28
n = 6; m = 4; sc = 'c';
case 29
n = 6; m = 4; sc = 's';
case 30
n = 7; m = 3; sc = 'c';
case 31
n = 7; m = 3; sc = 's';
case 32
n = 8; m = 2; sc = 'c';
case 33
n = 8; m = 2; sc = 's';
case 34
n = 9; m = 1; sc = 'c';
case 35
n = 9; m = 1; sc = 's';
case 36
n = 10; m = 0; sc = ' ';
case 37
n = 12; m = 0; sc = ' ';
end % switch j
end % function fringe
function [n, m, sc] = bw(j)
% sc = 's' = sin(), sc = 'c' = cos(), sc = ' ' for neither.
% sc = 1 = sin(), sc = 2 = cos().
% calculate the n value
n1 = (-1 + sqrt(1 + 8 * j)) / 2;
n = floor(n1);
if n1 == n
n = n - 1;
end % if statement
% calculate the m value
k = (n+1)*(n+2)/2;
d = k - j;
m1 = n - 2*d;
m = abs(m1);
% calculate the sc value
sc = ' ';
if m1 > 0
sc = 's';
end % if statement
if m1 < 0
sc = 'c';
end % if statement
end % function bw
function [n, m, sc] = mahajan(j)
% sc = 's' = sin(), sc = 'c' = cos(), sc = ' ' for neither.
% sc = 1 = sin(), sc = 2 = cos().
% calculate the n value
n1 = (-1 + sqrt(1 + 8 * j)) / 2;
n = floor(n1);
if n1 == n
n = n - 1;
end % if statement
% calculate the m value
k = (n+1)*(n+2)/2;
m = n - 2 * floor((k - j)/2);
% calculate the sc value
sc = ' ';
if (m ~= 0) && (mod(j,2) ~= 0)
% m ~= 0 and j odd
sc = 's';
end % if statement
if (m ~= 0) && (mod(j,2) == 0)
% m ~= 0 and j even
sc = 'c';
end % if statement
end % function mahajan
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ok = validateInput(ZernikeList, Zdata, mask, ...
ZernikeDef, ShapeParameter, ...
centerRow, centerCol, radiusInPixels, MMESorC)
%
% Validates the input.
% See above function definition for definition of input.
%
% OUTPUT
% ok = true.
%
ok = true;
% Check validity of input parameters.
theIOCSize = size(Zdata);
theMaskSize = size(mask);
if (theIOCSize(1,2) > 1) && (sum(theIOCSize == theMaskSize) ~= 2)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Surface data and mask matrices are not the same size.');
throwAsCaller(ME);
end % if statement
if centerRow < 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: centerRow must be positive.');
throwAsCaller(ME);
end % if statement
if centerCol < 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: centerCol must be positive.');
throwAsCaller(ME);
end % if statement
if radiusInPixels < 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: radiusInPixels must be positive.');
throwAsCaller(ME);
end % if statement
hlda = mask == 0;
hldb = mask == 1;
if sum(sum(hlda + hldb)) ~= (theMaskSize(1,1)*theMaskSize(1,2))
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: mask matrix must contain 0 or 1 only.');
throwAsCaller(ME);
end % if statement
% Now for the fun stuff: Validating the Zernike ordering.
switch ZernikeDef
case {'FRINGE', 'ISO', 'WYANT' 'MAHAJAN', 'NOLL', 'B&W', 'STANDARD', ...
'HEXAGON', 'HEXAGONR30', 'ELLIPSE', 'RECTANGLE', ...
'SQUARE', 'ANNULUS'}
% These are the valid values.
otherwise
% ZernikeType is not valid
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ZernikeDef is not valid.');
throwAsCaller(ME);
end % switch statement
theSize = size(ZernikeList);
rows = theSize(1,1);
cols = theSize(1,2);
hld1 = sum(sum(zeros(rows, cols) == (abs(ZernikeList) - floor(abs(ZernikeList)))));
if hld1 ~= rows*cols
% a number in ZernikeList is not an integer
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ZernikeList must contain only integers.');
throwAsCaller(ME);
end % if statement
if rows == 1
% j ordering
if sum(abs(ZernikeList(1, :)) ~= ZernikeList(1, :)) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ZernikeList j values must be positive or 0.');
throwAsCaller(ME);
end % if statement
switch ZernikeDef
case 'FRINGE'
if sum(ZernikeList(1, :) < 1) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: FRINGE order value j must be greater than 0.');
throwAsCaller(ME);
end % if statement
if sum(ZernikeList(1, :) > 37) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: FRINGE order value j must be less than 38.');
throwAsCaller(ME);
end % if statement
case 'ISO'
if sum(ZernikeList(1, :) > 35) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ISO order value j out of range.');
throwAsCaller(ME);
end % if statement
case 'WYANT'
if sum(ZernikeList(1, :) > 36) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: WYANT order value j out of range.');
throwAsCaller(ME);
end % if statement
case {'MAHAJAN', 'NOLL', 'B&W', 'STANDARD', ...
'HEXAGON', 'HEXAGONR30', 'ELLIPSE', 'RECTANGLE', ...
'SQUARE', 'ANNULUS'}
if sum(ZernikeList(1, :) == 0) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Zernike order value j can not be zero.');
throwAsCaller(ME);
end % if statement
end % switch statement
else
% (n, m) specification
if sum(abs(ZernikeList(1, :)) ~= ZernikeList(1, :)) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ZernikeList n values must be positive or 0.');
throwAsCaller(ME);
end % if statement
if sum(abs(ZernikeList(2, :)) > ZernikeList(1, :)) ~= 0
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: abs(m) values must not exceed n values.');
throwAsCaller(ME);
end % if statement
switch ZernikeDef
case 'FRINGE'
if rows > 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: FRINGE can only be specified with j ordering.');
throwAsCaller(ME);
end % if statement
case 'ISO'
if rows > 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ISO can only be specified with j ordering.');
throwAsCaller(ME);
end % if statement
case 'WYANT'
if rows > 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: WYANT can only be specified with j ordering.');
throwAsCaller(ME);
end % if statement
case {'HEXAGON', 'HEXAGONR30', 'ELLIPSE', 'RECTANGLE', ...
'SQUARE', 'ANNULUS'}
if rows > 1
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: ZernikeDef value can only be specified with j ordering.');
throwAsCaller(ME);
end % if statement
otherwise
if (sum(mod(ZernikeList(1, :) - abs(ZernikeList(2, :)), 2)) ~= 0)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: n-abs(m) must be even.');
throwAsCaller(ME);
end % if statement
end % switch
switch MMESorC
case {'-m=sin', '-m=cos'}
otherwise
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Sign convention value for (n, m) not valid.');
throwAsCaller(ME);
end % switch
end % (n,m) order specified
% Validate ZernikeShape parameter
switch ZernikeDef
case {'ELLIPSE', 'RECTANGLE','ANNULUS'}
% Check that the ShapeParameter is valid
if (ShapeParameter <= 0) || (ShapeParameter >= 1)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: The ShapeParameter is out of range: (0,1).');
throwAsCaller(ME);
end % if statement
end % switch ZernikeDef
end % function validateInput
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This section is for implementing non-circular Zernike polynomials.
% This immplements the tables in Mahajan's paper "Orthonomal polynomials
% in wavefront analysis: analytical soultion," J. Opt. Soc. Am. A, 24(9)
% pp. 2994-3016 2007
%
function result = Z(j, rho, theta)
%
% Caluclates the Mahajan Zernike polynomial value.
%
[n, m, sc] = convertj(j, 'MAHAJAN');
result = calculateZernike(n, m, sc, rho, theta, 'MAHAJAN');
end % function Z
function result = ZHexagon(j, rho, theta)
% Orthonormal hexagonal polynomials
if (j < 1) || (j > 45)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Hexagon order j out of range.');
throwAsCaller(ME);
end % if statement
switch j
case 1
result = Z(1,rho,theta);
case 2
result = sqrt(6/5)*Z(2,rho,theta);
case 3
result = sqrt(6/5)*Z(3,rho,theta);
case 4
result = sqrt(5/43)*Z(1,rho,theta)+(2*sqrt(15/43))*Z(4,rho,theta);
case 5
result = sqrt(10/7)*Z(5,rho,theta);
case 6
result = sqrt(10/7)*Z(6,rho,theta);
case 7
result = 16*sqrt(14/11055)*Z(3,rho,theta)+10*sqrt(35/2211)*Z(7,rho,theta);
case 8
result = 16*sqrt(14/11055)*Z(2,rho,theta)+10*sqrt(35/2211)*Z(8,rho,theta);
case 9
result = (2*sqrt(5)/3)*Z(9,rho,theta);
case 10
result = 2*sqrt(35/103)*Z(10,rho,theta);
case 11
result = (521/sqrt(1072205))*Z(1,rho,theta) ...
+88*sqrt(15/214441)*Z(4,rho,theta) ...
+14*sqrt(43/4987)*Z(11,rho,theta);
case 12
result = 225*sqrt(6/492583)*Z(6,rho,theta)+42*sqrt(70/70369)*Z(12,rho,theta);
case 13
result = 225*sqrt(6/492583)*Z(5,rho,theta)+42*sqrt(70/70369)*Z(13,rho,theta);
case 14
result = -2525*sqrt(14/297774543)*Z(6,rho,theta) ...
-(1495*sqrt(70/99258181)/3)*Z(12,rho,theta) ...
+(sqrt(378910/18337)/3)*Z(14,rho,theta);
case 15
result = 2525*sqrt(14/297774543)*Z(5,rho,theta) ...
+(1495*sqrt(70/99258181)/3)*Z(13,rho,theta) ...
+(sqrt(378910/18337)/3)*Z(15,rho,theta);
case 16
result = 30857*sqrt(2/3268147641)*Z(2,rho,theta) ...
+(49168/sqrt(3268147641))*Z(8,rho,theta) ...
+42*sqrt(1474/1478131)*Z(16,rho,theta);
case 17
result = 30857*sqrt(2/3268147641)*Z(3,rho,theta) ...
+(49168/sqrt(3268147641))*Z(7,rho,theta) ...
+42*sqrt(1474/1478131)*Z(17,rho,theta);
case 18
result = 386*sqrt(770/295894589)*Z(10,rho,theta) ...
+6*sqrt(118965/2872763)*Z(18,rho,theta);
case 19
result = 6*sqrt(10/97)*Z(9,rho,theta) ...
+14*sqrt(5/291)*Z(19,rho,theta);
case 20
result = -0.71499593*Z(2,rho,theta) ...
-0.72488884*Z(8,rho,theta)-0.46636441*Z(16,rho,theta) ...
+1.72029850*Z(20,rho,theta);
case 21
result = 0.71499594*Z(3,rho,theta)+0.72488884*Z(7,rho,theta) ...
+0.46636441*Z(17,rho,theta)+1.72029850*Z(21,rho,theta);
case 22
result = 0.58113135*Z(1,rho,theta)+0.89024136*Z(4,rho,theta) ...
+0.89044507*Z(11,rho,theta)+1.32320623*Z(22,rho,theta);
case 23
result = 1.15667686*Z(5,rho,theta)+1.10775599*Z(13,rho,theta) ...
+0.43375081*Z(15,rho,theta)+1.39889072*Z(23,rho,theta);
case 24
result = 1.15667686*Z(6,rho,theta)+1.10775599*Z(12,rho,theta) ...
-0.43375081*Z(14,rho,theta)+1.39889072*Z(24,rho,theta);
case 25
result = 1.31832566*Z(5,rho,theta)+1.14465174*Z(13,rho,theta) ...
+1.94724032*Z(15,rho,theta)+0.67629133*Z(23,rho,theta) ...
+1.75496998*Z(25,rho,theta);
case 26
result = -1.31832566*Z(6,rho,theta)-1.14465174*Z(12,rho,theta) ...
+1.94724032*Z(14,rho,theta)-0.67629133*Z(24,rho,theta) ...
+1.75496998*Z(26,rho,theta);
case 27
result = 2*sqrt(77/93)*Z(27,rho,theta);
case 28
result = -1.07362889*Z(1,rho,theta)-1.52546162*Z(4,rho,theta) ...
-1.28216588*Z(11,rho,theta)-0.70446308*Z(22,rho,theta) ...
+2.09532473*Z(28,rho,theta);
case 29
result = 0.97998834*Z(3,rho,theta)+1.16162002*Z(7,rho,theta) ...
+1.04573775*Z(17,rho,theta)+0.40808953*Z(21,rho,theta) ...
+1.36410394*Z(29,rho,theta);
case 30
result = 0.97998834*Z(2,rho,theta)+1.16162002*Z(8,rho,theta) ...
+1.04573775*Z(16,rho,theta)-0.40808953*Z(20,rho,theta) ...
+1.36410394*Z(30,rho,theta);
case 31
result = 3.63513758*Z(9,rho,theta)+2.92084414*Z(19,rho,theta) ...
+2.11189625*Z(31,rho,theta);
case 32
result = 0.69734874*Z(10,rho,theta)+0.67589740*Z(18,rho,theta) ...
+1.22484055*Z(32,rho,theta);
case 33
result = 1.56189763*Z(3,rho,theta)+1.69985309*Z(7,rho,theta) ...
+1.29338869*Z(17,rho,theta)+2.57680871*Z(21,rho,theta) ...
+0.67653220*Z(29,rho,theta)+1.95719339*Z(33,rho,theta);
case 34
result = -1.56189763*Z(2,rho,theta)-1.69985309*Z(8,rho,theta) ...
-1.29338869*Z(16,rho,theta)+2.57680871*Z(20,rho,theta) ...
-0.67653220*Z(30,rho,theta)+1.95719339*Z(34,rho,theta);
case 35
result = -1.63832594*Z(3,rho,theta)-1.74759886*Z(7,rho,theta) ...
-1.27572528*Z(17,rho,theta)-0.77446421*Z(21,rho,theta) ...
-0.60947360*Z(29,rho,theta)-0.36228537*Z(33,rho,theta) ...
+2.24453237*Z(35,rho,theta);
case 36
result = -1.63832594*Z(2,rho,theta)-1.74759886*Z(8,rho,theta) ...
-1.27572528*Z(16,rho,theta)+0.77446421*Z(20,rho,theta) ...
-0.60947360*Z(30,rho,theta)+0.36228537*Z(34,rho,theta) ...
+2.24453237*Z(36,rho,theta);
case 37
result = 0.82154671*Z(1,rho,theta)+1.27988084*Z(4,rho,theta) ...
+1.32912377*Z(11,rho,theta)+1.11636637*Z(22,rho,theta) ...
-0.54097038*Z(28,rho,theta)+1.37406534*Z(37,rho,theta);
case 38
result = 1.54526522*Z(6,rho,theta)+1.57785242*Z(12,rho,theta) ...
-0.89280081*Z(14,rho,theta)+1.28876176*Z(24,rho,theta) ...
-0.60514082*Z(26,rho,theta)+1.43097780*Z(38,rho,theta);
case 39
result = 1.54526522*Z(5,rho,theta)+1.57785242*Z(13,rho,theta) ...
+0.89280081*Z(15,rho,theta)+1.28876176*Z(23,rho,theta) ...
+0.60514082*Z(25,rho,theta)+1.43097780*Z(39,rho,theta);
case 40
result = -2.51783502*Z(6,rho,theta)-2.38279377*Z(12,rho,theta) ...
+3.42458933*Z(14,rho,theta)-1.69296616*Z(24,rho,theta) ...
+2.56612920*Z(26,rho,theta)-0.85703819*Z(38,rho,theta) ...
+1.89468756*Z(40,rho,theta);
case 41
result = 2.51783502*Z(5,rho,theta)+2.38279377*Z(13,rho,theta) ...
+3.42458933*Z(15,rho,theta)+1.69296616*Z(23,rho,theta) ...
+2.56612920*Z(25,rho,theta)+0.85703819*Z(39,rho,theta) ...
+1.89468756*Z(41,rho,theta);
case 42
result = -2.72919646*Z(1,rho,theta)-4.02313214*Z(4,rho,theta) ...
-3.69899239*Z(11,rho,theta)-2.49229315*Z(22,rho,theta) ...
+4.36717121*Z(28,rho,theta)-1.13485132*Z(37,rho,theta) ...
+2.52330106*Z(42,rho,theta);
case 43
result = 1362*sqrt(77/20334667)*Z(27,rho,theta)+(260/3) ...
*sqrt(341/655957)*Z(43,rho,theta);
case 44
result = -2.76678413*Z(6,rho,theta)-2.50005278*Z(12,rho,theta) ...
+1.48041348*Z(14,rho,theta)-1.62947374*Z(24,rho,theta) ...
+0.95864121*Z(26,rho,theta)-0.69034812*Z(38,rho,theta) ...
+0.40743941*Z(40,rho,theta)+2.56965299*Z(44,rho,theta);
case 45
result = -2.76678413*Z(5,rho,theta)-2.50005278*Z(13,rho,theta) ...
-1.48041348*Z(15,rho,theta)-1.62947374*Z(23,rho,theta) ...
-0.95864121*Z(25,rho,theta)-0.69034812*Z(39,rho,theta) ...
-0.40743941*Z(41,rho,theta)+2.56965299*Z(45,rho,theta);
end % switch j statement
end % function ZHexagon
function result = ZHexagonR30(j, rho, theta)
% Orthonormal hexagonal polynomials (hexagon rotated 30 degrees)
if (j < 1) || (j > 28)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Hexagon R30 order j out of range.');
throwAsCaller(ME);
end % if statement
switch j
case 1
result = Z(1,rho,theta);
case 2
result = sqrt(6/5)*Z(2, rho, theta);
case 3
result = sqrt(6/5)*Z(3, rho, theta);
case 4
result = sqrt(5/43)*Z(1, rho, theta)+2*sqrt(15/43)*Z(4, rho, theta);
case 5
result = sqrt(10/7)*Z(5,rho,theta);
case 6
result = sqrt(10/7)*Z(6,rho,theta);
case 7
result = 16*sqrt(14/11055)*Z(3,rho,theta)+10*sqrt(35/2211)*Z(7,rho,theta);
case 8
result = 16*sqrt(14/11055)*Z(2,rho,theta)+10*sqrt(35/2211)*Z(8,rho,theta);
case 9
result = 2*sqrt(35/103)*Z(9,rho,theta);
case 10
result = (2*sqrt(5)/3)*Z(10,rho,theta);
case 11
result = (521/sqrt(1072205))*Z(1,rho,theta)+88*sqrt(15/214441)*Z(4,rho,theta) ...
+ 14*sqrt(43/4987)*Z(11,rho,theta);
case 12
result = 225*sqrt(6/492583)*Z(6,rho,theta)+42*sqrt(70/70369)*Z(12,rho,theta);
case 13
result = 225*sqrt(6/492583)*Z(5,rho,theta)+42*sqrt(70/70369)*Z(13,rho,theta);
case 14
result = 2525*sqrt(14/297774543)*Z(6,rho,theta)+(1495*sqrt(70/99258181)/3)*Z(12,rho,theta) ...
+ (sqrt(378910/18337)/3)*Z(14,rho,theta);
case 15
result = -2525*sqrt(14/297774543)*Z(5,rho,theta)-(1495*sqrt(70/99258181)/3)*Z(13,rho,theta) ...
+ (sqrt(378910/18337)/3)*Z(15,rho,theta);
case 16
result = 30857*sqrt(2/3268147641)*Z(2,rho,theta)+49168/sqrt(3268147641)*Z(8,rho,theta) ...
+ 42*sqrt(1474/1478131)*Z(16,rho,theta);
case 17
result = 30857*sqrt(2/3268147641)*Z(3,rho,theta)+(49168/sqrt(3268147641))*Z(7,rho,theta) ...
+ 42*sqrt(1474/1478131)*Z(17,rho,theta);
case 18
result = 6*sqrt(10/97)*Z(10,rho,theta)+14*sqrt(5/291)*Z(18,rho,theta);
case 19
result = 386*sqrt(770/295894589)*Z(9,rho,theta)+6*sqrt(118965/2872763)*Z(19,rho,theta);
case 20
result = 0.71499593*Z(2,rho,theta)+0.72488884*Z(8,rho,theta) ...
+ 0.46636441*Z(16,rho,theta)+1.72029850*Z(20,rho,theta);
case 21
result = -0.71499593*Z(3,rho,theta)-0.72488884*Z(7,rho,theta) ...
- 0.46636441*Z(17,rho,theta)+1.72029850*Z(21,rho,theta);
case 22
result = 0.58113135*Z(1,rho,theta)+0.89024136*Z(4,rho,theta) ...
+ 0.89044507*Z(11,rho,theta)+1.32320623*Z(22,rho,theta);
case 23
result = 1.15667686*Z(5,rho,theta)+1.10775599*Z(13,rho,theta) ...
- 0.43375081*Z(15,rho,theta)+1.39889072*Z(23,rho,theta);
case 24
result = 1.15667686*Z(6,rho,theta)+1.10775599*Z(12,rho,theta) ...
+ 0.43375081*Z(14,rho,theta)+1.39889072*Z(24,rho,theta);
case 25
result = -1.31832566*Z(5,rho,theta)-1.14465174*Z(13,rho,theta) ...
+ 1.94724032*Z(15,rho,theta)-0.67629133*Z(23,rho,theta)+1.75496998*Z(25,rho,theta);
case 26
result = 1.31832566*Z(6,rho,theta)+1.14465174*Z(12,rho,theta) ...
+ 1.94724032*Z(14,rho,theta)+0.67629133*Z(24,rho,theta)+1.75496998*Z(26,rho,theta);
case 27
result = 1.81984283*Z(27,rho,theta);
case 28
result = 1.07362889*Z(1,rho,theta)+1.52546162*Z(4,rho,theta) ...
+ 1.28216588*Z(11,rho,theta)+0.70446308*Z(22,rho,theta)+2.09532473*Z(28,rho,theta);
end % switch j
end % function ZHexagonR30
function result = ZEllipse(j, rho, theta, b)
% Orthonormal elliptical polynomials
if (j < 1) || (j > 15)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Elliptical order j out of range.');
throwAsCaller(ME);
end % if statement
alpha = sqrt(45-60*b^2+94*b^4-60*b^6+45*b^8);
beta = sqrt(1575 - 4800*b^2 + 12020*b^4 - 17280*b^6 + 21066*b^8 - 17280*b^10 ...
+ 12020*b^12 - 4800*b^14 + 1575*b^16);
gamma = sqrt(35*b^8 - 60*b^6 + 114*b^4 - 60*b^2 + 35);
delta = sqrt(5 - 6*b^2 + 5*b^4);
switch j
case 1
result = Z(1,rho,theta);
case 2
result = Z(2,rho,theta);
case 3
result = Z(3,rho,theta)/b;
case 4
result = (sqrt(3)*(1-b^2)/sqrt(3-2*b^2+3*b^4))*Z(1,rho,theta) ...
+(2/sqrt(3-2*b^2+3*b^4))*Z(4,rho,theta);
case 5
result = Z(5,rho,theta)/b;
case 6
result = -(sqrt(3)*(3-4*b^2+b^4)/(2*b^2*sqrt(6-4*b^2+6*b^4)))*Z(1,rho,theta) ...
-3*((1-b^4)/(2*b^2*sqrt(6-4*b^2+6*b^4)))*Z(4,rho,theta) ...
+(sqrt(3-2*b^2+3*b^4)/(2*b^2))*Z(6,rho,theta);
case 7
result = (6*(1-b^2)/(b*sqrt(5-6*b^2+9*b^4)))*Z(3,rho,theta) ...
+(2*sqrt(2)/(b*sqrt(5-6*b^2+9*b^4)))*Z(7,rho,theta);
case 8
result = (2*(1-b^2)/sqrt(9-6*b^2+5*b^4))*Z(2,rho,theta) ...
+(2*sqrt(2)/sqrt(9-6*b^2+5*b^4))*Z(8,rho,theta);
case 9
result = -((5-8*b^2+3*b^4)/(b^3*sqrt(5-6*b^2+9*b^4)))*Z(3,rho,theta) ...
-((5-2*b^2-3*b^4)/(2*sqrt(2)*b^3*sqrt(5-6*b^2+9*b^4)))*Z(7,rho,theta) ...
+(sqrt(5-6*b^2+9*b^4)/(2*sqrt(2)*b^3))*Z(9,rho,theta);
case 10
result = -((3-4*b^2+b^4)/(b^2*sqrt(9-6*b^2+5*b^4)))*Z(2,rho,theta) ...
-((3+2*b^2-5*b^4)/(2*sqrt(2)*b^2*sqrt(9-6*b^2+5*b^4)))*Z(8,rho,theta) ...
+(sqrt(9-6*b^2+5*b^4)/(2*sqrt(2)*b^2))*Z(10,rho,theta);
case 11
result = sqrt(5)*(7 - 10*b^2 + 3*b^4)*alpha^(-1)*Z(1,rho,theta) ...
+ 4*sqrt(15)*(1 - b^2)*alpha^(-1)*Z(4,rho,theta) ...
-2*sqrt(30)*(1 - b^2)*alpha^(-1)*Z(6,rho,theta) ...
+8*alpha^(-1)*Z(11,rho,theta);
case 12
result = (b^2-1)*(sqrt(10)/4)*alpha^(-1)*gamma^(-1)*b^(-2)* ...
(195 - 280*b^2+278*b^4-144*b^6+15*b^8)*Z(1,rho,theta) ...
+(b^2-1)*(sqrt(30)/4)*b^(-2)*alpha^(-1)*gamma^(-1)* ...
(105-100*b^2+94*b^4-20*b^6-15*b^8)*Z(4,rho,theta) ...
-(b^2-1)*(sqrt(15)/2)*b^(-2)*alpha^(-1)*gamma^(-1)* ...
(75-80*b^2+94*b^4-40*b^6+15*b^8)*Z(6,rho,theta) ...
-10*sqrt(2)*b^(-2)*alpha^(-1)*gamma^(-1)*(3-2*b^2+2*b^6-3*b^8)*Z(11,rho,theta) ...
+ alpha*b^(-2)*gamma^(-1)*Z(12,rho,theta);
case 13
result = sqrt(15)*(1 - b^2)/(b*sqrt(5 - 6*b^2 + 5*b^4))*Z(5,rho,theta) ...
+2/(b*sqrt(5 - 6*b^2 + 5*b^4))*Z(13,rho,theta);
case 14
result = (b^2 - 1)^2*(sqrt(10)/8)*(35-10*b^2-b^4)*b^(-4)*gamma^(-1)*Z(1,rho,theta) ...
+(b^2 - 1)^2*5*(sqrt(30)/16)*(7+2*b^2-b^4)*b^(-4)*gamma^(-1)*Z(4,rho,theta) ...
-(sqrt(15)/8)*(35-70*b^2+56*b^4-26*b^6+5*b^8)*b^(-4)*gamma^(-1)*Z(6,rho,theta) ...
+(b^2 - 1)^2*5*(sqrt(2)/16)*(7+10*b^2+7*b^4)*b^(-4)*gamma^(-1)*Z(11,rho,theta) ...
-(5/8)*(7-6*b^2+6*b^6-7*b^8)*b^(-4)*gamma^(-1)*Z(12,rho,theta) ...
+(1/8)*gamma*b^(-4)*Z(14,rho,theta);
case 15
result = -(sqrt(15)/4)*(5-8*b^2+3*b^4)*b^(-3)*delta^(-1)*Z(5,rho,theta) ...
+(b^4-1)*(5/4)*b^(-3)*delta^(-1)*Z(13,rho,theta) ...
+(1/2)*b^(-3)*delta*Z(15,rho,theta);
end % switch statement
end % function ZElliptical
function result = ZRectangle(j, rho, theta, a)
% Orthonormal rectangle polynomials
if (j < 1) || (j > 15)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Rectangle order j out of range.');
throwAsCaller(ME);
end % if statement
mu = sqrt(9-36*a^2+103*a^4-134*a^6+67*a^8);
nu = sqrt(49-196*a^2+330*a^4-268*a^6+134*a^8);
tau = 1/(128*nu*a^4*(1-a^2)^2);
eta = 9 - 45*a^2 + 139*a^4 - 237*a^6 + 210*a^8 - 67*a^10;
switch j
case 1
result = Z(1,rho,theta);
case 2
result = (sqrt(3)/(2*a))*Z(2,rho,theta);
case 3
result = (sqrt(3)/(2*sqrt(1-a^2)))*Z(3,rho,theta);
case 4
result = (sqrt(5)/(4*sqrt(1-2*a^2+2*a^4)))*(Z(1,rho,theta)+sqrt(3)*Z(4,rho,theta));
case 5
result = (sqrt(3/2)/(2*a*sqrt(1-a^2)))*Z(5,rho,theta);
case 6
result = (sqrt(5)/(8*a^2*(1-a^2)*sqrt(1-2*a^2+2*a^4)))*((3-10*a^2+12*a^4-8*a^6)*Z(1,rho,theta) ...
+sqrt(3)*(1-2*a^2)*Z(4,rho,theta) + sqrt(6)*(1-2*a^2+2*a^4)*Z(6,rho,theta));
case 7
result = (sqrt(21)/(4*sqrt(2)*sqrt(27-81*a^2+116*a^4-62*a^6))) ...
* (sqrt(2)*(1+4*a^2)*Z(3,rho,theta)+5*Z(7,rho,theta));
case 8
result = (sqrt(21)/(4*sqrt(2)*a*sqrt(35-70*a^2+62*a^4))) ...
* (sqrt(2)*(5-4*a^2)*Z(2,rho,theta)+5*Z(8,rho,theta));
case 9
result = (sqrt(5/2)*sqrt((27-54*a^2+62*a^4)/(1-a^2)) ...
/ (16*a^2*(27-81*a^2+116*a^4-62*a^6))) ...
* (2*sqrt(2)*(9-36*a^2+52*a^4-60*a^6)*Z(3,rho,theta) ...
+ (9-18*a^2-26*a^4)*Z(7,rho,theta) ...
+ (27-54*a^2+62*a^4)*Z(9,rho,theta));
case 10
result = (sqrt(5/2)/(16*a^3*(1-a^2)*sqrt(35-70*a^2+62*a^4))) ...
* (2*sqrt(2)*(35-112*a^2+128*a^4-60*a^6)*Z(2,rho,theta) ...
+ (35-70*a^2+26*a^4)*Z(8,rho,theta)+(35-70*a^2+62*a^4)*Z(10,rho,theta));
case 11
result = (1/(16*mu))*(8*(3+4*a^2-4*a^4)*Z(1,rho,theta)+25*sqrt(3)*Z(4,rho,theta) ...
+ 10*sqrt(6)*(1-2*a^2)*Z(6,rho,theta) ...
+ 21*sqrt(5)*Z(11,rho,theta));
case 12
alpha = (134*a^8 - 268*a^6 + 330*a^4 - 196*a^2 + 49)^(1/2);
beta = (67*a^8 - 134*a^6 + 103*a^4 - 36*a^2 + 9)^(1/2);
result = ((3234*a^10 - 8085*a^8 + 8508*a^6 - 4677*a^4 + 1650*a^2 - 315)/((16*a^4 - 16*a^2)*beta*alpha))* ...
Z(1,rho,theta) + ...
-(15*3^(1/2)*(14-74*a^2+205*a^4-360*a^6+335*a^8-134*a^10))/(16*a^2*(a^2 - 1)*beta*alpha)* ...
Z(4,rho,theta) + ...
-(6^(1/2)*(2340-3975*a^2+3975*a^4+(525)/(a^2*(a^2-1)))/(64*alpha*beta)) * ...
Z(6,rho,theta) + ...
-(63*sqrt(5)*(1-4*a^2+6*a^4-4*a^6)/(16*a^2*(a^2-1)*alpha*beta))* ...
Z(11,rho,theta) + ...
-(21*sqrt(10)*beta/(64*a^2*(a^2-1)*alpha))* ...
Z(12,rho,theta);
case 13
result = (sqrt(21)/(16*sqrt(2)*a*sqrt(1-3*a^2+4*a^4-2*a^6))) ...
* (sqrt(3)*Z(5,rho,theta) + sqrt(5)*Z(13,rho,theta));
case 14
result = tau*(6*(245-1400*a^2+3378*a^4-4452*a^6+3466*a^8-1488*a^10+496*a^12)*Z(1,rho,theta) ...
+ 15*sqrt(3)*(49-252*a^2+522*a^4-540*a^6+270*a^8)*Z(4,rho,theta) ...
+ 15*sqrt(6)*(49-252*a^2+534*a^4-596*a^6+360*a^8-144*a^10)*Z(6,rho,theta) ...
+ 3*sqrt(5)*(49-196*a^2+282*a^4-172*a^6+86*a^8)*Z(11,rho,theta) ...
+ 147*sqrt(10)*(1-4*a^2+6*a^4-4*a^6)*Z(12,rho,theta) ...
+ 3*sqrt(10)*nu^2*Z(14,rho,theta));
case 15
result = (1/(32*a^3*(1-a^2)*sqrt(1-3*a^2+4*a^4-2*a^6))) ...
* (3*sqrt(7/2)*(5-18*a^2+24*a^4-16*a^6)*Z(5,rho,theta) ...
+ sqrt(105/2)*(1-2*a^2)*Z(13,rho,theta) ...
+ sqrt(210)*(1-2*a^2+2*a^4)*Z(15,rho,theta));
end % switch statement
end % function ZRectangle
function result = ZSquare(j, rho, theta)
% Orthonormal square polynomials
if (j < 1) || (j > 45)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Square order j out of range.');
throwAsCaller(ME);
end % if statement
switch j
case 1
result = Z(1,rho,theta);
case 2
result = sqrt(3/2)*Z(2,rho,theta);
case 3
result = sqrt(3/2)*Z(3,rho,theta);
case 4
result = (sqrt(5/2)/2)*Z(1,rho,theta)+(sqrt(15/2)/2)*Z(4,rho,theta);
case 5
result = sqrt(3/2)*Z(5,rho,theta);
case 6
result = (sqrt(15)/2)*Z(6,rho,theta);
case 7
result = (3*sqrt(21/31)/2)*Z(3,rho,theta)+(5*sqrt(21/62)/2)*Z(7,rho,theta);
case 8
result = (3*sqrt(21/31)/2)*Z(2,rho,theta)+(5*sqrt(21/62)/2)*Z(8,rho,theta);
case 9
result = -(7*sqrt(5/31)/2)*Z(3,rho,theta)-(13*sqrt(5/62)/4)*Z(7,rho,theta) ...
+(sqrt(155/2)/4)*Z(9,rho,theta);
case 10
result = (7*sqrt(5/31)/2)*Z(2,rho,theta)+(13*sqrt(5/62)/4)*Z(8,rho,theta) ...
+(sqrt(155/2)/4)*Z(10,rho,theta);
case 11
result = (8/sqrt(67))*Z(1,rho,theta)+(25*sqrt(3/67)/4)*Z(4,rho,theta) ...
+(21*sqrt(5/67)/4)*Z(11,rho,theta);
case 12
result = (45*sqrt(3)/16)*Z(6,rho,theta)+(21*sqrt(5)/16)*Z(12,rho,theta);
case 13
result = (3*sqrt(7)/8)*Z(5,rho,theta)+(sqrt(105)/8)*Z(13,rho,theta);
case 14
result = 261/(8*sqrt(134))*Z(1,rho,theta)+(345*sqrt(3/134)/16)*Z(4,rho,theta) ...
+(129*sqrt(5/134)/16)*Z(11,rho,theta)+(3*sqrt(335)/16)*Z(14,rho,theta);
case 15
result = (sqrt(105)/4)*Z(15,rho,theta);
case 16
result = ((41*sqrt(55/1966))/4)*Z(2,rho,theta)+((29*sqrt(55/983))/4)*Z(8,rho,theta) ...
+((11*sqrt(55/983))/4)*Z(10,rho,theta)+((21*sqrt(165/1966))/4)*Z(16,rho,theta);
case 17
result = ((41*sqrt(55/1966))/4)*Z(3,rho,theta)+((29*sqrt(55/983))/4)*Z(7,rho,theta) ...
-((11*sqrt(55/983))/4)*Z(9,rho,theta)+((21*sqrt(165/1966))/4)*Z(17,rho,theta);
case 18
result = ((34843*sqrt(3/844397))/16)*Z(2,rho,theta)+((20761*sqrt(3/1688794))/8)*Z(8,rho,theta) ...
+((32077*sqrt(3/1688794))/8)*Z(10,rho,theta) ...
+(22323/(16*sqrt(844397)))*Z(16,rho,theta)+((21*sqrt(983/859))/8)*Z(18,rho,theta);
case 19
result = -((34843*sqrt(3/844397))/16)*Z(3,rho,theta)-((20761*sqrt(3/1688794))/8)*Z(7,rho,theta) ...
+((32077*sqrt(3/1688794))/8)*Z(9,rho,theta) ...
-(22323/(16*sqrt(844397)))*Z(17,rho,theta)+((21*sqrt(983/859))/8)*Z(19,rho,theta);
case 20
result = ((1975*sqrt(7/859))/32)*Z(2,rho,theta)+((557*sqrt(7/1718))/8)*Z(8,rho,theta) ...
+((377*sqrt(7/1718))/8)*Z(10,rho,theta)+((349*sqrt(21/859))/32)*Z(16,rho,theta) ...
+((239*sqrt(21/859))/32)*Z(18,rho,theta)+(sqrt(18039)/32)*Z(20,rho,theta);
case 21
result = ((1975*sqrt(7/859))/32)*Z(3,rho,theta)+((557*sqrt(7/1718))/8)*Z(7,rho,theta) ...
-((377*sqrt(7/1718))/8)*Z(9,rho,theta)+((349*sqrt(21/859))/32)*Z(17,rho,theta) ...
-((239*sqrt(21/859))/32)*Z(19,rho,theta)+(sqrt(18039)/32)*Z(21,rho,theta);
case 22
result = ((77*sqrt(65/849))/16)*Z(1,rho,theta)+((65*sqrt(65/283))/16)*Z(4,rho,theta) ...
+((75*sqrt(39/283))/16)*Z(11,rho,theta)+((5*sqrt(39/566))/2)*Z(14,rho,theta) ...
+((11*sqrt(1365/283))/16)*Z(22,rho,theta);
case 23
result = ((51*sqrt(11/7846))/2)*Z(5,rho,theta)+(7*sqrt(165/7846))*Z(13,rho,theta) ...
+((15*sqrt(231/7846))/2)*Z(23,rho,theta);
case 24
result = ((147*sqrt(195/2698))/4)*Z(6,rho,theta)+(105*sqrt(13/2698))*Z(12,rho,theta) ...
+((33*sqrt(455/2698))/4)*Z(24,rho,theta);
case 25
result = ((7*sqrt(165))/16)*Z(15,rho,theta)+((3*sqrt(231))/16)*Z(25,rho,theta);
case 26
result = (20525/(64*sqrt(849)))*Z(1,rho,theta)+(15077/(64*sqrt(283)))*Z(4,rho,theta) ...
+((2565*sqrt(15/283))/64)*Z(11,rho,theta) ...
+((2665*sqrt(15/566))/32)*Z(14,rho,theta)+((749*sqrt(21/283))/64)*Z(22,rho,theta) ...
+((3*sqrt(5943/2))/32)*Z(26,rho,theta);
case 27
result = ((4911*sqrt(3/3923))/32)*Z(5,rho,theta)+((2429*sqrt(5/3923))/32)*Z(13,rho,theta) ...
+((641*sqrt(7/3923))/32)*Z(23,rho,theta)+(sqrt(27461)/32)*Z(27,rho,theta);
case 28
result = ((16877*sqrt(3/2698))/32)*Z(6,rho,theta)+((8295*sqrt(5/2698))/32)*Z(12,rho,theta) ...
+((2247*sqrt(7/2698))/32)*Z(24,rho,theta) ...
+((3*sqrt(9443/2))/32)*Z(28,rho,theta);
case 29
result = 2.42764289*Z(3,rho,theta)+2.69721906*Z(7,rho,theta) ...
-1.56598065*Z(9,rho,theta)+2.12208902*Z(17,rho,theta) ...
-0.93135654*Z(19,rho,theta)+0.25252773*Z(21,rho,theta) ...
+1.59017528*Z(29,rho,theta);
case 30
result = 2.42764289*Z(2,rho,theta)+2.69721906*Z(8,rho,theta)+1.56598064*Z(10,rho,theta) ...
+2.12208902*Z(16,rho,theta)+0.93135653*Z(18,rho,theta) ...
+0.25252773*Z(20,rho,theta)+1.59017528*Z(30,rho,theta);
case 31
result = -9.10300982*Z(3,rho,theta)-8.79978208*Z(7,rho,theta)+10.69381427*Z(9,rho,theta) ...
-5.37383386*Z(17,rho,theta)+7.01044701*Z(19,rho,theta) ...
-1.26347273*Z(21,rho,theta)-1.90131757*Z(29,rho,theta)+3.07960207*Z(31,rho,theta);
case 32
result = 9.10300982*Z(2,rho,theta)+8.79978208*Z(8,rho,theta)+10.69381427*Z(10,rho,theta) ...
+5.37383385*Z(16,rho,theta)+7.01044701*Z(18,rho,theta) ...
+1.26347272*Z(20,rho,theta)+1.90131756*Z(30,rho,theta)+3.07960207*Z(32,rho,theta);
case 33
result = 21.39630883*Z(3,rho,theta)+19.76696884*Z(7,rho,theta)-12.70550260*Z(9,rho,theta) ...
+11.05819453*Z(17,rho,theta)-7.02178757*Z(19,rho,theta) ...
+15.80286172*Z(21,rho,theta)+3.29259996*Z(29,rho,theta) ...
-2.07602716*Z(31,rho,theta)+5.40902889*Z(33,rho,theta);
case 34
result = 21.39630883*Z(2,rho,theta)+19.76696884*Z(8,rho,theta)+12.70550260*Z(10,rho,theta) ...
+11.05819453*Z(16,rho,theta)+7.02178756*Z(18,rho,theta) ...
+15.80286172*Z(20,rho,theta)+3.29259996*Z(30,rho,theta) ...
+2.07602718*Z(32,rho,theta)+5.40902889*Z(34,rho,theta);
case 35
result = -16.54454463*Z(3,rho,theta)-14.89205549*Z(7,rho,theta)+22.18054997*Z(9,rho,theta) ...
-7.94524850*Z(17,rho,theta)+11.85458952*Z(19,rho,theta) ...
-6.18963458*Z(21,rho,theta)-2.19431442*Z(29,rho,theta)+3.24324400*Z(31,rho,theta) ...
-1.72001172*Z(33,rho,theta)+8.16384008*Z(35,rho,theta);
case 36
result = 16.54454462*Z(2,rho,theta)+14.89205549*Z(8,rho,theta)+22.18054997*Z(10,rho,theta) ...
+7.94524849*Z(16,rho,theta)+11.85458952*Z(18,rho,theta) ...
+6.18963457*Z(20,rho,theta)+2.19431441*Z(30,rho,theta)+3.24324400*Z(32,rho,theta) ...
+1.72001172*Z(34,rho,theta)+8.16384008*Z(36,rho,theta);
case 37
result = 1.75238960*Z(1,rho,theta)+2.72870567*Z(4,rho,theta)+2.76530671*Z(11,rho,theta) ...
+1.43647360*Z(14,rho,theta)+2.12459170*Z(22,rho,theta) ...
+0.92450043*Z(26,rho,theta)+1.58545010*Z(37,rho,theta);
case 38
result = 19.24848143*Z(6,rho,theta)+16.41468913*Z(12,rho,theta)+9.76776798*Z(24,rho,theta) ...
+1.47438007*Z(28,rho,theta)+3.83118509*Z(38,rho,theta);
case 39
result = 0.46604820*Z(5,rho,theta)+0.84124290*Z(13,rho,theta)+1.00986774*Z(23,rho,theta) ...
-0.42520747*Z(27,rho,theta)+1.30579570*Z(39,rho,theta);
case 40
result = 28.18104531*Z(1,rho,theta)+38.52219208*Z(4,rho,theta)+30.18363661*Z(11,rho,theta) ...
+36.44278147*Z(14,rho,theta)+15.52577202*Z(22,rho,theta) ...
+19.21524879*Z(26,rho,theta)+4.44731721*Z(37,rho,theta)+6.00189814*Z(40,rho,theta);
case 41
result = 9.12899823*Z(15,rho,theta)+6.15821579*Z(25,rho,theta)+2.96653218*Z(41,rho,theta);
case 42
result = 85.33469748*Z(6,rho,theta)+64.01249391*Z(12,rho,theta)+30.59874671*Z(24,rho,theta) ...
+34.09158819*Z(28,rho,theta)+7.75796322*Z(38,rho,theta) ...
+9.37150432*Z(42,rho,theta);
case 43
result = 14.30642479*Z(5,rho,theta)+11.17404702*Z(13,rho,theta)...
+5.68231935*Z(23,rho,theta)+18.15306055*Z(27,rho,theta)...
+1.54919583*Z(39,rho,theta)+5.90178984*Z(43,rho,theta);
case 44
result = 36.12567424*Z(1,rho,theta)+47.95305224*Z(4,rho,theta)+35.30691679*Z(11,rho,theta) ...
+56.72014548*Z(14,rho,theta)+16.36470429*Z(22,rho,theta) ...
+26.32636277*Z(26,rho,theta)+3.95466397*Z(37,rho,theta)+6.33853092*Z(40,rho,theta) ...
+12.38056785*Z(44,rho,theta);
case 45
result = 21.45429746*Z(15,rho,theta)+9.94633083*Z(25,rho,theta) ...
+2.34632890*Z(41,rho,theta)+10.39130049*Z(45,rho,theta);
end % switch statement
end % function ZSquare
function result = ZAnnulus(j, n, m, rho, theta, e)
% Orthonormal annulus polynomials
if (j < 1) || (j > 37)
ME = MException('VerifyData:InvalidData', ...
'ZernikeCalc: Annulus order j out of range.');
throwAsCaller(ME);
end % if statement
switch j
case 1
% n=0, m=0
result = ones(size(rho));
case {2, 3}
% n=1, m=1
result = rho./sqrt(1+e^2);
case 4
% n=2, m=0
result = (2.*rho.^2-1-e^2)./(1-e^2);
case {5, 6}
% n=2, m=2
result = rho.^2./sqrt(1+e^2+e^4);
case {7, 8}
% n=3, m=1
result = (3*(1+e^2).*rho.^3-2*(1+e^2+e^4).*rho)./...
((1-e^2)*sqrt((1+e^2)*(1+4*e^2+e^4)));
case {9, 10}
% n=3, m=3
result = rho.^3./sqrt(1+e^2+e^4+e^6);
case 11
% n=4, m=0
result = (6.*rho.^4-6*(1+e^2).*rho.^2+1+4*e^2+e^4)/(1-e^2)^2;
case {12, 13}
% n=4, m=2
result = (4.*rho.^4-3*((1-e^8)/(1-e^6)).*rho.^2)./...
sqrt((1-e^2)^(-1)*(16*(1-e^10)-15*(1-e^8)^2/(1-e^6)));
case {14, 15}
% n=4, m=4
result = rho.^4./sqrt(1+e^2+e^4+e^6+e^8);
case {16, 17}
% n=5, m=1
result = (10*(1+4*e^2+e^4).*rho.^5-12*(1+4*e^2+4*e^4+e^6).*rho.^3 ...
+ 3*(1+4*e^2+10*e^4+4*e^6+e^8).*rho)./...
((1-e^2)^2*sqrt((1+4*e^2+e^4)*(1+9*e^2+9*e^4+e^6)));
case {18, 19}
% n=5, m=3
result = (5.*rho.^5-4*((1-e^10)/(1-e^8).*rho.^3))./...
sqrt(1/(1-e^2)*(25*(1-e^12)-24*(1-e^10)^2/(1-e^8)));
case {20, 21}
% n=5, m=5
result = rho.^5./sqrt(1+e^2+e^4+e^6+e^8+e^10);
case 22
% n=6, m=0
result = (20.*rho.^6-30*(1+e^2).*rho.^4+12*(1+3*e^2+e^4).*rho.^2- ...
(1+9*e^2+9*e^4+e^6))/(1-e^2)^3;
case {23, 24}
% n=6, m=2
result = (15*(1+4*e^2+10*e^4+4*e^6+e^8).*rho.^6-...
20*(1+4*e^2+10*e^4+10*e^6+4*e^8+e^10).*rho.^4+...
6*(1+4*e^2+10*e^4+20*e^6+10*e^8+4*e^10+e^12).*rho.^2)./...
((1-e^2)^2*sqrt((1+4*e^2+10*e^4+4*e^6+e^8)*...
(1+9*e^2+45*e^4+65*e^6+45*e^8+9*e^10+e^12)));
case {25, 26}
% n=6, m=4
result = (6.*rho.^6-5*((1-e^12)/(1-e^10)).*rho.^4)./...
sqrt((1/(1-e^2))*(36*(1-e^14)-35*(1-e^12)^2/(1-e^10)));
case {27, 28}
% n=6, m=6
result = rho.^6/sqrt(1+e^2+e^4+e^6+e^8+e^10+e^12);
case {29, 30}
% n=7, m=1
A710 = (1-e^2)^3 * sqrt(1 + 9*e^2 + 9*e^4 + e^6) * ...
sqrt(1 + 16*e^2 + 36*e^4 + 16*e^6 + e^8);
a71 = 35*(1 + 9*e^2 + 9*e^4 + e^6)/A710;
b71 = -60 * (1 + 9*e^2 + 15*e^4 + 9*e^6 + e^8)/A710;
c71 = 30*(1 + 9*e^2 + 25*e^4 + 25*e^6 + 9*e^8 + e^10)/A710;
d71 = -4*(1 + 9*e^2 + 45*e^4 + 65*e^6 + 45*e^8 + 9*e^10 + e^12)/A710;
result = a71.*rho.^7+b71.*rho.^5+c71.*rho.^3+d71.*rho;
case {31, 32}
% n=7, m=3
A730 = (1-e^2)^2 * sqrt(1 + 4*e^2 + 10*e^4 + 20*e^6 + 10*e^8 + 4*e^10 + e^12) * ...
sqrt(1 + 9*e^2 + 45*e^4 + 165*e^6 + 270*e^8 + 270*e^10 + ...
165*e^12 + 45*e^14 + 9*e^16 + e^18);
a73 = 21*(1 + 4*e^2 + 10*e^4 + 20*e^6 + 10*e^8 + 4*e^10 + e^12)/A730;
b73 = -30*(1 + 4*e^2 + 10*e^4 + 20*e^6 + 20*e^8 + 10*e^10 + 4*e^12 + e^14)/A730;
c73 = 10*(1 + 4*e^2 + 10*e^4 + 20*e^6 + 35*e^8 + 20*e^10 + 10*e^12 + 4*e^14 + e^16)/A730;
result = a73.*rho.^7+b73.*rho.^5+c73.*rho.^3;
case {33, 34}
% n=7, m=5
result = (7.*rho.^7-6*((1-e^14)/(1-e^12)).*rho.^5)/...
sqrt((1/(1-e^2))*(49*(1-e^16)-48*(1-e^14)^2/(1-e^12)));
case {35, 36}
% n=7, m=7
result = rho.^7/sqrt(1+e^2+e^4+e^6+e^8+e^10+e^12+e^14);
case 37
% n=8, m=0
result = (70*rho.^8-140*(1+e^2)*rho.^6+30*(3+8*e^2+3*e^4)*rho.^4 ...
-20*(1+6*e^2+6*e^4+e^6)*rho.^2+(1+16*e^2+36*e^4+16*e^6+e^8))./(1-e^2)^4;
end % switch j
% Calculate normalization factor.
if m == 0
factor = sqrt(n+1);
else
factor = sqrt(2*(n+1));
end % if statement
result = result * factor;
% Calculate sin(), cos() factor
if (m ~= 0) && (mod(j,2) == 0)
% j is even
result = result .* cos(m.*theta);
end % if statement
if (m ~= 0) && (mod(j,2) ~= 0)
result = result .* sin(m.*theta);
end % if statement
end % function ZAnnulus
function mask=makeDefaultMask(maskType, defaultRows, defaultCols, ShapeParameter)
% This function makes a default mask. Since it is a default mask, the
% mask matrix is to be square.
mask = zeros(defaultRows, defaultCols);
% the circle into which the mask shape is to fit
cr = (defaultRows+1)/2;
cc = (defaultCols+1)/2;
defaultRadiusInPixels = (min(defaultRows, defaultCols) - 1)/2;
switch maskType
case 'HEXAGON'
% make a Hexagon mask
for r=1:defaultRows
for c=1:defaultCols
x = (c-cc);
y = (r-cr);
rho = sqrt(x^2+y^2);
eTheta = atan2(y,x);
if (eTheta >= 0) && (eTheta <= (60*pi/180))
beta = 120*pi/180 - eTheta;
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if eTheta >= (120*pi/180)
beta = 120*pi/180 - (pi - eTheta);
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if eTheta <= -(120*pi/180)
beta = 120*pi/180 - (pi - abs(eTheta));
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if (eTheta <= 0) && (eTheta >= -(60*pi/180))
beta = 120*pi/180 - abs(eTheta);
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if (eTheta >= (60*pi/180)) && (eTheta <= (120*pi/180))
if abs(cr-r) <= ((sqrt(3)/2)*defaultRadiusInPixels)
mask(r,c) = 1;
end % if statement
end % if statement
if (eTheta <= (-60*pi/180)) && (eTheta >= (-120*pi/180))
if abs(cr-r) <= ((sqrt(3)/2)*defaultRadiusInPixels)
mask(r,c) = 1;
end % if statement
end % if statement
end % for c statement
end % for r statement
case 'HEXAGONR30'
% make a Hexagon mask rotated 30 degrees
for r=1:defaultRows
for c=1:defaultCols
x = (c-cc);
y = (r-cr);
rho = sqrt(x^2+y^2);
eTheta = atan2(y,x);
if (eTheta >= (30*pi/180)) && (eTheta <= (90*pi/180))
beta = 120*pi/180 - (eTheta - 30*pi/180);
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if (eTheta >= (90*pi/180)) && (eTheta <= (150*pi/180))
beta = 120*pi/180 - (pi - eTheta - 30*pi/180);
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if (eTheta <= -(90*pi/180)) && (eTheta >= -(150*pi/180))
beta = 120*pi/180 - (pi - abs(eTheta) - 30*pi/180);
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if (eTheta <= -30*pi/180) && (eTheta >= -(90*pi/180))
beta = 120*pi/180 - (abs(eTheta) - 30*pi/180);
R = defaultRadiusInPixels * sind(60) / sin(beta);
if rho < R
mask(r,c) = 1;
end % if R statement
end % if statement
if (eTheta >= (-30*pi/180)) && (eTheta <= (30*pi/180))
if abs(cc-c) <= ((sqrt(3)/2)*defaultRadiusInPixels)
mask(r,c) = 1;
end % if statement
end % if statement
if (eTheta >= (150*pi/180)) || (eTheta <= (-150*pi/180))
if abs(cc-c) <= ((sqrt(3)/2)*defaultRadiusInPixels)
mask(r,c) = 1;
end % if statement
end % if statement
end % for c statement
end % for r statement
case 'ELLIPSE'
% an ellipse mask is needed
a = defaultRadiusInPixels;
b = ShapeParameter * defaultRadiusInPixels;
for r=1:defaultRows
for c=1:defaultCols
rho = sqrt((cr-r)^2 + (cc-c)^2);
eTheta = atan2((cr-r),(cc-c));
er = a*b/sqrt((b*cos(eTheta))^2+(a*sin(eTheta))^2);
if rho <= er
mask(r,c) = 1;
end % if statement
end % for c statement
end % for r stateemnt
case 'RECTANGLE'
% a rectangular mask is needed
halfEdgec = ShapeParameter * defaultRadiusInPixels;
halfEdger = sqrt(1 - ShapeParameter^2) * defaultRadiusInPixels;
for r=1:defaultRows
for c=1:defaultCols
if (r > abs(cr-halfEdger)) && (r < (cr+halfEdger)) && ...
(c > abs(cr-halfEdgec)) && (c < (cr+halfEdgec))
mask(r,c) = 1;
end % if statement
end % for c statement
end % for r stateemnt
case 'SQUARE'
% a square mask is needed
halfEdge = (1/sqrt(2)) * defaultRadiusInPixels;
for r=1:defaultRows
for c=1:defaultCols
if (r > abs(cr-halfEdge)) && (r < (cr+halfEdge)) && ...
(c > abs(cr-halfEdge)) && (c < (cr+halfEdge))
mask(r,c) = 1;
end % if statement
end % for c statement
end % for r stateemnt
case 'ANNULUS'
% an annulus mask is needed
obscuration = ShapeParameter * defaultRadiusInPixels;
for r=1:defaultRows
for c=1:defaultCols
rho = sqrt((cr-r)^2 + (cc-c)^2);
if (rho <= defaultRadiusInPixels) && (rho > obscuration)
mask(r,c) = 1;
end % if statement
end % for c statement
end % for r stateemnt
otherwise
% a circle mask is needed
for r=1:defaultRows
for c=1:defaultCols
rho = sqrt((cr-r)^2 + (cc-c)^2);
if rho <= defaultRadiusInPixels
mask(r,c) = 1;
end % if statement
end % for c statement
end % for r stateemnt
end % switch statement
end % function makeDefaultMask
function validateZ()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate Circle
% fprintf(1,'\n\n CIRCLE VALIDATION \n');
%
% ZernikeType = 'MAHAJAN';
%
% for j1=1:37
%
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
% [n1,m1,sc1] = convertj(j1, ZernikeType);
% [n2,m2,sc2] = convertj(j2, ZernikeType);
%
% fun1 = @(rho,theta) calculateZernike(n1,m1,sc1,rho,theta,ZernikeType).*calculateZernike(n1,m1,sc1,rho,theta,ZernikeType).*rho;
% Q1 = quad2d(fun1, 0,1, 0,2*pi);
% diff1 = pi - Q1;
%
% fun2 = @(rho,theta) calculateZernike(n1,m1,sc1,rho,theta,ZernikeType).*calculateZernike(n2,m2,sc2,rho,theta,ZernikeType).*rho;
% Q2 = quad2d(fun2, 0,1, 0,2*pi);
% diff2 = 0.0 - Q2;
%
% fprintf(1,'j1=%d, n1=%d, m1=%d, result = %+20.18f Difference = %+20.18f \n',j1,n1,m1,Q1,diff1);
% fprintf(1,'j2=%d, n2=%d, m2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,n2,m2,Q2,diff2);
%
%
% end % for statement
%
% end Circle validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate Circle Fringe
% fprintf(1,'\n\n CIRCLE (FRINGE) VALIDATION \n');
%
% ZernikeType = 'FRINGE';
%
% for j1=1:37
%
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
% [n1,m1,sc1] = convertj(j1, ZernikeType);
% [n2,m2,sc2] = convertj(j2, ZernikeType);
%
% fun1 = @(rho,theta) calculateZernike(n1,m1,sc1,rho,theta,ZernikeType).*calculateZernike(n1,m1,sc1,rho,theta,ZernikeType).*rho;
% Q1 = quad2d(fun1, 0,1, 0,2*pi);
% % The normalization is different.
% em = 1;
% if m1 == 0
% em = 2;
% end % if statement
% diff1 = (em*pi)/(2*n1+2) - Q1;
%
% fun2 = @(rho,theta) calculateZernike(n1,m1,sc1,rho,theta,ZernikeType).*calculateZernike(n2,m2,sc2,rho,theta,ZernikeType).*rho;
% Q2 = quad2d(fun2, 0,1, 0,2*pi);
% diff2 = 0.0 - Q2;
%
% fprintf(1,'j1=%d, n1=%d, m1=%d, result = %+20.18f Difference = %+20.18f \n',j1,n1,m1,Q1,diff1);
% fprintf(1,'j2=%d, n2=%d, m2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,n2,m2,Q2,diff2);
%
%
% end % for statement
%
% end Circle Fringe validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate ZHexagon
% fprintf(1,'\n\n HEXAGON VALIDATION \n');
% for j1=1:45
%
% % Need another j value to compare orthogonality condition
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
%
% ymina = @(x) -(sqrt(3)*x + sqrt(3));
% ymaxa = @(x) sqrt(3)*x + sqrt(3);
%
% yminc = @(x) -(-sqrt(3)*x + sqrt(3));
% ymaxc = @(x) -sqrt(3)*x + sqrt(3);
%
% fun1 = @(x,y) ZHexagon(j1,sqrt(x.*x+y.*y),atan2(y,x)).*ZHexagon(j1,sqrt(x.*x+y.*y),atan2(y,x));
% Q1a = quad2d(fun1, -1,-0.5, ymina,ymaxa);
% Q1b = quad2d(fun1, -0.5,0.5, -sqrt(3)/2,sqrt(3)/2);
% Q1c = quad2d(fun1, 0.5,1, yminc,ymaxc);
% diff1 = 3*sqrt(3)/2 - (Q1a+Q1b+Q1c);
%
% fun2 = @(x,y) ZHexagon(j1,sqrt(x.*x+y.*y),atan2(y,x)).*ZHexagon(j2,sqrt(x.*x+y.*y),atan2(y,x));
% Q2a = quad2d(fun2, -1,-0.5, ymina,ymaxa);
% Q2b = quad2d(fun2, -0.5,0.5, -sqrt(3)/2,sqrt(3)/2);
% Q2c = quad2d(fun2, 0.5,1, yminc,ymaxc);
% diff2 = 0.0 - (Q2a+Q2b+Q2c);
%
% fprintf(1,'j1=%d, result = %+20.18f Difference = %+20.18f \n',j1,(Q1a+Q1b),diff1);
% fprintf(1,'j2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,(Q2a+Q2b),diff2);
%
%
% end % for statement
%
% end ZHexagon validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate ZHexagonR30
% fprintf(1,'\n\n HEXAGON ROTATED 30 DEGREES VALIDATION \n');
% for j1=1:28
%
% % Need another j value to compare orthogonality condition
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
%
% ymina = @(x) -(x/sqrt(3) + 1);
% ymaxa = @(x) x/sqrt(3) + 1;
%
% yminb = @(x) -(-x/sqrt(3) + 1);
% ymaxb = @(x) -x/sqrt(3) + 1;
%
% fun1 = @(x,y) ZHexagonR30(j1,sqrt(x.*x+y.*y),atan2(y,x)).*ZHexagonR30(j1,sqrt(x.*x+y.*y),atan2(y,x));
% Q1a = quad2d(fun1, -sqrt(3)/2,0, ymina,ymaxa);
% Q1b = quad2d(fun1, 0,sqrt(3)/2, yminb,ymaxb);
% diff1 = 3*sqrt(3)/2 - (Q1a+Q1b);
%
% fun2 = @(x,y) ZHexagonR30(j1,sqrt(x.*x+y.*y),atan2(y,x)).*ZHexagonR30(j2,sqrt(x.*x+y.*y),atan2(y,x));
% Q2a = quad2d(fun2, -sqrt(3)/2,0, ymina,ymaxa);
% Q2b = quad2d(fun2, 0,sqrt(3)/2, yminb,ymaxb);
% diff2 = 0.0 - (Q2a+Q2b);
%
% fprintf(1,'j1=%d, result = %+20.18f Difference = %+20.18f \n',j1,(Q1a+Q1b),diff1);
% fprintf(1,'j2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,(Q2a+Q2b),diff2);
%
%
% end % for statement
%
% end ZHexagonR30 validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate ZEllipse
%
% b = 0.8; % semi-minor axis
% %b = 1.0; % semi-minor axis
%
% fprintf(1,'\n\n ELLIPSE VALIDATION (b=%f)\n',b);
% for j1=1:15
%
% % Need another j value to compare orthogonality condition
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
%
% ymax = @(x) b*sqrt(1-x.*x);
% ymin = @(x) -b*sqrt(1-x.*x);
%
% fun1 = @(x,y) ZEllipse(j1,sqrt(x.*x+y.*y),atan2(y,x),b).*ZEllipse(j1,sqrt(x.*x+y.*y),atan2(y,x),b);
% Q1 = quad2d(fun1, -1,1, ymin,ymax);
% diff1 = pi*b - Q1;
%
% fun2 = @(x,y) ZEllipse(j1,sqrt(x.*x+y.*y),atan2(y,x),b).*ZEllipse(j2,sqrt(x.*x+y.*y),atan2(y,x),b);
% Q2 = quad2d(fun2, -1,1, ymin,ymax);
% diff2 = 0.0 - Q2;
%
% fprintf(1,'j1=%d, result = %+20.18f Difference = %+20.18f \n',j1,Q1,diff1);
% fprintf(1,'j2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,Q2,diff2);
%
%
% end % for statement
%
% end ZEllipse validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate ZRectangle
% a = 0.37; % half length of rectangle
% %a = 1/sqrt(2); % for a square
%
% fprintf(1,'\n\n RECTANGLE VALIDATION (a=%f)\n',a);
% for j1=1:15
%
% % Need another j value to compare orthogonality condition
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
% fun1 = @(x,y) ZRectangle(j1,sqrt(x.*x+y.*y),atan2(y,x),a).*ZRectangle(j1,sqrt(x.*x+y.*y),atan2(y,x),a);
% Q1 = quad2d(fun1, -a,a, -sqrt(1-a*a),sqrt(1-a*a));
% diff1 = (2*a * 2*sqrt(1-a*a)) - Q1;
%
% fun2 = @(x,y) ZRectangle(j1,sqrt(x.*x+y.*y),atan2(y,x),a).*ZRectangle(j2,sqrt(x.*x+y.*y),atan2(y,x),a);
% Q2 = quad2d(fun2, -a,a, -sqrt(1-a*a),sqrt(1-a*a));
% diff2 = 0.0 - Q2;
%
% fprintf(1,'j1=%d, result = %+20.18f Difference = %+20.18f \n',j1,Q1,diff1);
% fprintf(1,'j2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,Q2,diff2);
%
%
% end % for statement
%
% end ZRectangle validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate ZSquare
% fprintf(1,'\n\n SQUARE VALIDATION \n');
% for j1=1:45
%
% % Need another j value to compare orthogonality condition
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
% fun1 = @(x,y) ZSquare(j1,sqrt(x.*x+y.*y),atan2(y,x)).*ZSquare(j1,sqrt(x.*x+y.*y),atan2(y,x));
% Q1 = quad2d(fun1, -1/sqrt(2),1/sqrt(2), -1/sqrt(2),1/sqrt(2));
% diff1 = (sqrt(2)*sqrt(2)) - Q1;
%
% fun2 = @(x,y) ZSquare(j1,sqrt(x.*x+y.*y),atan2(y,x)).*ZSquare(j2,sqrt(x.*x+y.*y),atan2(y,x));
% Q2 = quad2d(fun2, -1/sqrt(2),1/sqrt(2), -1/sqrt(2),1/sqrt(2));
% diff2 = 0.0 - Q2;
%
% fprintf(1,'j1=%d, result = %+20.18f Difference = %+20.18f \n',j1,Q1,diff1);
% fprintf(1,'j2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,Q2,diff2);
%
%
% end % for statement
%
% end ZSquare validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Validate ZAnnulus
% fprintf(1,'\n\n ANNULUS VALIDATION \n');
% for j1=1:37
%
% j2 = j1 - 1;
% if j1 == 1
% j2 = 2;
% end % if statement
%
%
% [n1,m1,sc1] = convertj(j1, 'MAHAJAN');
% [n2,m2,sc2] = convertj(j2, 'MAHAJAN');
%
% ir = 0.3; % annulus inner radius
% fun1 = @(rho,theta) ZAnnulus(j1,n1,m1,rho,theta,ir).*ZAnnulus(j1,n1,m1,rho,theta,ir).*rho;
% Q1 = quad2d(fun1, ir,1, 0,2*pi);
% diff1 = (pi - pi*ir^2) - Q1;
%
% fun2 = @(rho,theta) ZAnnulus(j1,n1,m1,rho,theta,ir).*ZAnnulus(j2,n2,m2,rho,theta,ir).*rho;
% Q2 = quad2d(fun2, ir,1, 0,2*pi);
% diff2 = 0.0 - Q2;
%
% fprintf(1,'j1=%d, n1=%d, m1=%d, result = %+20.18f Difference = %+20.18f \n',j1,n1,m1,Q1,diff1);
% fprintf(1,'j2=%d, n2=%d, m2=%d, result = %+20.18f Difference = %+20.18f \n\n',j2,n2,m2,Q2,diff2);
%
% end % for statement
%
%
% end annulus validation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end % validateZ
543
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
