function [newimage]=calAltPer(phMap,beat,varargin)
%
% Calculate the temporal persistence of period-2 (alternans) instabilities in
% spaiotemporal data sets contaminated with noise, and return spatial maps of 
% their amplitude and phase
%
% 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
% emilia.entcheva@sunysb.edu
%    
% This program is designed to quantify the temporal persistence of
% alternans in traces of biological signals, exhibiting such instabilities 
% (ECG, action potentials, calcium transients etc), recorded from multiple spatial locations. 
% 
% 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 parameter Peak height Map (phMap) is the main data source for indentifying 
% alternans patterns. It is a 3D matrix containg the peak height 
% (or parameter of interest for each beat) for all events (transients) over time and 
% at each spatial location. 
%
% Input parameter 'beat' in format [beat_start beat_end] is the beat range for which
% one desires to apply the algorithm.
%
% Usage:
%
%   [Num_Alt_Per]=calAltPer(phMap,[beat_start beat_end],'param1', 'value1', ...)
%   sets optional parameters in ('param', value) pairs. Valid parameters
%   are:
%        'threshold' The output will only show the pixels whose alternans
%                   percentage is higher than a set threshold value for persistence. 
%                   Default value is .60 (60%)
%        'bicolor'  to show the image as two colors, the pixels whose
%                   alternans percentage is higher than threshold have one
%                   colors and all other pixels have another color. 
%                   Default value is 1.
%        'Alt_ratio'to show the extent of alternans. 
%                   It is represented by alternans ratio (AR): 
%		    the difference in magnitude of two neighbouring beats divided
%                   by the larger magnitude of two beats. 
%		    The default is 0; set to 1 if you want a map of AR to be outputed.  
%   Output:
%        'newimag'  is a 2D matrix, which can be binary image to show
%                   alternans area if 'bicolor' is set to 1. Or TP
%                   values if 'bicolor' is set to 0. Also can be alternans
%                   ratio if 'Alt_ratio' is set to 1.
%                            
%  Example:  [newimage]=calAltPer(phMap,[1 30],'threshold',0.6,'bicolor',1);


% Check for input beats values
if(nargin<2)
   beat_start=1;
   beat_end=size(phMap,3);
else
    beat_start=beat(1);
   beat_end=beat(2);
end

% set default values for optional parameters
% Temporal persistence threshold 
Threshold=0.60;
bicolor=1;
Alt_ratio=0;
atMap=[];

% get input optional parameters
narg=1;
while(narg<length(varargin))
    switch(varargin{narg})
        case 'threshold'
            Threshold=varargin{narg+1};
            narg=narg+2;
        case 'bicolor'
            bicolor=varargin{narg+1};
            narg=narg+2;
        case 'Alt_ratio'
            Alt_ratio=varargin{narg+1};
            narg=narg+2;
        case 'atMap'
            atMap=varargin{narg+1};
            narg=narg+2;
        otherwise
            warning('calAltPer:arg', 'Unknown argument %s', varargin{narg});
            narg=narg+2;
    end
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));
BlockMap=zeros(size(phMap,1),size(phMap,2));

%Calculated beats
beats=beat_start:beat_end;

%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,:)==0);size(phMap_Alt,3)]));
            % 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,:)==0);size(phMap_Alt,3)])==phMap_in(i,j),1);   
            break_end=break_start+1;
            % the index for continous alternans
            tem_array=[1;find(phMap_Alt(i,j,:)==0);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))./phMap(i,j,alt_in),3);
end
end

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

% To cut use grayimg
newimage=Num_Alt_Per>Threshold;
if ~bicolor
    newimage=newimage.*Num_Alt_Per;
end

if Alt_ratio
    newimage=newimage.*phMap_ma;
end