showms.new.m

Continuous Profile Models (CPM) Matlab Toolbox.

% function H = showMS(qmzD,merged,myWindow,lineMarkerPos,...
%    timeRatio,keepScans,plotLog,numDig,LOW,SHOW_DIFF,colorScheme)
%
% SHOW_DIFF will also show the difference between the two spectra
% if merged==1
%
% Display one figure per dd item as a 2D MS data, with
% color for the third dimension.
%
% If 'merged'=1, it can only take two samples, and
% superimposes them using two different color schemes,
% and transparency
%
% lineMarkPos is the position of line markers to put in to track
% the shifting traces

function H = showMS(qmzD,merged,myWindow,lineMarkerPos,...
    timeRatio,SHOW_DIFF,keepScans,plotLog,numDig,LOW,colorScheme)

fontSize=12;

if ~iscell(qmzD)
    qmzD={qmzD};
end

if ~exist('merged')|isempty(merged)
    merged=0;
end

if ~(exist('SHOW_DIFF')==1)
    SHOW_DIFF=0;
elseif ~merged
    SHOW_DIFF=0;
end

if length(qmzD)>2 & (merged)
    error('Not yet equipped to handle more than 2 images');
end

if ~exist('colorScheme')
    colorScheme=6;
end

if ~exist('myWindow')
    zoomIn=0;
elseif isempty(myWindow)
    zoomIn=0;
else
    zoomIn=1;
end

if ~exist('timeRatio')
    timeRatio=1;
end

if ~exist('keepScans')
    keepScans=1:length(qmzD{1});
end

if ~(exist('plotLog')==1)
    plotLog=0;
end

if ~(exist('numDig')==1)
    numDig=2;
end
if ~(exist('LOW')==1)
    LOW=400;
end

if merged &(length(qmzD)~=2)
    error('Can only use 2 merged images');
end

%% Load up colormaps
theseMaps{1} = {'blackToRed' 'blackToGreen'};
theseMaps{2} = {'blackToRedNowhite' 'blackToGreenNowhite'};
theseMaps{3} = {'whiteToRed' 'whiteToGreen'};
theseMaps{4} = {'blackToRed2' 'blackToGreen2'};
theseMaps{5} = {'whiteToRed2' 'whiteToGreen2'};
theseMaps{6} = {'whiteToRed3' 'whiteToGreen3'};
thisMap=theseMaps{colorScheme};
clmaps=getColorMaps(thisMap);

if merged
    %% concatenate them to get dual colormaps on same image
    clMap = [clmaps{1} ; clmaps{2}];
else
    clMap = clmaps{1};
end

LOGBASE=2;

for ii=1:length(qmzD)
    %%display the data in 2D
    plotData=qmzD{ii};
    if zoomIn
        %% Y - scan data,  X - M/Z
        rangeYLim = [myWindow(1),myWindow(2)];
        %rangeYLim = rangeYLim*timeRatio;
        rangeY=rangeYLim(1):rangeYLim(2);

        rangeXLim=[myWindow(3),myWindow(4)];
        % account for qmz quanitzation
        rangeXLim=(rangeXLim-(LOW-1))*numDig;
        rangeX=rangeXLim(1):rangeXLim(2);

        %[minx(rangeY) maxx(rangeY)]
        %[minx(rangeX) maxx(rangeX)]
        %whos plotData

        plotData=plotData(rangeY,rangeX);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Tricks with data to get good contrast given the huge
    %% dynamic range of the data and tons of garbage
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% scale the plot data to have min=0 and max=10^7, in a
    %weird way
    if (0)%newWay
        %% get rid outliers, and rescale the data
        %% so that they have same dynamic range
        tempDat = plotData(find(plotData));
        tempSort = sort(tempDat(:));
        tempNum = length(tempSort);
        mymin=0; mymax=1e+08;
        clipmax=tempSort(round(0.99*tempNum));
        clipmin=mymin;
        newPlotData=rescaledata(plotData,mymin,mymax,clipmin,clipmax);
        %% find out where the background is and set everything
        %below it to zero so that there isn't too much crap visible
        tempDat = newPlotData(find(newPlotData));
        tempSort = sort(tempDat(:));
        fracEmpty=0.98;
        backValue = tempSort(round(fracEmpty*tempNum));
        %figure,plot(tempSort,'-'); hold on, plot(backChangeInd,garb,'r+');
        %figure,plot(tempSort,'-'); hold on, plot(backChangeInd,tempSort(backChangeInd),'r+');
        numBelow=length(find(newPlotData<backValue));
        display(['Frac below threshold: ' num2str(numBelow/prod(size(newPlotData)))]);
        newPlotData(find(newPlotData<backValue))=0; %backValue;

        plotData=newPlotData;
        %length(find(plotData==0))/prod(size(plotData))

    else
        BLUR=1;
        if BLUR
            %% blur the data along m/z
            blurMat = [1/3 1/3 1/3];
            plotData2 = conv2(plotData,blurMat);
            %% remove extra padding that got from convolution
            plotData = plotData2(:,2:(end-1));
        end

        mymin=0; mymax=100; midway=(mymax+mymin)/2;
        %mymin=0; mymax=max(plotData(:));
        tempSort = sort(plotData(:));
        clipmin=0; clipmax=tempSort(end-10); % in case of outliers

        %plotData=rescaledata(plotData,mymin,mymax,clipmin,clipmax);
        %fracRemove=0;%.05;
        %plotData(plotData<fracRemove*mymax)=0;

        %mymin=''; mymax='';
        plotData=rescaledata(plotData,mymin,mymax,clipmin,clipmax);
        %% now push the bottom junk to 0
        plotDataRaw{ii}=rescaledata(plotData,mymin,mymax,0,100);
        %% apply a power transform to make stuff more visible      
        plotData = plotData.^0.55;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if ~merged
        colormap(clMap);
        %plotData(1:1,1:1)=2*max(max(plotDataKeep{1})); %HERE?
        H{ii}=image(plotData,'CDataMapping','scaled'); hold on;
    else
        if merged & (ii==2)
            %% Hack for adding two colour maps together on same figure
            %% Found this on the mathworks website!
            plotData = plotData + max(max(plotDataKeep{1}));
        end
        plotDataKeep{ii}=plotData;
        %imstats(plotData);
    end
    if ~merged & (ii>1)
        figure,
    end
end  %% end iterating over images


if merged
    %% play with the data so that when only one of the two
    %% merged images has non zero data, that that one gets
    %% full pixel intensity, but otherwise, they share
    minInk=0;
    for jj=1:2
        if jj==1
            hasInk{jj} = (plotDataKeep{jj}>minInk);
        elseif jj==2
            hasInk{jj} = (plotData>minInk);
        end
    end
    %% find out where there is no competion
    noComp{1} = hasInk{1} & ~hasInk{2};
    noComp{2} = hasInk{2} & ~hasInk{1};
    noCompDat = plotDataKeep{1};
    for jj=1:2
        tmpDat = plotDataKeep{jj};
        noCompDat(find(noComp{jj})) = tmpDat(find(noComp{jj}));
    end
    isComp = ~noComp{1} | ~noComp{2};
    compSecondImage = zeros(size(plotData));
    tmpDat = plotDataKeep{2};
    compSecondImage(find(isComp)) = tmpDat(find(isComp));
    H{1}=image(noCompDat,'CDataMapping','scaled'); hold on;
    H{2}=image(compSecondImage,'CDataMapping','scaled'); hold on;
    colormap(clMap);
    colorbar;
    
    imAlphaData = 0.5*ones(size(plotData));
    imAlphaData(noComp{1})=0;
    imAlphaData(noComp{2})=1;
    set(H{ii},'AlphaData',imAlphaData);

    %% add the colour axes together
    newcaxis=[mymin,length(qmzD)*mymax];
    %newcaxis=[mymin,length(qmzD)*clipmax];
    caxis(newcaxis);
end

ylabel('Time','FontSize',fontSize);
xlabel('Mass/Charge','FontSize',fontSize);
Xlabels=get(gca,'XTickLabel');
if zoomIn
    newXLabels=(str2num(Xlabels))/numDig + myWindow(3);
    Ylabels=get(gca,'YTickLabel');
    newYLabels=str2num(Ylabels);
    %newYLabels=newYLabels + myWindow(1)*timeRatio;
    newYLabels=(newYLabels + myWindow(1)*timeRatio)/timeRatio;
    set(gca,'YTickLabel',newYLabels);
    %end
else
    newXLabels=str2num(Xlabels)/numDig+(LOW-1);
end

%Xlabels
%str2num(Xlabels)/numDig
%str2num(Xlabels)/numDig +LOW

set(gca,'XTickLabel',newXLabels);
myTit='Time Vs. M/Z';
myTit = [myTit ' Experiment ' num2str(ii)];
title(myTit,'FontSize',fontSize);

if (exist('lineMarkerPos')==1)
    drawMarkers(lineMarkerPos,timeRatio,myWindow);        
end

if merged
    ax = findobj(gcf,'Type','axes');
    set(ax,'CLim', [minx(plotDataKeep{1}) maxx(plotDataKeep{2})]);
    %% something about ensuring proper usage of the color map
    set(gcf,'MinColormap',length(get(gcf,'ColorMap')));
end

set(gca,'FontSize',fontSize);

%% also show a map of the differences
if SHOW_DIFF
    %% get axes from previous graph
    Ylabels=get(gca,'YTickLabel');
    YTick=get(gca,'YTick');
    Xlabels=get(gca,'XTickLabel');
    XTick=get(gca,'XTick');
    
    tmpDat = plotDataRaw{1}-plotDataRaw{2};
    figure,image(tmpDat,'CDataMapping','scaled'); hold on;
    useMyColorMap('greenToWhiteToRed');
    trimFrac=0.9999;
    centerColorMap(tmpDat,trimFrac)
    colorbar;    
    
    ylabel('Time','FontSize',fontSize);
    xlabel('Mass/Charge','FontSize',fontSize);
    
    set(gca,'YTick',YTick);
    set(gca,'XTick',XTick);
    set(gca,'YTickLabel', Ylabels);
    set(gca,'XTickLabel', Xlabels);
    
    drawMarkers(lineMarkerPos,timeRatio,myWindow);    
    set(gca,'FontSize',fontSize);
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drawMarkers(lineMarkerPos,timeRatio,myWindow)

for ii=1:length(lineMarkerPos)
    % draw guidelines to track alignments
    if exist('lineMarkerPos')
        if ~isempty(lineMarkerPos)
            guideLabels=['A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I'];
            lsp={'-', '--', ':', '-.','-', '--', ':', '-.','-', '--', ':', '-.'};
            if timeRatio==1
                mycolours1 = {'k','k'};
                mycolours2 = {'k','k'};
            else
                mycolours1 = {'r','g'}; %lines
                mycolours2 = {'r','k'}; %matching text
            end
            ms=lineMarkerPos{ii};
            for jj=1:length(ms)
                thisSpec = [mycolours1{ii} lsp{jj}];
                tmpaxis=axis;
                myLineX=[ceil(tmpaxis(1)),floor(tmpaxis(2))];
                %myLineY=[ms(jj),ms(jj)] - myWindow(1)*timeRatio;
                myLineY=[ms(jj),ms(jj)] - myWindow(1);
                plot(myLineX,myLineY,thisSpec,'LineWidth',1);
                if (timeRatio>1)
                    %% this just offsets the two sets of labels in
                    %case they lie on top of each other
                    text(myLineX(2)+4+10*(ii-1),myLineY(1), guideLabels(jj),'Color', mycolours2{ii},'FontSize',12);
                else
                    text(myLineX(2)+4,myLineY(1), guideLabels(jj),'Color', mycolours2{ii},'FontSize',12);
                end
            end
        end
    end
end