% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%                                                        &
%                        2-D SIFT                        &
%                                                        &
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% last update, dec2814, ym
% last update, jan3016, ym

%% General instructions: 
% This run script guides you through the SIFT processing. In each section,
% edit arguments for your data, and execute the commands in the section.
% The "section" is deliniated by the symbol "%%" and horizontal
% lines in the script.

% There are number of ways to execute the commands. 
% 1. Place the cursor in a "section", and press the "Run Section" button in
% the EDITOR tab at the top of this script viewer. An equivalent shortcut
% is "Apple+return".

% 2. "Evaluate Selection" command is also very useful for executing only
% the highlighted lines. You can set up a shortcut for the "Evaluate
% Selection" at Preferences/Keyboard/Shortcus.

% 3. Above is same to simply copy and paste the relevant lines into the
% "Command Window".


% Usually, fold all the if-blocks in this run script by "Apple+'='". ;use
% "Apple+shift+'='" to unfold all the if-blocks to make internal changes.

% If if-blocks don't fold, go to Preferences/Editor/Code Folding, and make
% if-blocks foldable.

%%
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
% 0. Preparation using nmrPipe
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

% 0-1. 
% Convert your NUS data to the nmrPipe format, determine phase correction
% value and bright spectral region in the direct acquisition dimension.

% See the "Section 0" of the instruction manual for more detail.

%%
%\\\\\\\\\\\\\\\\\\\\
% 1. Set Parameters
%\\\\\\\\\\\\\\\\\\\\

%% 1-1. 
% Set the path to the SIFT script
    cd ~/SIFT2DinMATLAB/SIFT_Suite_2D

%% 1-2. 
% Read NUS table    
    load ~/SIFT2DinMATLAB/example/sched
	NUStab=sched;

%% 1-3
% Set other acquisition and processing parameters

nt1     =128;           % # of t1 pts, total t1 grid
nt2     =1024;          % # of t2 pts

nt1z    =256;           % # of t1 pts, after zero-fill
nt2z    =nt2;           % # of t1 pts, after zero-fill

cent1   =116.897;       % center ppm in F1
cent2   =4.804;         % center ppm in F2

obs1    =59.882;        % observ MHz in F1
obs2    =590.967;       % observ MHz in F2

sw1     =4000.0;        % SW for F2
sw2     =8012.8;        % SW for F1

ph1     =[0,0];         % standard phase or
ph2     =[0,0];         % that found during the nmrPipe processing

wnd1    =1;             % sineb function, on=1, off=0
wnd2    =1;             % sineb function, on=1, off=0

fctr1   =1.0;           % sine function ends at Pi*fctr1, << sine function reaches to zero at fctr1*(the length of t1-max), ym
fctr2   =1.0;           % sine function ends at Pi*fctr2

lab1    ='15N';         % axis label F1
lab2    ='1H';          % axis label F2

lowedge1 =(sw1/obs1)/2 + cent1; % calculate lowest field ppm
lowedge2 =(sw2/obs2)/2 + cent2; % calculate lowest field ppm

disp('======================')
disp('section 2 is executed!')
disp('======================')

%%
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
% 2. Import pre-SIFT time domain data
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

%% 2-1(a)
% If zeros are already inserted to unacquired t1 points
%--arguments--
    path= '~/SIFT2DinMATLAB/example/pre.fid'; 
%--commands--
    [read1,read2]=read2d(nt1,nt2,path);
    pre_c=read1;
    pre_s=read2;
%% 2-1(b)
% If not:
%--arguments--
    path= '~/SIFT2DinMATLAB/example/pre.fid'; 
%--commands--
    [read1,read2]=read2dnl(nt1,nt2,path,NUStab);
    pre_c=read1;
    pre_s=read2;

%% 2-2
% Check if the data is correctly imported and reshuffled.

% Set the data just read in, and execute the section.
%--arguments--
    in1=pre_c;
    in2=pre_s;
%--commands--
    if isempty(in1)==0 &&isempty(in2)==0
    figure()
    subplot(2,1,1)
    plot(real(in1(:,1,1)))

    subplot(2,1,2)
    temp2=sineb(in1,90,1,1);
    NFFT=2^nextpow2(size(temp2,1));
    temp2=fft(temp2(:,1),NFFT);
    temp2=mphasec(temp2,sw2,ph2,1);
    temp2=(real(temp2));
%     temp2=fftshift(real(temp2));
    plot(temp2)
    end
    if isempty(in1)==0 &&isempty(in2)==0
    clear temp
    [pref,junk,t1F2]=proc2d(in1,in2,nt1,nt2z,sw1,sw2,ph1,ph2,wnd1,wnd2,fctr1);
    temp2=pref;
    temp2(1:25,:)=0;
    temp2(nt2z-25:nt2z,:)=0;
    [a,b]=max(real(temp2));
    [c,d]=max(a);
    maxpick=b(d);

    figure()
    subplot(2,1,1)
    plot(real(t1F2(maxpick,:)),'.-')
    axis([0,nt1,min(real(t1F2(maxpick,:))),max(real(t1F2(maxpick,:)))])
    grid on
    xlabel('t1 (pts)')
    subplot(2,1,2)
    tab=zeros(1,nt1);
    for n=1:nt1
        if sum(NUStab==n)
            tab(n)=1;
        end
    end
    scatter(1:nt1,tab,'.')
    grid on
    axis([0,nt1,0.9 1.1])
    xlabel('t1')
    end
    
%% 2-3. 
% Axis swap or inversion

% Change flags (0=False, 1=True) according to the acquisition mode you
% used in the experiment, and execute this section.
%--arguments--
    in1=pre_c;
    in2=pre_s;
                        % in t1 -------
    alt1=0;             % States-TPPI/ Echo-Antiecho (axis left/right swap)
    neg1=0;             % Echo-Antiecho (axis inversion)
                        % in t2 -------
    alt2=0;             % States-TPPI/ Echo-Antiecho (axis left/right swap)
    neg2=0;             % Echo-Antiecho (axis inversion)
    
%--commands--
    [pre_c,pre_s]=acqmode2d(in1,in2,alt1,neg1,alt2,neg2);

%%    
%\\\\\\\\\\\\\\\\\\\\\\\\
% 3. Define bright region
%\\\\\\\\\\\\\\\\\\\\\\\\

% SIFT needs you to speficy bright spectral regions, or signal-containing
% regions, so that SIFT can use the dark spectral points to restore the NUS
% data.

%% 3-1.
% Process the pre-SIFT data

%--arguments--
in1=pre_c;
in2=pre_s;
con_min=0.4; % minimum contour level, ratio to the tallest peak
con_max=1.00; % maximum c
ncon=15;
%--commands--
[pre_ft2d,junk,t1F2pre]=proc2d(in1,in2,nt1z,nt2z,sw1,sw2,ph1,ph2,wnd1,wnd2,fctr1);
clear in1 in2
spec1=pre_ft2d;
if isempty(spec1)==0
    con1=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec1(25:nt2z-25,:)))));
    figure()
    spec=real(squeeze(spec1));
    mplp2d(squeeze(real(spec)),sw2,sw1,cent2,cent1,obs2,obs1,con1);
    legend('pre-SIFT')
    xlabel('F1 (ppm)','fontsize',16)
    ylabel('F2 (ppm)','fontsize',16)
end

%% 3-2
% Set the bright regions in one of three different ways

% Method A: 
% Look up the spectral limits from the spectrum plotted in the section 3-1,
% and by directly typing in the limits below.

%--arguments--
    brightLimits_F1=[100 132]; % ppm, [high-field limit, low-field limit]
    brightLimits_F2=[6 11]; 
%--commands--
if isempty(brightLimits_F1)==0 && isempty(brightLimits_F2)==0
    sigf1=[];
    sigf2=[];
    nr=length(brightLimits_F1(:,1));
    for n=1:nr
        [limit1]=ppm2pts(brightLimits_F1(n,:),sw1,nt1z,cent1,obs1,[]);
        [limit2]=ppm2pts(brightLimits_F2(n,:),sw2,nt2z,cent2,obs2,[]);
        
        if cent1>brightLimits_F1(n,1) && cent1<brightLimits_F1(n,2)...
                && cent2>brightLimits_F2(n,1) && cent2<brightLimits_F2(n,2)
            sigf1=[sigf1;
                1            limit1(1,2);
                1            limit1(1,2);
                limit1(1,1)  nt1z;
                limit1(1,1)  nt1z
                ];
            sigf2=[sigf2;
                1            limit2(1,2);
                limit2(1,1)  nt2z;
                1            limit2(1,2);
                limit2(1,1)  nt2z;
                ];
        elseif cent1>brightLimits_F1(n,1) && cent1<brightLimits_F1(n,2)
            sigf1=[sigf1;
                1            limit1(1,2);
                limit1(1,1)  nt1z
                ];
            sigf2=[sigf2;
                limit2(1,1)  limit2(1,2);
                limit2(1,1)  limit2(1,2)
                ];
        elseif cent2>brightLimits_F2(n,1) && cent2<brightLimits_F2(n,2)
            sigf1=[sigf1;
                limit1(1,1)  limit1(1,2);
                limit1(1,1)  limit1(1,2)
                ];
            sigf2=[sigf2;
                1            limit2(1,2);
                limit2(1,1)  nt2z
                ];
        else
            sigf1=[sigf1;
                limit1(1,1)  limit1(1,2)
                ];
            sigf2=[sigf2;
                limit2(1,1)  limit2(1,2)
                ];
        end
        n
    end
    pick=[];
end    

disp('======================')
disp('section 3-2 is executed!')
disp('======================')

%% 3-3
% Method B:
% Use the graphycal input tool 

% Activate the spectrum plotted in the section 3-1, zoom in to the relevant
% spectral region, then execute this section.

% A crosshair cursor will appear on the spectrum. Click on all the peaks
% you need, and hit "return" to exit the cursor mode. A peak list is
% generated, and bright regions are calculated.

%--arguments--
lw1=1;   % F1, ppm
lw2=0.1; % F2, ppm
%--commands--
F1F2ppm=ginput;
if isempty(F1F2ppm)==0
    clear peakList
    peakList(:,1)=F1F2ppm(:,1);
    peakList(:,2)=F1F2ppm(:,2);
    peakList(:,3)=0;
    peakList(:,4)=lw1*obs1; %40;
    peakList(:,5)=lw2*obs2; %20
    
    spec=pre_ft2d;
    list=peakList; 
    step=1;
    refn=1;
    % then, do the following:
    %~*~*~*~*~*~* make mask ~*~*~*~*~*~*~
    [sigf1 sigf2 vpls hpls intensity]=brightPeak_v3(spec,list,step,refn,...
        sw1,sw2,nt1z,nt2z,obs1,obs2,cent1,cent2);
    %~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~
    
    % ppm to pts
    F2=peakList(:,2);
    clear F2pts
    del2=(sw2/obs2)/(nt2z-1);
    for n=1:size(F2,1)
        if F2(n,1)>=cent2
            F2pts(1,n)=floor((F2(n,1)-cent2)/del2)+1;
        else
            F2pts(1,n)=nt2z-floor((cent2-F2(n,1))/del2);
        end
    end
        pick=F2pts    
end

disp('======================')
disp('section 3-3 is executed!')
disp('======================')

%% 3-4
% Method C:
% Use a peak list made in Sparky

% Copy and paste the part of the Sparky peak list as below, and execute the
% section.
%--arguments--
clear peakList
%   F1(ppm)     F2(ppm)    Int          F1,Lw(Hz)   F2,Lw(Hz)
peakList=[
    103.312      7.021     60.409184     40          20
    130.659     10.529     64.048592     40          20
    128.658      8.258     80.483104     40          20
    128.415      7.652     59.215568     40          20
    ];
list=peakList;
step=1; 
refn=1;
spec=pre_ft2d;

%--commands--
if isempty(spec)==0
[sigf1 sigf2 vpls hpls intensity]=brightPeak_v3(spec,list,step,refn,...
    sw1,sw2,nt1z,nt2z,obs1,obs2,cent1,cent2);

% ppm to pts
    F2=peakList(:,2);
    clear F2pts
    del2=(sw2/obs2)/(nt2z-1);
    for n=1:size(F2,1)
        if F2(n,1)>=cent2
            F2pts(1,n)=floor((F2(n,1)-cent2)/del2)+1;
        else
            F2pts(1,n)=nt2z-floor((cent2-F2(n,1))/del2);
        end
    end
        pick=F2pts

end

%% 3-5
% Plot the bright regions just made on the spectrum

%--arguments--
spec=pre_ft2d;
con_min=0.15; % ratio to the tallest peak
con_max=1.00;
ncon=15;
label=('pre-SIFT with masks');

%--commands--
if isempty(spec)==0
    figure()
    con=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec(25:nt2z-25,:)))));
    spec=real(squeeze(spec));
    mplp2d(squeeze(real(spec)),sw2,sw1,cent2,cent1,obs2,obs1,con);
    legend(label)
    xlabel('F1 (ppm)','fontsize',16)
    ylabel('F2 (ppm)','fontsize',16)
    hold on
    
    for n=1:length(sigf1(:,1))
        for m=1:size(sigf1,1)
            pick=sigf1(m,:);
            del1=(sw1/obs1)/(nt1z-1);
            if pick>nt1z/2
                pickppm=cent1-del1.*(nt1z-pick+1);
            else
                pickppm=cent1+del1.*(pick-1);
            end
            sigf1ppm(m,:)=pickppm;
            
            pick=sigf2(m,:);
            del2=(sw2/obs2)/(nt2z-1);
            if pick>nt2z/2
                pickppm=cent2-del2.*(nt2z-pick+1);
            else
                pickppm=cent2+del2.*(pick-1);
            end
            sigf2ppm(m,:)=pickppm;
        end
        
        temp=[sigf1ppm(n,1),sigf2ppm(n,1);
            sigf1ppm(n,1),sigf2ppm(n,2);
            sigf1ppm(n,2),sigf2ppm(n,2);
            sigf1ppm(n,2),sigf2ppm(n,1);
            sigf1ppm(n,1),sigf2ppm(n,1)];
        
        plot(temp(:,1),temp(:,2),'r-')
        
    end
    %~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~
end

%%
%\\\\\\\\\
% 4. SIFT
%\\\\\\\\\

%% 4-1
% Exploratory processing
%--arguments--
    in1=pre_c;
    in2=pre_s;
    ncyc=40;
%--commands--
    if isempty(in1)==0 && isempty(in2)==0
        sift=1;                 % =1 for SIFT, =0 for threshoulding.
        thre=10;                % threshould level, only used when sift=0.
        [j1,j2,s2n]=sift2dv1_2_peak(in1,in2,NUStab,ncyc,nt1z,nt2z,0,[],...
            sw1,sw2,ph1,ph2,wnd1,fctr1,sift,thre,...
            obs1,obs2,cent1,cent2,sigf1,sigf2,pick);
        clear j1 j2
    end
    if isempty(s2n)==0
        figure()
        subplot(3,1,1)
        slc1(s2n,[],[1])
        xlabel('Number of SIFT cycle','fontsize',16)
        title('nussft','fontsize',16)
        legend('S/N')

        subplot(3,1,2)
        slc1(s2n,[],[2])
        xlabel('Number of SIFT cycle','fontsize',16)
        legend('Max Peak')

        subplot(3,1,3)
        slc1(s2n,[],[3])
        xlabel('Number of SIFT cycle','fontsize',16)
        legend('Noise')
    end

%% 4-2
% Rerun with the optimal number of cycles

%--arguments--
    ncyc  = 10;         % # of SIFT cycle
%--commands--
    if isempty(in1)==0 && isempty(in2)==0
        sift=1;                 % =1 for SIFT, =0 for threshoulding.
        thre=10;                % threshould level, only used when sift=0.
        [j1,j2,s2n,sigf1,sigf2]=sift2dv1_2_peak(in1,in2,NUStab,ncyc,nt1z,nt2z,0,[],...
            sw1,sw2,ph1,ph2,wnd1,fctr1,sift,thre,...
            obs1,obs2,cent1,cent2,sigf1,sigf2,pick);
    end
    post_c=j1;
    post_s=j2;
    clear in1 in2 j1 j2

%% 
% Save all data in the Matlab working memory for later use.
    save ~/SIFT2DinMATLAB/example/workspace.mat
    
%% 
% To load saved set of data:
    load ~/SIFT2DinMATLAB/exmample/workspace.mat

%%
%\\\\\\\\\\\\\
% 5. Analysis
%\\\\\\\\\\\\\
% You can output the time domain data to process them in nmrPipe, or
% process directly in Matlab for a quick overview.

%% 5-1 
% Output time domain data to bruker format, for processing data with
% nmrPipe:

% pre-SIFT

%--arguments--
    out1=pre_c;
    out2=pre_s;
    path='~/SIFT2DinMATLAB/pre_reshuffled.bruker';
%--commands--
    if isempty(out1)==0 && isempty(out2)==0
%     scale=1e6/max(max(real(out1(25:nt2-25,:))));
    scale=1e6/max(max(real(out1)));
    in1=out1.*scale;
    in2=out2.*scale;
    brukout2dtimev1(in1,in2,nt1,nt2,path);
    end
%% 

% post-SIFT

%--arguments--
    out1=post_c;
    out2=post_s;
    path='~/SIFT2DinMATLAB/example/post.bruker';
%--commands--
    if isempty(out1)==0 && isempty(out2)==0
%     scale=1e6/max(max(real(out1(25:nt2-25,:))));
    scale=1e6/max(max(real(out1)));
    in1=out1.*scale;
    in2=out2.*scale;
    brukout2dtimev1(in1,in2,nt1,nt2,path);
    end
    
%% 5-2
% Process data in Matlab.

% pre-SIFT data:
    in1=pre_c;
    in2=pre_s;
    [pre_ft2d,junk,t1F2pre]=proc2d(in1,in2,nt1z,nt2z,sw1,sw2,ph1,ph2,wnd1,wnd2,fctr1);
    clear in1 in2 
%%
% post-SIFT data
    in1=post_c;
    in2=post_s;
    [post_ft2d,junk,t1F2post]=proc2d(in1,in2,nt1z,nt2z,sw1,sw2,ph1,ph2,wnd1,wnd2,fctr1);
    clear in1 in2 
%% 5-3
% Output spectra to Sparky format
    home='/Applications/Sparky.app/Contents/Resources/bin/matrix2ucsf';
    in1=pre_ft2d;
    in2=post_ft2d;
    ucsfout1= '~/SIFT2DinMATLAB/example/pre_ft2d';
    ucsfout2= '~/SIFT2DinMATLAB/example/post_ft2d';
    if isempty(in1)==0 && isempty(in2)==0
	binout2dfreq(in1,nt1z,nt2z,obs1,obs2,sw1,sw2,lowedge1,lowedge2,lab1,lab2,ucsfout1,home);
	binout2dfreq(in2,nt1z,nt2z,obs1,obs2,sw1,sw2,lowedge1,lowedge2,lab1,lab2,ucsfout2,home);
    end
    
%% 5-4
% Plot spectra in Matlab

% pre- and post-SIFT spectra side by side
    spec1=pre_ft2d;
    spec2=post_ft2d;
    con_min=0.10;
    con_max=1.00;
    ncon=15;
    if isempty(spec1)==0 && isempty(spec2)==0
    con1=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec1(25:nt2z-25,:)))));
    con2=logspace(log10(con_min),log10(con_max),ncon).*max(max(max(real(spec2(25:nt2z-25,:)))));
    figure()
    subplot(1,2,1)
    spec=real(squeeze(spec1));
    mplp2d(squeeze(real(spec)),sw2,sw1,cent2,cent1,obs2,obs1,con1);
    legend('pre-SIFT')
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F1 (ppm)','fontsize',16)
    subplot(1,2,2)
    spec=real(squeeze(spec2));
    mplp2d(squeeze(real(spec)),sw2,sw1,cent2,cent1,obs2,obs1,con2);
    legend('post-SIFT')
    xlabel('F2 (ppm)','fontsize',16)
    ylabel('F1 (ppm)','fontsize',16)
    end
%%
% 1D slices in the F1/F2 plane.

    F1=120.0;    % ppm
    F2=8.34;     % ppm
    if isempty(F1)==0 && isempty(F2)==0 
            %-- ppm to pts
            del1=(sw1/obs1)/(nt1z-1);
            if F1>=cent1
                F1pts=floor((F1-cent1)/del1)+1;
            else
                F1pts=nt1z-floor((cent1-F1)/del1);
            end
            F1pts
            %----
            del2=(sw2/obs2)/(nt2z-1);
            if F2>=cent2
                F2pts=floor((F2-cent2)/del2)+1;
            else
                F2pts=nt2z-floor((cent2-F2)/del2);
            end
            F2pts
    spec=real(spec1);
    figure()
    subplot(2,2,1)
    mplp2dslcv2(squeeze(spec),sw2,sw1,cent2,cent1,obs2,obs1,1,F2pts);
    title('pre-SIFT')
	xlabel('F2 (ppm)','fontsize',16)
    subplot(2,2,3)
    mplp2dslcv2(squeeze(spec),sw2,sw1,cent2,cent1,obs2,obs1,2,F1pts);
	xlabel('F1 (ppm)','fontsize',16)

    spec=real(spec2);
    subplot(2,2,2)
    mplp2dslcv2(squeeze(spec),sw2,sw1,cent2,cent1,obs2,obs1,1,F2pts);
    title('post-SIFT')
	xlabel('F2 (ppm)','fontsize',16)
    subplot(2,2,4)
    mplp2dslcv2(squeeze(spec),sw2,sw1,cent2,cent1,obs2,obs1,2,F1pts);
	xlabel('F1 (ppm)','fontsize',16)
    end
    
        % Change spectral limits in 1D slices
        subplot(2,2,1);
        xlim([lowedge1-sw1/obs1 lowedge1])
        ylim([-0.5e8 3e8])

        subplot(2,2,3);
        xlim([6 11])
        ylim([-0.5e8 3e8])

        subplot(2,2,2);
        xlim([lowedge1-sw1/obs1 lowedge1])
        ylim([-0.5e8 3e8])

        subplot(2,2,4);
        xlim([6 11])
        ylim([-0.5e8 3e8])

%%
% plot interferograms

spec_pre=t1F2pre;
spec_post=t1F2post;

figure()
plot(real(spec_pre(:,1)))
xlabel('F2 (Points)')

% In the spectrum above, chose F2 point to plot the corresponding t1
% interferogram.
itrfgmAtF2=162;
if isempty(itrfgmAtF2)==0
    figure()
    plot(1:nt1,real(spec_post(itrfgmAtF2,:)),1:nt1,real(spec_pre(itrfgmAtF2,:)))
    title('t1 interferograms')
    legend('post','pre')
end


% end of run-script.






    