dftcat.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include <string.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Catenate the TACs from the second TAC file to the first one, or to a new",
  "datafile. Sample times or correction for physical decay are not changed.",
  " ",
  "Deprecated. Please use taccat instead.",
  " ",
  "Usage: @P tacfile1 tacfile2 [catenated_file]",
  " ",
  "Options:",
  " -both | -first | -second | -cut=<time>",
  "     In case of overlapping samples, either samples from both (-both),",
  "     first (-first), or second (-second, default) TAC are saved in",
  "     combined file, or specified cut time is used.",
  "     Cut time can also be given in file with keys 'x' or 'time'.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P t455ap_pump.kbq t455ap_manual.kbq t455ap_combined.kbq",
  " ",
  "See also: taccat, taccut, tacadd, tactime, inpstart, taccross, tacframe",
  " ",
  "Keywords: TAC, input, blood, modelling, tool, simulation",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  int    ret, fi, voi=-1;
  int    useboth=0;
  int    overlap=0;
  int    more_reliable=1; /* In case of overlapping data, which one is 
         considered more reliable: 0=first, or 1=second curve */
  DFT    data1, data2, data3;
  char  *cptr, dfile1[FILENAME_MAX], dfile2[FILENAME_MAX], rfile[FILENAME_MAX];
  double cut_time=-1.0E+20;
  int    unit1, unit2;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile1[0]=dfile2[0]=rfile[0]=(char)0; 
  dftInit(&data1); dftInit(&data2); dftInit(&data3);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "BOTH", 1)==0) {
      useboth=1; cut_time=-1.0E+20; continue;
    } else if(strncasecmp(cptr, "FIRST", 2)==0) {
      more_reliable=0; useboth=0; continue;
    } else if(strncasecmp(cptr, "SECOND", 2)==0) {
      more_reliable=1; useboth=0; continue;
    } else if(strncasecmp(cptr, "CUT=", 4)==0) {
      cptr+=4;
      /* Try as a value first */
      if(atof_with_check(cptr, &cut_time)==0) {if(cut_time>0.0) continue;} 
      /* Then try as a IFT file */
      if(access(cptr, 0)!=-1) {
        IFT ift; iftInit(&ift); ret=iftRead(&ift, cptr, 1, 0);
        if(ret==0) {
          if(iftGetDoubleValue(&ift, 0, "x", &cut_time, 0)>=0 && cut_time>0.0) {
            iftEmpty(&ift); continue;}
          if(iftGetDoubleValue(&ift, 0, "time", &cut_time, 0)>=0 && cut_time>0.0) {
            iftEmpty(&ift); continue;}
          if(iftGetDoubleValue(&ift, 0, "cut", &cut_time, 0)>=0 && cut_time>0.0) {
            iftEmpty(&ift); continue;}
          iftEmpty(&ift);
        }
      }
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc; ai++) {
    if(!dfile1[0]) {strcpy(dfile1, argv[ai]); continue;}
    else if(!dfile2[0]) {strcpy(dfile2, argv[ai]); continue;}
    else if(!rfile[0]) {strcpy(rfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!dfile2[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n",
      argv[0]);
    return(1);
  }
  /* First file is also output file, if output file was not specified */
  if(!rfile[0]) strcpy(rfile, dfile1);
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile1 := %s\n", dfile1);
    printf("dfile2 := %s\n", dfile2);
    printf("rfile := %s\n", rfile);
    printf("useboth := %d\n", useboth);
    printf("more_reliable := %d\n", more_reliable);
    printf("cut_time := %g\n", cut_time);
  }
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", dfile1);
  if(dftRead(dfile1, &data1)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile1, dfterrmsg);
    return(2);
  }
  if(verbose>2) {
    printf("voiNr1 := %d\n", data1.voiNr);
    printf("frameNr1 := %d\n", data1.frameNr);
  }
  dftSortByFrame(&data1);
 
  if(verbose>1) printf("reading %s\n", dfile2);
  if(dftRead(dfile2, &data2)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile2, dfterrmsg);
    return(2);
  }
  if(verbose>2) {
    printf("voiNr2 := %d\n", data2.voiNr);
    printf("frameNr2 := %d\n", data2.frameNr);
  }
  dftSortByFrame(&data2);
 
  /* Check that files have equal nr of regions */
  if(data1.voiNr!=data2.voiNr) {
    fprintf(stderr, "Error: different nr of regions.\n");
    dftEmpty(&data1); dftEmpty(&data2); return(3);
  }
 
  /* Check that files have the same time unit */
  if(data1.timeunit==TUNIT_UNKNOWN || data2.timeunit==TUNIT_UNKNOWN) {
    fprintf(stderr, "Warning: unknown time units.\n");
    if(data1.timeunit!=TUNIT_UNKNOWN) data2.timeunit=data1.timeunit;
    else if(data2.timeunit!=TUNIT_UNKNOWN) data1.timeunit=data2.timeunit;
  } else {
    if(data1.timeunit!=data2.timeunit) {
      fprintf(stderr, "Warning: different time units.\n");
      if(verbose>0)
        fprintf(stdout, "  converting units '%s' to '%s'\n",
                petTunit(data2.timeunit), petTunit(data1.timeunit));
      ret=dftTimeunitConversion(&data2, data1.timeunit);
      if(ret!=0) {
        fprintf(stderr, "Error: cannot convert time units.\n");
        dftEmpty(&data1); dftEmpty(&data2); return(3);
      }
    }
  }
 
  /* Check that files have the same concentration unit */
  unit1=dftUnitId(data1.unit);
  unit2=dftUnitId(data2.unit);
  if(unit1==CUNIT_UNKNOWN || unit2==CUNIT_UNKNOWN) {
    fprintf(stderr, "Warning: unknown concentration units.\n");
    if(unit1!=CUNIT_UNKNOWN) strcpy(data2.unit, data1.unit);
    else if(unit2!=CUNIT_UNKNOWN) strcpy(data1.unit, data2.unit);
  } else {
    if(unit1!=unit2) {
      fprintf(stderr, "Warning: different concentration units.\n");
      if(verbose>0)
        fprintf(stdout, "  converting units '%s' to '%s'\n",
                petCunit(unit2), petCunit(unit1));
      ret=dftUnitConversion(&data2, unit1);
      if(ret!=0) {
        fprintf(stderr, "Error: cannot convert concentration units.\n");
        dftEmpty(&data1); dftEmpty(&data2); return(3);
      }
    }
  }
 
  /* Deal with possible overlap */
  /* Check that start times do not overlap, unless delt with options */
  ret=0;
  if(data1.timetype==DFT_TIME_STARTEND) {
    if(data1.x1[0]>data2.x[0] && useboth==0 && cut_time<0.0) ret=1;
  } else if(data2.timetype==DFT_TIME_STARTEND) {
    if(data1.x[0]>data2.x2[0] && useboth==0 && cut_time<0.0) ret=2;
  } else {
    if(data1.x[0]>data2.x[0] && useboth==0 && cut_time<0.0) ret=3;
  }
  if(ret) {
    fprintf(stderr, "Error: check datafile times and their order!\n");
    dftEmpty(&data1); dftEmpty(&data2); return(4);
  }
  if(cut_time>=0.0 && data1.x[data1.frameNr-1]>=data2.x[0]) {
    /* use user-given cut time to solve overlap */
    if(data1.x[data1.frameNr-1]>=cut_time && data2.x[0]<cut_time) {
      /* It is necessary to cut from both tacs */
      if(verbose>1) printf("cutting both TACs at %g\n", cut_time);
      ret=dftRemoveTimeRange(&data1, data1.x[0]-1.0, cut_time);
      if(ret==0) ret=dftRemoveTimeRange(&data2, cut_time,
                                        data2.x[data2.frameNr-1]+1.0);
      if(ret!=0) {
        fprintf(stderr, "Error: invalid cut time.\n");
        dftEmpty(&data1); dftEmpty(&data2); return(4);
      }
    } else if(data2.x[0]<cut_time) {
      /* 2nd tac will be cut later */
      more_reliable=0;
    } else if(data1.x[data1.frameNr-1]>=cut_time) {
      /* 1st tac will be cut later */
      more_reliable=1;
    }
  }
  if(data1.x[data1.frameNr-1]>=data2.x[0]) {
    overlap=1;
    if(verbose>0 && cut_time<0.0 && useboth==0)
      fprintf(stderr, "Warning: times in datafiles overlap.\n");
    if(useboth==0) {
      /* Overlapping data is not allowed */
      if(more_reliable==1) { // second tac is more reliable
        ret=dftRemoveTimeRange(&data1, data1.x[0]-1.0, data2.x[0]-1.0E-10);
        if(ret!=0 || data1.frameNr==0) {
          fprintf(stderr, "Error: check datafile times and their order!\n");
          dftEmpty(&data1); dftEmpty(&data2); return(4);
        }
      } else { // first tac is more reliable
        /* Remove the overlapping part from the second tac */
        ret=dftRemoveTimeRange(&data2, data1.x[data1.frameNr-1]+1.0E-10,
                                       data2.x[data2.frameNr-1]+1.0);
        if(ret!=0) {
          fprintf(stderr, "Error: check datafile times and their order!\n");
          dftEmpty(&data1); dftEmpty(&data2); return(4);
        }
        if(data2.frameNr<1)
          if(verbose>0) fprintf(stderr, "Warning: datafiles overlap fully.\n");
      }
    }
    if(verbose>2) {
      printf("first_tac_range := %g - %g\n", 
        data1.x[0], data1.x[data1.frameNr-1]);
      printf("second_tac_range := %g - %g\n", 
        data2.x[0], data2.x[data2.frameNr-1]);
    }
  }
 
 
  /*
   *  Make catenated data
   */
  /* Allocate memory for it */
  if(dftSetmem(&data3, data1.frameNr+data2.frameNr, data1.voiNr)){
    fprintf(stderr, "Error: cannot allocate memory.\n");
    dftEmpty(&data1); dftEmpty(&data2); return(5);
  }
  data3.frameNr=data1.frameNr+data2.frameNr;
  data3.voiNr=data1.voiNr;
  /* Set headers */
  (void)dftCopymainhdr2(&data1, &data3, 1);
  (void)dftCopymainhdr2(&data2, &data3, 0); // only missing fields
  for(voi=0; voi<data1.voiNr; voi++)
    (void)dftCopyvoihdr(&data1, voi, &data3, voi);
  if(strlen(data3.studynr)==0) strcpy(data3.studynr, data2.studynr);
  if(strlen(data3.unit)==0) strcpy(data3.unit, data2.unit);
  if(data3.timeunit==TUNIT_UNKNOWN) data3.timeunit=data2.timeunit;
  if(data3._type!=data2._type) {
    if(data2._type==DFT_FORMAT_STANDARD) data3._type=data2._type;
  }
  if(data2.timetype==DFT_TIME_MIDDLE) data3.timetype=DFT_TIME_MIDDLE;
  /* Copy data from first file */
  for(fi=0; fi<data1.frameNr; fi++) {
    data3.x[fi]=data1.x[fi]; 
    data3.x1[fi]=data1.x1[fi]; data3.x2[fi]=data1.x2[fi];
    for(voi=0; voi<data1.voiNr; voi++)
      data3.voi[voi].y[fi]=data1.voi[voi].y[fi];
  }
  /* Copy data from second file */
  for(fi=0; fi<data2.frameNr; fi++) {
    data3.x[data1.frameNr+fi]=data2.x[fi];
    data3.x1[data1.frameNr+fi]=data2.x1[fi];
    data3.x2[data1.frameNr+fi]=data2.x2[fi];
    for(voi=0; voi<data1.voiNr; voi++)
      data3.voi[voi].y[data1.frameNr+fi]=data2.voi[voi].y[fi];
  }
  /* Turn weights off */
  data3.isweight=0;
  /* Sort data by increasing sample time, if necessary */
  if(overlap && useboth) ret=dftSortByFrame(&data3);
 
  /* Add comment */
  dftSetComments(&data3);
  //strcpy(data3.comments, data.comments);
  //strcat(data3.comments, data2.comments);
  strcat(data3.comments, "# catenated\n");
 
 
  /* Save data */
  ret=dftWrite(&data3, rfile);
  if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, rfile, dfterrmsg);
    dftEmpty(&data1); dftEmpty(&data2); dftEmpty(&data3);
    return(11);
  }
  if(verbose>0)
    printf("  %d+%d=%d frames written in %s\n", data1.frameNr, data2.frameNr,
                                                data3.frameNr, rfile);
 
  /* Free memory */
  dftEmpty(&data1); dftEmpty(&data2); dftEmpty(&data3);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/