#! /usr/bin/perl
#
# Program to get QT values for the Computers in Cardiology Challenge 2006.
#
# Copyright 2006 by Raphael Schneider (rasch@librasch.org)
#
# All rights reserved. This program is free software; you can redistribute
# it and/or modify it under the same terms as Perl itself.
#
# This Perl script needs libRASCH version 0.8.0pre2 or greater from www.librasch.org
# and the Perl API bindings have to be installed also (see INSTALL file coming
# with the libRASCH distribution).
#
# Usage:
#
# perl get_qt_cinc_2006.pl --template=templ-file --file=out-file
# [--detect] [--ch=X] [--check-morphology]
# [--skip-powerline] <ecg-dir(s)>
#
# with --template=templ-file selecting the template of the entry-form
# --file=out-file selecting the name of the final entry-form
# --detect if given, a beat detection is performed, if not,
# the beat information (position and annotation) which
# is available is used
# --ch=X do the QT measurement using channel 'X'
# (channel-numbers start with 0 -> lead I is 0, lead II is 1, ...)
# --check-morphology if given, the morphologies of the QRS complexes and T-waves
# are forced to be the same for all beats (e.g. when the Q-wave
# on some beats is too small to be detected but seen on the
# majority of the beats, the threshold used to detect Q-waves
# is reduced so if small Q-waves will be detected)
# --skip-powerline if given, the powerline-noise filter is NOT applied to
# the raw ecg data
# <ecg-dir(s)> the directory(s) with the ECGs to process
#
# Example:
#
# perl get_qt_cinc_2006.pl --detect --ch=1
# --template=template_div3.txt --file=results_div3.txt
# ./patient*
#
#
# Version 0.0: rasch@librasch.org (2006-04-19)
# first verision
# Version 1.0: rasch@librasch.org (2006-09-11)
# final changes before sending it to PhysioNet
#
use strict;
use warnings;
$|++;
use RASCH;
use Getopt::Long;
my $debug = 1;
########## get and check options ##########
my $do_detect;
my $ch = 1; # ch-1 is lead II
my $check_morphology;
my $method = 0;
my $outfile; # default is print to stdout
my $template_file = './template.txt';
my $skip_powerline_filter;
GetOptions("detect" => \$do_detect,
"ch=s" => \$ch,
"check-morphology" => \$check_morphology,
"method=s" => \$method,
"file=s" => \$outfile,
"template=s" => \$template_file,
"skip-powerline" => \$skip_powerline_filter);
# if no arguments left -> no files/dirs given to process
exit 0 if ($#ARGV < 0);
# check if file-/dir-names are given on stdin
my @files_use;
if ($ARGV[0] eq '-')
{
while (<>)
{
my $fn = $_;
chomp($fn);
push(@files_use, $fn);
}
}
else
{
@files_use = @ARGV;
}
########## start processing the files/dirs ##########
my $ra = new RASCH or die "can't open libRASCH\n";
my $meas_cnt = 0;
my %results = ();
for my $curr (@files_use)
{
my $meas = $ra->open_meas($curr);
if (defined($meas))
{
$meas_cnt++;
print STDERR "process measurement $curr (#$meas_cnt)...";
process_meas($ra, $meas, $do_detect, $ch, $check_morphology, $method, $skip_powerline_filter, \%results);
print STDERR "\n";
}
else
{
# open was not sucessful, maybe it is a directory
my @sigs = @{$ra->find_meas($curr)};
if ($#sigs < 0)
{
print STDERR "can't handle file/dir $curr\n";
next;
}
for my $curr2 (@sigs)
{
my $fn = $curr2->filename();
$meas_cnt++;
print STDERR "process measurement $fn in directory $curr (#$meas_cnt)...";
$meas = $ra->open_meas($fn);
if (not defined($meas))
{
print STDERR "\nError opening $fn. This should not happen\n";
next;
}
process_meas($ra, $meas, $do_detect, $ch, $check_morphology, $method, $skip_powerline_filter, \%results);
print STDERR "\n";
}
}
}
########## write results ##########
write_results($template_file, $outfile, \%results);
exit 0;
######################### sub-routines #########################
sub process_meas
{
my ($ra, $meas, $do_detect, $ch, $check_morph, $method, $skip_powerline_filter, $results) = @_;
my $clh = undef;
if ($do_detect)
{
$clh = do_detect($ra, $meas, $ch, $check_morph, $skip_powerline_filter);
return undef if (not defined($clh));
}
else
{
my $eva = $meas->get_default_eval();
if (not defined($eva))
{
print STDERR "nok";
return undef;
}
my $all_clh = $eva->get_class('heartbeat');
$clh = $all_clh->[0];
if ((not defined($clh)) || ($clh == 0))
{
print STDERR "nok";
return undef;
}
}
my $rec = $meas->get_first_rec();
my ($value, $name, $desc) = $rec->get_info(info => 'rec_name');
my $rec_name = $value->value();
($value, $name, $desc) = $rec->get_info(info => 'ch_samplerate', ch => $ch);
my $samplerate = $value->value();
my $res;
if ($method == 0)
{
$res = qt_method0($clh, $ch);
}
elsif ($method == 1)
{
$res = qt_method1($clh, $ch);
}
elsif ($method == 2)
{
$res = qt_method2($clh, $ch);
}
else
{
print STDERR "Unknown method to get QT times. Please select a known one (0-2).\n";
exit -1;
}
$res->[2] = $samplerate;
$results->{$rec_name} = $res;
return $res;
} # sub process_meas()
#################### beat detection code ####################
sub do_detect
{
my ($raq, $meas, $ch, $check_morph, $skip_powerline_filter) = @_;
my $eva = $meas->add_eval("get_qt.pl", "evaluation for CinC Challenge 2006");
if (not defined($eva))
{
print STDERR "Error creating evaluation";
return 0;
}
do_beat_detect($ra, $meas, $eva, $ch, $skip_powerline_filter);
my $all_clh = $eva->get_class('heartbeat');
my $clh = $all_clh->[0];
if ((not defined($clh)) || ($clh == 0))
{
print STDERR "nok";
return undef;
}
# get most used morphologies for qrs and t
my ($qrs_morph_use, $t_morph_use) = get_max_morph($clh);
# now do beat-detection again forcing to use the
# maxed used morphologies
do_beat_detect($ra, $meas, $eva, $ch, $skip_powerline_filter, $clh, $qrs_morph_use, $t_morph_use);
$all_clh = $eva->get_class('heartbeat');
$clh = $all_clh->[0];
if ((not defined($clh)) || ($clh == 0))
{
print STDERR "nok";
return undef;
}
do_fiducial($ra, $meas, $clh);
do_template($ra, $meas, $clh);
if ($check_morph)
{
do_morphology_check($ra, $meas, $eva, $clh, $skip_powerline_filter);
}
do_ecg($ra, $meas, $clh);
$meas->save_eval();
return $clh;
} # sub do_detect()
sub do_beat_detect
{
my $ra = shift;
my $meas = shift;
my $eva = shift;
my $ch = shift;
my $skip_powerline_filter = shift;
my $clh = (defined(@_) ? shift : undef);
my $qrs_morph_force = (defined(@_) ? shift : undef);
my $t_morph_force = (defined(@_) ? shift : undef);
if (defined($clh))
{
$clh->delete();
}
my $pl = $ra->get_plugin_by_name("detect-ecg", 0);
if (not defined($pl))
{
print STDERR "Can not find 'detect-ecg' plugin";
return 0;
}
my $proc = $pl->get_process($meas->get_handle());
if (not defined($proc))
{
print STDERR "can't initialize processing";
return -1;
}
my $ret = $proc->set_option('save_in_eval', 'l', 1);
$ret = $proc->set_option('do_interpolation', 'l', 1);
$ret = $proc->set_option('eh', 'h', $eva->get_handle());
$ret = $proc->set_option('filter_powerline_noise', 'l', ($skip_powerline_filter ? 0 : 1));
$ret = $proc->set_option('force_qrs_type', 'l', $qrs_morph_force) if (defined($qrs_morph_force));
$ret = $proc->set_option('force_t_type', 'l', $t_morph_force) if (defined($t_morph_force));
my $rec = $meas->get_first_rec();
my $num_ch = $rec->get_info(info => 'rec_num_channel')->value();
my @chs = ();
if ($ch == -1)
{
for (my $i = 0; $i < $num_ch; $i++)
{
$chs[$i] = $i;
}
}
else
{
$chs[0] = $ch;
}
$ret = $proc->set_option('ch', 'pl', \@chs);
my $results = $proc->process();
return 0;
} # sub do_beat_detect()
sub get_max_morph
{
my ($clh) = @_;
my $events = $clh->get_events();
# get morphology types
my $prop_qrs_morph = $clh->get_prop('ecg-qrs-type');
if (not defined($prop_qrs_morph))
{
print STDERR "no 'ecg-qrs-type' event-property in the event-class";
return -1;
}
my $prop_t_morph = $clh->get_prop('ecg-t-type');
if (not defined($prop_t_morph))
{
print STDERR "no 'ecg-t-type' event-property in the event-class";
return -1;
}
my %qrs_morph = ();
my %t_morph = ();
my $cnt_used = 0;
for (@$events)
{
my $qrs_type = $prop_qrs_morph->get_value($_, $ch);
next if ($qrs_type <= 0);
my $t_type = $prop_t_morph->get_value($_, $ch);
next if ($t_type <= 0);
$qrs_morph{$qrs_type} = 0 if (not defined($qrs_morph{$qrs_type}));
$qrs_morph{$qrs_type}++;
$t_morph{$t_type} = 0 if (not defined($t_morph{$t_type}));
$t_morph{$t_type}++;
$cnt_used++;
}
my $qrs_morph_use;
my $num = 0;
for (keys(%qrs_morph))
{
if ($num < $qrs_morph{$_})
{
$qrs_morph_use = $_;
$num = $qrs_morph{$_};
}
}
my $t_morph_use;
$num = 0;
for (keys(%t_morph))
{
# if bi-phasic, use only when more than 75% of the beats where found
# having a bi-phasic T-wave
next if (($_ > 100) && (($t_morph{$_}/$cnt_used) < 0.75));
if ($num < $t_morph{$_})
{
$t_morph_use = $_;
$num = $t_morph{$_};
}
}
return ($qrs_morph_use, $t_morph_use);
} # sub get_max_morph()
sub do_fiducial
{
my ($ra, $meas, $clh) = @_;
my $pl = $ra->get_plugin_by_name("fiducial-point", 0);
if (not defined($pl))
{
print STDERR "Can not find 'fiducial-point' plugin";
return 0;
}
my $proc = $pl->get_process($meas->get_handle());
if (not defined($proc))
{
print STDERR "can't initialize processing";
return -1;
}
my $ret = $proc->set_option('use_class', 'l', 1);
$ret = $proc->set_option('clh', 'h', $clh->get_handle());
$ret = $proc->set_option('save_in_class', 'l', 1);
$ret = $proc->set_option('ch', 'l', 0);
my $results = $proc->process();
return 0;
} # sub do_fiducial()
sub do_template
{
my ($ra, $meas, $clh) = @_;
my $pl = $ra->get_plugin_by_name("template", 0);
if (not defined($pl))
{
print STDERR "Can not find 'template' plugin";
return 0;
}
my $proc = $pl->get_process($meas->get_handle());
if (not defined($proc))
{
print STDERR "can't initialize processing";
return -1;
}
my $ret = $proc->set_option('use_class', 'l', 1);
$ret = $proc->set_option('clh', 'h', $clh->get_handle());
$ret = $proc->set_option('save_in_class', 'l', 1);
$ret = $proc->set_option('pos_name', 'c', 'qrs-pos');
$ret = $proc->set_option('templ_name', 'c', 'qrs-template');
$ret = $proc->set_option('templ_corr', 'c', 'qrs-template-corr');
$ret = $proc->set_option('accept', 'd', 0.85);
$ret = $proc->set_option('slope_accept_low', 'd', 0.6);
$ret = $proc->set_option('slope_accept_high', 'd', 1.4);
$ret = $proc->set_option('win_before', 'd', 0.15);
$ret = $proc->set_option('win_after', 'd', 0.15);
my $rec = $meas->get_first_rec();
my $num_ch = $rec->get_info(info => 'rec_num_channel')->value();
my @chs = ();
for (my $i = 0; $i < $num_ch; $i++)
{
$chs[$i] = $i;
}
$ret = $proc->set_option('ch', 'pl', \@chs);
my $results = $proc->process();
return 0;
} # sub do_template()
sub do_morphology_check
{
my ($ra, $meas, $eva, $clh, $skip_powerline_filter) = @_;
my $pl = $ra->get_plugin_by_name("detect-ecg", 0);
if (not defined($pl))
{
print STDERR "Can not find 'detect-ecg' plugin";
return 0;
}
my $proc = $pl->get_process($meas->get_handle());
if (not defined($proc))
{
print STDERR "can't initialize processing";
return -1;
}
my $ret = $proc->set_option('save_in_eval', 'l', 1);
$ret = $proc->set_option('do_interpolation', 'l', 1);
$ret = $proc->set_option('eh', 'h', $eva->get_handle());
$ret = $proc->set_option('clh', 'h', $clh->get_handle());
$ret = $proc->set_option('filter_powerline_noise', 'l', ($skip_powerline_filter ? 0 : 1));
my $rec = $meas->get_first_rec();
my $num_ch = $rec->get_info(info => 'rec_num_channel')->value();
my @chs = ();
if ($ch == -1)
{
for (my $i = 0; $i < $num_ch; $i++)
{
$chs[$i] = $i;
}
}
else
{
$chs[0] = $ch;
}
$ret = $proc->set_option('ch', 'pl', \@chs);
# my $events = get_max_template_events($clh);
my $events = $clh->get_events();
$ret = $proc->set_option('events', 'pl', $events);
# get morphology types
my $prop_qrs_morph = $clh->get_prop('ecg-qrs-type');
if (not defined($prop_qrs_morph))
{
print STDERR "no 'ecg-qrs-type' event-property in the event-class";
return -1;
}
my $prop_t_morph = $clh->get_prop('ecg-t-type');
if (not defined($prop_t_morph))
{
print STDERR "no 'ecg-t-type' event-property in the event-class";
return -1;
}
my %qrs_morph = ();
my %t_morph = ();
for (@$events)
{
my $qrs_type = $prop_qrs_morph->get_value($_, $ch);
next if ($qrs_type <= 0);
my $t_type = $prop_t_morph->get_value($_, $ch);
next if ($t_type <= 0);
$qrs_morph{$qrs_type} = 0 if (not defined($qrs_morph{$qrs_type}));
$qrs_morph{$qrs_type}++;
$t_morph{$t_type} = 0 if (not defined($t_morph{$t_type}));
$t_morph{$t_type}++;
}
my $qrs_morph_use;
my $num = 0;
for (keys(%qrs_morph))
{
if ($num < $qrs_morph{$_})
{
$qrs_morph_use = $_;
$num = $qrs_morph{$_};
}
}
my $t_morph_use;
$num = 0;
for (keys(%t_morph))
{
if ($num < $t_morph{$_})
{
$t_morph_use = $_;
$num = $t_morph{$_};
}
}
print STDERR "qrs-force=$qrs_morph_use t-force=$t_morph_use " if ($debug);
$ret = $proc->set_option('force_qrs_type', 'l', $qrs_morph_use);
$ret = $proc->set_option('force_t_type', 'l', $t_morph_use);
my $results = $proc->process();
} # sub do_morphology_check()
sub do_ecg
{
my ($ra, $meas, $clh) = @_;
my $pl = $ra->get_plugin_by_name("ecg", 0);
if (not defined($pl))
{
print STDERR "Can not find 'ecg' plugin";
return 0;
}
my $proc = $pl->get_process($meas->get_handle());
if (not defined($proc))
{
print STDERR "can't initialize processing";
return -1;
}
my $ret = $proc->set_option('use_class', 'l', 1);
$ret = $proc->set_option('clh', 'h', $clh->get_handle());
$ret = $proc->set_option('save_in_class', 'l', 1);
my $results = $proc->process();
return 0;
} # sub do_ecg()
#################### QT time assessment ####################
#
# (1) select template containing the most beats
# (2) use beat with lowest noise
#
sub qt_method0
{
my ($clh, $ch) = @_;
my $prop_noise = $clh->get_prop('ecg-noise');
if (not defined($prop_noise))
{
print STDERR "no 'ecg-noise' event-property in the event-class";
return undef;
}
my $prop_qt = $clh->get_prop('ecg-qt');
if (not defined($prop_qt))
{
print STDERR "no 'ecg-qt' event-property in the event-class";
return undef;
}
my $events = get_max_template_events($clh);
# get event-id with lowest noise
my $ev_id_min_noise;
my $min_noise = 10000000;
for (@{$events})
{
# skip first beat because of possible measurement problems for it
next if ($_ == 1);
my $value = $prop_qt->get_value($_, $ch);
next if ($value == -1);
$value = $prop_noise->get_value($_, $ch);
next if ($value == -1);
if ($value < $min_noise)
{
$ev_id_min_noise = $_;
$min_noise = $value;
}
}
return undef if (not defined($ev_id_min_noise));
print STDERR "min-noise=$min_noise ev-id=$ev_id_min_noise " if ($debug);
# get QT related values for export and return it
return get_qt_values($clh, $ev_id_min_noise, $ch);
} # sub qt_method0()
#
# (1) select template containing the most beats
# (2) use beat with median QT value in the template-class
#
sub qt_method1
{
my ($clh, $ch) = @_;
my $prop_qt = $clh->get_prop('ecg-qt');
if (not defined($prop_qt))
{
print STDERR "no 'ecg-qt' event-property in the event-class";
return undef;
}
my $events = get_max_template_events($clh);
# get event-id with median QT value
my @data = ();
for (@{$events})
{
my $value = $prop_qt->get_value($_, $ch);
next if ($value == -1);
my @curr = ($value, $_);
push(@data, \@curr);
}
my @data_sort = sort { $a->[0] <=> $b->[0] } @data;
my $len = @data_sort;
my $ev_use = int($len/2);
my $ev_id = $data_sort[$ev_use]->[1];
# skip first beat because of possible measurement problems for it
if ($ev_id == 1)
{
$ev_use++;
$ev_id = $data_sort[$ev_use]->[1];
}
print STDERR "median-qt=$data_sort[$ev_use]->[0] ev-id=$ev_id " if ($debug);
# get QT related values for export and return it
return get_qt_values($clh, $ev_id, $ch);
} # sub qt_method1()
#
# (1) use beat with median QT value over all beats
#
sub qt_method2
{
my ($clh, $ch) = @_;
my $prop_qt = $clh->get_prop('ecg-qt');
if (not defined($prop_qt))
{
print STDERR "no 'ecg-qt' event-property in the event-class";
return undef;
}
my $events = $clh->get_events();
# get event-id with median QT value
my @data = ();
for (@{$events})
{
my $value = $prop_qt->get_value($_, $ch);
next if ($value == -1);
my @curr = ($value, $_);
push(@data, \@curr);
}
my @data_sort = sort { $a->[0] <=> $b->[0] } @data;
my $len = @data_sort;
my $ev_use = int($len/2);
my $ev_id = $data_sort[$len]->[1];
# skip first beat because of possible measurement problems for it
if ($ev_id == 1)
{
$ev_use++;
$ev_id = $data_sort[$ev_use]->[1];
}
print STDERR "median-qt=$data_sort[$ev_use]->[0] ev-id=$ev_id " if ($debug);
# get QT related values for export and return it
return get_qt_values($clh, $ev_id, $ch);
} # sub qt_method2()
########## QT time assessment helper functions ##########
sub get_max_template_events
{
my ($clh) = @_;
my $summaries = $clh->get_sum('template');
if (not defined($summaries))
{
print STDERR "no 'template' summaries in the event-class";
return undef;
}
my $templ = $summaries->[0];
my $parts = $templ->get_part_all();
# get template (and the event-id's) containing the most beats
my @events_all = ();
my $num_ev_max = 0;
my $ev_idx_max;
my $cnt = 0;
for (@{$parts})
{
my $ev = $templ->get_part_events($_);
push(@events_all, $ev);
if (@{$ev} > $num_ev_max)
{
$num_ev_max = @{$ev};
$ev_idx_max = $cnt;
}
$cnt++;
}
print STDERR "num-ev=$num_ev_max " if ($debug);
return $events_all[$ev_idx_max];
} # sub get_max_template_events()
sub get_qt_values
{
my ($clh, $ev_id, $ch) = @_;
my $prop_qrs_pos = $clh->get_prop('qrs-pos');
if (not defined($prop_qrs_pos))
{
print STDERR "no 'qrs-pos' event-property in the event-class";
return undef;
}
my $prop_q = $clh->get_prop('ecg-qrs-start');
if (not defined($prop_q))
{
print STDERR "no 'ecg-qrs-start' event-property in the event-class";
return undef;
}
my $prop_t = $clh->get_prop('ecg-t-end');
if (not defined($prop_t))
{
print STDERR "no 'ecg-t-end' event-property in the event-class";
return undef;
}
my @res = ();
my $qrs_pos = $prop_qrs_pos->get_value($ev_id, $ch);
$res[0] = $qrs_pos + $prop_q->get_value($ev_id, $ch);
$res[1] = $qrs_pos + $prop_t->get_value($ev_id, $ch);
return \@res;
} # sub get_qt_values()
#################### writing results ####################
sub write_results
{
my ($template_file, $outfile, $results) = @_;
my $f = *STDOUT;
if (defined($outfile))
{
open OUT_FILE, ">$outfile" or die "can't create file $outfile: $!\n";
$f = *OUT_FILE;
}
open TEMPL_FILE, "<$template_file" or die "can't open template-file $template_file: $!\n";
while (<TEMPL_FILE>)
{
my $curr_line = $_;
chomp($curr_line);
$curr_line =~ m|^patient[\w]+\/([\w]+)$|;
if (defined($1))
{
my $values = $results->{$1};
if ((not defined($values)) || ($values->[0] == 0))
{
print $f $curr_line . "\n";
}
else
{
my $scale = 1000 / $values->[2];
$values->[0] *= $scale;
$values->[1] *= $scale;
print $f $curr_line . "\t$values->[0]\t$values->[1]\n";
if ($debug)
{
my $qt = $values->[1] - $values->[0];
print $curr_line . ": $qt\n";
}
}
}
else
{
print $f $curr_line . "\n";
}
}
close(TEMPL_FILE);
close($f) if (defined($outfile));
} # sub write_results()