function [L1,L2,L3,L4,L5,L6,L7,L8,L9] = ...
source_grid(sensor_r,source_maxr,source_minr,Factor,STATS,VERBOSE);
%SOURCE_GRID - generate variable grid of possible source locations
% function [L1,L2,L3,L4,L5,L6,L7,L8,L9] = ...
% source_grid(sensor_r,source_maxr,source_minr,Factor,STATS,VERBOSE);
% Given the scalar radius to the sensors sensor_r,
% the maximum scalar radius to the sources source_rmax,
% the minimum scalar radius to the sources source_rmin,
% then generate possible dipole source locations L in the upper hemisphere.
%
% Optional Factor is the fraction considered significant, default 0.1;
% Optional STATS flag means prepend stats onto the first rows of L
% Optional VERBOSE flag means describe the results.
% Optional [L1,L2,L3,...,L9] = source_grid(...
% generate higher order expansion points.
%
% If(STATS), then in each location matrix,
% L(1) give the number of surfaces used,
% L(2,:) to L(L(1)+1,:) are statistics. The first column of the statistics
% rows gives the radial distance of the surface, the second column gives
% the spacing used on that surface, and the third gives the number of
% source points assigned on that surface. Thus delete L(1:(L(1)+1),:)
% after viewing or extracting the statistics.
% ---------------------- 26-May-2004 11:34:30 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\sensor_ring.m
%
% At Check-in: $Author: Mosher $ $Revision: 12 $ $Date: 5/26/04 10:02a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
% University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
% Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
% CNRS, Hopital de la Salpetriere, Paris, France
%
% See BrainStorm website at http://neuroimage.usc.edu for further information.
%
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
% ------------------------ 26-May-2004 11:34:30 -----------------------
% John C. Mosher, Ph.D.
% See Copyright.m file for information
% $Date: 5/26/04 10:02a $ $Revision: 12 $
%
% February 8, 1995 author
% JCM 2/27/95 cleanup from test version to production
% JCM 7/26/95 move first gridding surface down somewhat below the source_maxr
% uses: sensor_ring
if(exist('Factor') ~= 1), % user gave no Factor
Factor = 0.1; % what Factor is considered insignificant
end
if(exist('STATS') ~= 1), % user did not give
STATS = 0; % do not prepend statistics
end
if(exist('VERBOSE') ~= 1), % user did not give
VERBOSE = 0; % silent running
end
for order = 1:nargout, % foreach requested expansion order
% Theory is that an uncovered source is min_d from the nearest sensor. We
% want to be b_offset below this source, such that b_offset + min_d provide
% sufficient coverage. The radius of coverage is such that two
% intersecting circles intersect such that the radius by distance is still
% equivalent to the Factor. So radius of coverage = Factor *
% (min_d+offset). Spacing = sqrt(2)*radius of coverage. Interlace
% sensor_ring routine by 1/2 the spacing on every other surface.
next_source_r = source_maxr; % initialize first source surface
Rdip = []; % source locations
Rstat = []; % statistics
isurfaces = 0; % number of surfaces used
while(next_source_r > source_minr), % for appropriate spherical surfaces
% minimum distance between source and sensor
min_d = sensor_r - next_source_r;
% calculate radius of coverage of these expansions
if(isurfaces == 0), % first surface only
b_offset = min_d/((1/Factor)*sqrt(2)-1);
else
b_offset = min_d/((1/Factor)-1);
end
radius_of_coverage = Factor*(min_d+b_offset);
% set our expansion surface
next_expan_r = next_source_r - b_offset;
% keep sensible
next_expan_r = max(next_expan_r,0);
% create equally spaced expansion points on this surface,
% spaced sqrt(2)*radius_of_coverage apart
isurfaces = isurfaces + 1;
% spacing between sensors on this expansion surface
spacing = sqrt(2)*radius_of_coverage; % (equal to the b_offset*2)
if(spacing > (2*next_expan_r)), % spacing is bigger than radius
Rn = [0 0 next_expan_r]; % single sensor point
else % radius bigger than spacing
theta_inc = acos(1 - 0.5*(spacing/next_expan_r)^2); % sensor increment
if(rem(isurfaces,2)), % if an odd numbered surface
Rn = sensor_ring(next_expan_r,theta_inc,'half');
theta_inc_old = theta_inc; % drop into old bucke
else % interleave from previous theta_inc
Rn = sensor_ring(next_expan_r,[theta_inc_old/2 theta_inc],'half');
end
end % generate spacing on surface
Rdip = [Rdip;Rn]; % add these source locations
% what was the radius and spacing and how many
Rstat = [Rstat; [next_expan_r spacing size(Rn,1)]]; % statistics
% next surface of uncovered sources is at least b_offset below this
% surface.
next_source_r = next_expan_r - b_offset;
% then repeat for this source radius
end
if(STATS), % user wants statistics as well
eval(sprintf('L%.0f = [isurfaces 0 0;Rstat;Rdip];',order));
else % use doesn't want stats
eval(sprintf('L%.0f = [Rdip];',order));
end
if(VERBOSE),
disp(' ')
disp(sprintf('Order %.0f, %.0f surfaces, %.0f sources.',...
order,isurfaces,size(Rdip,1)));
disp(['Surface distances' sprintf(' %.1f',Rstat(:,1)*1000) ' (mm)'])
disp(['Surface spacings ' sprintf(' %.1f',Rstat(:,2)*1000) ' (mm)'])
end
end
if(STATS & VERBOSE), % wanted statistics and verbosity
disp(' ')
disp(...
'Remember to delete the first "surfaces + 1" statistical rows of each L')
end
return