zernike拟合

来源：华佗养生网

Zernike_main.m

% ------------------------------------------------------------------------- % Copyright C 2012 Amir Tahmasbi % The University of Texas at Dallas % a.tahmasbi@utdallas.edu

% http://www.utdallas.edu/~a.tahmasbi/research.html %

% -------------------------------------------------------------------------

% Licence Agreement: You are free to use this code in your scientific

% research but you should cite the following papers: %

% [1] - A. Tahmasbi, F. Saki, S. B. Shokouhi,

% Classification of Benign and Malignant Masses Based on Zernike Moments, % J. Computers in Biology and Medicine, vol. 41, no. 8, pp. 726-735, 2011. %

% [2] - A. Tahmasbi, F. Saki, H. Aghapanah, S. B. Shokouhi,

% A Novel Breast Mass Diagnosis System based on Zernike Moments as Shape and Density Descriptors,

% in Proc. IEEE, 18th Iranian Conf. on Biomedical Engineering (ICBME'2011), % Tehran, Iran, 2011, pp. 100-104. %

% ------------------------------------------------------------------------- %%

% A demo of how to use the Zernike moment function. %

% Example:

% 1- calculate the Zernike moment (n,m) for an oval shape, % 2- rotate the oval shape around its centeroid, % 3- calculate the Zernike moment (n,m) again,

% 4- the amplitude of the moment (A) should be the same for both images % 5- the phase (Phi) should be equal to the angle of rotation % ------------------------------------------------------------------------- clc;

n = 4; m = 2; % Define the order and the iteration of the moment

disp('------------------------------------------------');

disp(['Calculating Zernike moments ..., n = ' num2str(n) ', m = ' num2str(m)]);

p = rgb2gray(imread('Oval_H.png')); figure(1);subplot(1,3,1);imshow(p); title('Horizantal oval');

p = logical(not(p)); N = size(p,1); tic

[ZOH AOH PhiOH] = Zernikmoment(p,n,m); % Call Zernikemoment fuction Elapsed_time = toc;

xlabel({['A = ' num2str(AOH)]; ['\\phi = ' num2str(PhiOH)]});

p = rgb2gray(imread('Oval_V.png')); figure(1);subplot(1,3,2);imshow(p); title('Vertical oval'); p = logical(not(p));

[ZOV AOV PhiOV] = Zernikmoment(p,n,m); % Call Zernikemoment fuction xlabel({['A = ' num2str(AOV)]; ['\\phi = ' num2str(PhiOV)]});

p = rgb2gray(imread('Rectangular_H.png')); figure(1);subplot(1,3,3);imshow(p); title('Vertical rectangular'); p = logical(not(p));

[ZRH ARH PhiRH] = Zernikmoment(p,n,m); % Call Zernikemoment fuction xlabel({['A = ' num2str(ARH)]; ['\\phi = ' num2str(PhiRH)]});

disp('Calculation is complete.');

disp(['The elapsed time per image is ' num2str(Elapsed_time) ' seconds']);

Zernikmoment(p,n,m)

% ------------------------------------------------------------------------- % Copyright C 2012 Amir Tahmasbi % The University of Texas at Dallas % a.tahmasbi@utdallas.edu

% http://www.utdallas.edu/~a.tahmasbi/research.html %

% -------------------------------------------------------------------------

% Licence Agreement: You are free to use this code in your scientific

% research but you should cite the following paper: %

% [1] - A. Tahmasbi, F. Saki, S. B. Shokouhi,

% Classification of Benign and Malignant Masses Based on Zernike Moments, % J. Computers in Biology and Medicine, vol. 41, no. 8, pp. 726-735, 2011. %

% [2] - A. Tahmasbi, F. Saki, H. Aghapanah, S. B. Shokouhi,

% A Novel Breast Mass Diagnosis System based on Zernike Moments as Shape and Density Descriptors,

% in Proc. IEEE, 18th Iranian Conf. on Biomedical Engineering (ICBME'2011), % Tehran, Iran, 2011, pp. 100-104. %

% ------------------------------------------------------------------------- %%

% Function to find the Zernike moments for an N x N binary ROI %

% [Z A Phi] = Zernikmoment(p,n,m) % where

% p = input image N x N (N should be an even number) % n = The order of Zernike moment

% m = The iterration number of Zernike moment % and

% Z = Complex Zernike moment % A = Amplitude of the moment

% Phi = phase (angle) of the mement (in degrees) %

% Example:

% 1- calculate the Zernike moment (n,m) for an oval shape, % 2- rotate the oval shape around its centeroid, % 3- calculate the Zernike moment (n,m) again,

function [Z A Phi] = Zernikmoment(p,n,m) N = size(p,1); x = 1:N; y = x;

[X,Y] = meshgrid(x,y);

R = sqrt((2.*X-N-1).^2+(2.*Y-N-1).^2)/N; Theta = atan2((N-1-2.*Y+2),(2.*X-N+1-2)); R = (R<=1).*R;

Rad = radialpoly(R,n,m); % get the radial polynomial

Product = p(x,y).*Rad.*exp(-1i*m*Theta);

Z = sum(Product(:)); % calculate the moments

cnt = nnz(R)+1; % count the number of pixels inside the unit circle Z = (n+1)*Z/cnt; % normalize the amplitude of moments A = abs(Z); % calculate the amplitude of the moment

Phi = angle(Z)*180/pi; % calculate the phase of the mement (in degrees)

% -------------------------------------------------------------------------

radialpoly(r,n,m)

% ------------------------------------------------------------------------- % Copyright C 2012 Amir Tahmasbi % The University of Texas at Dallas % a.tahmasbi@utdallas.edu

% http://www.utdallas.edu/~a.tahmasbi/research.html %

% -------------------------------------------------------------------------

% Licence Agreement: You are free to use this code in your scientific

% research but you should cite the following paper: %

% [1] - A. Tahmasbi, F. Saki, S. B. Shokouhi,

% Classification of Benign and Malignant Masses Based on Zernike Moments, % J. Computers in Biology and Medicine, vol. 41, no. 8, pp. 726-735, 2011. %

% [2] - A. Tahmasbi, F. Saki, H. Aghapanah, S. B. Shokouhi,

% A Novel Breast Mass Diagnosis System based on Zernike Moments as Shape and Density Descriptors,

% in Proc. IEEE, 18th Iranian Conf. on Biomedical Engineering (ICBME'2011), % Tehran, Iran, 2011, pp. 100-104. %

% ------------------------------------------------------------------------- %%

% Function to compute Zernike Polynomials: %

% f = radialpoly(r,n,m) % where

% r = radius

% n = the order of Zernike polynomial

% m = the iteration Number of Zernike moment

% -------------------------------------------------------------------------

function rad = radialpoly(r,n,m)

rad = zeros(size(r)); % Initilization for s = 0:(n-abs(m))/2

c = (-1)^s*factorial(n-s)/(factorial(s)*factorial((n+abs(m))/2-s)*factorial((n-abs(m))/2-s)); rad = rad + c*r.^(n-2*s); end

% -------------------------------------------------------------------------

zernikes.m

% Zernike polynomial calculator % Christina Dunn

% 7 Sept 2010

% christina.r.dunn@gmail.com

% This calculator displays Zernike polynomials for several types of % apertures, including circular, annular, rectangular, hexagonal and % elliptical apertures.

% The circular Zernikes displayed are the FRINGE set, as defined in the % Code V reference Manual, R. C. Juergens, ed. (Optical Research % Associates, Pasadena, CA, 1998), Vol 1, version 8.30.

% The annular Zernikes are from V. N. Mahajan, \"Zernike annular polynomials % for imaging systems with annular pupils,\" J. Opt. Soc. Am., Vol. 71, No.

% 1, pg 75-85 (1981). All other Zernikes are from V. N. Mahajan and G. M. Dai, % \"Orthonormal polynomials in wavefront analysis: analytical solution,\" % Vol. 24, No. 9, pg 2994-3016 (Sept. 2007).

% This calculator makes use of the POLAR3D function written by J M De Freitas. % POLAR3D: A 3-Dimensional Polar Plot Function

% in Matlab. Version 1.2. QinetiQ Ltd, Winfrith Technology Centre, Winfrith, % Dorchester DT2 8XJ. UK. 2 June 2005. %%

function varargout = zernikes(varargin) % ZERNIKES M-file for zernikes.fig

% ZERNIKES, by itself, creates a new ZERNIKES or raises the existing % singleton*. %

% H = ZERNIKES returns the handle to a new ZERNIKES or the handle to % the existing singleton*. %

% ZERNIKES('CALLBACK',hObject,eventData,handles,...) calls the local

% function named CALLBACK in ZERNIKES.M with the given input arguments. %

% ZERNIKES('Property','Value',...) creates a new ZERNIKES or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before zernikes_OpeningFcn gets called. An

% unrecognized property name or invalid value makes property application % stop. All inputs are passed to zernikes_OpeningFcn via varargin. %

% *See GUI Options on GUIDE's Tools menu. Choose \"GUI allows only one % instance to run (singleton)\". %

% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help zernikes

% Last Modified by GUIDE v2.5 11-Aug-2010 11:25:57

% Begin initialization code - DO NOT EDIT gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename, ...

'gui_Singleton', gui_Singleton, ...

'gui_OpeningFcn', @zernikes_OpeningFcn, ... 'gui_OutputFcn', @zernikes_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1})

gui_State.gui_Callback = str2func(varargin{1}); end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else

gui_mainfcn(gui_State, varargin{:}); end

% End initialization code - DO NOT EDIT

%% --- Executes just before zernikes is made visible.

function zernikes_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn. % hObject handle to figure

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % varargin command line arguments to zernikes (see VARARGIN)

warning off all

% Choose default command line output for zernikes handles.output = hObject;

% Add tool bar to figure

set(hObject, 'toolbar', 'figure');

% Colors for tabs

handles.inactiveColor = [0.9 0.9 0.9]; handles.activeColor = [0.80 0.80 0.95];

% Plot style

% Options are 'SURFACE', 'MESH', 'FRINGES', 'PROFILE' handles.plotStyle = 'SURFACE';

set(handles.plotModePanel,'SelectionChangeFcn',@plotModePanel_SelectionChangeFcn);

% The application opens in FRINGE Zernike mode % Possible modes are:

% 'FRINGE','ANNULAR', 'RECT', 'HEX', 'ELLIPSE'

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.activeColor); set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor); handles.currentZernikeMode = 'FRINGE';

% No special inputs for the FRINGE set, so these are set to be invisible set(handles.obscurationSlider, 'Visible', 'off'); set(handles.obscurationLabel, 'Visible', 'off'); set(handles.obscurationMin, 'Visible', 'off'); set(handles.obscurationMax, 'Visible', 'off'); set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'off'); set(handles.rectHalfWidthBox, 'Visible', 'off'); set(handles.rectHalfWidthSlider, 'Visible', 'off'); set(handles.rectHalfWidthMin, 'Visible', 'off'); set(handles.rectHalfWidthMax, 'Visible', 'off');

set(handles.ellipseAspectRatioBox, 'Visible', 'off'); set(handles.ellipseAspectRatioLabel, 'Visible', 'off'); set(handles.ellipseAspectSlider, 'Visible', 'off'); set(handles.ellipseAspectMin, 'Visible', 'off'); set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'on');

% Default values for apertures handles.resolution = 100; handles.obscuration = 0.5; handles.rectHalfWidth = 0.5; handles.ellipseAspectRatio = 0.5;

% Coefficients of all the Zernike terms start out as zero handles.table = zeros(5,6);

handles.table(:,1) = [1:5]; handles.table(:,3) = [6:10]; handles.table(:,5) = [11:15];

set(handles.uitable, 'data', handles.table); handles.terms = zeros(1,15);

% Update handles structure guidata(hObject, handles);

% UIWAIT makes zernikes wait for user response (see UIRESUME) % uiwait(handles.figure1);

%% --- Outputs from this function are returned to the command line. function varargout = zernikes_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT); % hObject handle to figure

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure varargout{1} = handles.output;

%% --- Executes on button press in FRINGEZernikeButton.

function FRINGEZernikeButton_Callback(hObject, eventdata, handles) % hObject handle to FRINGEZernikeButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

handles.currentZernikeMode = 'FRINGE';

set(handles.noVariablesText, 'Visible', 'on');

guidata(hObject, handles);

%% --- Executes on button press in annularZernikeButton.

function annularZernikeButton_Callback(hObject, eventdata, handles) % hObject handle to annularZernikeButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) handles.currentZernikeMode = 'ANNULAR';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.annularZernikeButton, 'BackgroundColor', handles.activeColor); set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.obscurationSlider, 'Visible', 'on'); set(handles.obscurationLabel, 'Visible', 'on'); set(handles.obscurationMin, 'Visible', 'on'); set(handles.obscurationMax, 'Visible', 'on'); set(handles.obscurationRatioBox, 'Visible', 'on');

set(handles.ellipseAspectRatioBox, 'Visible', 'off'); set(handles.ellipseAspectRatioLabel, 'Visible', 'off');

set(handles.ellipseAspectSlider, 'Visible', 'off'); set(handles.ellipseAspectMin, 'Visible', 'off'); set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'off');

guidata(hObject, handles);

%% --- Executes on button press in rectZernikeButton.

% The Rectangular Zernikes tab sets the GUI to draw Zernikes with a % rectangular aperture.

function rectZernikeButton_Callback(hObject, ~, handles) % hObject handle to rectZernikeButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) handles.currentZernikeMode = 'RECT';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.rectZernikeButton, 'BackgroundColor', handles.activeColor); set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.obscurationSlider, 'Visible', 'off'); set(handles.obscurationLabel, 'Visible', 'off'); set(handles.obscurationMin, 'Visible', 'off'); set(handles.obscurationMax, 'Visible', 'off'); set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'on'); set(handles.rectHalfWidthBox, 'Visible', 'on'); set(handles.rectHalfWidthSlider, 'Visible', 'on'); set(handles.rectHalfWidthMin, 'Visible', 'on'); set(handles.rectHalfWidthMax, 'Visible', 'on');

set(handles.noVariablesText, 'Visible', 'off');

guidata(hObject, handles);

%% --- Executes on button press in hexZernikeButton.

function hexZernikeButton_Callback(hObject, eventdata, handles) % hObject handle to hexZernikeButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) handles.currentZernikeMode = 'HEX';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.hexZernikeButton, 'BackgroundColor', handles.activeColor);

set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.noVariablesText, 'Visible', 'on');

guidata(hObject, handles);

%% --- Executes on button press in ellipticalZernikeButton.

function ellipticalZernikeButton_Callback(hObject, eventdata, handles) % hObject handle to ellipticalZernikeButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) handles.currentZernikeMode = 'ELLIPSE';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor); set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.activeColor);

set(handles.ellipseAspectRatioBox, 'Visible', 'on'); set(handles.ellipseAspectRatioLabel, 'Visible', 'on'); set(handles.ellipseAspectSlider, 'Visible', 'on'); set(handles.ellipseAspectMin, 'Visible', 'on'); set(handles.ellipseAspectMax, 'Visible', 'on');

set(handles.noVariablesText, 'Visible', 'off');

guidata(hObject, handles);

%% --- Executes on slider movement.

% Sets the obscuration ratio for the annular Zernike

function obscurationSlider_Callback(hObject, eventdata, handles) % hObject handle to obscurationSlider (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider

% get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles.obscuration = get(hObject, 'Value');

%set the text box

set(handles.obscurationRatioBox, 'String', handles.obscuration);

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties. function obscurationSlider_CreateFcn(hObject, eventdata, handles) % hObject handle to obscurationSlider (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.

if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end %%

function rectHalfWidthBox_Callback(hObject, eventdata, handles) % hObject handle to rectHalfWidthBox (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of rectHalfWidthBox as text

% str2double(get(hObject,'String')) returns contents of rectHalfWidthBox as a double

a = str2double(get(hObject, 'String'));

if a <= 0.9 && a >= 0.1

handles.rectHalfWidth = a;

set(handles.rectHalfWidthSlider, 'Value', a); else

errordlg('The rectangle half width must be between 0 and 1.') end

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties. function rectHalfWidthBox_CreateFcn(hObject, eventdata, handles) % hObject handle to rectHalfWidthBox (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER.

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end

function ellipseAspectRatioBox_Callback(hObject, eventdata, handles) % hObject handle to ellipseAspectRatioBox (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ellipseAspectRatioBox as text

% str2double(get(hObject,'String')) returns contents of ellipseAspectRatioBox as a double

a = str2double(get(hObject, 'String'));

if a <= 1 && a >= 0

handles.ellipseAspectRatio = a;

set(handles.ellipseAspectSlider, 'Value', a); else

errordlg('The aspect ratio must be between 0 and 1.') end

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.

function ellipseAspectRatioBox_CreateFcn(hObject, eventdata, handles) % hObject handle to ellipseAspectRatioBox (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER.

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end %%

function obscurationRatioBox_Callback(hObject, eventdata, handles) % hObject handle to obscurationRatioBox (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of obscurationRatioBox as text

% str2double(get(hObject,'String')) returns contents of obscurationRatioBox as a double

a = str2double(get(hObject, 'String'));

if a <= 1 && a >= 0

handles.obscuration = a;

set(handles.obscurationSlider, 'Value', a); else

errordlg('The obscuration ratio must be between 0 and 1.') end

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties. function obscurationRatioBox_CreateFcn(hObject, eventdata, handles) % hObject handle to obscurationRatioBox (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows. % See ISPC and COMPUTER.

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor','white'); end

%% --- Executes on slider movement.

function rectHalfWidthSlider_Callback(hObject, eventdata, handles) % hObject handle to rectHalfWidthSlider (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider

% get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles.rectHalfWidth = get(hObject, 'Value');

%set the text box

set(handles.rectHalfWidthBox, 'String', handles.rectHalfWidth);

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties. function rectHalfWidthSlider_CreateFcn(hObject, eventdata, handles) % hObject handle to rectHalfWidthSlider (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.

if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end

%% --- Executes on slider movement.

function ellipseAspectSlider_Callback(hObject, eventdata, handles) % hObject handle to ellipseAspectSlider (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider

% get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles.ellipseAspectRatio = get(hObject, 'Value');

%set the text box

set(handles.ellipseAspectRatioBox, 'String', handles.ellipseAspectRatio);

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties. function ellipseAspectSlider_CreateFcn(hObject, eventdata, handles) % hObject handle to ellipseAspectSlider (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.

if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) set(hObject,'BackgroundColor',[.9 .9 .9]); end

% --- Executes on button press in displayButton.

function displayButton_Callback(hObject, eventdata, handles) % hObject handle to displayButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) cla reset

switch handles.currentZernikeMode case 'FRINGE'

drawFringeZernikes(hObject, handles);

case 'ANNULAR'

drawAnnularZernikes(hObject, handles); case 'RECT'

drawRectZernikes(hObject, handles); case 'HEX'

drawHexZernikes(hObject, handles); case 'ELLIPSE'

drawEllipseZernikes(hObject, handles); end

function drawFringeZernikes(hObject, handles)

resolution = handles.resolution;

resolution = handles.resolution; r = linspace(1, 0, resolution)';

a = linspace(0, 2 * pi + 2*pi/resolution, resolution);

F = zeros(resolution, resolution, 15);

F(:,:,1) = ones(resolution, resolution); F(:,:,2) = r * cos(a); F(:,:,3) = r * sin(a);

F(:,:,4) = (2*r.^2 - ones(size(r))) * ones(size(a)); F(:,:,5) = r.^2 * cos(2*a); F(:,:,6) = r.^2 * sin(2*a);

F(:,:,7) = (3*r.^3 - 2*r) * cos(a); F(:,:,8) = (3*r.^3 - 2*r) * sin(a);

F(:,:,9) = (6*r.^4 - 6*r.^2 + 1) * ones(size(a)); F(:,:,10) = r.^3 * cos(3*a); F(:,:,11) = r.^3 * sin(3*a);

F(:,:,12) = (4*r.^4 -3*r.^2) * cos(2*a); F(:,:,13) = (4*r.^4 -3*r.^2) * sin(2*a);

F(:,:,14) = (10*r.^5 - 12*r.^3 + 3*r) * cos(a); F(:,:,15) = (10*r.^5 - 12*r.^3 + 3*r) * sin(a);

result = zeros(resolution, resolution); for n = 1:15

result = result + handles.terms(n) * F(:,:,n); end

[x, y, z] = polar3d(result, 0, 2*pi, 0, 1, 1);

plotZernike(x, y, z, handles);

function drawAnnularZernikes(hObject, handles)

e = handles.obscuration;

resolution = handles.resolution;

r = linspace(1, e, resolution)';

a = linspace(0, 2 * pi + 2*pi/resolution, resolution);

A = zeros(resolution, resolution, 15);

A(:,:,1) = ones(resolution, resolution); A(:,:,2) = 2 * (r/sqrt(1 + e^2)) * cos(a); A(:,:,3) = 2 * (r/sqrt(1 + e^2)) * sin(a);

A(:,:,4) = sqrt(3) *((2*r.^2 - 1 - e^2)/(1 - e^2))*ones(size(a)); A(:,:,5) = sqrt(6) * (r.^2 / sqrt(1 + e^2 + e^4)) * sin(2*a); A(:,:,6) = sqrt(6) * (r.^2 / sqrt(1 + e^2 + e^4)) * cos(2*a); A(:,:,7) = sqrt(8) * ((3*(1+e^2)*r.^3 - 2*(1+ e^2 + e^4)*r)/ ... ((1-e^2)*sqrt((1+4*e^2+e^4)))) * sin(a);

A(:,:,8) = sqrt(8) * ((3*(1+e^2)*r.^3 - 2*(1+ e^2 + e^4)*r)/... ((1-e^2)*sqrt((1+4*e^2+e^4)))) * cos(a);

A(:,:,9) = sqrt(8) * (r.^3 / sqrt(1 + e^2 + e^4 + e^6)) * sin(3*a); A(:,:,10) = sqrt(8) * (r.^3 / sqrt(1 + e^2 + e^4 + e^6)) * cos(3*a);

A(:,:,11) = sqrt(5) * (6*r.^4 - 6*(1+e^2)*r.^2 + 1+ 4*e^2 + e^4)/(1-e^2)^2 ... * ones(size(a));

A(:,:,12) = sqrt(10) * (4 * r.^2 - 3*((1-e^8)/(1-e^6))*r.^2) / ...

sqrt((1-e^2)^-1 * (16*(1-e^10) - 15*(1-e^8)^2/(1-e^6))) * cos(2*a); A(:,:,13) = sqrt(10) * (4 * r.^2 - 3*((1-e^8)/(1-e^6))*r.^2) / ...

sqrt((1-e^2)^-1 * (16*(1-e^10) - 15*(1-e^8)^2/(1-e^6))) * sin(2*a); A(:,:,14) = sqrt(10) * r.^4/sqrt(1 + e^2 + e^4 + e^6 + e^8) * cos(4*a); A(:,:,15) = sqrt(10) * r.^4/sqrt(1 + e^2 + e^4 + e^6 + e^8) * sin(4*a);

result = zeros(resolution, resolution); for n = 1:15

result = result + handles.terms(n) * A(:,:,n); end

[x, y, z] = polar3d(result, 0, 2*pi, e, 1, 1);

plotZernike(x, y, z, handles); %%

function drawRectZernikes(hObject, handles)

a = handles.rectHalfWidth;

resolution = handles.resolution; xRes = linspace(-a, a, resolution);

yRes = linspace(-sqrt(1-a^2), sqrt(1-a^2), resolution);

[x, y] = meshgrid(xRes, yRes); r = sqrt(x.^2 + y.^2);

u = sqrt(9 - 36*a^2 + 103*a^4 - 134*a^6 + 67*a^8); v = sqrt(49 - 196*a^2 + 330*a^4 - 268*a^6 + 134*a^8); t = 1/(128*v*a^4*(1-a^2)^2);

n = 9 - 45*a^2 + 139*a^4 - 237*a^6 + 210*a^8 - 67*a^10;

R = zeros(resolution, resolution,15); R(:,:,1) = ones(size(x)); R(:,:,2) = (sqrt(3)/a) * x; R(:,:,3) = sqrt(3/(1-a^2)) * y;

R(:,:,4) = (sqrt(5)/(2 * sqrt(1 - 2*a^2 + 2*a^4)))*(3*(x.^2 + y.^2) - 1); R(:,:,5) = (3/(a*sqrt(1-a^2))) * x .* y;

R(:,:,6) = (sqrt(5)*(2*a^2*(1-a^2)*sqrt(1-2*a^2+2*a^4))) * ... (3*(1-a^2)^2*x.^2 - 3*a^4*y.^2 - a^2*(1 - 3*a^2 + 2*a^4));

R(:,:,7) = (sqrt(21)/(2*sqrt(27-81*a^2+116*a^4 - 62*a^6)))*(15*(x.^2 + y.^2) - 9 + 4*a^2) .* y; R(:,:,8) = (sqrt(21)/(2*a*sqrt(35-70*a^2+62*a^4)))*(15*r.^2 - 5 - 4*a^2)*x;

R(:,:,9) = (sqrt(5)*sqrt((27-54*a^2+62*a^4)/(1-a^2))/(2*a^2*(27-81*a^2+116*a^4-62*a^6))) ... * (27*(1-a^2)^2*x.^2 - 35*a^4*y.^2 - a^2*(9 - 39*a^2 + 30*a^4)) .* y; R(:,:,10) = (sqrt(5)/(2*a^3*(1-a^2)*sqrt(35-70*a^2+62*a^4)))*...

(35*(1-a^2)^2*x.^2 - 27*a^4*y.^2 - a^2*(21 - 51*a^2 + 30*a^4))*x;

R(:,:,11) = (1/(8*u))*(315*r.^4 - 30*(7+2*a^2)*x.^2-30*(9-2*a^2)*y.^2 + 27 + 16*a^2 - 16*a^4); R(:,:,12) = (3*u/(8*a^2*v*n))*(35*(1-a^2)^2*(18-36*a^2+67*a^4)*x.^4 + ... 630*(1-2*a^2)*(1-2*a^2+2*a^4)*x.^2*y.^2 ... - 35*a^4*(49-98*a^2+67*a^4)*y.^4 - 30*(1-a^2)*(7-10*a^2-12*a^4+75*a^6-67*a^8)*x.^2 ...

-30*a^2*(7-77*a^2+1*a^4-193*a^6+67*a^8)*y.^2 + ... a^2*(1-a^2)*(1-2*a^2)*(70-233*a^2+233*a^4));

R(:,:,13) = (sqrt(21)/(2*a*sqrt(1-3*a^2+4*a^4-2*a^6)))*(5*r.^2-3)* x .* y;

R(:,:,14) = 16*t*(735*(1-a^2)^4*x.^4 - 540*a^4*(1-a^2)^2*x.^2 .* y.^2 + 735*a^8*y.^4 ... - 90*a^2*(1-a^2)^3*(7-9*a^2)*x.^2 ...

+ 90 * a^6 *(1-a^2)*(2-9*a^2)*y.^2 + 3*a^4*(1-a^2)^2*(21-62*a^2+62*a^4)); R(:,:,15) = (sqrt(21)/(2*a^3*(1-a^2)*sqrt(1-3*a^2+4*a^4-2*a^6))) * ... (5*(1-a^2)^2*x.^2 - 5*a^4*y.^2 - a^2*(3 - 9*a^2 + 6*a^4)) .* x .* y;

z = zeros(resolution, resolution); for n = 1:15

z = z + handles.terms(n) * R(:,:,n);

end

plotZernike(x, y, z, handles);

guidata(hObject, handles); %%

function drawHexZernikes(hObject, handles)

resolution = 2*handles.resolution; r = linspace(1, 0, resolution)';

a = linspace(0, 2 * pi + 2*pi/resolution, resolution);

H = zeros(resolution, resolution,15); H(:,:,1) = ones(resolution);

H(:,:,2) = 2 * sqrt(6/5) * r * cos(a); H(:,:,3) = 2 * sqrt(6/5) * r * sin(a);

H(:,:,4) = sqrt(5/43) * (-5 + 12 * r.^2) * ones(size(a)); H(:,:,5) = 2 * sqrt(15/7) * r.^2 * sin(2*a); H(:,:,6) = 2 * sqrt(15/7) * r.^2 * cos(2*a);

H(:,:,7) = 4 * sqrt(42/3685) * (-14 * r + 25 * r.^3) * sin(a); H(:,:,8) = 4 * sqrt(42/3685) * (-14 * r + 25 * r.^3) * cos(a); H(:,:,9) = (4*sqrt(10)/3) * r.^3 * sin(3*a); H(:,:,10) = 4 * sqrt(70/103) * r.^3 * cos(3*a);

H(:,:,11) = (3/sqrt(1072205)) * (737 - 5140 * r.^2 + 6020 * r.^4) * ones(size(a)); H(:,:,12) = (30 / sqrt(492583)) * (-249 * r.^2 + 392 * r.^4) * cos(2*a); H(:,:,13) = (30 / sqrt(492583)) * (-249 * r.^2 + 392 * r.^4) * sin(2*a);

H(:,:,14) = (10/3) * sqrt(7/99258181) * (10 * (297 - 598 * r.^2) .* r.^2 * cos(2*a) + 5413 * r.^4 * cos(4*a));

H(:,:,15) = (10/3) * sqrt(7/99258181) * (-10 * (297 - 598 * r.^2) .* r.^2 * sin(2*a) + 5413 * r.^4 * sin(4*a));

result = zeros(resolution, resolution); for n = 1:15

result = result + handles.terms(n) * squeeze(H(:,:,n)); end

[x, y, z] = polar3d(result, 0, 2*pi, 0, 1, 1);

xHex = [1/2; 1; 1/2; -1/2; -1; -1/2; 1/2];

yHex = [sqrt(3)/2; 0 ; -sqrt(3)/2; -sqrt(3)/2; 0; sqrt(3)/2; sqrt(3)/2]; in = inpolygon(x, y, xHex, yHex); x(~in) = NaN; y(~in) = NaN;

z(~in) = NaN;

plotZernike(x, y, z, handles); %%

function drawEllipseZernikes(hObject, handles)

resolution = 2*handles.resolution; r = linspace(1, 0, resolution)';

a = linspace(0, 2 * pi + 2*pi/resolution, resolution); b = handles.ellipseAspectRatio;

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 - 60*b^2 + 114*b^4 - 60*b^6 + 35*b^8); delta = sqrt(5 - 6*b^2 + 5*b^4);

E = zeros(resolution, resolution,15); E(:,:,1) = ones(resolution); E(:,:,2) = 2 * r * cos(a); E(:,:,3) = (2 * r * sin(a))/b;

E(:,:,4) = (sqrt(3) / sqrt(3 - 2*b^2 + 3*b^4)) * (-1 - b^2 + 4*r.^2) * ones(size(r))'; E(:,:,5) = (sqrt(6)/b) * r.^2 * sin(2*a); E(:,:,6) = (1/(2*b)) * r.^2 * sin(2*a);

E(:,:,7) = (4/(b*sqrt(5-6*b^2+9*b^4)))*(-1*(1+3*b^2)*r + 6*r.^3) * sin(a); E(:,:,8) = (4/sqrt(9-6*b^2+5*b^4))*(-1*(3+b^2)*r + 6*r.^3)*cos(a); E(:,:,9) = (1/(b^3*sqrt(5-6*b^2+9*b^4)))* ...

(3*(4*b^2*(1-b^2)*r -(5-2*b^2-3*b^4)*r.^3)*sin(a) + ... (5-6*b^2+9*b^4)*r.^3 * sin(3*a));

E(:,:,10) = (1/(b^2*sqrt(9-6*b^2+5*b^4))*(3*(4*b^2*(1-b^2)*r-(3+2*b^2-5*b^4)*r.^3)*... cos(a) + (9-6*b^2+5*b^4)*r.^3*cos(3*a)));

E(:,:,11) = sqrt(5)*(3 + 2*b^2 + 3*b^4 - 24*(1 + b^2)*r.^2 *ones(size(a)) + ... 48*r.^4*ones(size(a)) - 12*(1-b^2)*r.^2 * cos(2*a)) * alpha;

E(:,:,12) = (sqrt(10)*alpha/(gamma*b^2))*(-3*r.^2+4*r.^4)*cos(2*a) + (sqrt(5/2)/(2*b^2*beta)) * ...

(-12*b^2*(5-2*b^2+2*b^6-5*b^8) + 6*(15+125*b^2-194*b^4 + 194*b^6 - 125*b^8 - 15*b^10)*r.^2*ones(size(a)) + ...

240*(-3 + 2*b^2 - 2*b^6 + 3*b^8)*r.^4*ones(size(a)) ...

+ 6*(75-155*b^2 + 174*b^4 - 134*b^6 + 55*b^8 - 15*b^10)*r.^2 * cos(2*a)); E(:,:,13) = (sqrt(10)/b)*(1/delta)*(-3*(1+b^2)*r.^2 + 8*r.^4)*sin(2*a); E(:,:,14) = (sqrt(10)/(8*b^4*gamma))*(3*(1-b^2)^2 * (8*b^4-40*b^2*(1+b^2)*r.^2+5*(7+10*b^2+7*b^4)*r.^4)*ones(size(a)) ...

+4*(6*b^2*(5-7*b^2+7*b^4-5*b^6)-5*(7-6*b^2+6*b^6-7*b^8)*r.^2).*r.^2*cos(2*a) + (35 - 60*b^2 + 114*b^4 ...

-60*b^6 + 35*b^8)*r.^4*cos(4*a));

E(:,:,15) = (sqrt(10)/b^3)*(1/delta)*((6*b^2*(1-b^2)-5*(1-b^4)*r.^2).*r.^2*sin(2*a) + ... ((5 - 6*b^2+5*b^4)/2)*r.^4*sin(4*a));

result = zeros(resolution, resolution); for n = 1:15

result = result + handles.terms(n) * squeeze(E(:,:,n)); end

[x, y, z] = polar3d(result, 0, 2*pi, 0, 1, 1); theta = linspace(0, 2*pi); xEllipse = cos(theta); yEllipse = b * sin(theta);

in = inpolygon(x, y, xEllipse, yEllipse); x(~in) = NaN; y(~in) = NaN; z(~in) = NaN;

plotZernike(x, y, z, handles);

% --- Executes when entered data in editable cell(s) in uitable. function uitable_CellEditCallback(hObject, eventdata, handles) % hObject handle to uitable (see GCBO)

% eventdata structure with the following fields (see UITABLE) % Indices: row and column indices of the cell(s) edited % PreviousData: previous data for the cell(s) edited % EditData: string(s) entered by the user

% NewData: EditData or its converted form set on the Data property. Empty if Data was not changed

% Error: error string when failed to convert EditData to appropriate value for Data % handles structure with handles and user data (see GUIDATA)

handles.table(eventdata.Indices(1), eventdata.Indices(2)) = eventdata.NewData; handles.terms = [handles.table(:, 2)' handles.table(:,4)' handles.table(:,6)'];

guidata(hObject, handles);

% --- Executes on button press in clearButton.

function clearButton_Callback(hObject, eventdata, handles) % hObject handle to clearButton (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Coefficients of all the Zernike terms start out as zero handles.table = zeros(5,6); handles.table(:,1) = [1:5]; handles.table(:,3) = [6:10]; handles.table(:,5) = [11:15];

set(handles.uitable, 'data', handles.table);

handles.terms = zeros(1,15);

axes(handles.profile1) cla reset

axes(handles.profile2) cla reset

axes(handles.profile3) cla reset

axes(handles.display) cla reset

set(handles.display, 'Visible', 'on'); set(handles.profile1, 'Visible', 'off'); set(handles.profile2, 'Visible', 'off'); set(handles.profile3, 'Visible', 'off');

guidata(hObject, handles);

function plotModePanel_SelectionChangeFcn(hObject, eventdata)

%retrieve GUI data, i.e. the handles structure handles = guidata(hObject);

switch get(eventdata.NewValue,'Tag') % Get Tag of selected object case 'surfButton'

handles.plotStyle = 'SURFACE'; case 'meshButton'

handles.plotStyle = 'MESH'; case 'phaseButton'

handles.plotStyle = 'FRINGES'; case 'profileButton'

handles.plotStyle = 'PROFILE'; end

%updates the handles structure

guidata(hObject, handles); %%

function plotZernike(x, y, z, handles)

axes(handles.profile1) cla reset

axes(handles.profile2) cla reset

axes(handles.profile3) cla reset

axes(handles.display) cla reset

switch handles.plotStyle case 'SURFACE';

set(handles.display, 'Visible', 'on'); set(handles.profile1, 'Visible', 'off'); set(handles.profile2, 'Visible', 'off'); set(handles.profile3, 'Visible', 'off');

axes(handles.display) rotate3d on;

axis(gca,'normal'); colormap('default') surf(x, y, z) lighting phong shading interp light axis equal xlim([-1 1]) ylim([-1 1]) view(2) case 'MESH'

set(handles.display, 'Visible', 'on'); set(handles.profile1, 'Visible', 'off'); set(handles.profile2, 'Visible', 'off'); set(handles.profile3, 'Visible', 'off');

axes(handles.display) rotate3d on;

colormap('default') meshc(x, y, z) axis equal

xlim([-1 1]) ylim([-1 1]) view(2) case 'FRINGES'

set(handles.display, 'Visible', 'on'); set(handles.profile1, 'Visible', 'off'); set(handles.profile2, 'Visible', 'off'); set(handles.profile3, 'Visible', 'off');

axes(handles.display) rotate3d off a = [1 1 1];

bFringe = sin(linspace(0, pi, 50)); b = [];

for t = 1:sum(handles.terms) b = [b bFringe]; end

if ~isempty(b)

map = ((a'*b)'); colormap(map)

surf(x, y, zeros(size(x)), z, 'LineStyle', 'none') view([0 90]) axis equal xlim([-1 1]) ylim([-1 1]) end

case 'PROFILE'

set(handles.display, 'Visible', 'off'); set(handles.profile1, 'Visible', 'on'); set(handles.profile2, 'Visible', 'on'); set(handles.profile3, 'Visible', 'on');

axes(handles.profile1) cla reset rotate3d off;

colormap('default') surf(x, y, z) lighting phong shading interp light hold on

zmax = max(max(z));

plot3([-1 1], [0 0], [zmax zmax], 'r'); plot3([0 0], [-1 1], [zmax zmax], 'g'); view([0 90])

dims = size(z);

axes(handles.profile2);

x(isnan(x)) = []; y(isnan(y)) = []; z(isnan(z)) = [];

xprofile = linspace(-1, 1, 40)'; yprofile = linspace(0, 0, 40)';

zprofile = griddata(x(1:10:end,1:10:end), y(1:10:end,1:10:end), z(1:10:end,1:10:end), xprofile, yprofile)';

plot(xprofile, zprofile, 'r'); ylim([-zmax zmax])

axes(handles.profile3)

xprofile = linspace(0, 0, 40)'; yprofile = linspace(-1, 1, 40)';

zprofile = griddata(x(1:10:end,1:10:end), y(1:10:end,1:10:end), z(1:10:end,1:10:end), xprofile, yprofile)';

plot(yprofile, zprofile, 'g'); ylim([-zmax zmax]) end %%

%Following code written by:

%J M De Freitas. POLAR3D: A 3-Dimensional Polar Plot Function

% in Matlab. Version 1.2. QinetiQ Ltd, Winfrith Technology Centre, Winfrith, % Dorchester DT2 8XJ. UK. 2 June 2005.

function [Xout, Yout, Zout] = polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,varargin) %

% POLAR3D Plots a 3D polar surface. %

% POLAR3D(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale) plots

% the profiles as a mesh plot for Zin between radii Rho_min and

% Rho_max for polar angles equally spaced between theta_min and theta_max. %

% POLAR3D(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,plotspec) % plots the profiles Zin between radii Rho_min and Rho_max for % polar angles between theta_min and theta_max with a plot type % specification. If plotspec = 'surf' a standard Matlab surface % is plotted,whereas 'mesh', 'surfc' or 'meshc' will plot mesh, % surface with countour, or mesh with contour, respectively.

% The size of the squares on the mesh or surf plots is determined % by meshscale. The default plot is a mesh plot. % %

% [Xout,Yout,Zout] = POLAR3D(Zin,theta_min,theta_max,Rho_min, % Rho_max, meshscale)returns Zout values corresponding to the % Cartesian positions (Xout,Yout). %

% SYNTAX polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale)

% polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,plotspec) % polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,interpspec)

% polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale, plotspec,interpspec)

% [Xout,Yout,Zout] = polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,...) %

% INPUT Zin input magnitude profiles where each column in Zin is % assumed to represent the radial information in the % plot i.e. each column represents the information % along a radial line defined by theta, where theta % varies between theta_min and theta_max. %

% Zin is a (M x N) matrix, where M and N are not % necessarily equal. If M is not equal to N then the % data are interpolated to make them equal. The final % size is determined by the larger value of (M,N). %

% The N columns of Zin are assumed to be equally

% spaced measurements of the radial components, where % Zin(1,1) corresponds to (theta_min,Rho_max) and

% Zin(M,1) corresponds to (theta_min,Rho_min), and so on. % Zin(1,N) corresponds to (theta_max,Rho_max) and % Zin(M,N) corresponds to (theta_max,Rho_min). Theta % increases in the anticlockwise direction. %

% theta_min the lower value in radians of the angular range

% over which the data is defined. Theta_min is a % scalar quantity. %

% theta_max the upper value in radians of the angular range % over which the data is defined. Theta_max is a % scalar quantity. %

% Rho_min the lower value of the radial component for which % the data is defined. Rho_min is a scalar quantity. %

% Rho_max % %

% meshscale % % % % % % % % %

% plotspec % % % % % % %

% % % % % % %

% % % % % %

the upper value of the radial component for which the data is defined. Rho_max is a scalar quantity. a scalar that determines the size of the squares on the mesh or surf plots. If meshscale is 1, the mesh remains unchanged relative to the input grid; if meshscale is 2, the size of the squares is doubled, if 3.2, say, it is scaled accordingly. Moreover, if meshscale is less than 1, e.g. 0.2, the size of the squares is reduced into a finer mesh. The dimensions of Xout, Yout and Zout are reduced or enlarged by meshscale. = 'surf' produces a surface plot.

= 'surfc' produces a surface plot with contour. = 'mesh' produces a mesh plot.

= 'meshc' produces a mesh plot with countour. = 'contour' produces a 2D contour plot. = 'off' disengages plot function.

= 'meshl' produces a mesh plot with labelled axes. It may be necessary for you to scale Zin when this option is used because the DataAspectRatio is 1:1:1. This is done by replacing 'Zin' with 'n*Zin', e.g. polar3d(7.12*Zin,...) where the scaling factor in this case is 7.12. Note also that when this option is

invoked, the data does not undergo any interpolation; the same size matrix is returned and plotted. 'meshscale' is ignored. % interpspec = 'linear' bilinear interpolation on Zin. % = 'spline' spline interpolation on Zin.

% = 'nearest' nearest neighbour interpolation on Zin. % = 'cubic' bicubic interpolation on Zin. %

% OUTPUT Zout output magnitude profiles defined by Zin at % positions (Xout,Yout). %

% Zout is square with dimensions determined by the % maximum dimension of the input matrix Zin. The

% dimensions of Zout are reduced or enlarged by meshscale. %

% If 'meshl' option is used Zout is same size as Zin. %

% Xout output X-positions corresponding to polar positions % (rho,theta). Xout is square with dimensions

% determined by the maximum dimension of the input % matrix Zin. The dimensions of Xout are reduced or % enlarged by meshscale. %

% If 'meshl' option is used Xout is same size as Zin %

% Yout output Y-positions corresponding to polar positions % (rho,theta). Yout is square with dimensions

% determined by the maximum dimension of the input % matrix Zin. The dimensions of Yout are reduced or % enlarged by meshscale. %

% If 'meshl' option is used Yout is same size as Zin %

% See also POLAR, SPHERE3D, CYL3D, POL2CART and INTERP2

% Written by JM DeFreitas, QinetiQ Ltd, Winfrith Technology % Centre, Dorchester DT2 8XJ, UK. jdefreitas@qinetiq.com. %

% Revision 1.0 4 April 2005.

% Released 5 April 2005. (Beta Release). %

% Revision 1.1 17 April 2005

% 1. Inroduced new check on Zin to cover case of dimensions (M x 2) % or (2 x N) matrix. These are not allowed.

% 2. Changed L2 so that when meshscale > 1 and theta ranges over % 360 degrees the data wraps around without spaces. % 3. Removed Xout1, Yout1 and Zout1.

% 4. Changed 'theta(j,:) = ones([(L2/n) 1])*angl(j);' to % 'theta(j,:) = ones([1 (L2/n)])*angl(j)'; so that it is % compatible with Matlab6 R12.1.

% 5. Reorganised meshgrid so that interpolation now works with % meshscale < 1.

% 6. Changed error traps from '((p ~= 1)&(q ~= 1))' to % '((p ~= 1)|(q ~= 1))' where used. %

% Revision 1.2 2 March 2006

% 1. Changed the 'contourf' option to 'meshl' option. This new % option allows the user to plot a standard mesh plot with % labelled polar axes.

% 2. The 'meshl' option plots with a DataAspectRatio of 1:1:1 and % as such does not automatically scale Zin. It is important for % the user to suitably scale Zin in this eventuality. % 3. The angles are labelled in degrees. % %

% Terms and Conditions of Use %

% 1. This function is made available to Matlab users under the % terms and conditions set out in the Matlab Exchange by % The Mathworks, Inc.

% 2. Where the use of POLAR3D is to be cited, the following is recommended: % J M De Freitas. POLAR3D: A 3-Dimensional Polar Plot Function

% in Matlab. Version 1.2. QinetiQ Ltd, Winfrith Technology Centre, Winfrith, % Dorchester DT2 8XJ. UK. 2 June 2005.

% 3. No offer of warranty is made or implied by the author and % use of this work implies that the user has agreed to take full % responsibility for its use. %

if (nargin < 6)

disp('Polar3d Error: Too few input arguments.'); return

elseif (nargin > 8)

disp('Polar3d Error: Too many input arguments.'); return end

[p,q] = size(theta_min);

if (((p ~= 1)|(q ~= 1))|~isreal(theta_min))|ischar(theta_min)

disp('Polar3d Error: theta_min must be scalar and real.'); return end

[p,q] = size(theta_max);

if (((p ~= 1)|(q ~= 1))|~isreal(theta_max))|ischar(theta_max) disp('Polar3d Error: theta_max must be scalar and real.'); return end

if theta_max <= theta_min

disp('Polar3d Error: theta_max less than or equal theta_min.'); return end

if abs(theta_max - theta_min) > 2*pi

disp('Polar3d Error: range of theta greater than 2pi.'); return end

[p,q] = size(Rho_max);

if (((p ~= 1)|(q ~= 1))|~isreal(Rho_max))|(ischar(Rho_max)|Rho_max < 0) disp('Polar3d Error: Rho_max must be scalar, positive and real.'); return end

[p,q] = size(Rho_min);

if (((p ~= 1)|(q ~= 1))|~isreal(Rho_min))|(ischar(Rho_min)|Rho_min < 0) disp('Polar3d Error: Rho_min must be scalar, positive and real.'); return end

if Rho_max <= Rho_min

disp('Polar3d Error: Rho_max less than or equal Rho_min.'); return end

[p,q] = size(meshscale);

if (((p ~= 1)|(q ~= 1))|~isreal(meshscale))|ischar(meshscale)

disp('Polar3d Warning: mesh scale must be scalar and real.'); meshscale = 1; end

if (meshscale <= 0)

disp('Polar3d Warning: mesh scale must be scalar and positive.'); meshscale = 1; end

% Set up default plot and interpolation specifications. str1 = 'mesh';

str2 = 'linear';

if length(varargin) == 2

% Sort out plot and interpolation specification if both strings given. str1 = [varargin{1}(:)]'; str2 = [varargin{2}(:)]';

g1 = (~isequal(str1,'mesh')&~isequal(str1,'surf'))&~isequal(str1,'off'); g2 = (~isequal(str1,'meshc')&~isequal(str1,'surfc')); g5 = (~isequal(str1,'contour')&~isequal(str1,'meshl')); g3 = (~isequal(str2,'cubic')&~isequal(str2,'linear')); g4 = (~isequal(str2,'spline')&~isequal(str2,'nearest')); if (g1&g2&g5)

disp('Polar3d Warning: Incorrect plot specification. Default to mesh plot.'); str1 = 'mesh'; end

if (g3&g4)

disp('Polar3d Warning: Incorrect interpolation specification.'); disp('Default to linear interpolation.'); str2 = 'linear'; end

elseif length(varargin) == 1

% Sort out plot or interpolation specification from single string input. str1 = [varargin{1}(:)]';

g4 = (~isequal(str1,'spline')&~isequal(str1,'nearest')); if (g1&g2)&(g3&g4&g5)

disp('Polar3d Error: Incorrect plot or interpolation specification.'); return

elseif isequal(str1,'cubic') str2 = str1; str1 = 'mesh';

elseif isequal(str1,'linear') str2 = str1; str1 = 'mesh';

elseif isequal(str1,'spline') str2 = str1; str1 = 'mesh';

elseif isequal(str1,'nearest') str2 = str1; str1 = 'mesh'; elseif isequal(str1,'off') str2 = 'linear';

end end;

% Check if dimensions of input data are acceptable. [r,c] = size(Zin'); if (r < 5)&(c < 5)

disp('Polar3d Error: Input matrix dimensions must be greater than (4 x 4).'); return end

% Check if input data has two rows or columns or less. if (r < 3)|(c < 3)

disp('Polar3d Error: One or more input matrix dimensions too small.'); return end

% Transpose and setup input magnitude matrix. temp = Zin'; for j = 1:c

P(:,j) = temp(:,c-j+1); % swap columns over end Zin = P;

[r,c] = size(Zin);

% Check if meshscale is compatible with dimensions of input data. scalefactor = round(max(r,c)/meshscale); if scalefactor < 3

disp('Polar3d Error: mesh scale incompatible with dimensions of input data.'); return end

% Set up meshgrid corresponding to larger matrix dimension of Zin % for interpolation if required. if ~isequal(str1,'meshl') n = meshscale; if r > c L = r;

L2 = fix(L/n)*n; step = r/(c-1);

[X1,Y1] = meshgrid(0:step:r,1:r); if n < 1

[X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n)); else

[X,Y] = meshgrid(1:n:L2,1:n:L2); end

T = interp2(X1,Y1,Zin,X,Y,str2); elseif c > r L = c;

L2 = fix(L/n)*n; step = c/(r-1);

[X1,Y1] = meshgrid(1:c,0:step:c); if n < 1

[X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n)); else

[X,Y] = meshgrid(1:n:L2,1:n:L2); end

T = interp2(X1,Y1,Zin,X,Y,str2); else

L = r;

L2 = fix(L/n)*n;

[X1,Y1] = meshgrid(1:r,1:r); if n < 1

[X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n)); else

[X,Y] = meshgrid(1:n:L2,1:n:L2); end

T = interp2(X1,Y1,Zin,X,Y,str2); end

[p,q] = size(T); L2 = max(p,q);

% Set up angles

angl = theta_min:abs(theta_max-theta_min)/(L2-1):theta_max; for j = 1:L2

theta(j,:) = ones([1 L2])*angl(j); end

% Set up radial components

Rho = Rho_min:abs(Rho_max-Rho_min)/(L2-1):Rho_max;

% Convert to Cartesian coordinates for j = 1:L2

[Xout(j,:) Yout(j,:) Zout(j,:)] = pol2cart(theta(j,:),Rho,T(j,:)); end

end %if str1 ~= 'meshl'

% Plot Cartesian surface switch str1;

case 'mesh'

colormap([0 0 0]); mesh(Xout,Yout,Zout); axis off; grid off; case 'meshc'

colormap([0 0 0]);

meshc(Xout,Yout,Zout); axis off; grid off; case 'surf'

surf(Xout,Yout,Zout); axis off; grid off; case 'surfc'

surfc(Xout,Yout,Zout); axis off; grid off; hold off case 'contour' axis equal;

h = polar([theta_min theta_max], [Rho_min Rho_max]); delete(h) hold on

contour(Xout,Yout,Zout,20); hold off colorbar; case 'meshl'

% Set up mesh plot with polar axes and labels

angl = theta_min:abs(theta_max-theta_min)/(r-1):theta_max; Rho = ones([1 c])*Rho_max*1.005; X = Rho'*cos(angl); Y = Rho'*sin(angl); Z = zeros(size(X));

% set up output data - this is exactly the same as the input Rho = Rho_min:abs(Rho_max-Rho_min)/(c-1):Rho_max; Xout = Rho'*cos(angl); Yout = Rho'*sin(angl); Zout = Zin'; % plot the data axis equal;

mesh(Xout,Yout,Zout); hold on

% plot the axis

mesh(X,Y,Z,'edgecolor',[0 0 0]); hold on;

% set up tic mark and labels

ticangle = round(theta_min*18/pi)*pi/18:pi/18:round(theta_max*18/pi)*pi/18; ticlength = [Rho_max*1.005 Rho_max*1.03]; Xtic = ticlength'*cos(ticangle); Ytic = ticlength'*sin(ticangle); Ztic = zeros(size(Xtic));

Xlbl = Rho_max*1.1*cos(ticangle); Ylbl = Rho_max*1.1*sin(ticangle); Zlbl = zeros(size(Xlbl));

line(Xtic,Ytic,Ztic,'Color',[0 0 0]); if abs(theta_min-theta_max)==2*pi Ntext = round(length(ticangle)/2)-1; else

Ntext = round(length(ticangle)/2); end

for i = 1:Ntext

text(Xlbl(2*i-1),Ylbl(2*i-1),Zlbl(2*i-1),... num2str(ticangle(2*i-1)*180/pi),... 'horizontalalignment','center') end

set(gca,'DataAspectRatio',[1 1 1]) axis off; grid off; end return

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文