clear
clc
close all
warning('OFF');


% Description: calculate origin point from XYZ lead and generate image of XYZ 
% median beat based on the modified baseline.
% SYNTAX : origin_point_algorithm(file_name,source_folder,image_folder,destination_folder)
% INPUT  : file_name          - file name for the matfile containing XYZ 
%                               median beat detected R peak indexed location.
%          source_folder      - the folder location where the file is placed.
%          image_folder       - the destination folder to place generated image.
%          destination_folder - where the matfile will be placed after 
%                               calculation is done.

% Origin Point Calculation Code V.1.1
% Authors:
%       Erick Andres Perez Alday, PhD, <perezald@ohsu.edu>
%       Annabel Li-Pershing, BS, <annabel2@pdx.edu>
% Principle Investigator:
%       Larisa Tereshchenko, MD, PhD, < tereshch@ohsu.edu > 
% Last update: August 2nd, 2018



% sampling frequency
prompt='What is the sampling frequency in Hz';
dlg_title='Sample Frequency';
fs_d = inputdlg(prompt,dlg_title,1);
fs=str2num(fs_d{1});


% load mat files
[file_name,path_name] = uigetfile('*','Select file for Median beat');

%% file import
file_ID          = strsplit(file_name,'.');
file_path        = fullfile(path_name,file_name);

A= load(file_path);

%% Create Results folder
Results_folder = [path_name, '/', 'Results', '/'];
if (exist(Results_folder,'file') ~= 7)
    mkdir (Results_folder);
end



    %% Load data from matfile
%   A=load(strcat(source_folder,'/',file_name));
    XYZ_M=A.XYZ_M_T;
    Rx=A.Rpeaks(:,1);	
    disp(file_name)
    
    %% Calculate various windows analysis
    w=floor((mean(diff(Rx))*.8)-(372.69*(fs/1000)));

    %% Is the RR interval < 600ms (0.6*fs)
    if mean(diff(Rx))<(600*(fs/1000)) && (A.Label_beat~='S')
	cent=(600-260)*(fs/1000);
    else
	cent=(600-320)*(fs/1000);
    end


    w1_start=cent-w;
    w1_end=cent+w;

    w2_start=cent-(160*(fs/1000));
    w2_end=cent+(160*(fs/1000));

    if w1_start < w2_start
       w1_start=w2_start;
    end

    if w1_end > w2_end
        w1_end=w2_end;
    end

    %% cut xyz median beat based on the window
    xyz_w1 = XYZ_M(w1_start:w1_end,:);
    xyz_w2 = XYZ_M(w2_start:w2_end,:);

    %% calculate origin point with clustering algorithm ori_clustering
    [cluster_pt,~]=ori_clustering(xyz_w1);

    if isnan(cluster_pt)
        [cluster_pt,~]=ori_clustering(xyz_w2);
        Ori_pt(1) = cluster_pt+w2_start;
    else 
        Ori_pt(1) = cluster_pt+w1_start;  % convert the calculated index to real time window
    end


    % calculate the first origin point and the modified vector magnitude
    for ii=1:length(XYZ_M(:,1))
        VecMag_M1(ii)=norm([XYZ_M(Ori_pt(1),1)-XYZ_M(ii,1),XYZ_M(Ori_pt(1),2)-XYZ_M(ii,2),XYZ_M(Ori_pt(1),3)-XYZ_M(ii,3)]);
    end


    % calculate the sum of the absolute gradient on XYZ within w1
    [abs_grad_ori,~]=abs_grad_fun(xyz_w1);

    % calcaulate the second origin point
    Ori_pt2(1)=abs_grad_ori+w1_start;

    % calculate the second origin point and the modified vector magnitude
    VecMag_M2 = nan(1,length(XYZ_M(:,1)));
    for ii=1:length(XYZ_M(:,1))
        VecMag_M2(ii)=norm([XYZ_M(Ori_pt2(1),1)-XYZ_M(ii,1),XYZ_M(Ori_pt2(1),2)-XYZ_M(ii,2),XYZ_M(Ori_pt2(1),3)-XYZ_M(ii,3)]);
    end

    % calculate the area under the curve of the given angle
    Area_m1=trapz(VecMag_M1(w1_start:w1_end));
    Area_m2=trapz(VecMag_M2(w1_start:w1_end));

    % modify the XYZ median beat with the two origon point calculated
    XYZ_M1 = repmat(XYZ_M(Ori_pt(1),:),size(XYZ_M,1),1) - XYZ_M;
    XYZ_M2 = repmat(XYZ_M(Ori_pt2(1),:),size(XYZ_M,1),1) - XYZ_M;

    % calculate the sum of the absolute gradient on modified XYZ median beats
    data_grad_M1 = data_processing(XYZ_M1,'gradsum',5);
    data_grad_M2 = data_processing(XYZ_M2,'gradsum',5);

    % calcualte the average of the sum of the absolute gradient about the
    % origin points with 3 points on either side
    Sum_M1=mean(data_grad_M1(Ori_pt-3:Ori_pt+3));
    Sum_M2=mean(data_grad_M2(Ori_pt2-3:Ori_pt2+3));

    % decide which origin point to use based on the area under the curve and
    % the average gradient about the origin points
    if Area_m1 < Area_m2 && Sum_M1 < Sum_M2 
        dat_ori=Ori_pt;
        VecMag_ML=VecMag_M1;
    else
        dat_ori=Ori_pt2;
        VecMag_ML=VecMag_M2;
    end

    %% plotting
    VCG_fig = figure('visible','off');

    subplot(211)
    hold on
    plot(VecMag_ML,'b','Linewidth',1.5);
    plot([0 600],[0 0],'k:');
    hold off



    subplot(212)
    hold on;
    plot(XYZ_M(:,1),'Linewidth',1);
    plot(XYZ_M(:,2),'Linewidth',1);
    plot(XYZ_M(:,3),'Linewidth',1);
    ylabel('XYZ');
    hold off

    %% write variable to file and move file to destination
    saveas(VCG_fig,strcat(Results_folder,file_name(1:length(file_name)-4),'_origin.jpg'));
    copyfile(file_path,Results_folder);
    save(strcat(Results_folder,file_name(1:length(file_name)-4),'.mat'),'-append','VecMag_ML','dat_ori','-mat');




%% INCLUDED FUNCTIONS



function [ori_pt_idx,ori_val]=abs_grad_fun(xyz_median)

    % calculate the sum of the absolute gradient on XYZ within w1
    diff_w1 = data_processing(xyz_median,'gradsum',2);

    iter_num=0;
    [~,min_idx]=min(diff_w1);

    % find the point at which has the least angle variation at both side of the
    % window
    for ii=1:length(diff_w1)-1
        if abs(diff_w1(length(diff_w1)-ii)-diff_w1(min_idx))<0.1
            iter_num=length(diff_w1)-ii;

            break
        end

    end

    ori_pt_idx=iter_num;
    ori_val=xyz_median(iter_num);	

end


function [ori_pt_idx,ori_pt_cluster] = ori_clustering(xyz_median)

    %% Find the clusters of possible origin point candidate 
    % SYNTAX :  [ori_pt_idx,ori_pt_cluster] = ori_clustering(xyz_median)
    % INPUT  :  xyz_median  - XYZ median beat (segmented to a window) 
    %                         size(M x 3).
    % OUTPUT :  ori_pt_idx - the index of the calculated origin point
    %           ori_pt_cluster - origin point candidate clusters index

    %% Process data
    processed_data = data_processing(xyz_median,'norm_mvvar',10);
    [N,~,bin_label] = histcounts(processed_data,'BinWidth',10);
    xyz_norm = nan(size(xyz_median,1),1);

    for ii=1:length(xyz_median(:,1))
        xyz_norm(ii)=norm([xyz_median(ii,1),xyz_median(ii,2),xyz_median(ii,3)]);
    end

    %% Find origin point candidate clusters
    idx = find(N == max(N));
    candidate_idx = find(bin_label==idx(1));
    idx_diff = diff(candidate_idx);
    diff_idx = [0; find(idx_diff>=4)];

    if length(diff_idx) == 1
        % when there's only one point found
        ori_pt_cluster = candidate_idx;
        ori_pt_idx = round(median(candidate_idx));
    else
        % look for two largest clusters
        cluster_ele_count = zeros(length(diff_idx),1);

        for ii = 1:length(diff_idx)
            if ii == length(diff_idx)
                cluster_ele_count(ii) = sum(idx_diff(diff_idx(ii)+1:...
                                            length(idx_diff)));
            else
                cluster_ele_count(ii) = sum(idx_diff(diff_idx(ii)+1:...
                                            diff_idx(ii+1)-1));
            end
        end

        [~,I]=sort(cluster_ele_count,'descend');

        if I(1) == max(I)
            cluster_group1 = candidate_idx(diff_idx(I(1))+1:end);
        else
            cluster_group1 = candidate_idx(diff_idx(I(1))+1:...
                                           diff_idx(I(1)+1)-1);
        end

        if I(2) == max(I)
                    cluster_group2 = candidate_idx(diff_idx(I(2))+1:end);
        else
            cluster_group2 = candidate_idx(diff_idx(I(2))+1:...
                                           diff_idx(I(2)+1)-1);
        end

        %% find the best cluster by total average gradient of the clusters
        data_gradient_sum = data_processing(xyz_median,'gradsum',5);
        cluster_avg_slope(1) = mean(data_gradient_sum(cluster_group1));
        cluster_avg_slope(2) = mean(data_gradient_sum(cluster_group2));

        [~,bestgroup] = min(cluster_avg_slope);

        if bestgroup == 1
            ori_pt_cluster = cluster_group1;
            ori_pt_idx = round(median(cluster_group1));
        else
            ori_pt_cluster = cluster_group2;
            ori_pt_idx = round(median(cluster_group2));
        end
    end
end

function processed_data = data_processing(data,process,process_window)
    %% find the clusters of possible origin point candidate 
    % SYNTAX : processed_data = data_processing(data,process,process_window)
    % INPUT  : data            - data matrix with size(M x 3) [X,Y,X leads].
    %          process         - the method of processing data.
    %                            'norm_var': norn of the moving variance 
    %                            with given window.
    %                            'gradsum': sum of the absolute gradient
    %                            of the smoothed XYZ signals with given window.
    %          process_window  - window for specific process
    % OUTPUT : processed_data  - data processed based on specified process
    %                             and process window
    switch process
        case 'norm_mvvar'
            data_norm = nan(size(data,1),1);
            for ii=1:length(data(:,1))
                data_norm(ii)=norm([data(ii,1),data(ii,2),data(ii,3)]);
            end
            processed_data = movvar(data_norm,process_window);
        case 'gradsum'
            processed_data = abs(gradient(smooth(data(:,1),process_window)))...
                            +abs(gradient(smooth(data(:,2),process_window)))...
                            +abs(gradient(smooth(data(:,3),process_window)));

    end
end



