function msheatmap(MZ,Y,varargin)
%MSHEATMAP creates a heat map image of a set of spectra
%
%   MSHEATMAP(MZ,Y) shows a heat map, or pseudocolor image, of the
%   spectrograms in Y. Y is a matrix with several spectra arranged rowwise,
%   all sharing the same MZ scale.
%
%   MSHEATMAP(...,'MARKERS',M) sets a list of markers. Positions are marked
%   along the top axis. Default is M = []. 
%
%   MSHEATMAP(...,'LIMITS',L) uses a 1-by-2 vector (L) with the mass-charge
%   range for the desired heat map.
%
%   MSHEATMAP(...,'GROUP',G) sets the class label for every spectrogram
%   used to group the rows of the heat map. G can be a numeric vector or a
%   cell array of strings with the same number of elements as spectrograms
%   in Y.
%  
%   Examples:
%
%       load sample_lo_res
%       P = [3991.4 4598 7964 9160];
%         
%       % Show a heat map of a set of spectrograms limiting the
%       % viewing range and marking the location of the peaks.
%       msheatmap(MZ_lo_res,Y_lo_res,'markers',P,'limit',[3000 10000])
%
%       % Show a heat map of a set of spectrograms grouping them
%       % into two sets.
%       groups = [1 1 2 2 1 1 2 2];
%       msheatmap(MZ_lo_res,Y_lo_res,'markers',P,'group',groups)
%  
%   See also MSALIGN, MSBACKADJ, MSLOWESS, MSNORM, MSRESAMPLE, MSSGOLAY,
%   MSVIEWER. 

%   Copyright 2003-2004 The MathWorks, Inc.
%   $Revision: 1.1.10.1 $  $Date: 2004/12/24 20:43:26 $

% validate MZ and Y

if ~isnumeric(Y) || ~isreal(Y)
   error('Bioinfo:msheatmap:IntensityNotNumericAndReal',...
         'Ion intensities must be numeric and real.') 
end

if ~isnumeric(MZ) || ~isreal(MZ) || ~isvector(MZ)
   error('Bioinfo:msheatmap:MZNotNumericAndReal',...
         'MZ scale must be a numeric and real vector.') 
end

if size(MZ,1) ~= size(Y,1)
   error('Bioinfo:msheatmap:NotEqualNumberOfSamples',...
         ['The ion intensity vector(s) must have the same number'...
          '\nof samples as the MZ scale.'])
end

[numSamples, numSpectrograms] = size(Y);  %#ok

% set defaults
B = [];
L = [0 inf];
classIdsProvided = false;

% get optional input arguments

if numel(varargin)
    if rem(numel(varargin),2)
        error('Bioinfo:msheatmap:IncorrectNumberOfArguments',...
            'Incorrect number of arguments to %s.',mfilename);
    end
    okargs = {'markers','group','limits'};
    for j=1:2:numel(varargin)
        pname = varargin{j};
        pval = varargin{j+1};
        k = strmatch(lower(pname), okargs);  %#ok
        if isempty(k)
            error('Bioinfo:msheatmap:UnknownParameterName',...
                'Unknown parameter name: %s.',pname);
        elseif length(k)>1
            error('Bioinfo:msheatmap:AmbiguousParameterName',...
                'Ambiguous parameter name: %s.',pname);
        else
            switch(k)
                case 1 % markers
                    B = pval;
                    if ~isnumeric(B) || ~isreal(B) || ~isvector(B)
                        error('Bioinfo:msheatmap:MarkersNotNumericAndReal',...
                              'Markers must be numeric and real.') 
                    end
                    B=B(:);
                    if max(B)>max(MZ) || min(B)<min(MZ)
                        error('Bioinfo:msheatmap:MrkersOutOfRange',...
                              'Markers must be within the MZ scale.')
                    end
                case 2 % class id
                    ID = pval(:);
                    classIdsProvided = true;
                    if iscell(ID)
                        if any(~cellfun('isclass',ID,'char'))
                            error('Bioinfo:msheatmap:InvalidCellForID',...
                            'When GROUP is a cell, all its elements must be strings.')
                        end
                        if ismember('',ID)
                            error('Bioinfo:msheatmap:InvalidCellForID',...
                            'Empty strings are not allowed in the grouping vector.')
                        end
                       [ID,validGroups] = grp2idx(ID);
                   elseif isnumeric(ID) || islogical(ID)
                       if islogical(ID)
                           ID = double(ID);
                       end
                       if ~isreal(ID) || any(isnan(ID)) || any(rem(ID,1)) 
                       error('Bioinfo:msheatmap:InvalidNumericForID',...
                            'The grouping input GROUP is a not valid numeric vector.')
                       end
                       [ID,validGroups] = grp2idx(ID);
                   else
                       error('Bioinfo:msheatmap:InvalidTypeForID',...
                       'The grouping input GROUP must be a cell of strings or a numeric vector.')
                   end
                   %check ID has the correct size
                   if numel(ID)~=numSpectrograms
                        error('Bioinfo:msheatmap:InvalidSizeForID',...
                       'Class ID must have the same number of elements as spectrograms in Y.')
                   end 
               case 3 % limits
                    L = pval;
                    if numel(L)~=2 || diff(L)<=0 
                        error('Bioinfo:msheatmap:badRange',...
                              'The limits must be provided in the form [MZmin, MZmax]')
                    end
                    L(1) = max(L(1),0);
            end
        end
    end
end

% new mz vector (for the image axis)
newMZ = (round(max(L(1),max(0,min(MZ)))):round(min(L(2),max(MZ))))';

% remove any negative value
Y = max(0,Y);

% interpolate Y to the new mz values
Y = interp1q(MZ,Y,newMZ);

% reorder Y if class ids were provided
if classIdsProvided 
    [ID,hID]=sort(ID);
    Y = Y(:,hID);
end

% create image
fig = figure;
imagesc(newMZ,1:numSpectrograms,Y')
xlabel('Mass/Charge (M/Z)')
if classIdsProvided 
    ylabel('Spectrogram Groups')
    q = find([diff(ID);1]);
    set(gca,'Ytick',q-[q(1);diff(q)]/2+0.5)
    set(gca,'YtickLabel',validGroups)
    hold on
    for i = 1:length(q)-1
        plot([0 newMZ(end)],q([i i])+0.5,'r','Linewidth',2)
    end
else
    ylabel('Spectrogram Indices')
end
xlabel('Mass/Charge (M/Z)')

% for better exploration of the image ...
toolsmenufcn(fig,'PanX')         % set zoom mode to horizontal constraining
toolsmenufcn(fig,'ZoomX')        % set pan  mode to horizontal constraining

% draw markers
hold on
for i = 1:numel(B)
  plot(B(i),.5-numSpectrograms/100,'Vr','clipping','off',...
       'markersize',5,'linewidth',2,'Tag','biomarker');
end
hold off
setupMarkerListeners; % listeners to update markers when zooming and panning

colormap(privateColorMap(1));

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function setupMarkerListeners  
% helper function to setsup listeners for the markers, so we can detect if
% we would need to change the relocate
hgp     = findpackage('hg');
axesC   = findclass(hgp,'axes');
figureC = findclass(hgp,'figure');  %#ok
% listens when the Ylim of axes has changed
YLimListener = handle.listener(gca,axesC.findprop('YLim'),...
               'PropertyPostSet',{@markersListener,gcf,gca});
% listens when the Xlim of axes has changed
XLimListener = handle.listener(gca,axesC.findprop('XLim'),...
               'PropertyPostSet',{@markersListener,gcf,gca});           
% store the listeners
setappdata(gcf,'msheatmapListeners',[YLimListener, XLimListener]);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function markersListener(hSrc,event,hf,ha) %#ok
% Auto sizes the ylabels
xlim = get(ha,'Xlim');
ylim = get(ha,'Ylim');
hMarkers = findobj(ha,'Tag','biomarker');

for i = 1:numel(hMarkers)
    x = get(hMarkers(i),'Xdata');
    if x>=min(xlim) && x<=max(xlim)
        set(hMarkers(i),'Visible','on')
        set(hMarkers(i),'Ydata',min(ylim)-diff(ylim)/100)
    else
        set(hMarkers(i),'Visible','off')
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pcmap = privateColorMap(selection)
%PRIVATECOLORMAP returns a custom color map
switch selection
    case 1, pts = [0 0 .5 4;
            0 .1 .8 9;
            0 .9 .6 10;
            .9 1 .9 24;
            1 1 0 40;
            1 0 0 40;
            .4 0 0 0];
    otherwise, pts = [0 0 0 128; 1 1 1 0];
end
xcl=1;
for i=1:size(pts,1)-1
    xcl=[xcl,i+1/pts(i,4):1/pts(i,4):i+1];
end
pcmap = interp1(pts(:,1:3),xcl);