tesselate

(Toolbox/tesselate.m in BrainStorm 2.0 (Alpha))

Function Synopsis

[x,y,z,R,geo,tri_num] = tesselate(shell,spacing,coverage,R,sensnum);

Help Text

TESSELATE - tesselate based on the sensor_ring program
 function [x,y,z,R,geo,tri_num] = tesselate(shell,spacing,coverage,R,sensnum);
 function [x,y,z,R,geo,tri_num] = tesselate(shell,spacing,coverage);
 or
 function [x,y,z,R,geo,tri_num] = tesselate([],[],[],R,sensnum);
 Input: shell is the radius at which to tesselate
        spacing is the nominal length of one side of the triangle
        coverage is the theta angle (radians), measured from the
        z-axis, over which to generate the triangles
        Coverage may also be 'half' or 'full' for pi/2 or pi coverage.
 Optionally, enter garbage for the first three, then enter your sensor
  locations and the number of sensors per ring.  This assumes that all
  locations in R are on a sphere, in circular rings similar to what
  sensor_ring generates, and sensnum might be for example [1 6 12 18] for
  the BTi 37 channel system.
 Output: x,y,z are suitable for fill3(x,y,z,z).  All are 3 x # triangles
  Each column of x represents the x-coordinates of the ith triangle,
  similarly for y and z.
  Optionally, use instead:
  R: is the coordinates of the vertices, one xyz location per row.
  geo: is the geometry matrix, 3 x # triangles.  Each column has
   the integer numbers representing the index from R that forms the
   triangle
   Ordering from vertice 1 to 2 to 3 is such that the vector from 1 to
   2  cross the vector from 1 to 3 is "outward" from the sphere.
  tri_num: the number of triangles in each ring, such that 
   sum(tri_num) is the total number of triangles.

Cross-Reference Information

colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m
sensor_spacing                     C:\BrainStorm_2001\Toolbox\sensor_spacing.m

overlapping_sphere                 C:\BrainStorm_2001\Toolbox\overlapping_sphere.m
source_tess                        C:\BrainStorm_2001\Toolbox\source_tess.m

Listing of function C:\BrainStorm_2001\Toolbox\tesselate.m

function [x,y,z,R,geo,tri_num] = tesselate(shell,spacing,coverage,R,sensnum);
%TESSELATE - tesselate based on the sensor_ring program
% function [x,y,z,R,geo,tri_num] = tesselate(shell,spacing,coverage,R,sensnum);
% function [x,y,z,R,geo,tri_num] = tesselate(shell,spacing,coverage);
% or
% function [x,y,z,R,geo,tri_num] = tesselate([],[],[],R,sensnum);
% Input: shell is the radius at which to tesselate
%        spacing is the nominal length of one side of the triangle
%        coverage is the theta angle (radians), measured from the
%        z-axis, over which to generate the triangles
%        Coverage may also be 'half' or 'full' for pi/2 or pi coverage.
% Optionally, enter garbage for the first three, then enter your sensor
%  locations and the number of sensors per ring.  This assumes that all
%  locations in R are on a sphere, in circular rings similar to what
%  sensor_ring generates, and sensnum might be for example [1 6 12 18] for
%  the BTi 37 channel system.
% Output: x,y,z are suitable for fill3(x,y,z,z).  All are 3 x # triangles
%  Each column of x represents the x-coordinates of the ith triangle,
%  similarly for y and z.
%  Optionally, use instead:
%  R: is the coordinates of the vertices, one xyz location per row.
%  geo: is the geometry matrix, 3 x # triangles.  Each column has
%   the integer numbers representing the index from R that forms the
%   triangle
%   Ordering from vertice 1 to 2 to 3 is such that the vector from 1 to
%   2  cross the vector from 1 to 3 is "outward" from the sphere.
%  tri_num: the number of triangles in each ring, such that 
%   sum(tri_num) is the total number of triangles.

%<autobegin> ---------------------- 26-May-2004 11:34:36 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\sensor_spacing.m
%
% Subfunctions in this file, in order of occurrence in file:
%   r = vnorm(x);
%
% 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.
%<autoend> ------------------------ 26-May-2004 11:34:36 -----------------------


% uses: sensor_spacing, sensor_ring

% John C. Mosher, Ph.D.
% See Copyright.m file for information
% $Date: 5/26/04 10:02a $ $Revision: 12 $
%
% 3/1/95 author
% 3/24/95 pulled xyz out of the loop, added reording for outward normals
% 3/5/96 added option for user to give R and sensnum
%  21 Dec 2001 R12 has problem with norm of infinity. Reported to Mathworks 21 Dec 2001
%  JCM here adds a local norm function that handles infinity just fine.


if(exist('sensnum')~= 1),    % user did not give
  % get vertices
  [R,s,sensnum,theta,no_rings] = sensor_spacing(shell,spacing,coverage);
  clear s theta         % don't need
else                 % user gave R and sensnum
  sensnum = sensnum(:)';    % make sure row
  no_rings = length(sensnum);
end

% R is the vertices, 
% sensnum is vector the number of sensors per ring
% theta is the theta increment
% no_rings is the number of rings (length of sensnum)

% Number of vertices on a closed surface is half the triangles + 2
% Our surface might be possibly open (based on coverage), so I'll
%  just shoot for overkill and reserve a space of three times
%  the vertices, then trim at the end.  Also need three rows per triangle

%xyz = zeros(9*size(R,1),3);    % reserved space for the triangles
geo = zeros(3,3*size(R,1));    % the indices of the triangles
xyzi = 0;             % indexer to xyz

csensnum = cumsum(sensnum);    % cumulative indexer
csensnum = [0 csensnum];    % to align first ring
tri_num = zeros(1,length(sensnum)-1); % number of triangles per ring

for i = 1:(no_rings-1),        % foreach ring of vertices
  pndx = [1:sensnum(i)] + csensnum(i);     % index of previous ring
  nndx = [1:sensnum(i+1)] + csensnum(i+1); % index to next ring
  lp = length(pndx);     % length of previous ring
  ln = length(nndx);     % length of next ring
  if(lp > 1),            % all but the single point ring
    pndx = [pndx pndx(1)]; % wrap around the ring
    lp = lp + 1;
  end
  if(ln > 1),            % all but the single point ring
    nndx = [nndx nndx(1)]; % wrap around the ring
    ln = ln + 1;
  end
  previ = 1;             % indexer to previous ring
  nexti = 1;            % indexer to next ring
  
  while((previ < lp) | (nexti < ln)), % while we walk around the ring
    
    if((previ + 1) <= lp),    % we're not all the way around yet
      test_prev = R(pndx(previ+1),:); % next test vertice for previous ring
    else
      test_prev = Inf;        % no more left on this ring
    end

    if((nexti + 1) <= ln),     % we're not all the way around yet
      test_next = R(nndx(nexti+1),:); % next test vertice for next ring
    else
      test_next = Inf;        % no more left on this ring
    end
    
    % Of the quadrilateral (possibly triangle) prev prev+1 next+1 next
    %  which is the shorter diagonal
    if(vnorm(R(pndx(previ),:)-test_next) < ...
      vnorm(R(nndx(nexti),:)-test_prev)),
      % diagonal from top left to bottom right is winner
%      xyz([1:3]+xyzi,:) = R([pndx(previ) nndx([nexti nexti+1])],:); % vertices
      xyzi = xyzi + 3;        % increment
      geo(:,xyzi/3) = [pndx(previ) nndx([nexti nexti+1])]'; % indices
      nexti = nexti + 1;    % increment
    else
      % winner is bottom left to top right
      % vertices:
%      xyz([1:3]+xyzi,:) = R([pndx(previ) nndx(nexti) pndx(previ+1)],:); 
      xyzi = xyzi + 3;        % increment
      geo(:,xyzi/3) = [pndx(previ) nndx(nexti) pndx(previ+1)]'; % indices
      previ = previ + 1;    % increment
    end                % which diagonal one
    tri_num(i) = tri_num(i) + 1; % number of triangles this ring
  end                % while we are on these rings
    
end                % for all rings

% xyzi represents three times number of triangles
% xyz = xyz(1:xyzi,:);        % trim off the blanks, JCM, pulled out of loop
if(xyzi > 0),            % then we have triangles to consider
  
  geo = geo(:,1:xyzi/3);     % trim off the blanks
  geo = geo([1 3 2],:);     % reverse ordering for outward direction


  %x = reshape(xyz(:,1),3,xyzi/3); % the x vertices, JCM, form directly below
  %y = reshape(xyz(:,2),3,xyzi/3); % the y vertices
  %z = reshape(xyz(:,3),3,xyzi/3); % the z vertices
  x = reshape(R(geo(:),1),3,xyzi/3); % the x vertices
  y = reshape(R(geo(:),2),3,xyzi/3); % the x vertices
  z = reshape(R(geo(:),3),3,xyzi/3); % the x vertices
end

while(0)            % other ideas on picture
%%%%%%%%%%%%%
h = fill3(x,y,z,z);         % return handles of patches
for i = 1:length(h),
  set(h(i),'edgecolor','interp') % no lines
  set(h(i),'facecolor','none')    % mesh grid
end
lightpoint = [-.5 -.5 .25]';        % light source point
C1 = lightpoint*ones(1,size(x,2)) - [x(1,:);y(1,:);z(1,:)];
C1 = colnorm(C1);
C2 = lightpoint*ones(1,size(x,2)) - [x(2,:);y(2,:);z(2,:)];
C2 = colnorm(C2);
C3 = lightpoint*ones(1,size(x,2)) - [x(3,:);y(3,:);z(3,:)];
C3 = colnorm(C3);
C = [C1;C2;C3];            % distances to vertices
%Cl = 1../C;             % lighting inverse to distance
Cl = 1../(C.*C);         % inverse distance squared
h = fill3(x,y,z,Cl,'face','none','edge','interp');
%%%%%%%%%%%%%
end

return

function r = vnorm(x);
% take the norm of a vector, including vectors containing Inf
% for some unknown reason, Matlab R12.1 does not handle this well

r = sqrt(sum(x.^2));
return