/* file: gqpost.c		G. Moody	16 November 2006
				Last revised:	15 January 2019
-------------------------------------------------------------------------------
gqpost: A post-processor for gqrs
Copyright (C) 2006-2012 George B. Moody

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, see <http://www.gnu.org/licenses/>.

You may contact the author by e-mail (wfdb@physionet.org) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA).  For updates to this software,
please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________

The companion to this program, gqrs, is a QRS detector that is optimized for
sensitivity.  To reduce the number of false positives, gqpost can read the
output annotation file produced by gqrs and can look for interpolated beat
detections in it. Such events are identified by looking for pairs or longer
groups of inter-beat intervals that sum to values close to the predicted
interbeat interval (in other words, an interpolated event is one that occurs
between a pair of beats that occur at their expected times).  Although
true beats can be interpolated, such beats are relatively rare in comparison
with false detections that are almost always interpolated, so a high proportion
of interpolated events is false, and a simple strategy for improving positive
predictivity is to reject all interpolated events.  To avoid rejecting
(at least some of) the true interpolated beats, gqpost also looks at the size
of each interpolated event as recorded by gqrs in the annotation 'num' field,
and it rejects only those events that are smaller than a threshold value
(which may be adjusted using the -m option).

gqpost is most effective in the context of regular rhythms.  During highly
irregular rhythms, such as atrial fibrillation, the predicted time of the
next beat is highly uncertain, and gqpost rejects very few if any events in
such cases.

gqpost is very fast since it operates entirely on the input annotation
file, without reference to the raw signals.  The output of gqpost is a
complete copy of its input, except that the anntyp fields of any rejected
events occurring between the start and end times (specified on the command
line using the -f and -t options) are set to ARFCT.  gqpost can accept its
own output as input, which allows one to vary the rejection threshold if
necessary by use of multiple post-processing passes.
*/

#include <stdio.h>
#include <stddef.h>
#include <wfdb/wfdb.h>
#include <wfdb/ecgmap.h>
#define PBLEN	12	/* size of predictor array */

/* The `getconf' macro is used to check a line of input (already in `buf') for
   the string named by getconf's first argument.  If the string is found, the
   value following the string (and an optional `:' or `=') is converted using
   sscanf and the format specifier supplied as getconf's second argument, and
   stored in the variable named by the first argument. */
#define getconf(a, fmt)	if (p=strstr(buf,#a)){sscanf(p,#a "%*[=: \t]" fmt,&a);}

void help(void);
char *prog_name(char *p);
void *gcalloc(size_t nmemb, size_t size);
void cleanup(int status);

char *pname;			/* name of this program, used in messages */
char *record = NULL;		/* name of input record */

main(int argc, char **argv)
{
    FILE *config = NULL;
    int dt, i, meanrr, meanrrd, minrrd, state = 0, thresh = 20;
    static double HR, RR, RRdelta;
    WFDB_Anninfo a[2];
    WFDB_Annotation annot;
    WFDB_Frequency sps;
    WFDB_Time from = (WFDB_Time)0, to = (WFDB_Time)0;

    a[0].name = "qrs"; a[0].stat = WFDB_READ;
    a[1].name = "gqp"; a[1].stat = WFDB_WRITE;
    pname = prog_name(argv[0]);

    for (i = 1; i < argc; i++) {
	if (*argv[i] == '-') switch (*(argv[i]+1)) {
	  case 'a':	/* input annotator name */
	    if (++i >= argc) {
		(void)fprintf(stderr,
			      "%s: input annotator name must follow -a\n",
			      pname);
		cleanup(1);
	    }
	    a[0].name = argv[i];
	    break;
	  case 'c':	/* configuration file */
	    if (++i >= argc) {
		(void)fprintf(stderr,
		     "%s: name of configuration file must follow -c\n", pname);
		cleanup(1);
	    }
	    if ((config = fopen(argv[i], "rt")) == NULL) {
		(void)fprintf(stderr,
		     "%s: can't read configuration file %s\n", pname, argv[i]);
		cleanup(1);
	    }
	    break;
	  case 'f':	/* starting time */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: time must follow -f\n", pname);
		cleanup(1);
	    }
	    from = i;
	    break;
	  case 'h':	/* help requested */
	    help();
	    cleanup(0);
	    break;
	  case 'm':	/* threshold */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: threshold must follow -m\n", pname);
		cleanup(1);
	    }
	    thresh = 10*atof(argv[i]) + 0.5;
	    break;
	  case 'o':	/* output annotator name */
	    if (++i >= argc) {
		(void)fprintf(stderr,
			      "%s: output annotator name must follow -o\n",
			      pname);
		cleanup(1);
	    }
	    a[1].name = argv[i];
	    break;
	  case 'r':	/* record name */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: input record name must follow -r\n",
			      pname);
		cleanup(1);
	    }
	    record = argv[i];
	    break;
	  case 't':	/* end time */
	    if (++i >= argc) {
		(void)fprintf(stderr, "%s: time must follow -t\n", pname);
		cleanup(1);
	    }
	    to = i;
	    break;
	  default:
	    (void)fprintf(stderr, "%s: unrecognized option %s\n", pname,
			  argv[i]);
	    cleanup(1);
	}
	else {
	    (void)fprintf(stderr, "%s: unrecognized argument %s\n", pname,
			  argv[i]);
	    cleanup(1);
	}
    }
    if (record == NULL) {
	help();
	cleanup(1);
    }

    /* Read a priori physiologic parameters from the configuration file if
       available. They can be adjusted in the configuration file for pediatric,
       fetal, or animal ECGs. */
    if (config) {
	char buf[256], *p;

	/* Read the configuration file a line at a time. */
	while (fgets(buf, sizeof(buf), config)) {
	    /* Skip comments (empty lines and lines beginning with `#'). */
	    if (buf[0] == '#' || buf[0] == '\n') continue;
	    /* Set parameters.  Each `getconf' below is executed once for
	       each non-comment line in the configuration file. */
	    getconf(HR, "%lf");
	    getconf(RR, "%lf");
	    getconf(RRdelta, "%lf");
	}
	fclose(config);
    }

    /* If any a priori parameters were not specified in the configuration file,
       initialize them here (using values chosen for adult human ECGs). */
    if (HR != 0.0) RR = 60.0/HR;
    if (RR == 0.0) RR = 0.8;
    if (RRdelta == 0.0) RRdelta = RR/4;

    sps = sampfreq(record);
    if (annopen(record, a, 2) < 0) cleanup(2);
    if (getiaorigfreq(0) > 0) {
	setifreq(sps = getiaorigfreq(0));
	setiafreq(0, sps);
    }
    meanrr = RR * sps;
    meanrrd = RRdelta * sps;
    if ((minrrd = meanrrd/2) < 4) minrrd = 4;
    if ((dt = meanrr/80) < 1) dt = 1;

    if (from) {
	from = strtim(argv[from]);
	if (from < 0) from = -from;
    }
    if (to) {
	to = strtim(argv[to]);
	if (to < 0) to = -to;
    }

    while (getann(0, &annot) == 0) {
	switch (state) {
	  case 0:
	    if (annot.time < from) break;
	    else state++;
	    /* fall through to case 1 */
  	  case 1:
	    if (to > (WFDB_Time)0 && annot.time > to) state++;
	    else {	/* reject sub-threshold interpolated events */
		static WFDB_Time tpred, tpredmin, tprev;
		static int rr, rrd;

		if (annot.time < tpredmin &&
		    isqrs(annot.anntyp) &&
		    annot.num < thresh) {   /* event is small and early */
		    WFDB_Annotation nextann;
		    if (getann(0, &nextann) == 0) {	/* peek ahead */
			if (nextann.time <= tpred)
			    annot.anntyp = ARFCT;	/* reject event */
			ungetann(0, &nextann);
		    }
		}
		if (isqrs(annot.anntyp)) {
		    rr = annot.time - tprev;
		    if ((rrd = rr - meanrr) < 0) rrd = -rrd;
		    if (rr > meanrr + dt) meanrr += dt;
		    else if (rr > meanrr - dt) meanrr = rr;
		    else meanrr -= dt;
		    if (rrd > meanrrd + dt) meanrrd += dt;
		    else if (rrd > meanrrd - dt) meanrrd = rrd;
		    else meanrrd -= dt;
		    if (meanrrd < minrrd) meanrrd = minrrd;
		    tpred = annot.time + meanrr + minrrd;
		    tpredmin = annot.time + meanrr - meanrrd;
		    tprev = annot.time;
		}
	    }
	    break;
	  case 2:
	    /* do nothing */
	    break;
	}
	putann(0, &annot);
    }
    cleanup(0);
}


/* prog_name() extracts this program's name from argv[0], for use in error and
   warning messages. */

char *prog_name(char *s)
{
    char *p = s + strlen(s);

#ifdef MSDOS
    while (p >= s && *p != '\\' && *p != ':') {
	if (*p == '.')
	    *p = '\0';		/* strip off extension */
	if ('A' <= *p && *p <= 'Z')
	    *p += 'a' - 'A';	/* convert to lower case */
	p--;
    }
#else
    while (p >= s && *p != '/')
	p--;
#endif
    return (p+1);
}

/* help() prints a (very) concise summary of how to use this program.
   A more detailed summary is in the man page (gqpost.1). */

static char *help_strings[] = {
 "usage: %s -r RECORD [OPTIONS ...]\n",
 "where RECORD is the name of the record to be analyzed, and OPTIONS may",
 "include any of:",
 " -a ANNOTATOR read annotations from the specified ANNOTATOR (default: qrs)",
 " -c FILE     initialize parameters from the specified configuration FILE",
 " -f TIME     begin processing at specified time",
 " -h          print this usage summary",
 " -H          read multifrequency signals in high resolution mode",
 " -m THRESH   set interpolated event acceptance threshold to THRESH",
 "              (default: 1)",
 " -o ANNOTATOR write annotations to the specified ANNOTATOR (default: gqp)",
 " -t TIME     stop processing at specified time",
 "If too many true beats are rejected, decrease THRESH;  if too many false",
 "detections are accepted, increase THRESH.",
 "Note that the output is a complete copy of the input (with rejected events",
 "flagged as ARFCT).  The -f and -t options only limit the interval during",
 "which events may be rejected.",
NULL
};

void help()
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++)
	(void)fprintf(stderr, "%s\n", help_strings[i]);
}

/* gcalloc() is a wrapper for calloc() that handles errors and maintains
   a list of allocated buffers for automatic on-exit deallocation via
   cleanup(). */

size_t gcn = 0, gcnmax = 0;
void **gclist = NULL;

void *gcalloc(size_t nmemb, size_t size)
{
    void *p = calloc(nmemb, size);

    if ((p == NULL) ||
	((gcn >= gcnmax) &&
	 (gclist = realloc(gclist, (gcnmax += 32)*(sizeof(void *)))) == NULL))
	cleanup(3);	/* could not allocate requested memory */
    return (gclist[gcn++] = p);
}

void cleanup(int status)	/* close files and free allocated memory */
{
    if (status == 3)
	fprintf(stderr, "%s: insufficient memory (%d)\n", pname, status);
    if (record) wfdbquit();
    while (gcn-- > 0)
	if (gclist[gcn]) free(gclist[gcn]);
    exit(status);
}