% QT_gradient_detection.m - Find QT intervals based on the gradient rules
%        Copyright (C) 2006  Vaclav Chudacek
%        This software is released under the terms of the GNU General
%        Public License (http://www.gnu.org/copyleft/gpl.html).
     
function QT_detection_using_Gradient_Rules_CinC_2006(path,plot_enabled)

if (nargin~=2),
 disp(sprintf('This function requires two arguments. Setting to default'));
 path = 'C:\!Vasek\Work\!CINC 2006\!Data\ptbdb';
 plot_enabled = 0;
end

recordFile = 'RECORDS';           % Record names
load 'CINC_2006_Reference_results.txt';     % Official results for comparison  as of 07.09.2006
CINC_2006_results = zeros(549,2); % Our results
CINC_2006_QT_lengths = zeros(1,549);
CINC_2006_differences = zeros(1,549);     % Difference of our results to the official
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
recordFileFull = fullfile(path, recordFile);
fid_record = fopen(recordFileFull,'r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin_at = 119;
% for i = 1:begin_at,
%     fileName = fgetl(fid_record);
% end
for counter = 1:549,
    fileName = fgetl(fid_record);
    headerFile = fullfile(path,[fileName '.hea']);
%     vcgFile = fullfile(path,[fileName '.xyz']);
    ecgFile = fullfile(path,[fileName '.dat']);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %------     LOAD HEADER     -------------------------------------------
    fprintf(strcat('Loading new record #',num2str(counter),': ',fileName));
    fprintf('\n');
    fid_header = fopen(headerFile,'r');
    z = fgetl(fid_header);
    fclose(fid_header);
    %------     LOAD DATA      --------------------------------------------
    fid_ecg = fopen(ecgFile,'r');
    ECG = fread(fid_ecg, [12, inf], 'int16');
    fclose(fid_ecg);
%     fid_vcg = fopen(vcgFile,'r');
%     VCG = fread(fid_vcg, [3, inf], 'int16');
%     fclose(fid_vcg);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%% Very, very basic filtration - just to get rid of the worst...
    for i = 1:12
        herz = 1000/length(ECG(i,:));
        A = fft(ECG(i,:));
        A(1:round(0.66/herz)) = 0;
        A(end-round(0.66/herz):end)=0;
        A((round(50/herz)-5):(round(50/herz)+5)) = 0;
        A(end-(round(50/herz)+5):end-(round(50/herz)-5)) = 0;
        ECGf(i,:) = (real(ifft(A)));
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%% R-peak detection based on PT_Hamilton algorithm %%%%%%%
    % Source code distributed under GNU license can be found on
    % http://www.eplimited.com/
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%% START %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [nbr_heart_beats, r_peak_sample_list] = PanTomkinsHamilton(ECGf(5,:),1000);
    [nbr_heart_beatsI, r_peak_sample_listI] = PanTomkinsHamilton(ECGf(1,:),1000);
    [nbr_heart_beatsV4, r_peak_sample_listV4] = PanTomkinsHamilton(ECGf(10,:),1000);
    r_peak_sample_list = r_peak_sample_list*5;
    posR = r_peak_sample_list;
    posRI = r_peak_sample_listI;
    posRV4 = r_peak_sample_listV4;

    RR_dist = zeros(1,length(min(100,length(posR)-1)));
    selected_beat = 5; next_beat = 6;
    for i = 1:min(100,length(posR)-1),
        RR_dist(i) = (posR(i+1) - posR(i));
    end

    %%% Correction of detected peaks - R-T mess, wrong detection
    if ((RR_dist(selected_beat))>(1.2*RR_dist(selected_beat-1))||(RR_dist(selected_beat))<(0.8*RR_dist(selected_beat-1)))||...
            ((RR_dist(selected_beat))>(1.2*median(RR_dist))||(RR_dist(selected_beat))), %%% added because of file #543
        if ((RR_dist(selected_beat))>=1500), %%% added because of files #543, #315, #415
            selected_beat = selected_beat - 1; next_beat = next_beat +1;
        elseif (RR_dist(selected_beat+2))<=(1.2*RR_dist(selected_beat))&&...
                (RR_dist(selected_beat+2))>=(0.8*RR_dist(selected_beat))&&...
                (RR_dist(selected_beat))<425,
            if RR_dist(selected_beat)<320,
                if RR_dist(selected_beat)<500
                    next_beat = next_beat+1;
                else
                    selected_beat = selected_beat+1; next_beat = next_beat+2;
                end
            else
                selected_beat = selected_beat-1;
            end
        elseif (RR_dist(selected_beat+1))<=(1.2*RR_dist(selected_beat))&&(RR_dist(selected_beat+1))>=(0.8*RR_dist(selected_beat))
            selected_beat = selected_beat+1; next_beat = next_beat +1;
        end
    end
    posRpeak_out = posR(selected_beat);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%% R-peak detection by PanTomkins&Hamilton algortihm END %%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%% Q-on computation START %%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Q_window_span_left = 150; Q_window_span_right = 50;
    R_window_span_left = 70; R_window_span_right = 10;
    T_window_span_right_1 = 200; T_window_span_right_2 = 200;
    Q_comps = 0; mistakeQ = 0; T_comps = 0; mistakeT = 0;
    Q_on_abs = zeros(1,12);

    for i = 1:12,
        %%% Correction of R-peak for each lead
        [temp correctionR_temp] = max((ECGf(i,posR(selected_beat)-R_window_span_left:posR(selected_beat)+R_window_span_right)));
        correctionR = correctionR_temp - R_window_span_left;
        posR_adjusted = posR(selected_beat) + correctionR;

        %%% 1st step of Q-on computation - see function: FIND_Q_ON
        try
            Q_on_temp = find_Q_on(ECGf(i,posR_adjusted - Q_window_span_left:posR_adjusted + Q_window_span_right));
        catch
            disp(strcat(['!!!!! MISTAKE of: Q-on computation on lead:', (num2str(i)),'!!!!!']));
            mistakeQ = 1;
            Q_on_temp = 0;
        end

        if mistakeQ == 0,
            Q_comps = Q_comps+1;
            Q_on_rel(Q_comps) = Q_on_temp + correctionR + 1;
            Q_on_abs(i) = posR(selected_beat) - Q_on_rel(Q_comps);
        else
            mistakeQ = 0;
        end
    end

    % All possible Qons are sorted and the smallest are deleted
    if min(Q_on_rel)<10,
        Q_on_rel_tmp = sort(Q_on_rel,'ascend');
        clear('Q_on_rel');
        Q_on_rel = Q_on_rel_tmp(max(find(Q_on_rel_tmp<10))+1:end);
    else
        Q_on_rel = sort(Q_on_rel,'ascend');
    end

    % 2nd step of Q-on computations - see function PRESCISE_Q
    for i = 2:2
        Q_on_final_temp(i,:) = precise_Q_on(ECGf(i,posR(selected_beat) - 150 : posR(selected_beat) + Q_on_rel(end)+25 -150),Q_on_rel);
    end

    % Maximal score is computed to select the best fitting Q-on
    maxScoreQ = max(sum(Q_on_final_temp(:,:)));
    if maxScoreQ == 0,
        Q_on_rel_final = median(Q_on_rel);
    else
        maxScoreSelected = find(sum(Q_on_final_temp(:,:))>=maxScoreQ);
        Q_on_rel_final = Q_on_rel(maxScoreSelected(1));
    end
    if (median(Q_on_rel) - Q_on_rel_final) < -25,
        Q_on_rel_final = round(median(Q_on_rel));
    end

    Q_on_rel_final = Q_on_rel_final - 150;
    Q_on_abs_final = Q_on_rel_final + posR(selected_beat);

    % 3rd step of Q-on computations - see function: CORRECT_Q
    correctionQ = [];
    leads = [1 2 8 9];
    for i = 1:length(leads);
        correctionQ =  [correctionQ correct_Q_on(ECGf(leads(i),:),Q_on_abs_final,Q_on_abs_final-20,Q_on_abs_final+30)];
    end
    correctionQ = median(correctionQ);

    Q_on_rel_final = Q_on_rel_final + correctionQ;
    Q_on_abs_final = Q_on_abs_final + correctionQ;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%% Q-on computation END %%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    R_pos_in_summedECG = 75; nextR_pos_summedECG_distance = 200;

    summedECG = zeros(1,751);
     for i = 1:12,
        % Computation of summed and filtered ECG
        if (i == 8)||(i == 9),
            if(i == 8),
                summedECG = slidingavg(slidingavg(slidingavg(ECGf(12,(posR(selected_beat)-R_pos_in_summedECG):(posR(next_beat)-nextR_pos_summedECG_distance)),50),10),3);
            else
                summedECG = summedECG + slidingavg(slidingavg(slidingavg(ECGf(12,(posR(selected_beat)-R_pos_in_summedECG):(posR(next_beat)-nextR_pos_summedECG_distance)),50),10),3);
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Inverting signal when negative
    if abs(min(summedECG))>max(summedECG),
        summedECG = -1*summedECG;
    end
    %Correction of R-peak on summedECG in +-20 sample radius
    [unused correctionR] = max(summedECG(R_pos_in_summedECG-20:R_pos_in_summedECG+20));
    R_pos_in_summedECG = R_pos_in_summedECG - (20-correctionR+1);
    Q_on_for_summedECG = R_pos_in_summedECG - Q_on_rel_final;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%% T-off computation START %%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    T_error = 0;
    %Searching for first and second maxima/minima of the T-wave
    [unused maxT] = max(abs(summedECG(R_pos_in_summedECG+min(150,round(RR_dist(selected_beat)/2)):(min(length(summedECG(R_pos_in_summedECG:end)),400+R_pos_in_summedECG))))-summedECG(Q_on_for_summedECG));
    if ~exist('maxT','var'),
        maxT = 0;
    end
    if summedECG(R_pos_in_summedECG+200+maxT)<0,
        zaporne_T = 1;
        maxT = maxT + R_pos_in_summedECG + 200;
        [velikostT2 maxT2] = max(summedECG(maxT:end - min(length(summedECG)-maxT,100))-summedECG(maxT));
        maxT2 = maxT2 -1;
        if ((velikostT2<(summedECG(end - min(length(summedECG)-maxT,100))-summedECG(maxT))+100)&&(velikostT2>(summedECG(end - min(length(summedECG)-maxT,100))-summedECG(maxT))-100)),
            maxT2=0;
        end
    else
        maxT2 = 0;
        maxT = maxT + R_pos_in_summedECG + 200;
        zaporne_T =0;
    end
    if maxT>length(summedECG)-50,
        maxT = length(summedECG)-200;
    end
    maxT2_rel = maxT2 - maxT;

    %Signal used for T-off detection is shifted by value of Q-on 
    Tecg = summedECG(maxT:end-50)-summedECG(Q_on_for_summedECG);

    %%% 1st step of T-off computation - see function: FIND_T_OFF
    [Toff_temp, T_error] = find_T_off(Tecg,maxT, maxT2, maxT2_rel,zaporne_T);
    
    %   posTmax_out = maxT - R_pos_in_summedECG + posR(selected_beat);
    T_off_rel_final = Toff_temp + maxT;
    T_off_abs_final = Toff_temp + maxT - R_pos_in_summedECG + posR(selected_beat);

    % 2nd step of T-off computations - see function: CORRECT_T_OFF
    correctionT = [];
    leads = [1 2 8 9 10];
    for i = 1:length(leads);
        correctionT(i) =  correct_T_off(ECGf(leads(i),T_off_abs_final-200:T_off_abs_final));
    end

    correctionTfinal = min(correctionT);
    T_off_abs_final = T_off_abs_final - correctionTfinal;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%% T-off computation START %%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Plotting of the results START 
    if plot_enabled,
        %%% Draws detected R(selected_beat):R(next_beat)
        % figure(counter*100);plot(ECGf(1,posR(selected_beat):posR(next_beat)));
        % Draws preliminary and corrected results
        figure(counter);
        plot(ECGf(2,Q_on_abs_final -100:T_off_abs_final+100),'r');
        hold on;
        plot(ECGf(1,Q_on_abs_final -100:T_off_abs_final+100));
        plot(ECGf(8,Q_on_abs_final -100:T_off_abs_final+100));
        plot(ECGf(9,Q_on_abs_final -100:T_off_abs_final+100));
        plot(100,-500:500,'k');
        plot(100 + correctionQ, -500:500,'g');
        %%%%%%%%%%%
        plot(length(ECGf(i,Q_on_abs_final -100:T_off_abs_final+100))-100+correctionTfinal,-500:500,'k');
        plot(length(ECGf(i,Q_on_abs_final -100:T_off_abs_final+100))-100 - correctionTfinal,-500:500,'g');
    end
    
    % Saving final results
    CINC_2006_results(counter,:) = [Q_on_abs_final,T_off_abs_final];
    CINC_2006_QT_lengths(counter) = (T_off_abs_final-Q_on_abs_final);
    CINC_2006_differences(counter) = abs(CINC_2006_Reference_results(counter)-CINC_2006_QT_lengths(counter));

    clear ECG*; clear Q*; clear R*; clear T*; clear cor*; clear m*; clear n*; clear pos*; clear r*; clear z*; clear s*; clear t*; clear u*;
end
fclose(fid_record);
% Final score computation based on final reference QT measurements (not included) 
omitted_results = find((CINC_2006_QT_lengths > 500) | (CINC_2006_QT_lengths < 250));
if length(omitted_results)>27
    disp('Sorry you are omitting to many results. More than 27.')
end
square = (CINC_2006_differences.^2);
square(omitted_results) = -1;

RESULT = sqrt(mean(square(find(square ~= -1))));
RESULT = RESULT*((549-length(omitted_results))/549)
save('QTresults','CINC_2006_results')