elec_fit_sphere

 elec_fit_sphere - find radius of sphere and closest spherical points to XYZ

 Usage: [r,x,y,z] = elec_fit_sphere(X,Y,Z,xo,yo,zo)

 Notes:    The general formula for a sphere, with radius r is given by:

           (x - xo)^2  +  (y - yo)^2  +  (z - zo)^2  =  r^2

           This function takes arguments for cartesian co-ordinates
           of X,Y,Z (assume Z > 0) and the center of the sphere (xo,yo,zo).
           If (xo,yo,zo) are not provided, they are assumed (0,0,0).

           Returns values are the radius 'r' and the (x,y,z) Cartesian
           co-ordinates for the fitted sphere.

           See also, elec_sphere_project.m

0001 function [r,x,y,z] = elec_fit_sphere(X,Y,Z,xo,yo,zo,estim)
0002 
0003 % elec_fit_sphere - find radius of sphere and closest spherical points to XYZ
0004 %
0005 % Usage: [r,x,y,z] = elec_fit_sphere(X,Y,Z,xo,yo,zo)
0006 %
0007 % Notes:    The general formula for a sphere, with radius r is given by:
0008 %
0009 %           (x - xo)^2  +  (y - yo)^2  +  (z - zo)^2  =  r^2
0010 %
0011 %           This function takes arguments for cartesian co-ordinates
0012 %           of X,Y,Z (assume Z > 0) and the center of the sphere (xo,yo,zo).
0013 %           If (xo,yo,zo) are not provided, they are assumed (0,0,0).
0014 %
0015 %           Returns values are the radius 'r' and the (x,y,z) Cartesian
0016 %           co-ordinates for the fitted sphere.
0017 %
0018 %           See also, elec_sphere_project.m
0019 %
0020 
0021 % $Revision: 1.2 $ $Date: 2005/07/12 21:26:24 $
0022 
0023 % Licence:  GNU GPL, no implied or express warranties
0024 % History:  06/01   Darren.Weber_at_radiology.ucsf.edu
0025 %
0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 
0028 % initialise centroid, unless input parameters defined
0029 if ~exist('xo','var') xo = 0; else, if isempty(xo) xo = 0; end, end
0030 if ~exist('yo','var') yo = 0; else, if isempty(yo) yo = 0; end, end
0031 if ~exist('zo','var') zo = 0; else, if isempty(zo) zo = 0; end, end
0032 
0033 % The data is assumed a rough hemisphere.
0034 % To fit to a sphere, we duplicate first, with negative Z.
0035 
0036 Zmin = min(Z); Z = Z + ( 0 - Zmin ); % translate Z so min(Z) = zero
0037 % also used below to recalc Z
0038 X2 = [X;X];
0039 Y2 = [Y;Y];
0040 Z2 = [Z;-Z];
0041 
0042 % Initialise r0 as a rough guess at the sphere radius
0043 rX = (max(X2) - min(X2)) / 2;
0044 rY = (max(Y2) - min(Y2)) / 2;
0045 rZ = (max(Z2) - min(Z2)) / 2;
0046 r0 = mean([ rX rY rZ ]);
0047 
0048 if ~exist('estim','var') estim = 1; 
0049 else, if isempty(estim) estim = 1; end, end
0050 
0051 if estim ==1, 
0052   fprintf('\n%s%f\n', 'Initial spherical radius   = ', r0); end
0053 
0054 % perform least squares estimate of spherical radius (r)
0055 options = optimset('fminsearch');
0056 [r,fval,exitflag,output] = fminsearch('elec_fit_sphere_optim',...
0057   r0, options, X2, Y2, Z2, xo, yo, zo);
0058 %fprintf('\n%s%f\n', 'Iterations = ', output.iterations);
0059 %fprintf('\n%s%d\n', 'Exit = ', exitflag);
0060 
0061 if estim ==1,
0062   fprintf('%s%f\n',   'Estimated spherical radius = ', r); end
0063 
0064 
0065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0066 %   Fit all X,Y,Z to sphere
0067 %
0068 %   It'd be better to find the intersection of the line from (xo,yo,zo)
0069 %   to (x,y,z) intersecting the sphere:
0070 %
0071 %   (sX,sY,sZ) where s = 1 / (X/r)^2 + (Y/r)^2 + (Z/r)^2
0072 %
0073 %   Instead, this routine generates points on the sphere
0074 %   and then locates the closest of these to each (x,y,z).
0075 
0076 [Xs,Ys,Zs] = elec_sphere_points(40,80,r);
0077 
0078 % For each given X,Y,Z point, locate the nearest Xs,Ys,Zs.
0079 
0080 x = zeros(size(X)); y = zeros(size(Y)); z = zeros(size(Z));
0081 
0082 for s = 1:length(X)
0083   
0084   xs = X(s);  ys = Y(s);  zs = Z(s);  % Get original values
0085   
0086   % Generate matrix of linear distances between point (X,Y,Z)
0087   % and all points on the sphere (Xs,Ys,Zs).  Could be arc
0088   % length, but doesn't matter here.
0089   d = sqrt( (Xs-xs).^2 + (Ys-ys).^2 + (Zs-zs).^2 );
0090   m = min(min(d));
0091   
0092   index = find (d <= m); i = index(1);
0093   
0094   if isfinite(Xs(i)) x(s) = Xs(i); 
0095   else fprintf('\n\aWarning: Xs(i) is NaN: %d', Xs(i)); end
0096   if isfinite(Ys(i)) y(s) = Ys(i); 
0097   else fprintf('\n\aWarning: Ys(i) is NaN: %d', Ys(i)); end
0098   if isfinite(Zs(i)) z(s) = Zs(i); 
0099   else fprintf('\n\aWarning: Zs(i) is NaN: %d', Zs(i)); end
0100   
0101 end
0102 z = z + Zmin;    % restore z to original height
0103 Zs = Zs + Zmin;