利用前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 interp

zernike函數

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.') endif length(n)~=length(m)error('zernfun:NMlength','N and M must be the same length.') endn=n(:); m=m(:); if any(mod(n-m,2))error('zernfun:NMmultiplesof2', ...'All N and M must differ by multiples of 2 (including 0).') endif any(m>n)error('zernfun:MlessthanN', ...'Each M must be less than or equal to its corresponding N.') endif any(r>1|r<0)error('zernfun:Rlessthan1','All R must be between 0 and 1.') endif (~any(size(r)==1))||(~any(size(theta)==1))error('zernfun:RTHvector','R and THETA must be vectors.') endr=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 ~isnormerror('zernfun:normalization','Unrecognized normalization flag.')end elseisnorm=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)==0rpowern=arrayfun(@(p)r.^p,rpowers(2:end),'UniformOutput',false);rpowern=cat(2,rpowern{:});rpowern=[ones(length_r,1) rpowern]; elserpowern=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:1p=(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);endif isnormy(:,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)');elsez(:,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)');elsez(:,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 statementif nargin <= 3 || isempty(ZernikeDef)ZernikeDef = 'FRINGE'; end % if statement % Convert to upper caseZernikeDef = upper(ZernikeDef);if nargin <= 4 || isempty(ShapeParameter)ShapeParameter = 0; end % if statementif 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 statementmask = makeDefaultMask(ZernikeDef, defaultRows, defaultCols, ShapeParameter);end % if no mask statementsm = 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 nccase 1% only the radiusInPixels has been specifiedradiusInPixels = 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-centerColcm = (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 specifiedcenterRow = unitCircle(1,1);centerCol = unitCircle(1,2);radiusInPixels = unitCircle(1,3);otherwise% error.end % switch statementif 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 MMESorCn(k) = ZernikeList(1, k);m(k) = abs(ZernikeList(2, k));sc(k) = ' ';switch MMESorCcase '-m=sin'if ZernikeList(2, k) < 0sc(k) = 's';end % if statementif ZernikeList(2, k) > 0sc(k) = 'c'; end % if statementcase '-m=cos'if ZernikeList(2, k) < 0sc(k) = 'c';end % if statementif ZernikeList(2, k) > 0sc(k) = 's'; end % if statement end % switch statementend % for k statementend % 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-centerColcm = 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 statementfor 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 figuressizeZD = 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 differencetoPlot(:,:,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 >= 1if nargout < 2% This means nargout == 1.if interferogramInput% return calculated coefficients onlyoutput1 = 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) == 1output1 = sum(zm,3);output2 = zm;return; end % if statementoutput1 = zm;output2 = zcoeff;end % ZernikeCalcfunction 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 inputif 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 = []; endsizedata = 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; endif isfield(datastruct,'mask'), mask(:,:,i) = datastruct(i).mask; endif isfield(datastruct,'psf'), data(:,:,i) = datastruct(i).psf; endend 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 matrixsqrtval = sqrt(length(data(:,1)));if mod(length(data(:,1)),10000) == 0 %if data is integer lengths of 100x100disp('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 == 1data = reshape(data,100,100,size(data,2)); %3d dataelsedata = reshape(data,100,100,num100x100,size(data,2)); %4d dataendelseif mod(length(data(:,1)), sqrtval) == 0 % if data squaredisp('Data appears to be square in column format, reshaping into 2d.');data = reshape(data,sqrtval,sqrtval,size(data,2));elsesqrtval = sqrt(length(data(:,2)));if mod(length(data(:,1)), sqrtval) == 0 % if data squaredisp('Data appears to be square in row format, reshaping into 2d.');data = reshape(data',sqrtval,sqrtval,size(data',2));endend 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; elsenumplots = 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*numplotcolsnumfigs = ceil(numplots/(numplotrows*numplotcols)); %number of figures elsenumfigs = 1; endif ~iscell(titles) %array of text givenfor i=1:size(titles,1)t{i} = titles(i,:);endtitles = 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:numfigsscrsz = 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);endset(gcf,'Name',inputname(1));for i=1:numplotrowsfor j=1:numplotcolsplotnum = (i-1)*numplotcols+j+(k-1)*numplotrows*numplotcols; %data to plotsubplotnum = (i-1)*numplotcols+j; %subplot location to plotif plotnum<numplots+1,plotdata = data(:,:,plotnum)*1; % *1 converts logical to doublesif size(mask,3)>1, maskdata = mask(:,:,plotnum); else maskdata = mask(:,:,1); endsubplot(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<1title(int2str(plotnum));elseif plotnum>numtitles %do nothing if partial title list givenelse title(titles{plotnum}); end put_xlabel(plotdata,maskdata,labeloption,labelprecision); %subfunction included belowendif nargin>3, if length(titles)>=plotnum, title(titles(plotnum)); end; endendend end end % function figs% X Label for each plot function put_xlabel(plotdata,maskdata,labeloption,labelprecision) if nargin<2maskdata = ~isnan(maskdata) & data~=0;datavals = plotdata(maskdata); elseif ~islogical(maskdata), maskdata = logical(maskdata); enddatavals = plotdata(maskdata); end if ~isempty(datavals) %only analyze statistics on a real data setstdval = 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))); elsestdval = 0; pk2val = 0; maxval = 0; minval = 0; meanval = 0; p95 = 0; end if labeloption==1xlabel(['RMS = ' num2str(stdval,labelprecision)]); elseif labeloption==2xlabel({['AVG = ' num2str(meanval,labelprecision)];['RMS = ' num2str(stdval,labelprecision)];['PV = ' num2str(pk2val,labelprecision)]}); elseif labeloption==3xlabel(['RMS = ' num2str(stdval*1e9,labelprecision) ' nm']); endend % function put_xlabelfunction 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 RnmRnm = zeros(size(rho));for s=0:(n-m)/2numerator = (-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 sccase ' '% 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 ZernikeDefcase {'MAHAJAN', 'NOLL', ...'HEXAGON', 'HEXAGONR30', 'RECTANGLE', 'SQUARE', 'ANNULUS'}result = theFactor * result; end % switchend % 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 ztypecase '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 ztypeend % function convertjfunction [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 jcase 1n = 0; m = 0; sc = ' ';case 2n = 1; m = 1; sc = 'c'; case 3n = 1; m = 1; sc = 's';case 4n = 2; m = 0; sc = ' ';case 5n = 2; m = 2; sc = 'c';case 6n = 2; m = 2; sc = 's';case 7n = 3; m = 1; sc = 'c';case 8n = 3; m = 1; sc = 's';case 9n = 4; m = 0; sc = ' ';case 10n = 3; m = 3; sc = 'c';case 11n = 3; m = 3; sc = 's';case 12n = 4; m = 2; sc = 'c';case 13n = 4; m = 2; sc = 's';case 14n = 5; m = 1; sc = 'c';case 15n = 5; m = 1; sc = 's';case 16n = 6; m = 0; sc = ' ';case 17n = 4; m = 4; sc = 'c';case 18n = 4; m = 4; sc = 's';case 19n = 5; m = 3; sc = 'c';case 20n = 5; m = 3; sc = 's';case 21n = 6; m = 2; sc = 'c';case 22n = 6; m = 2; sc = 's';case 23n = 7; m = 1; sc = 'c';case 24n = 7; m = 1; sc = 's';case 25n = 8; m = 0; sc = ' ';case 26n = 5; m = 5; sc = 'c';case 27n = 5; m = 5; sc = 's';case 28n = 6; m = 4; sc = 'c';case 29n = 6; m = 4; sc = 's';case 30n = 7; m = 3; sc = 'c';case 31n = 7; m = 3; sc = 's';case 32n = 8; m = 2; sc = 'c';case 33n = 8; m = 2; sc = 's';case 34n = 9; m = 1; sc = 'c';case 35n = 9; m = 1; sc = 's';case 36n = 10; m = 0; sc = ' ';case 37n = 12; m = 0; sc = ' ';end % switch jend % function fringefunction [n, m, sc] = bw(j)% sc = 's' = sin(), sc = 'c' = cos(), sc = ' ' for neither.% sc = 1 = sin(), sc = 2 = cos().% calculate the n valuen1 = (-1 + sqrt(1 + 8 * j)) / 2;n = floor(n1);if n1 == nn = n - 1; end % if statement% calculate the m valuek = (n+1)*(n+2)/2;d = k - j;m1 = n - 2*d;m = abs(m1);% calculate the sc valuesc = ' ';if m1 > 0 sc = 's';end % if statementif m1 < 0 sc = 'c';end % if statementend % function bwfunction [n, m, sc] = mahajan(j)% sc = 's' = sin(), sc = 'c' = cos(), sc = ' ' for neither.% sc = 1 = sin(), sc = 2 = cos().% calculate the n valuen1 = (-1 + sqrt(1 + 8 * j)) / 2;n = floor(n1);if n1 == nn = n - 1; end % if statement% calculate the m valuek = (n+1)*(n+2)/2;m = n - 2 * floor((k - j)/2);% calculate the sc valuesc = ' ';if (m ~= 0) && (mod(j,2) ~= 0)% m ~= 0 and j odd sc = 's';end % if statementif (m ~= 0) && (mod(j,2) == 0)% m ~= 0 and j even sc = 'c';end % if statementend % 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 < 1ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: centerRow must be positive.');throwAsCaller(ME); end % if statement if centerCol < 1ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: centerCol must be positive.');throwAsCaller(ME); end % if statementif radiusInPixels < 1ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: radiusInPixels must be positive.');throwAsCaller(ME); end % if statementhlda = 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 ZernikeDefcase {'FRINGE', 'ISO', 'WYANT' 'MAHAJAN', 'NOLL', 'B&W', 'STANDARD', ...'HEXAGON', 'HEXAGONR30', 'ELLIPSE', 'RECTANGLE', ...'SQUARE', 'ANNULUS'}% These are the valid values. otherwise% ZernikeType is not validME = 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 integerME = MException('VerifyData:InvalidData', ...'ZernikeCalc: ZernikeList must contain only integers.');throwAsCaller(ME); end % if statementif rows == 1% j orderingif sum(abs(ZernikeList(1, :)) ~= ZernikeList(1, :)) ~= 0ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: ZernikeList j values must be positive or 0.');throwAsCaller(ME); end % if statement switch ZernikeDefcase 'FRINGE'if sum(ZernikeList(1, :) < 1) ~= 0ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: FRINGE order value j must be greater than 0.');throwAsCaller(ME); end % if statement if sum(ZernikeList(1, :) > 37) ~= 0ME = 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) ~= 0ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: ISO order value j out of range.');throwAsCaller(ME); end % if statement case 'WYANT' if sum(ZernikeList(1, :) > 36) ~= 0ME = 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) ~= 0ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Zernike order value j can not be zero.');throwAsCaller(ME); end % if statement end % switch statement else% (n, m) specificationif sum(abs(ZernikeList(1, :)) ~= ZernikeList(1, :)) ~= 0ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: ZernikeList n values must be positive or 0.');throwAsCaller(ME); end % if statementif sum(abs(ZernikeList(2, :)) > ZernikeList(1, :)) ~= 0ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: abs(m) values must not exceed n values.');throwAsCaller(ME); end % if statement switch ZernikeDefcase 'FRINGE'if rows > 1ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: FRINGE can only be specified with j ordering.');throwAsCaller(ME); end % if statementcase 'ISO'if rows > 1ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: ISO can only be specified with j ordering.');throwAsCaller(ME); end % if statementcase 'WYANT'if rows > 1ME = 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 > 1ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: ZernikeDef value can only be specified with j ordering.');throwAsCaller(ME); end % if statement otherwiseif (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 % switchswitch MMESorCcase {'-m=sin', '-m=cos'}otherwiseME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Sign convention value for (n, m) not valid.');throwAsCaller(ME); end % switchend % (n,m) order specified % Validate ZernikeShape parameterswitch ZernikeDef case {'ELLIPSE', 'RECTANGLE','ANNULUS'}% Check that the ShapeParameter is validif (ShapeParameter <= 0) || (ShapeParameter >= 1)ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: The ShapeParameter is out of range: (0,1).');throwAsCaller(ME);end % if statementend % switch ZernikeDefend % 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 Zfunction result = ZHexagon(j, rho, theta)% Orthonormal hexagonal polynomialsif (j < 1) || (j > 45)ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Hexagon order j out of range.');throwAsCaller(ME); end % if statementswitch jcase 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 statementend % 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 statementswitch jcase 1result = Z(1,rho,theta); case 2result = sqrt(6/5)*Z(2, rho, theta); case 3result = sqrt(6/5)*Z(3, rho, theta); case 4result = sqrt(5/43)*Z(1, rho, theta)+2*sqrt(15/43)*Z(4, rho, theta); case 5result = sqrt(10/7)*Z(5,rho,theta); case 6result = sqrt(10/7)*Z(6,rho,theta); case 7result = 16*sqrt(14/11055)*Z(3,rho,theta)+10*sqrt(35/2211)*Z(7,rho,theta); case 8result = 16*sqrt(14/11055)*Z(2,rho,theta)+10*sqrt(35/2211)*Z(8,rho,theta); case 9result = 2*sqrt(35/103)*Z(9,rho,theta); case 10result = (2*sqrt(5)/3)*Z(10,rho,theta); case 11result = (521/sqrt(1072205))*Z(1,rho,theta)+88*sqrt(15/214441)*Z(4,rho,theta) ...+ 14*sqrt(43/4987)*Z(11,rho,theta); case 12result = 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 14result = 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 15result = -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 16result = 30857*sqrt(2/3268147641)*Z(2,rho,theta)+49168/sqrt(3268147641)*Z(8,rho,theta) ...+ 42*sqrt(1474/1478131)*Z(16,rho,theta); case 17result = 30857*sqrt(2/3268147641)*Z(3,rho,theta)+(49168/sqrt(3268147641))*Z(7,rho,theta) ...+ 42*sqrt(1474/1478131)*Z(17,rho,theta); case 18result = 6*sqrt(10/97)*Z(10,rho,theta)+14*sqrt(5/291)*Z(18,rho,theta); case 19result = 386*sqrt(770/295894589)*Z(9,rho,theta)+6*sqrt(118965/2872763)*Z(19,rho,theta); case 20result = 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 22result = 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 23result = 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 24result = 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 25result = -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 26result = 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 27result = 1.81984283*Z(27,rho,theta); case 28result = 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 polynomialsif (j < 1) || (j > 15)ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Elliptical order j out of range.');throwAsCaller(ME); end % if statementalpha = 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 jcase 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 10result = -((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 statementend % function ZEllipticalfunction result = ZRectangle(j, rho, theta, a)% Orthonormal rectangle polynomialsif (j < 1) || (j > 15)ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Rectangle order j out of range.');throwAsCaller(ME); end % if statementmu = 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 jcase 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 statementend % function ZRectanglefunction result = ZSquare(j, rho, theta)% Orthonormal square polynomialsif (j < 1) || (j > 45)ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Square order j out of range.');throwAsCaller(ME); end % if statementswitch jcase 1result = Z(1,rho,theta);case 2result = sqrt(3/2)*Z(2,rho,theta);case 3result = sqrt(3/2)*Z(3,rho,theta);case 4result = (sqrt(5/2)/2)*Z(1,rho,theta)+(sqrt(15/2)/2)*Z(4,rho,theta);case 5result = sqrt(3/2)*Z(5,rho,theta);case 6result = (sqrt(15)/2)*Z(6,rho,theta);case 7result = (3*sqrt(21/31)/2)*Z(3,rho,theta)+(5*sqrt(21/62)/2)*Z(7,rho,theta);case 8result = (3*sqrt(21/31)/2)*Z(2,rho,theta)+(5*sqrt(21/62)/2)*Z(8,rho,theta);case 9result = -(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 10result = (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 11result = (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 12result = (45*sqrt(3)/16)*Z(6,rho,theta)+(21*sqrt(5)/16)*Z(12,rho,theta);case 13result = (3*sqrt(7)/8)*Z(5,rho,theta)+(sqrt(105)/8)*Z(13,rho,theta);case 14result = 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 15result = (sqrt(105)/4)*Z(15,rho,theta);case 16result = ((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 17result = ((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 18result = ((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 19result = -((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 20result = ((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 21result = ((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 22result = ((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 23result = ((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 24result = ((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 25result = ((7*sqrt(165))/16)*Z(15,rho,theta)+((3*sqrt(231))/16)*Z(25,rho,theta); case 26result = (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 27result = ((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 28result = ((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 29result = 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 30result = 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 31result = -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 32result = 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 33result = 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 34result = 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 35result = -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 36result = 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 37result = 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 38result = 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 39result = 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 40result = 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 41result = 9.12899823*Z(15,rho,theta)+6.15821579*Z(25,rho,theta)+2.96653218*Z(41,rho,theta);case 42result = 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 43result = 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 44result = 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 45result = 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 statementend % function ZSquarefunction result = ZAnnulus(j, n, m, rho, theta, e)% Orthonormal annulus polynomialsif (j < 1) || (j > 37) ME = MException('VerifyData:InvalidData', ...'ZernikeCalc: Annulus order j out of range.');throwAsCaller(ME); end % if statementswitch jcase 1% n=0, m=0result = 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=1A710 = (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 == 0factor = sqrt(n+1);elsefactor = sqrt(2*(n+1)); end % if statementresult = result * factor;% Calculate sin(), cos() factorif (m ~= 0) && (mod(j,2) == 0)% j is even result = result .* cos(m.*theta);end % if statementif (m ~= 0) && (mod(j,2) ~= 0) result = result .* sin(m.*theta);end % if statementend % function ZAnnulusfunction 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 fitcr = (defaultRows+1)/2; cc = (defaultCols+1)/2;defaultRadiusInPixels = (min(defaultRows, defaultCols) - 1)/2;switch maskTypecase 'HEXAGON'% make a Hexagon maskfor r=1:defaultRowsfor c=1:defaultColsx = (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 statementend % 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 statementend % 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 statementend % 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 statementend % 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 statementend % if statementif (eTheta <= (-60*pi/180)) && (eTheta >= (-120*pi/180))if abs(cr-r) <= ((sqrt(3)/2)*defaultRadiusInPixels)mask(r,c) = 1; end % if statementend % if statementend % for c statementend % for r statementcase 'HEXAGONR30' % make a Hexagon mask rotated 30 degreesfor r=1:defaultRowsfor c=1:defaultColsx = (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 statementend % 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 statementend % 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 statementend % 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 statementend % 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 statementend % if statementif (eTheta >= (150*pi/180)) || (eTheta <= (-150*pi/180))if abs(cc-c) <= ((sqrt(3)/2)*defaultRadiusInPixels)mask(r,c) = 1; end % if statementend % if statementend % for c statementend % for r statement case 'ELLIPSE'% an ellipse mask is neededa = defaultRadiusInPixels;b = ShapeParameter * defaultRadiusInPixels;for r=1:defaultRowsfor c=1:defaultColsrho = 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 <= ermask(r,c) = 1;end % if statement end % for c statementend % for r stateemntcase 'RECTANGLE'% a rectangular mask is neededhalfEdgec = ShapeParameter * defaultRadiusInPixels;halfEdger = sqrt(1 - ShapeParameter^2) * defaultRadiusInPixels;for r=1:defaultRowsfor c=1:defaultColsif (r > abs(cr-halfEdger)) && (r < (cr+halfEdger)) && ...(c > abs(cr-halfEdgec)) && (c < (cr+halfEdgec)) mask(r,c) = 1;end % if statement end % for c statementend % for r stateemnt case 'SQUARE'% a square mask is neededhalfEdge = (1/sqrt(2)) * defaultRadiusInPixels;for r=1:defaultRowsfor c=1:defaultColsif (r > abs(cr-halfEdge)) && (r < (cr+halfEdge)) && ...(c > abs(cr-halfEdge)) && (c < (cr+halfEdge)) mask(r,c) = 1;end % if statement end % for c statementend % for r stateemntcase 'ANNULUS'% an annulus mask is neededobscuration = ShapeParameter * defaultRadiusInPixels;for r=1:defaultRowsfor c=1:defaultColsrho = sqrt((cr-r)^2 + (cc-c)^2);if (rho <= defaultRadiusInPixels) && (rho > obscuration)mask(r,c) = 1;end % if statement end % for c statementend % for r stateemntotherwise% a circle mask is needed for r=1:defaultRowsfor c=1:defaultColsrho = sqrt((cr-r)^2 + (cc-c)^2);if rho <= defaultRadiusInPixelsmask(r,c) = 1;end % if statement end % for c statementend % for r stateemntend % switch statementend % function makeDefaultMaskfunction 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

?

總結

以上是生活随笔為你收集整理的Zernike函数拟合曲面--MATLAB实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。