Daniel@0: function Coords = som_unit_coords(topol,lattice,shape)
Daniel@0: 
Daniel@0: %SOM_UNIT_COORDS Locations of units on the SOM grid. 
Daniel@0: %
Daniel@0: % Co = som_unit_coords(topol, [lattice], [shape])
Daniel@0: % 
Daniel@0: %  Co = som_unit_coords(sMap);
Daniel@0: %  Co = som_unit_coords(sMap.topol);
Daniel@0: %  Co = som_unit_coords(msize, 'hexa', 'cyl');
Daniel@0: %  Co = som_unit_coords([10 4 4], 'rect', 'toroid');
Daniel@0: %
Daniel@0: %  Input and output arguments ([]'s are optional): 
Daniel@0: %   topol              topology of the SOM grid
Daniel@0: %             (struct) topology or map struct
Daniel@0: %             (vector) the 'msize' field of topology struct
Daniel@0: %   [lattice] (string) map lattice, 'rect' by default
Daniel@0: %   [shape]   (string) map shape, 'sheet' by default
Daniel@0: %
Daniel@0: %   Co        (matrix, size [munits k]) coordinates for each map unit    
Daniel@0: %
Daniel@0: % For more help, try 'type som_unit_coords' or check out online documentation.
Daniel@0: % See also SOM_UNIT_DISTS, SOM_UNIT_NEIGHS.
Daniel@0: 
Daniel@0: %%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: %
Daniel@0: % som_unit_coords
Daniel@0: %
Daniel@0: % PURPOSE
Daniel@0: %
Daniel@0: % Returns map grid coordinates for the units of a Self-Organizing Map.
Daniel@0: %
Daniel@0: % SYNTAX
Daniel@0: %
Daniel@0: %  Co = som_unit_coords(sTopol);
Daniel@0: %  Co = som_unit_coords(sM.topol);
Daniel@0: %  Co = som_unit_coords(msize);
Daniel@0: %  Co = som_unit_coords(msize,'hexa');
Daniel@0: %  Co = som_unit_coords(msize,'rect','toroid');
Daniel@0: %
Daniel@0: % DESCRIPTION
Daniel@0: %
Daniel@0: % Calculates the map grid coordinates of the units of a SOM based on 
Daniel@0: % the given topology. The coordinates are such that they can be used to
Daniel@0: % position map units in space. In case of 'sheet' shape they can be 
Daniel@0: % (and are) used to measure interunit distances. 
Daniel@0: %
Daniel@0: % NOTE: for 'hexa' lattice, the x-coordinates of every other row are shifted 
Daniel@0: % by +0.5, and the y-coordinates are multiplied by sqrt(0.75). This is done 
Daniel@0: % to make distances of a unit to all its six neighbors equal. It is not 
Daniel@0: % possible to use 'hexa' lattice with higher than 2-dimensional map grids.
Daniel@0: %
Daniel@0: % 'cyl' and 'toroid' shapes: the coordinates are initially determined as 
Daniel@0: % in case of 'sheet' shape, but are then bended around the x- or the 
Daniel@0: % x- and then y-axes to get the desired shape. 
Daniel@0: % 
Daniel@0: % POSSIBLE BUGS
Daniel@0: %
Daniel@0: % I don't know if the bending operation works ok for high-dimensional
Daniel@0: % map grids. Anyway, if anyone wants to make a 4-dimensional
Daniel@0: % toroid map, (s)he deserves it.
Daniel@0: %
Daniel@0: % REQUIRED INPUT ARGUMENTS
Daniel@0: % 
Daniel@0: %  topol          Map grid dimensions.
Daniel@0: %        (struct) topology struct or map struct, the topology 
Daniel@0: %                 (msize, lattice, shape) of the map is taken from 
Daniel@0: %                 the appropriate fields (see e.g. SOM_SET)
Daniel@0: %        (vector) the vector which gives the size of the map grid
Daniel@0: %                 (msize-field of the topology struct).
Daniel@0: %  
Daniel@0: % OPTIONAL INPUT ARGUMENTS 
Daniel@0: % 
Daniel@0: %  lattice (string) The map lattice, either 'rect' or 'hexa'. Default
Daniel@0: %                   is 'rect'. 'hexa' can only be used with 1- or 
Daniel@0: %                   2-dimensional map grids.
Daniel@0: %  shape   (string) The map shape, either 'sheet', 'cyl' or 'toroid'. 
Daniel@0: %                   Default is 'sheet'. 
Daniel@0: %
Daniel@0: % OUTPUT ARGUMENTS
Daniel@0: %
Daniel@0: %  Co   (matrix) coordinates for each map units, size is [munits k] 
Daniel@0: %                where k is 2, or more if the map grid is higher
Daniel@0: %                dimensional or the shape is 'cyl' or 'toroid'
Daniel@0: %
Daniel@0: % EXAMPLES
Daniel@0: %
Daniel@0: % Simplest case:
Daniel@0: %  Co = som_unit_coords(sTopol);
Daniel@0: %  Co = som_unit_coords(sMap.topol);
Daniel@0: %  Co = som_unit_coords(msize);
Daniel@0: %  Co = som_unit_coords([10 10]);
Daniel@0: %
Daniel@0: % If topology is given as vector, lattice is 'rect' and shape is 'sheet'
Daniel@0: % by default. To change these, you can use the optional arguments:
Daniel@0: %  Co = som_unit_coords(msize, 'hexa', 'toroid');
Daniel@0: %
Daniel@0: % The coordinates can also be calculated for high-dimensional grids:
Daniel@0: %  Co = som_unit_coords([4 4 4 4 4 4]);
Daniel@0: %
Daniel@0: % SEE ALSO
Daniel@0: % 
Daniel@0: %  som_unit_dists    Calculate interunit distance along the map grid.
Daniel@0: %  som_unit_neighs   Calculate neighborhoods of map units.
Daniel@0: 
Daniel@0: % Copyright (c) 1997-2000 by the SOM toolbox programming team.
Daniel@0: % http://www.cis.hut.fi/projects/somtoolbox/
Daniel@0: 
Daniel@0: % Version 1.0beta juuso 110997
Daniel@0: % Version 2.0beta juuso 101199 070600
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: %% Check arguments 
Daniel@0: 
Daniel@0: error(nargchk(1, 3, nargin));
Daniel@0: 
Daniel@0: % default values
Daniel@0: sTopol = som_set('som_topol','lattice','rect');
Daniel@0: 
Daniel@0: % topol
Daniel@0: if isstruct(topol), 
Daniel@0:   switch topol.type, 
Daniel@0:   case 'som_map', sTopol = topol.topol;
Daniel@0:   case 'som_topol', sTopol = topol;
Daniel@0:   end
Daniel@0: elseif iscell(topol), 
Daniel@0:   for i=1:length(topol), 
Daniel@0:     if isnumeric(topol{i}), sTopol.msize = topol{i}; 
Daniel@0:     elseif ischar(topol{i}),  
Daniel@0:       switch topol{i}, 
Daniel@0:       case {'rect','hexa'}, sTopol.lattice = topol{i}; 
Daniel@0:       case {'sheet','cyl','toroid'}, sTopol.shape = topol{i}; 
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0: else
Daniel@0:   sTopol.msize = topol;
Daniel@0: end
Daniel@0: if prod(sTopol.msize)==0, error('Map size is 0.'); end
Daniel@0: 
Daniel@0: % lattice
Daniel@0: if nargin>1 & ~isempty(lattice) & ~isnan(lattice), sTopol.lattice = lattice; end
Daniel@0: 
Daniel@0: % shape 
Daniel@0: if nargin>2 & ~isempty(shape) & ~isnan(shape), sTopol.shape = shape; end
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: %% Action
Daniel@0: 
Daniel@0: msize = sTopol.msize;
Daniel@0: lattice = sTopol.lattice;
Daniel@0: shape = sTopol.shape;
Daniel@0: 
Daniel@0: % init variables
Daniel@0: 
Daniel@0: if length(msize)==1, msize = [msize 1]; end
Daniel@0: munits = prod(msize);
Daniel@0: mdim = length(msize);
Daniel@0: Coords = zeros(munits,mdim);
Daniel@0: 
Daniel@0: % initial coordinates for each map unit ('rect' lattice, 'sheet' shape)
Daniel@0: k = [1 cumprod(msize(1:end-1))]; 
Daniel@0: inds = [0:(munits-1)]';
Daniel@0: for i = mdim:-1:1, 
Daniel@0:   Coords(:,i) = floor(inds/k(i)); % these are subscripts in matrix-notation
Daniel@0:   inds = rem(inds,k(i)); 
Daniel@0: end
Daniel@0: % change subscripts to coordinates (move from (ij)-notation to (xy)-notation)
Daniel@0: Coords(:,[1 2]) = fliplr(Coords(:,[1 2])); 
Daniel@0: 
Daniel@0: % 'hexa' lattice
Daniel@0: if strcmp(lattice,'hexa'), 
Daniel@0:   % check
Daniel@0:   if mdim > 2, 
Daniel@0:     error('You can only use hexa lattice with 1- or 2-dimensional maps.');
Daniel@0:   end
Daniel@0:   % offset x-coordinates of every other row 
Daniel@0:   inds_for_row = (cumsum(ones(msize(2),1))-1)*msize(1); 
Daniel@0:   for i=2:2:msize(1), 
Daniel@0:     Coords(i+inds_for_row,1) = Coords(i+inds_for_row,1) + 0.5; 
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: % shapes
Daniel@0: switch shape, 
Daniel@0: case 'sheet', 
Daniel@0:   if strcmp(lattice,'hexa'), 
Daniel@0:     % this correction is made to make distances to all 
Daniel@0:     % neighboring units equal
Daniel@0:     Coords(:,2) = Coords(:,2)*sqrt(0.75); 
Daniel@0:   end
Daniel@0: 
Daniel@0: case 'cyl', 
Daniel@0:   % to make cylinder the coordinates must lie in 3D space, at least
Daniel@0:   if mdim<3, Coords = [Coords ones(munits,1)]; mdim = 3; end
Daniel@0: 
Daniel@0:   % Bend the coordinates to a circle in the plane formed by x- and 
Daniel@0:   % and z-axis. Notice that the angle to which the last coordinates
Daniel@0:   % are bended is _not_ 360 degrees, because that would be equal to 
Daniel@0:   % the angle of the first coordinates (0 degrees).
Daniel@0: 
Daniel@0:   Coords(:,1)     = Coords(:,1)/max(Coords(:,1));
Daniel@0:   Coords(:,1)     = 2*pi * Coords(:,1) * msize(2)/(msize(2)+1);
Daniel@0:   Coords(:,[1 3]) = [cos(Coords(:,1)) sin(Coords(:,1))];
Daniel@0:                     
Daniel@0: case 'toroid', 
Daniel@0: 
Daniel@0:   % NOTE: if lattice is 'hexa', the msize(1) should be even, otherwise 
Daniel@0:   % the bending the upper and lower edges of the map do not match 
Daniel@0:   % to each other
Daniel@0:   if strcmp(lattice,'hexa') & rem(msize(1),2)==1, 
Daniel@0:     warning('Map size along y-coordinate is not even.');
Daniel@0:   end
Daniel@0: 
Daniel@0:   % to make toroid the coordinates must lie in 3D space, at least
Daniel@0:   if mdim<3, Coords = [Coords ones(munits,1)]; mdim = 3; end
Daniel@0: 
Daniel@0:   % First bend the coordinates to a circle in the plane formed
Daniel@0:   % by x- and z-axis. Then bend in the plane formed by y- and
Daniel@0:   % z-axis. (See also the notes in 'cyl').
Daniel@0: 
Daniel@0:   Coords(:,1)     = Coords(:,1)/max(Coords(:,1));
Daniel@0:   Coords(:,1)     = 2*pi * Coords(:,1) * msize(2)/(msize(2)+1);
Daniel@0:   Coords(:,[1 3]) = [cos(Coords(:,1)) sin(Coords(:,1))];
Daniel@0: 
Daniel@0:   Coords(:,2)     = Coords(:,2)/max(Coords(:,2));
Daniel@0:   Coords(:,2)     = 2*pi * Coords(:,2) * msize(1)/(msize(1)+1);
Daniel@0:   Coords(:,3)     = Coords(:,3) - min(Coords(:,3)) + 1;
Daniel@0:   Coords(:,[2 3]) = Coords(:,[3 3]) .* [cos(Coords(:,2)) sin(Coords(:,2))];
Daniel@0:   
Daniel@0: end
Daniel@0: 
Daniel@0: return;
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: %% subfunctions
Daniel@0: 
Daniel@0: function C = bend(cx,cy,angle,xishexa)
Daniel@0: 
Daniel@0:   dx = max(cx) - min(cx);
Daniel@0:   if dx ~= 0, 
Daniel@0:     % in case of hexagonal lattice it must be taken into account that
Daniel@0:     % coordinates of every second row are +0.5 off to the right
Daniel@0:     if xishexa, dx = dx-0.5; end
Daniel@0:     cx = angle*(cx - min(cx))/dx; 
Daniel@0:   end    
Daniel@0:   C(:,1) = (cy - min(cy)+1) .* cos(cx);
Daniel@0:   C(:,2) = (cy - min(cy)+1) .* sin(cx);
Daniel@0: 
Daniel@0: % end of bend
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: