elec_ellipse_fit

 elec_ellipse_fit - Find radii of ellipse and elliptical points for XYZ

 Usage: [r,x,y,z,Xe,Ye,Ze] = elec_ellipse_fit(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

0001 function [r,x,y,z,Xe,Ye,Ze] = elec_ellipse_fit(X,Y,Z,xo,yo,zo,gridSize,plot)
0002 
0003 % elec_ellipse_fit - Find radii of ellipse and elliptical points for XYZ
0004 %
0005 % Usage: [r,x,y,z,Xe,Ye,Ze] = elec_ellipse_fit(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 22:16:48 $
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 eegversion = '$Revision: 1.2 $';
0045 fprintf('ELEC_ELLIPSE_FIT [v %s]\n',eegversion(11:15));
0046 tic;
0047 
0048 % The data is assumed a rough hemisphere.
0049 % To fit to a sphere, we duplicate first, with negative Z.
0050 
0051 Zmin = min(Z); Z = Z + ( 0 - Zmin ); % translate Z so min(Z) = zero
0052 % also used below to recalc Z
0053 X2 = [X;X];
0054 Y2 = [Y;Y];
0055 Z2 = [Z;-Z];
0056 
0057 
0058 % Initialise Ro(x,y,z) as a rough guess at axis radius in X,Y,Z
0059 rX = (max(X2) - min(X2)) / 2;
0060 rY = (max(Y2) - min(Y2)) / 2;
0061 rZ = (max(Z2) - min(Z2)) / 2;
0062 r0 = [ rX rY rZ ];
0063 
0064 fprintf('...initial elliptical radii   (X,Y,Z) = %6.4f, %6.4f, %6.4f\n', r0);
0065 
0066 if r0 == 0,
0067   msg = sprintf('...initial elliptical radii are zero!');
0068   error(msg);
0069 end
0070 
0071 % r(X,Y,Z) = radius in X,Y,Z of fitted ellipse
0072 options = optimset('fminsearch');
0073 [r,fval,exitflag,output] = fminsearch('elec_ellipse_fit_optim', r0, options, X2, Y2, Z2, xo, yo, zo);
0074 ax = r(1);  by = r(2);  cz = r(3);
0075 
0076 fprintf('...estimated elliptical radii (X,Y,Z) = %6.4f, %6.4f, %6.4f\n', r);
0077 
0078 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0079 %
0080 %   Fit all X,Y,Z to ellipsoid
0081 %
0082 %   It'd be better to find the intersection of the line from (xo,yo,zo)
0083 %   to (x,y,z) intersecting the ellipsoid:
0084 %
0085 %   (sX,sY,sZ) where s = 1 / (X/a)^2 + (Y/b)^2 + (Z/c)^2
0086 %
0087 %   Instead, this routine generates points on the ellipsoid surface
0088 %   and then locates the closest of these to each (x,y,z).
0089 
0090 if ~exist('gridSize','var')  gridSize = 50;   end
0091 
0092 [Xe,Ye,Ze] = ellipsoid(xo,yo,zo,ax,by,cz,(gridSize-1));
0093 
0094 Ze = max(Ze,0);    % Force into a hemi-ellipse.
0095 
0096 % For each given X,Y,Z point, locate the nearest Xe,Ye,Ze.
0097 
0098 x = zeros(size(X)); y = zeros(size(Y)); z = zeros(size(Z));
0099 
0100 for e = 1:length(X)
0101   
0102   xe = X(e);  ye = Y(e);  ze = Z(e);
0103   
0104   % Generate matrix of linear distances between point (xe,ye,ze)
0105   % and all points on the ellipsoid (Xe,Ye,Ze).  Could be arc
0106   % length, but doesn't matter here.
0107   d = (Xe-xe).^2 + (Ye-ye).^2 + (Ze-ze).^2; d = sqrt(d);
0108   m = min(min(d));
0109   
0110   index = find (d <= m); i = index(1);
0111   
0112   if isfinite(Xe(i)),
0113     x(e) = Xe(i);
0114   else
0115     fprintf('ELEC_ELLIPSE_FIT: Warning: Xe(i) is NaN: %d\n', Xe(i));
0116   end
0117   if isfinite(Ye(i)),
0118     y(e) = Ye(i);
0119   else
0120     fprintf('ELEC_ELLIPSE_FIT: Warning: Ye(i) is NaN: %d\n', Ye(i));
0121   end
0122   if isfinite(Ze(i)),
0123     z(e) = Ze(i);
0124   else
0125     fprintf('ELEC_ELLIPSE_FIT: Warning: Ze(i) is NaN: %d\n', Ze(i)); 
0126   end
0127   
0128 end
0129 
0130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0131 % Plot the input & projected electrode positions on a sphere
0132 if isequal(plot,1),
0133     figure('NumberTitle','off','Name','Electrode Placements');
0134     set(gca,'Projection','perspective');
0135     set(gca,'DataAspectRatio',[1 1 1]);
0136 
0137     hold on
0138 
0139     plot3(X,Y,Z,'b.');
0140     plot3(x,y,z,'ro');
0141     legend('input xyz','fitted ellipse','Location','BestOutside');
0142 
0143     surf(Xe,Ye,Ze,...
0144         'EdgeColor','none',...
0145         'FaceColor',[0.7 0.7 0.7],...
0146         'FaceAlpha',0.4);
0147     view(2);
0148     rotate3d on;
0149 
0150     axis tight; 
0151     hold off;
0152 end
0153 
0154 z = z + Zmin;                       % restore z to original height
0155 Ze = Ze + Zmin;
0156 
0157 t = toc; fprintf('...done (%6.2f sec).\n\n',t);
0158 
0159 return