zernike拟合.doc

Zernike_main.m
% -------------------------------------------------------------------------
% Copyright C 2012 Amir Tahmasbi
% The University of Texas at Dallas
% [e-mail]
% 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
% [e-mail]
% 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,
% 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
% -------------------------------------------------------------------------

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
% [e-mail]
% 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@

% 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.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);

% 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');

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
% han