function CCXY=STOCSY(X,ppm,peakID,pthr,CCmetric,peakPPM,twoD)
    % Description: this function performs the STOCSY analysis on a selected
    % peak intensity variable defined by the PeakID chemical shift

    % Use the Statistical Total Correlation SpectroscopY (STOCSY) method to
    % generate a contour map of correlations between spectral variables in
    % a 1D NMR dataset    
    
    %Input:  X =any data matrix of dimension observation-by-variables
    %           ppm = a vector of NMR chemical shifts
    %           peakID - the chemical shift value of an intensity variable for STOCSY 
    %           CCmetric – the selection of correlation coefficient measure {'pearson'  or 'spearman', by default ‘pearson’}
    %           pthr - the threshold p-value for testing the hypothesis of no correlation (by default 0.1)
    %           twoD - confirmation of 2D plot
    
    %Output: CC - correlation coefficients between NMR peakID
    %        intensity variable and all other variables
    
    
    % Author: Kirill Veselkov, Imperial College 2009
    
    [nsmpls, nvrbls]=size(X);
    
    if nsmpls>nvrbls
        warning(['a data set needs to be arranged in a input matrix observation-by-variables'])
    end
    
    if nargin<5
        CCmetric='pearson';
    end
    
    if nargin<4||isempty(pthr)
        pthr=0.1;
    end
    
    if ppm(1)>ppm(2)
        ppm=sort(ppm);
    end

    meanSp = mean(X);
    Xmean  = X-(meanSp(ones(1,nsmpls),:));    
    twoD = 1;

    if twoD == 1
        h=waitbar(0,'Calculation the Correlation Matrix for 2D STOCSY plots');
        for i = 1:size(ppm,2)
            h=waitbar(i/size(ppm,2));
            Y=X(:,i);
            
            [CCXY,pXY]     = corrcoeffs(X,Y,CCmetric);
            CCXY(pXY>pthr) = 0;            
            CovXY(i,:)     = Xmean'*(Y-mean(Y))./(nsmpls-1); % calcuation of the covariance matrix
            
%             if size(ppm,1)>size(ppm,2)
%                 ppm=ppm';
%             end

        end
        close(h)

        contlev = 11;
        clev=contlev*2; % Generate an option to be adjusted in the future
        clvl=2.^(0:clev); 
        clvl=clvl(1:clev); 
        clvl=1./clvl;
        maxc=max(max(max(abs(squeeze(CovXY(:,:,1))))));
        NmrData.clvl=clvl.*maxc;
        NmrData.clvl=NmrData.clvl(1:contlev);
        NmrData.highestValue=0.95;    
        NmrData.threshold=-0.95;

        [X,Y]=meshgrid(ppm', ppm');
        squeezed=squeeze(CovXY(:,:,1));

%         squeezed=squeezed.*(squeezed>0); % Filter for positive values only

        hSTOCSYFigure = figure(...
            'Units','normalized',...
            'MenuBar','none',...
            'NumberTitle','Off',...
            'Toolbar','Figure',...
            'OuterPosition',[0.0 0.0 0.8 0.80],...
            'Tag','STOCSY',...
            'Visible','on');
        hAxes = axes('Parent',hSTOCSYFigure,...
            'Units','normalized',...
            'Tag','Axes',...
            'Position',[0.065 0.08 0.90 0.90],...
            'XLim',[0 10],...
            'Box','on');

%         [~,hCB] = plotCCValues(ppm,CovXY',abs(CCXY'),0,1);


%         hcontour = contourfm(topo60c,topo60cR,levels);

        possqueezed=squeezed'.*(squeezed>0);
        [C ]=contour (hAxes,X,Y,possqueezed,NmrData.clvl,'b','tag','positive');
        [C ]=contourcbar (hAxes,X,Y,possqueezed,NmrData.clvl,'b','tag','positive');
        set(hC,'HitTest','off');
        set(hAxes,'Xdir','Reverse');
        set(hAxes,'Ydir','Reverse');
        hold(hAxes, 'on')
        negsqueezed=squeezed.*(squeezed<0);
        [C ]=contour(hAxes,X,Y,negsqueezed,-NmrData.clvl,'r','tag','negative');
        hold(hAxes, 'off')

        

%         fullscreen = get(0,'ScreenSize');
%         Main = figure('Position',fullscreen);
%         [~,hCB] = plotCCValues(ppm,CovXY',abs(CCXY'),0,1);
%         xlabel('\delta,ppm','FontSize',18);
%         set(get(hCB,'ylabel'),'String', 'abs(correlation coefficient)','FontSize',18);
%         ylabel(['Covariance(X, \delta = ',num2str(peakPPM),')'],'FontSize',18);
%         set(gca,'XDir','reverse');
%         title(['STOCSY',' (', '\delta = ',num2str(peakPPM),')'],'FontSize',24);
%         set(gca,'FontSize',16);        

    else   
    %     peakIndex=find(ppm>=peakID,1,'first');
    %     Y=X(:,peakIndex);
        Y=X(:,peakID);
        
        [CCXY,pXY]     = corrcoeffs(X,Y,CCmetric);
        CCXY(CCXY<pthr) = 0;
        CovXY          = Xmean'*(Y-mean(Y))./(nsmpls-1);
        
        if size(ppm,1)>size(ppm,2)
            ppm=ppm';
        end
    
%         fullscreen = get(0,'ScreenSize');
        Main = figure('Color','w');
        [~,hCB] = plotCCValues(ppm,CovXY',abs(CCXY'),0,1);
        xlabel('\delta,ppm','FontSize',18);
        set(gca,'Box','on');
        axis tight
        set(get(hCB,'ylabel'),'String', 'abs(correlation coefficient)','FontSize',18);
        ylabel(['Covariance(X, \delta = ',num2str(peakPPM),')'],'FontSize',18);
        set(gca,'XDir','reverse');
        title(['STOCSY',' (', '\delta = ',num2str(peakPPM),')'],'FontSize',24);
        set(gca,'FontSize',16);
    end
end