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, % Annabel Li-Pershing, BS, % 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