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 xyz0001 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