classdef ECGtask_QRS_detection < ECGtask % ECGtask for ECGwrapper (for Matlab) % --------------------------------- % % Description: % % Abstract class for defining ECGtask interface % % Adding user-defined QRS detectors: % A QRS detector that has the following interface can be added to the task: % % [positions_single_lead, position_multilead] = your_QRS_detector( ECG_matrix, ECG_header, progress_handle, payload_in); % where the arguments are: % + ECG_matrix, is a matrix size [ECG_header.nsamp ECG_header.nsig] % + ECG_header, is a struct with info about the ECG signal, such as: % .freq, the sampling frequency % .desc, description about the signals. % + progress_handle, is a handle to a waitbar object, that can be used % to track the progress within your function. See the % documentation about this class in this kit. % + payload_in, is a user data variable allowed to be sent each call to % your function. It is sent, via the payload property of this % class, for example: % % this_ECG_wrappers.ECGtaskHandle.payload = your_variable; % this_ECG_wrappers.ECGtaskHandle.payload = {your_var1 your_var2}; % this_ECG_wrappers.ECGtaskHandle.payload = load(cached_filenames); % % the output of your function must be: % + positions_single_lead, a cell array size ECG_header.nsig with the % QRS sample locations found in each lead. % + position_multilead, a numeric vector with the QRS locations % calculated using multilead rules. % % Author: Mariano Llamedo Soria (llamedom at {electron.frba.utn.edu.ar; unizar.es} % Version: 0.1 beta % Birthdate : 18/2/2013 % Last update: 18/2/2013 properties(GetAccess = public, Constant) name = 'QRS_detection'; target_units = 'ADCu'; doPayload = true; end properties( GetAccess = public, SetAccess = private) % if user = memory; % memory_constant is the fraction respect to user.MaxPossibleArrayBytes % which determines the maximum input data size. memory_constant = 0.3; started = false; end properties( Access = private, Constant) cQRSdetectors = {'all-detectors' 'wavedet' 'pantom' 'aristotle' 'gqrs' 'sqrs' 'wqrs' 'ecgpuwave' 'epltdqrs1' 'epltdqrs2' }; WFDBMaxSampRate = 260; % Hz end properties( Access = private ) detectors2do bWFDBdetectors WFDBtarget_samp_rate = 250; %Hz WFDBResampled = false; tmp_path_local WFDB_bin_path WFDB_cmd_prefix_str lead_idx = []; lead_names = []; wavedet_config end properties progress_handle tmp_path detectors = 'all-detectors'; only_ECG_leads = false; gqrs_config_filename = []; payload CalculatePerformance = false; detection_threshold = 1; end methods function obj = ECGtask_QRS_detection(obj) obj.wavedet_config.setup.wavedet.QRS_detection_only = true; end function Start(obj, ECG_header, ECG_annotations) if( obj.only_ECG_leads ) obj.lead_idx = get_ECG_idx_from_header(ECG_header); else % 'all-leads' obj.lead_idx = 1:ECG_header.nsig; end if( isempty(obj.lead_idx) ) return end if( any(strcmpi('all-detectors', obj.detectors)) ) obj.detectors2do = obj.cQRSdetectors(2:end); else if( ischar(obj.detectors) ) obj.detectors2do = cellstr(obj.detectors); else obj.detectors2do = obj.detectors; end end % lead names desambiguation str_aux = regexprep(cellstr(ECG_header.desc), '\W*(\w+)\W*', '$1'); obj.lead_names = regexprep(str_aux, '\W', '_'); [str_aux2, ~ , aux_idx] = unique( obj.lead_names ); aux_val = length(str_aux2); if( aux_val ~= ECG_header.nsig ) for ii = 1:aux_val bAux = aux_idx==ii; aux_matches = sum(bAux); if( sum(bAux) > 1 ) obj.lead_names(bAux) = strcat( obj.lead_names(bAux), repmat({'v'}, aux_matches,1), cellstr(num2str((1:aux_matches)')) ); end end end obj.lead_names = regexprep(obj.lead_names, '\W*(\w+)\W*', '$1'); obj.lead_names = regexprep(obj.lead_names, '\W', '_'); obj.bWFDBdetectors = ~isempty(intersect({'aristotle' 'gqrs' 'sqrs' 'wqrs' 'ecgpuwave' 'epltdqrs1' 'epltdqrs2'}, obj.detectors2do)); % local path required to avoid network bottlenecks in distributed filesystems if( isunix() && exist('/scratch/', 'dir') ) str_username = getenv('USER'); obj.tmp_path_local = ['/scratch/' str_username filesep]; if( ~exist(obj.tmp_path_local, 'dir') ) if(~mkdir(obj.tmp_path_local)) obj.tmp_path_local = '/scratch/'; end end obj.tmp_path = tempdir; else if( isempty(obj.tmp_path) ) obj.tmp_path = tempdir; end obj.tmp_path_local = obj.tmp_path; end obj.WFDB_bin_path = init_WFDB_library(obj.tmp_path_local); obj.WFDB_cmd_prefix_str = WFDB_command_prefix(obj.tmp_path_local); obj.started = true; end function payload_out = Process(obj, ECG, ECG_start_offset, ECG_sample_start_end_idx, ECG_header, ECG_annotations, ECG_annotations_start_end_idx ) payload_out = []; if( ~obj.started ) obj.Start(ECG_header); if( ~obj.started ) cprintf('*[1,0.5,0]', 'Task %s unable to be started for %s.\n', obj.name, ECG_header.recname); return end end if(obj.bWFDBdetectors) % MIT conversion is needed for WFDB detectors. % Resampled to obj.WFDBtarget_samp_rate Hz for performance ECG_header_WFDB = ECG_header; ECG_header_WFDB.recname = regexprep(ECG_header_WFDB.recname, '\W*(\w+)\W*', '$1'); ECG_header_WFDB.recname = regexprep(ECG_header_WFDB.recname, '\W', '_'); MIT_filename = [ECG_header_WFDB.recname '_' num2str(ECG_start_offset) '_' num2str(ECG_header_WFDB.nsamp+ECG_start_offset-1) ]; ECG_header_WFDB.recname = MIT_filename; if(ECG_header_WFDB.freq > obj.WFDBMaxSampRate ) [p, q] = rat(obj.WFDBtarget_samp_rate/ECG_header_WFDB.freq); obj.WFDBResampled = true; ECG_header_WFDB.freq = round(ECG_header_WFDB.freq * p / q); else obj.WFDBResampled = false; end MIT_filename = [obj.tmp_path_local MIT_filename '.dat']; fidECG = fopen(MIT_filename, 'w'); try if( obj.WFDBResampled ) fwrite(fidECG, int16(round(resample(double(ECG), p, q)))', 'int16', 0 ); else fwrite(fidECG, ECG', 'int16', 0 ); end fclose(fidECG); catch MEE fclose(fidECG); rethrow(MEE); end writeheader(obj.tmp_path_local, ECG_header_WFDB); if( isempty(obj.gqrs_config_filename) ) aux_str = obj.WFDB_bin_path; if( aux_str(end) == filesep ) obj.gqrs_config_filename = [aux_str 'gqrs.conf' ]; else obj.gqrs_config_filename = [aux_str filesep 'gqrs.conf' ]; end end end cant_QRSdetectors = length(obj.detectors2do); for ii = 1:cant_QRSdetectors this_detector = obj.detectors2do{ii}; [this_detector, this_detector_name] = strtok(this_detector, ':'); if( isempty(this_detector_name) ) this_detector_name = this_detector; else this_detector_name = this_detector_name(2:end); end cprintf( 'Blue', [ 'Processing QRS detector ' this_detector_name '\n' ] ); %% perform QRS detection switch( this_detector ) case 'wavedet' %% Wavedet delineation try obj.wavedet_config.setup.wavedet.QRS_detection_thr = repmat( obj.detection_threshold, 5, 1); [position_multilead, positions_single_lead] = wavedet_interface(ECG, ECG_header, [], obj.lead_idx, obj.wavedet_config, ECG_sample_start_end_idx, ECG_start_offset, obj.progress_handle); for jj = rowvec(obj.lead_idx) % QRS detections in milliseconds payload_out.(['wavedet_' obj.lead_names{jj} ]).time = colvec(positions_single_lead(jj).qrs); end % QRS detections in milliseconds payload_out.wavedet_multilead.time = colvec(position_multilead.qrs); catch aux_ME strAux = sprintf('Wavedet failed in recording %s\n', ECG_header.recname); strAux = sprintf('%s\n', strAux); report = getReport(aux_ME); fprintf(2, '%s\nError report:\n%s', strAux, report); end case 'pantom' %% Pan-Tompkins detector for jj = rowvec(obj.lead_idx) obj.progress_handle.checkpoint(['Lead ' ECG_header.desc(jj,:)]) try aux_val = PeakDetection2(double(ECG(:,jj)), ECG_header.freq, [], [], [], obj.detection_threshold * 0.2, [], obj.progress_handle); % filter and offset aux_val = aux_val( aux_val >= ECG_sample_start_end_idx(1) & aux_val <= ECG_sample_start_end_idx(2) ) + ECG_start_offset - 1; % QRS detections in milliseconds payload_out.(['pantom_' obj.lead_names{jj} ]).time = colvec(aux_val); catch aux_ME strAux = sprintf('Pan & Tompkins failed in recording %s lead %s\n', ECG_header.recname, ECG_header.desc(jj,:) ); strAux = sprintf('%s\n', strAux); report = getReport(aux_ME); fprintf(2, 'Error report:\n%s', report); end end case { 'aristotle' 'gqrs' 'sqrs' 'wqrs' 'ecgpuwave' 'epltdqrs1' 'epltdqrs2'} %% WFDB_comp_interface (sqrs, wqrs, aristotle, ecgpuwave) [status, ~] = system([ obj.WFDB_cmd_prefix_str this_detector ' -h' ]); if( status ~= 0 ) disp_string_framed(2, sprintf('Could not execute "%s" QRS detector.', this_detector)); else % QRS detection. for jj = rowvec(obj.lead_idx) if( any(strcmpi( 'sqrs', this_detector ) ) ) file_name_orig = [obj.tmp_path_local ECG_header_WFDB.recname '.qrs' ]; this_thrs = round(500 * obj.detection_threshold); elseif( any(strcmpi( 'wqrs' , this_detector ) ) ) file_name_orig = [obj.tmp_path_local ECG_header_WFDB.recname '.wqrs' ]; this_thrs = round(100 * obj.detection_threshold); elseif( any(strcmpi( {'epltdqrs1' 'epltdqrs2'} , this_detector ) ) ) file_name_orig = [obj.tmp_path_local ECG_header_WFDB.recname '.epl' ]; this_thrs = obj.detection_threshold; else file_name_orig = [obj.tmp_path_local ECG_header_WFDB.recname '.' this_detector num2str(jj) ]; this_thrs = obj.detection_threshold; end try if( any(strcmpi( {'sqrs' 'wqrs' 'epltdqrs1' 'epltdqrs2'} , this_detector ) ) ) % wqrs tiene una interface diferente [status, ~] = system([ obj.WFDB_cmd_prefix_str this_detector ' -r ' ECG_header_WFDB.recname ' -s ' num2str(jj-1) ' -m ' num2str(this_thrs) ]); if( status ~= 0 ); disp_string_framed(2, sprintf('%s failed in recording %s lead %s', this_detector, ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end % elseif( any(strcmpi( {'epltdqrs1' 'epltdqrs2'}, this_detector ) ) ) % [status, ~] = system([ obj.WFDB_cmd_prefix_str this_detector ' ' ECG_header_WFDB.recname ' ' num2str(jj-1)]); % if( status ~= 0 ); disp_string_framed(2, sprintf('%s failed in recording %s lead %s', this_detector, ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end elseif( any(strcmpi( 'ecgpuwave', this_detector ) ) ) [status, ~] = system([ obj.WFDB_cmd_prefix_str this_detector ' -r ' ECG_header_WFDB.recname ' -a ' this_detector num2str(jj) ' -s ' num2str(jj-1)]); if( status ~= 0 ); disp_string_framed(2, sprintf('%s failed in recording %s lead %s', this_detector, ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end elseif( any(strcmpi( 'aristotle', this_detector ) ) ) [status, ~] = system([ obj.WFDB_cmd_prefix_str this_detector ' -r ' ECG_header_WFDB.recname ' -o ' this_detector num2str(jj) ' -s ' num2str(jj-1) ' -G ' num2str(this_thrs) ]); if( status ~= 0 ); disp_string_framed(2, sprintf('%s failed in recording %s lead %s', this_detector, ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end else if( strcmpi( 'gqrs', this_detector ) && exist(obj.gqrs_config_filename, 'file') ) % using the configuration file to % post-process [status, ~] = system([ obj.WFDB_cmd_prefix_str 'gqrs -c ' obj.gqrs_config_filename ' -r ' ECG_header_WFDB.recname ' -s ' num2str(jj-1) ' -m ' num2str(this_thrs) ]); if( status == 0 ) [status, ~] = system([ obj.WFDB_cmd_prefix_str 'gqpost -c ' obj.gqrs_config_filename ' -r ' ECG_header_WFDB.recname ' -o ' this_detector num2str(jj) ]); if( status ~= 0 ) disp_string_framed(2, sprintf('gqpost failed in recording %s lead %s', ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end else disp_string_framed(2, sprintf('gqrs failed in recording %s lead %s', ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end delete([obj.tmp_path_local ECG_header_WFDB.recname '.qrs' ]); else % run only WFDB compatible detector [status, ~] = system([ obj.WFDB_cmd_prefix_str this_detector ' -r ' ECG_header_WFDB.recname ' -o ' this_detector num2str(jj) ' -s ' num2str(jj-1)]); if( status ~= 0 ); disp_string_framed(2, sprintf('%s failed in recording %s lead %s', this_detector, ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ) ); end end end catch aux_ME strAux = sprintf( '%s failed in recording %s lead %d\n', this_detector, ECG_header_WFDB.recname, jj); report = getReport(aux_ME); fprintf(2, '%s\nError report:\n%s', strAux, report); end if( status == 0 ) if( exist(file_name_orig, 'file') ) % reference comparison anns_test = []; try anns_test = readannot(file_name_orig); if( isempty(anns_test) ) payload_out.([this_detector '_' obj.lead_names{jj} ]).time = []; else anns_test = AnnotationFilterConvert(anns_test, 'MIT', 'AAMI'); if( obj.WFDBResampled ) %take sample locations %to its original %sampling rate anns_test.time = round(anns_test.time * ECG_header.freq / ECG_header_WFDB.freq); end % filter and offset anns_test.time = anns_test.time(anns_test.time >= ECG_sample_start_end_idx(1) & anns_test.time <= ECG_sample_start_end_idx(2)) + ECG_start_offset - 1; payload_out.([this_detector '_' obj.lead_names{jj} ]).time = colvec(anns_test.time); end catch aux_ME if( strcmpi(aux_ME.identifier, 'MATLAB:nomem') ) payload_out.([this_detector '_' obj.lead_names{jj} ]).time = []; else strAux = sprintf( '%s failed in recording %s lead %s\n', this_detector, ECG_header_WFDB.recname, ECG_header_WFDB.desc(jj,:) ); report = getReport(aux_ME); error('ECGtask_QRS_detection:WFDB', '%s\nError report:\n%s', strAux, report); end end delete(file_name_orig); else payload_out.([this_detector '_' obj.lead_names{jj} ]).time = []; end else payload_out.([this_detector '_' obj.lead_names{jj} ]).time = []; if( exist(file_name_orig, 'file') ) delete(file_name_orig); end end end end case 'user' %% user-defined QRS detector try if( exist(this_detector_name) == 2 ) % ud_func_pointer = eval(['@' this_detector_name]); ud_func_pointer = str2func(this_detector_name); obj.progress_handle.checkpoint([ 'User defined function: ' this_detector_name]) ECG_header_aux = trim_ECG_header(ECG_header, obj.lead_idx); [positions_single_lead, position_multilead] = ud_func_pointer( double(ECG(:,obj.lead_idx)), ECG_header_aux, obj.progress_handle, obj.payload); for jj = 1:length(obj.lead_idx) aux_val = positions_single_lead{jj}; % filter and offset aux_val = aux_val( aux_val >= ECG_sample_start_end_idx(1) & aux_val <= ECG_sample_start_end_idx(2) ) + ECG_start_offset - 1; % QRS detections in milliseconds payload_out.([ this_detector_name '_' obj.lead_names{jj} ]).time = colvec(aux_val); end if( ~isempty(position_multilead) ) % filter and offset position_multilead = position_multilead( position_multilead >= ECG_sample_start_end_idx(1) & position_multilead <= ECG_sample_start_end_idx(2) ) + ECG_start_offset - 1; payload_out.([ this_detector_name '_multilead' ]).time = colvec(position_multilead); end else disp_string_framed(2, sprintf('Function "%s" is not reachable in path.', this_detector_name)); fprintf(1, 'Make sure that exist(%s) == 2\n',this_detector_name); end catch aux_ME disp_string_framed(2, sprintf('Detector "%s" failed in recording %s lead %s', this_detector_name, ECG_header.recname, ECG_header.desc(jj,:) ) ); report = getReport(aux_ME); fprintf(2, 'Error report:\n%s', report); end end end % Add QRS detections quality metrics, Names, etc. payload_out = calculateSeriesQuality(payload_out, ECG_header, [1 ECG_header.nsamp] + ECG_start_offset - 1 ); % delete intermediate tmp files if(obj.bWFDBdetectors) delete([obj.tmp_path_local ECG_header.recname '.*' ]); end % calculate performance if( obj.CalculatePerformance ) AnnNames = payload_out.series_quality.AnnNames(:,1); cant_lead_name = size(AnnNames,1); payload_out.series_performance.conf_mat = zeros(2,2,cant_lead_name); if(isempty(ECG_annotations)) disp_string_framed(2, sprintf('Trusted references not found for %s', ECG_header.recname) ); else % offset refs, produced anns were already shifted ECG_annotations.time = ECG_annotations.time + ECG_start_offset - 1; for kk = 1:cant_lead_name payload_out.series_performance.conf_mat(:,:,kk) = bxb(ECG_annotations, payload_out.(AnnNames{kk}).time, ECG_header ); end end end end function payload = Finish(obj, payload, ECG_header) if( isfield(payload, 'series_quality') && isfield(payload.series_quality, 'ratios') ) payload.series_quality.ratios = mean(payload.series_quality.ratios, 2); end end function payload = Concatenate(obj, plA, plB) payload = ConcatenateQRSdetectionPayloads(obj, plA, plB); end %% property restriction functions function set.detectors(obj,x) if( (ischar(x) || iscellstr(x)) ) x = cellstr(x); aux_val = colvec(intersect(obj.cQRSdetectors, x)); aux_idx = find(cellfun(@(a)(~isempty(strfind(a, 'user:'))), x)); if( isempty(aux_idx) && isempty(aux_val) ) warning('ECGtask_QRS_detection:BadArg', 'Invalid detectors.'); else obj.detectors = [aux_val; colvec(x(aux_idx)) ]; end else warning('ECGtask_QRS_detection:BadArg', 'Invalid detectors.'); end end function set.only_ECG_leads(obj,x) if( islogical(x) ) obj.only_ECG_leads = x; else warning('ECGtask_QRS_detection:BadArg', 'Invalid lead configuration.'); end end function set.gqrs_config_filename(obj,x) if( ischar(x) ) if(exist(x, 'file')) obj.gqrs_config_filename = x; else warning('ECGtask_QRS_detection:BadArg', 'obj.gqrs_config_filename does not exist.'); end else warning('ECGtask_QRS_detection:BadArg', 'obj.gqrs_config_filename must be a string.'); end end function set.tmp_path_local(obj,x) if( ischar(x) ) if(exist(x, 'dir')) obj.tmp_path_local = x; else if(mkdir(x)) obj.tmp_path_local = x; else warning('ECGtask_QRS_detection:BadArg', ['Could not create ' x ]); end end else warning('ECGtask_QRS_detection:BadArg', 'tmp_path_local must be a string.'); end end function set.tmp_path(obj,x) if( ischar(x) ) if(exist(x, 'dir')) obj.tmp_path = x; else if(mkdir(x)) obj.tmp_path = x; else warning('ECGtask_QRS_detection:BadArg', ['Could not create ' x ]); end end else warning('ECGtask_QRS_detection:BadArg', 'tmp_path_local must be a string.'); end end end methods ( Access = private ) end end