function [SDA Alt]=calfSDA(phMap,varargin)
%
% Function calfSDA is an advanced version of the algorithm to identify period-2 (alternans) type of 
% instabilities in spatio-temporal data sets. 
% It analyses data traces from multiple locations, obtained under dynamic pacing - 
% successive increasing/decreasing pacing frequencies.
% It performs analysis and statistics on the identified alternans over space and time, 
% revealing evolution patterns as pacing frequency changes, keeping track of spatial concordance/discordance.
% 
% For simpler version - see function calcAltPer.
%
% Algorithm published in:
%
% IEEE Trans Biomed Eng. 2010 Feb;57(2):316-24.
% Detecting space-time alternating biological signals close to the bifurcation point.
% Jia Z, Bien H, Entcheva E.
% http://www.ncbi.nlm.nih.gov.ezproxy.hsclib.sunysb.edu/pubmed/19695992
% Contact:``emilia.entcheva@sunysb.edu
%    
% The program accepts pre-processed data, after beat detection and extraction of 
% a parameter of interest - peak height, duration etc.
%
% For the custom-developed pre-processing software please contact: 
% emilia.entcheva@sunysb.edu
%
% Input: 
%    'phMap' is essential. Both temporal persistence and phase will be
%            calculated using peak height (or other per-beat identified parameter).   
%
%    'master_act' is optional. The maps.master_act is used to 
%            calculate the frequencies which can be shown on the title of figures. 
%
%    'beat_range' [beat_start1 beat_end1 beat_start2 beat_end2...] is a 2*n-element 
%            array. It specifies the start beat and end beat number of each calculated 
%            temporal interval. By default, beat_start is 1, and beat_end is the 
%            length of phMap's third dimension.
%
%    'Tem_thres' a value from 0 to 1. It represents the temporal persistence  
%            of chosen period. By default, it is 0.60 (60%). 
%
%    'spat_thres' it represents the maxium distance (pixels) between the 
%            centroids of two spatially-discordant (SDA) regions. Recommended value is 10 for fine-SDA.
%
%    'scale' Specifies the number of unit distances represented by one
%            spatial sample, i.e. spatial resolution. Default is 0.088mm
%            per pixel (2 by 2 binning on top of 2x magnification) for our imaging system.
%
%    'clim'  Specifies the color axis limits for image scaling. Default
%            is [-15 15], typical for alternans ratio.
% 
%    'plot_flat' if its value is set to 1, visual results will be
%            generated. such as alternans and SDA map for each pacing frequency. 
%            The default valur is 0, no figures.
%
%    'block_check' The required input is atMap $ maps.t_start. the activation time 
%            map of all pixels. If there is a block, the activation time of the two
%            beats will be same, the just count how many zeros in the
%            derivative of atMap and divided by the total beats number, we
%            can get the percentage of blockage.
%
% Output: 
%    Alt.img, the two color image of alternans whose temporal persistence
%          is higher than the threshold.
%
%    Alt.num_po, the number of alternans regions with positive phase.
%
%    Alt.fre, the pacing frequencies, if master_act is given.
%
%    Alt.area_po, the number of alternans pixels with positive phase.
%
%    Alt.spa, the spatial percentage of alternans regions
%
%    Alt.RA, the area ratio for two phase alternans regions
%
%    Alt.NA, the number ratio for two phase alternans regions
%
%    Alt.size, the size of Alternans image
%
%    Alt.BlockMap, the image which shows the extent of blockage
%
%
%    SDA.img, the two color image of SDA whose temporal persistence
%          is higher than the threshold.
%
%    SDA.num_po, the number of SDA regions with positive phase.
%
%    SDA.fre, the pacing frequencies, if master_act is given. Same as
%             Alt.fre
%
%    SDA.area_po, the number of SDA pixels with positive phase.
%
%    SDA.spa, the spatial percentage of SDA regions
%
%    SDA.RA, the area ratio for two phase SDA regions
%
%    SDA.NA, the number ratio for two phase SDA regions
%
%    SDA.size, the size of SDA image
%
%    SDA.BlockMap, the image which shows the extent of blockage%
%
%
%
%  EXAMPLE: 
%        [SDA Alt]=calfSDA(maps.maxPeak,'master_act',maps.master_act,'beat_range',[1 29
%        31 59 62 88 94 119],'plot_flag',1,'clim',[-20 20],'block_check',maps.t_start)
%
% See also function [newimage]=calAltPer(phMap,beat,varargin);
%

% set up default parameters

% The switch to visualize data or not
    plot_flag=0;
% The spatial threshold to identify SDA regions
    spat_thres=6;
% The temporal threshold to identify alternans regions
    Tem_thres=0.6;
% The beats range
    beat_range=[1 size(phMap,3)];
% The limits for alternans ratio
    clim=[-15 15];
% The scale
    scale=[0.044 0.044];
% The activation map, the default is empty matrix
    atMap=zeros(0);
% The activation time of the master point
    master_act=ones(size(phMap,3));
    
narg=1;
while(narg<length(varargin))
    switch(lower(varargin{narg}))
        case 'Tem_thres'
            Tem_thres=varargin{narg+1};
            narg=narg+2;
        case 'spat_thres'
            spat_thres=varargin{narg+1};
            narg=narg+2;     
        case 'plot_flag'
            plot_flag=varargin{narg+1};
            narg=narg+2;     
        case 'beat_range'
            beat_range=varargin{narg+1};
            narg=narg+2;
        case 'clim'
            clim=varargin{narg+1};
            narg=narg+2;
        case 'scale'
            scale=varargin{narg+1};
            narg=narg+2;
        case 'master_act'
            master_act=varargin{narg+1};
            narg=narg+2;
        case 'block_check'
            atMap=varargin{narg+1};
            narg=narg+2;
            BlockMap=zeros(size(phMap,1),size(phMap,2));
        otherwise
            warning('calfSDA:arg', 'Unknown argument %s', varargin{narg});
            narg=narg+2;
    end
end



% Validate input parameters
if mod(size(beat_range,2),2)==1;
   errordlg('The input beat numbers should be even','Input Error','modal');
   return;
end

% pre allocate phMap_in to accerlerate computation
phMap_in=zeros(size(phMap,1),size(phMap,2));
phMap_ma=zeros(size(phMap,1),size(phMap,2));

%calculate for different frequency interval
for q=1:size(beat_range,2)/2
    beats=beat_range(q*2-1):beat_range(q*2);
    %phase beat
    phase_beat=round((beats(1)+beats(end))/2)-beats(1);
    %derivative
    phMap_der=diff(phMap(:,:,beats),1,3);
    %get the sign of phMap_n_de
    phMap_s=sign(phMap_der);
    %Suppose one positive and one negative is one alternate, calculate the Num
    phMap_Alt=diff(phMap_s,1,3);
    %to get rid of all other patterns except alternans
    phMap_Alt=(phMap_Alt==2)|(phMap_Alt==-2);

    for i=1:size(phMap_Alt,1)
        for j=1:size(phMap_Alt,2)
            % try to find block areas where the the diff of activation time
            % is zero
            if(~isempty(atMap))
               BlockMap(i,j)=size((find(diff(atMap(i,j,beats),1,3)==0)),1)/size(beats,2)*100;
            end
            phMap_in(i,j)=max(diff([0;find(phMap_Alt(i,j,:)~=1);size(phMap_Alt,3)+1]));
            % Get the start and end of max continous alternans with index of
            % zeros from diff([1;find(phMap_Alt(i,j,:)==0);size(phMap_Alt,3)])
            break_start=find(diff([0;find(phMap_Alt(i,j,:)~=1);size(phMap_Alt,3)+1])==phMap_in(i,j),1);   
            break_end=break_start+1;
            % the index for continous alternans
            tem_array=[1;find(phMap_Alt(i,j,:)~=1);size(phMap_Alt,3)];
            alt_in=tem_array(break_start):tem_array(break_end)+1;
            % calculate the mean magnitude of alternans
            phMap_ma(i,j)=100*mean(abs(phMap_der(i,j,alt_in))./max(phMap(i,j,alt_in),phMap(i,j,alt_in+1)),3);
        end
    end

    %final result of temporal persistence
    Num_Alt_Per=phMap_in./(size(phMap_Alt,3)+1);

    

    % Save the original alternans ratio image
    Alt.TP(:,:,q)=Num_Alt_Per;
    
    % Save the alternans ratio image
    Alt.AR(:,:,q)=phMap_s(:,:,phase_beat);
    
    % To cut use grayimg
    newimage=Num_Alt_Per>Tem_thres;

    %calculate phase biimage
    phase_biimage=newimage.*phMap_s(:,:,phase_beat);

    % Add two phase alternans image to output
    Alt.img(:,:,q)=phase_biimage.*phMap_ma;

        
    % Add block map
if(~isempty(atMap))
    Alt.BlockMap(:,:,q)=BlockMap;
end
    %  calculate the properties of positive phase image
    L = bwlabel(phase_biimage>0);
    Prop_po= regionprops(L, 'PixelList');
    Alt.num_po(q)=size(Prop_po,1);
    Alt.area_po(q)=sum(sum(phase_biimage>0,2),1);

    %  calculate the properties of positive phase image
    L = bwlabel(phase_biimage<0);
    Prop_ne= regionprops(L,'PixelList');
    Alt.num_ne(q)=size(Prop_ne,1);
    Alt.area_ne(q)=sum(sum(phase_biimage<0,2),1);

    % plot the alternans cluster-size profile
    L = bwlabel(newimage);
    s= regionprops(L, 'area' );
    a=cat(1,s.Area);
    N=hist(a,max(a));
if(~isempty(N))
    Alt.hist(q,1:size(N,2))=squeeze(N);
end
    % Calculate the average pacing frequency
    Alt.fre(q)=mean(200/mean(diff(master_act(beats))));
    % Calculate the CV
    if(~isempty(atMap))
        for CV_row=3:(size(atMap,1)-2);
           for CV_beat=beats(3):beats(end-2)
            [P,S] = polyfit(1:size(atMap,2),atMap(CV_row,:,CV_beat),1);
            CV(CV_row,CV_beat)=1./P(1)*scale(1)*200;
           %CV(CV_row)=P(1)*scale(1)*200;

           end;
        end
        Alt.CV(q)=nanmean(CV(:));
    end
%%    
    % SDA regions identification 
    % The idea is for every alternans-exhibiting pixel, to search in its
    % neigbourhood (a square) to check if there is an opposite phase pixel. If
    % such is found, the whole area wiil be considered SDA area. In order to
    % avoid the situation that the positive and negative phase regions are
    % directly connected, so SDA regions are identified separately.

    SDAmap=zeros(size(newimage));
    %Search for every area in positive phase areas
    for i=1:size(Prop_po,1)
            % Search every pixel in one area
            for n=1:size(Prop_po(i,1).PixelList,1)
                % Get the coordinates of  
                C_x=Prop_po(i,1).PixelList(n,1);
                C_y=Prop_po(i,1).PixelList(n,2);   
                % Difine the search region
                Search_x=max(1,(C_x-spat_thres)):min(size(newimage,2),(C_x+spat_thres));
                Search_y=max(1,(C_y-spat_thres/2)):min(size(newimage,1),(C_y+spat_thres/2));
                % To check if there are out of phase areas
                if(~isempty(find(phase_biimage(Search_y,Search_x)==-1,1)))
                    % if there are opposite phase region, set the whole area as
                    % SDA region
                    for m=1:size(Prop_po(i,1).PixelList,1)    
                        SDAmap(Prop_po(i,1).PixelList(m,2),Prop_po(i,1).PixelList(m,1))=1;      
                    end

                        break
                end
           end
    end

    %Search for every area in negative phase areas
    for i=1:size(Prop_ne,1)
            for n=1:size(Prop_ne(i,1).PixelList,1)
                C_x=Prop_ne(i,1).PixelList(n,1);
                C_y=Prop_ne(i,1).PixelList(n,2);     
                Search_x=max(1,(C_x-spat_thres)):min(size(newimage,2),(C_x+spat_thres));
                Search_y=max(1,(C_y-spat_thres/2)):min(size(newimage,1),(C_y+spat_thres/2));
                if(~isempty(find(phase_biimage(Search_y,Search_x)==1,1)))
                    for m=1:size(Prop_ne(i,1).PixelList,1)    
                        SDAmap(Prop_ne(i,1).PixelList(m,2),Prop_ne(i,1).PixelList(m,1))=-1;      
                    end

                        break
                end
           end
    end

    % Add two phase SDA image to output
    SDA.img(:,:,q)=SDAmap.*phMap_ma;

    % Add block map
    if(~isempty(atMap))
        SDA.BlockMap(:,:,q)=BlockMap;
    end
        
    
    %  calculate the properties of positive phase in SDA image
    L = bwlabel(SDAmap>0);
    Prop_po= regionprops(L, 'PixelList');
    SDA.num_po(q)=size(Prop_po,1);
    SDA.area_po(q)=sum(sum(SDAmap>0,2),1);

    %  calculate the properties of positive phase in SDA image
    L = bwlabel(SDAmap<0);
    Prop_ne= regionprops(L,'PixelList');
    SDA.num_ne(q)=size(Prop_ne,1);
    SDA.area_ne(q)=sum(sum(SDAmap<0,2),1);


        
    % calculate cluster-size profile
        L = bwlabel((SDAmap==-1)|(SDAmap==1));
        s= regionprops(L, 'area' );
        a=cat(1,s.Area);
        N=hist(a,max(a));
        if isempty(N)
             SDA.hist(q,1)=0;
        else
            SDA.hist(q,1:size(N,2))=squeeze(N);
        end
    
    %Calculate the average pacing frequency
    SDA.fre(q)=mean(200/mean(diff(master_act(beats))));
end

%% Calculate the ratios for all pacing frequecies
    % just plug the number of size of whole image for later compute the spatial
    % percentage of alternans or SDA
       Alt.ImSize=size(newimage);
       SDA.ImSize=size(newimage);

    % calculate the spatial percentage of alternans and SDA regions
        Alt.spa=100*(Alt.area_po+Alt.area_ne)./(Alt.ImSize(1)*Alt.ImSize(2));
        SDA.spa=100*(SDA.area_po+SDA.area_ne)./(SDA.ImSize(1)*SDA.ImSize(2));

    % calculate the spatila ratio for the two phases in SDA and alternans
    % regions
        RA=Alt.area_po./Alt.area_ne;
        Alt.RA=(RA<=1).*RA+(1./RA<1).*(1./RA);
        RA=SDA.area_po./SDA.area_ne;
        SDA.RA=(RA<=1).*RA+(1./RA<1).*(1./RA);   

    % calculate the number ratio for the two phases in SDA and alternans
    % regions
        RN=Alt.num_po./Alt.num_ne;
        Alt.RN=(RN<=1).*RN+(1./RN<1).*(1./RN);
        RN=SDA.num_po./SDA.num_ne;
        SDA.RN=(RN<=1).*RN+(1./RN<1).*(1./RN);  
%% Data visualization
if plot_flag==1
    for f=1:size(Alt.img,3)
        figure;

if(~isempty(atMap))  
        %plot alternan image
        subplot(3,1,1);
        imagesc((1:Alt.ImSize(2)).*scale(2),(1:Alt.ImSize(1)).*scale(1),Alt.img(:,:,f),clim);
        axis image
        colorbar
        set(gca,'FontSize',14);
        title(sprintf('%.2f Hz Alternans',Alt.fre(f)),'FontSize',16,'HorizontalAlignment','center','FontWeight','bold');
        
        %plot SDA image
        subplot(3,1,2);
        imagesc((1:SDA.ImSize(2)).*scale(2),(1:SDA.ImSize(1)).*scale(1),SDA.img(:,:,f),clim);
        axis image
        colorbar
        set(gca,'FontSize',14);
        title(sprintf('%.2f Hz SDA',SDA.fre(f)),'FontSize',16,'HorizontalAlignment','center','FontWeight','bold');
        
        % plot the image which can represent the blockage
        subplot(3,1,3);
        imagesc((1:SDA.ImSize(2)).*scale(2),(1:SDA.ImSize(1)).*scale(1),Alt.BlockMap(:,:,f),[0 100]);
        axis image
        colorbar
        set(gca,'FontSize',14);
        title(sprintf('%.2f Hz Blockage(%%)',SDA.fre(f)),'FontSize',16,'HorizontalAlignment','center','FontWeight','bold');
% if atMap is unavailable, it is no way to calculate the blockage, so
% simply only display the alternans and SDA images
else
        subplot(2,1,1);
        imagesc((1:Alt.ImSize(2)).*scale(2),(1:Alt.ImSize(1)).*scale(1),Alt.img(:,:,f),clim);
        axis image
        colorbar
        set(gca,'FontSize',14);
        title(sprintf('%.2f Hz Alternans',Alt.fre(f)),'FontSize',16,'HorizontalAlignment','center','FontWeight','bold');
        subplot(2,1,2);
        imagesc((1:SDA.ImSize(2)).*scale(2),(1:SDA.ImSize(1)).*scale(1),SDA.img(:,:,f),clim);
        axis image
        colorbar
        set(gca,'FontSize',14);
        title(sprintf('%.2f Hz SDA',SDA.fre(f)),'FontSize',16,'HorizontalAlignment','center','FontWeight','bold');
        
end
    end

    % plot the bar figure of different frequencies
    
    %SDA spatial areas
    figure;
    bar(SDA.fre,100*[SDA.area_po(:),SDA.area_ne(:)]./(size(newimage,1)*size(newimage,2)),'stacked')
    set(gca,'FontSize',14)
    Xlabel('Pacing Frequency (Hz)','FontSize',14)
    Ylabel('Spatial Percentage (%)','FontSize',14)
    legend('phase+','phase-','Location','NorthWest');
    title('SDA','FontSize',14);
    text(SDA.fre,100*[SDA.area_po(:)+SDA.area_ne(:)]./(size(newimage,1)*size(newimage,2)),num2str(SDA.RA','%3.2f'),'FontSize',14,'HorizontalAlignment','center','VerticalAlignment','bottom');
    
    
    %Alternans spatial areas
    figure;
    bar(Alt.fre,100*[Alt.area_po(:),Alt.area_ne(:)]./(size(newimage,1)*size(newimage,2)),'stacked')
    set(gca,'FontSize',14)
    Xlabel('Pacing Frequency (Hz)','FontSize',14)
    Ylabel('Spatial Percentage (%)','FontSize',14)
    legend('phase+','phase-','Location','NorthWest');
    title('Alternans','FontSize',14);
    text(Alt.fre,100*[Alt.area_po(:)+Alt.area_ne(:)]./(size(newimage,1)*size(newimage,2)),num2str(Alt.RA','%3.2f'),'FontSize',14,'HorizontalAlignment','center','VerticalAlignment','bottom');
    
    %SDA num of areas
    figure;
    bar(SDA.fre,[SDA.num_po(:),SDA.num_ne(:)],'stacked')
    set(gca,'FontSize',14)
    Xlabel('Pacing Frequency (Hz)','FontSize',14)
    Ylabel('Number of regions','FontSize',14)
    legend('phase+','phase-','Location','NorthWest');
    title('SDA','FontSize',14);
    text(SDA.fre,SDA.num_po(:)+SDA.num_ne(:),num2str(SDA.RN','%3.2f'),'FontSize',14,'HorizontalAlignment','center','VerticalAlignment','bottom');

    %Alternans num of areas
    figure;
    bar(Alt.fre,[Alt.num_po(:),Alt.num_ne(:)],'stacked')
    set(gca,'FontSize',14)
    Xlabel('Pacing Frequency (Hz)','FontSize',14)
    Ylabel('Number of regions','FontSize',14)
    legend('phase+','phase-','Location','NorthWest');
    title('Alternans','FontSize',14);
    text(Alt.fre,Alt.num_po(:)+Alt.num_ne(:),num2str(Alt.RN','%3.2f'),'FontSize',14,'HorizontalAlignment','center','VerticalAlignment','bottom');

end