wolffd@0: function [color,X]=som_fuzzycolor(sM,T,R,mode,initRGB,S)
wolffd@0: 
wolffd@0: % SOM_FUZZYCOLOR Heuristic contraction projection/soft cluster color coding for SOM 
wolffd@0: % 
wolffd@0: % function [color,X]=som_fuzzycolor(map,[T],[R],[mode],[initRGB],[S])
wolffd@0: %
wolffd@0: %  sM        (map struct)
wolffd@0: %  [T]       (scalar) parameter that defines the speed of contraction 
wolffd@0: %              T<1: slow contraction, T>1: fast contraction. Default: 1
wolffd@0: %  [R]       (scalar) number of rounds, default: 30
wolffd@0: %  [mode]    (string) 'lin' or 'exp', default: 'lin'  
wolffd@0: %  [initRGB] (string) Strings accepted by SOM_COLORCODE,  default: 'rgb2'
wolffd@0: %  [S]       (matrix) MxM matrix a precalculated similarity matrix 
wolffd@0: %  color     (matrix) of size MxRx3 resulting color codes at each step 
wolffd@0: %  X         (matrix) of size MxRx2 coordiantes for projected unit weight vectors 
wolffd@0: %             at each step of iteration. (Color code C is calculated using this
wolffd@0: %             projection.)
wolffd@0: %
wolffd@0: % The idea of the projection is to use a naive contraction model which
wolffd@0: % pulls the units together. Units that are close to each other in the
wolffd@0: % output space (clusters) contract faster into the same point in the
wolffd@0: % projection. The original position for each unit is its location in
wolffd@0: % the topological grid.
wolffd@0: % 
wolffd@0: % This is an explorative tool to color code the map units so that
wolffd@0: % similar units (in the sense of euclidean norm) have similar coloring
wolffd@0: % (See also SOM_KMEANSCOLOR) The tool gives a series of color codings
wolffd@0: % which start from an initial color coding (see SOM_COLORCODE) and
wolffd@0: % show the how the fuzzy clustering process evolves.
wolffd@0: %
wolffd@0: % The speed of contraction is controlled by the input parameter T. If
wolffd@0: % it is high the projection contracts more slowly and reveals more
wolffd@0: % intermediate stages (hierarchy).  A good value for T must be
wolffd@0: % searched manually. It is probable that the default values do not
wolffd@0: % yield good results.
wolffd@0: %
wolffd@0: % The conatrction process may be slow. In this case the mode can be
wolffd@0: % set to 'exp' instead of 'lin', however, then the computing becomes
wolffd@0: % heavier.
wolffd@0: %
wolffd@0: % EXAMPLE
wolffd@0: %
wolffd@0: %  load iris; % or any other map struct sM 
wolffd@0: %  [color]=som_fuzzycolor(sM,'lin',10);
wolffd@0: %  som_show(sM,'color',color);
wolffd@0: %
wolffd@0: % See also SOM_KMEANSCOLOR, SOM_COLORCODE, SOM_CLUSTERCOLOR
wolffd@0: %
wolffd@0: % REFERENCES
wolffd@0: % 
wolffd@0: % Johan Himberg, "A SOM Based Cluster Visualization and Its
wolffd@0: % Application for False Coloring", in Proceedings of International
wolffd@0: % Joint Conference on Neural Networks (IJCNN2000)},
wolffd@0: % pp. 587--592,Vol. 3, 2000
wolffd@0: % 
wolffd@0: % Esa Alhoniemi, Johan Himberg, and Juha Vesanto, Probabilistic
wolffd@0: % Measures for Responses of Self-Organizing Map Units, pp. 286--290,
wolffd@0: % in Proceedings of the International ICSC Congress on Computational
wolffd@0: % Intelligence Methods and Applications (CIMA '99)}, ICSC Academic
wolffd@0: % Press}, 1999
wolffd@0: %
wolffd@0: % Outline of the heuristic
wolffd@0: %
wolffd@0: % First a matrix D of squared pairwise euclidean distances
wolffd@0: % D(i,j)=d(i,j)^2 between map weight vectors is calculated. This
wolffd@0: % matrix is transformed into a similarity matrix S,
wolffd@0: % s(i,j)=exp(-(D(i,j)/(T.^2*v)), where T is a free input parameter and
wolffd@0: % v the variance of all elements of D v=var(D(:)). The matrix is
wolffd@0: % further normalized so that all rows sum to one. The original
wolffd@0: % topological coordinates X=som_unit_coords(sM) are successively
wolffd@0: % averaged using this matrix. X(:,:,i)=S^i*X(:,:,1); As the process is
wolffd@0: % actually a series of successive weighted averagings of the initial
wolffd@0: % coordinates, all projected points eventually contract into one
wolffd@0: % point.  T is a user defined parameter that defines how fast the
wolffd@0: % projection contracts into this center point. If T is too small, the
wolffd@0: % process will end into the center point at once.
wolffd@0: % 
wolffd@0: % In practise, we don't calculate powers of S, but compute
wolffd@0: % 
wolffd@0: %  X(:,:,i)=S.*X(:,:,i-1); % mode: 'lin'
wolffd@0: %
wolffd@0: % The contraction process may be slow if T is selected to be large,
wolffd@0: % then for each step the similarity matrix is squared
wolffd@0: %
wolffd@0: %  X(:,:,i)=S*X(:,:,1); S=S*S % mode: 'exp'
wolffd@0: %
wolffd@0: % The coloring is done using the function SOM_COLORCODE according to
wolffd@0: % the projections in X, The coordinates are rescaled in order to
wolffd@0: % achieve maximum color resolution.
wolffd@0: 
wolffd@0: % Contributed to SOM Toolbox vs2, 2000 by Johan Himberg
wolffd@0: % Copyright (c) by Johan Himberg
wolffd@0: % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0: 
wolffd@0: % Previously rownorm function normalized the rows of S erroneously
wolffd@0: % into unit length, this major bug was corrected 14042003. Now the
wolffd@0: % rownorm normalizes the rows to have unit sum as it should johan 14042003
wolffd@0: 
wolffd@0: %% Check input arguments
wolffd@0: 
wolffd@0: if isstruct(sM), 
wolffd@0:    if ~isfield(sM,'topol')
wolffd@0:       error('Topology field missing.');
wolffd@0:    end
wolffd@0:    M=size(sM.codebook,1);
wolffd@0: else
wolffd@0:    error('Requires a map struct.');
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin<2 | isempty(T),
wolffd@0:    T=1;
wolffd@0: end
wolffd@0: if ~vis_valuetype(T,{'1x1'})
wolffd@0:    error('Input for T must be a scalar.');
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin<3 | isempty(R),
wolffd@0:    R=30;
wolffd@0: end
wolffd@0: if ~vis_valuetype(R,{'1x1'})
wolffd@0:    error('Input for R must be a scalar.');
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 4 | isempty(mode),
wolffd@0:    mode='lin';
wolffd@0: end
wolffd@0: if ~ischar(mode),
wolffd@0:    error('String input expected for mode.');
wolffd@0: else
wolffd@0:    mode=lower(mode);
wolffd@0:    switch mode
wolffd@0:    case {'lin','exp'}
wolffd@0:       ;   
wolffd@0:    otherwise
wolffd@0:       error('Input for mode must be ''lin'' or ''exp''.');
wolffd@0:    end
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 5 | isempty(initRGB)
wolffd@0:    initRGB='rgb2';
wolffd@0: end
wolffd@0: 
wolffd@0: if ischar(initRGB),   
wolffd@0:    try
wolffd@0:       dummy=som_colorcode(sM,initRGB);
wolffd@0:    catch
wolffd@0:       error(['Color code ''' initRGB ''' not known, see SOM_COLORCODE.']);
wolffd@0:    end
wolffd@0: else
wolffd@0:    error('Invalid color code string');   
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin<6 | isempty(S),
wolffd@0:    S=fuzzysimilarity(sM,1./T);
wolffd@0: end
wolffd@0: 
wolffd@0: if ~vis_valuetype(S,{[M M]}),
wolffd@0:    error('Similarity matrix must be a MunitsxMunits matrix.')
wolffd@0: end
wolffd@0: 
wolffd@0: x = maxnorm(som_unit_coords(sM.topol.msize,sM.topol.lattice,'sheet'));
wolffd@0: 
wolffd@0: x = x-repmat(mean(x),size(x,1),1);
wolffd@0: 
wolffd@0: X(:,:,1)=x; 
wolffd@0: color(:,:,1)=som_colorcode(x,'rgb2',1);
wolffd@0: 
wolffd@0: %%% Actions
wolffd@0: 
wolffd@0: for i=1:R,
wolffd@0:    switch mode
wolffd@0:    case 'exp'
wolffd@0:       S=rownorm(S*S);
wolffd@0:       tmpX=S*X(:,:,1);
wolffd@0:    case 'lin'
wolffd@0:       tmpX=S*X(:,:,i);
wolffd@0:    end
wolffd@0:    X(:,:,i+1)=tmpX;
wolffd@0:    color(:,:,i+1)=som_colorcode(X(:,:,i+1),initRGB);
wolffd@0: end
wolffd@0: 
wolffd@0: color(isnan(color))=0;
wolffd@0: 
wolffd@0: function r=fuzzysimilarity(sM,p)
wolffd@0:   % Calculate a "fuzzy response" similarity matrix
wolffd@0:   % sM: map
wolffd@0:   % p: sharpness factor
wolffd@0:   d=som_eucdist2(sM,sM);
wolffd@0:   v=std(sqrt(d(:))).^2;
wolffd@0:   r=rownorm(exp(-p^2*(d./v)));
wolffd@0:   r(~isfinite(r))=0;
wolffd@0:   return;
wolffd@0: 
wolffd@0: 
wolffd@0: function X = rownorm(X)
wolffd@0: 
wolffd@0:   r = sum(X,2);
wolffd@0:   X = X ./ r(:,ones(size(X,2),1)); 
wolffd@0:   return;
wolffd@0: 
wolffd@0: 
wolffd@0: function X = maxnorm(X)
wolffd@0: 
wolffd@0:   for i=1:size(X,2), r = (max(X(:,i))-min(X(:,i))); if r, X(:,i) = X(:,i) / r; end, end
wolffd@0:   return;