Home > bioelectromagnetism > elec_fit_ellipse.m

elec_fit_ellipse

PURPOSE ^

elec_fit_ellipse - Find radii of ellipse and elliptical points for XYZ

SYNOPSIS ^

function [r,x,y,z,Xe,Ye,Ze] = elec_fit_ellipse(X,Y,Z,xo,yo,zo,gridSize,plot)

DESCRIPTION ^

 elec_fit_ellipse - Find radii of ellipse and elliptical points for XYZ

 Usage: [r,x,y,z,Xe,Ye,Ze] = elec_fit_ellipse(X,Y,Z,xo,yo,zo,gridSize,plot)

 Notes:    The general formula for a 3D ellipse, with center (xo,yo,zo) and
           major axis length 2a in X, 2b in Y, and 2c in Z, is given by:

           (x-xo/a)^2  +  (y-yo/b)^2  +  (z-zo/c)^2  =  1

           This function takes arguments for cartesian co-ordinates
           of X,Y,Z and returns the radius 'r(x,y,z)' and the 'x,y,z'
           co-ordinates for the estimated elliptical surface points
           corresponding to X,Y,Z (Z > 0 assumed).

           Also returns Xe,Ye,Ze - the matrix values for an ellipsoid
           of dimensions r(x,y,z).  The size of Xe,Ye,Ze depends on
           'gridSize' (default 50).  Note that x,y,z are a subset of
           Xe,Ye,Ze.

           'plot' is a boolean option for plotting of the input/fitted xyz

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [r,x,y,z,Xe,Ye,Ze] = elec_fit_ellipse(X,Y,Z,xo,yo,zo,gridSize,plot)
0002 
0003 % elec_fit_ellipse - Find radii of ellipse and elliptical points for XYZ
0004 %
0005 % Usage: [r,x,y,z,Xe,Ye,Ze] = elec_fit_ellipse(X,Y,Z,xo,yo,zo,gridSize,plot)
0006 %
0007 % Notes:    The general formula for a 3D ellipse, with center (xo,yo,zo) and
0008 %           major axis length 2a in X, 2b in Y, and 2c in Z, is given by:
0009 %
0010 %           (x-xo/a)^2  +  (y-yo/b)^2  +  (z-zo/c)^2  =  1
0011 %
0012 %           This function takes arguments for cartesian co-ordinates
0013 %           of X,Y,Z and returns the radius 'r(x,y,z)' and the 'x,y,z'
0014 %           co-ordinates for the estimated elliptical surface points
0015 %           corresponding to X,Y,Z (Z > 0 assumed).
0016 %
0017 %           Also returns Xe,Ye,Ze - the matrix values for an ellipsoid
0018 %           of dimensions r(x,y,z).  The size of Xe,Ye,Ze depends on
0019 %           'gridSize' (default 50).  Note that x,y,z are a subset of
0020 %           Xe,Ye,Ze.
0021 %
0022 %           'plot' is a boolean option for plotting of the input/fitted xyz
0023 %
0024 
0025 % $Revision: 1.2 $ $Date: 2005/07/12 21:26:24 $
0026 
0027 % Licence:  GNU GPL no implied or express warranties
0028 % History:  08/99   Chris Harvey
0029 %                   function needs to swap X,Y input for some reason?
0030 %           06/01   Darren.Weber_at_radiology.ucsf.edu
0031 %                   replaced deprecated fmins with fminsearch function;
0032 %                   rectified swapping of X,Y input & scaling of Z;
0033 %                   included input options for xo,yo,zo;
0034 %                   included output of fitted radius estimates
0035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0036 
0037 % initialise centroid, unless input parameters defined
0038 if ~exist('xo','var') xo = 0; else, if isempty(xo) xo = 0; end, end
0039 if ~exist('yo','var') yo = 0; else, if isempty(yo) yo = 0; end, end
0040 if ~exist('zo','var') zo = 0; else, if isempty(zo) zo = 0; end, end
0041 
0042 if ~exist('plot','var') plot = 0; else, if isempty(plot) plot = 0; end, end
0043 
0044 tic; fprintf('\nELEC_FIT_ELLIPSE...\n');
0045 
0046 % The data is assumed a rough hemisphere.
0047 % To fit to a sphere, we duplicate first, with negative Z.
0048 
0049 Zmin = min(Z); Z = Z + ( 0 - Zmin ); % translate Z so min(Z) = zero
0050 % also used below to recalc Z
0051 X2 = [X;X];
0052 Y2 = [Y;Y];
0053 Z2 = [Z;-Z];
0054 
0055 
0056 % Initialise Ro(x,y,z) as a rough guess at axis radius in X,Y,Z
0057 rX = (max(X2) - min(X2)) / 2;
0058 rY = (max(Y2) - min(Y2)) / 2;
0059 rZ = (max(Z2) - min(Z2)) / 2;
0060 r0 = [ rX rY rZ ];
0061 
0062 fprintf('...initial elliptical radii   (X,Y,Z) = %6.4f, %6.4f, %6.4f\n', r0);
0063 
0064 if r0 == 0,
0065   msg = sprintf('...initial elliptical radii are zero!');
0066   error(msg);
0067 end
0068 
0069 % r(X,Y,Z) = radius in X,Y,Z of fitted ellipse
0070 options = optimset('fminsearch');
0071 [r,fval,exitflag,output] = fminsearch('elec_fit_ellipse_optim', r0, options, X2, Y2, Z2, xo, yo, zo);
0072 ax = r(1);  by = r(2);  cz = r(3);
0073 
0074 fprintf('...estimated elliptical radii (X,Y,Z) = %6.4f, %6.4f, %6.4f\n', r);
0075 
0076 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0077 %
0078 %   Fit all X,Y,Z to ellipsoid
0079 %
0080 %   It'd be better to find the intersection of the line from (xo,yo,zo)
0081 %   to (x,y,z) intersecting the ellipsoid:
0082 %
0083 %   (sX,sY,sZ) where s = 1 / (X/a)^2 + (Y/b)^2 + (Z/c)^2
0084 %
0085 %   Instead, this routine generates points on the ellipsoid surface
0086 %   and then locates the closest of these to each (x,y,z).
0087 
0088 if ~exist('gridSize','var')  gridSize = 50;   end
0089 
0090 [Xe,Ye,Ze] = ellipsoid(xo,yo,zo,ax,by,cz,(gridSize-1));
0091 
0092 Ze = max(Ze,0);    % Force into a hemi-ellipse.
0093 
0094 % For each given X,Y,Z point, locate the nearest Xe,Ye,Ze.
0095 
0096 x = zeros(size(X)); y = zeros(size(Y)); z = zeros(size(Z));
0097 
0098 for e = 1:length(X)
0099   
0100   xe = X(e);  ye = Y(e);  ze = Z(e);
0101   
0102   % Generate matrix of linear distances between point (xe,ye,ze)
0103   % and all points on the ellipsoid (Xe,Ye,Ze).  Could be arc
0104   % length, but doesn't matter here.
0105   d = (Xe-xe).^2 + (Ye-ye).^2 + (Ze-ze).^2; d = sqrt(d);
0106   m = min(min(d));
0107   
0108   index = find (d <= m); i = index(1);
0109   
0110   if isfinite(Xe(i)),
0111     x(e) = Xe(i);
0112   else
0113     fprintf('ELEC_FIT_ELLIPSE: Warning: Xe(i) is NaN: %d\n', Xe(i));
0114   end
0115   if isfinite(Ye(i)),
0116     y(e) = Ye(i);
0117   else
0118     fprintf('ELEC_FIT_ELLIPSE: Warning: Ye(i) is NaN: %d\n', Ye(i));
0119   end
0120   if isfinite(Ze(i)),
0121     z(e) = Ze(i);
0122   else
0123     fprintf('ELEC_FIT_ELLIPSE: Warning: Ze(i) is NaN: %d\n', Ze(i)); 
0124   end
0125   
0126 end
0127 
0128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0129 % Plot the input & projected electrode positions on a sphere
0130 if isequal(plot,1),
0131     figure('NumberTitle','off','Name','Electrode Placements');
0132     set(gca,'Projection','perspective');
0133     set(gca,'DataAspectRatio',[1 1 1]);
0134 
0135     hold on
0136 
0137     plot3(X,Y,Z,'b.');
0138     plot3(x,y,z,'ro');
0139     legend('input xyz','fitted ellipse','Location','BestOutside');
0140 
0141     surf(Xe,Ye,Ze,...
0142         'EdgeColor','none',...
0143         'FaceColor',[0.7 0.7 0.7],...
0144         'FaceAlpha',0.4);
0145     view(2);
0146     rotate3d on;
0147 
0148     axis tight; 
0149     hold off;
0150 end
0151 
0152 z = z + Zmin;                       % restore z to original height
0153 Ze = Ze + Zmin;
0154 
0155 t = toc; fprintf('...done (%6.2f sec).\n\n',t);
0156 
0157 return

Generated on Mon 15-Aug-2005 15:36:19 by m2html © 2003