inpstart.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
int dftFixPeak(DFT *dft, int verbose);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Verifies that peak is not missed in an input TAC file.",
  "Return code is zero, if no problems were encountered.",
  "This may be used in scripts for a rough test if venous plasma TAC can",
  "be used as input in Patlak analysis, but it must not replace visual",
  "inspection.",
  " ",
  "Usage: @P [Options] <TAC file(s)>",
  " ",
  "Options:",
  " -fix",
  "     Minor problems are corrected; always check the result before use.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: check that specified data files are suitable as input TACs",
  "     @P *vp.kbq",
  " ",
  "See also: tacframe, tactime, avgbolus, taccat, interpol",
  " ",
  "Keywords: TAC, input, plasma, peak, zero",
  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    ri, fileNr=0, ffi=0, ret, errorCount=0, doFix=0, fixed, checked=0;
  char  *cptr, tmp[256], dftfile[FILENAME_MAX];
  DFT    dft;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftfile[0]=(char)0;
  dftInit(&dft);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strcasecmp(cptr, "FIX")==0) {
      doFix=1; continue;
    }
    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(ffi==0) ffi=ai;
    fileNr++;
  }
 
  /* Is something missing? */
  if(fileNr==0) {
    fprintf(stderr, "Error: missing file name for input data.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("fileNr := %d\n", fileNr);
    printf("doFix := %d\n", doFix);
  }
 
 
  /*
   *  Do one file at a time
   */
  for(ai=ffi; ai<argc; ai++) {
    fixed=0;
    strlcpy(dftfile, argv[ai], FILENAME_MAX);
    if(verbose>0 && fileNr>1) fprintf(stdout, "%s:\n", dftfile);
 
    /* Read data */
    if(verbose>3) printf("reading %s\n", dftfile);
    ret=dftRead(dftfile, &dft);
    if(ret) {
      fprintf(stderr, "Error in reading '%s': %s\n", dftfile, dfterrmsg);
      return(2);
    }
 
    /* Make sure data is sorted by increasing sample time */
    dftSortByFrame(&dft);
 
    /* Check the data for missing samples */
    if(dft_nr_of_NA(&dft)>0) {
      if(verbose>0) printf("  missing sample(s).\n");
      if(doFix==0) {dftEmpty(&dft); errorCount++; continue;}
      /* Try to fix NaNs */
      if(verbose>1) printf("  trying to fix...\n");
      ret=dftNAfill(&dft);
      if(ret==0) {
        if(verbose>1) printf("  corrected.\n");
        fixed++;
      } else {
        fprintf(stderr,
          "Error: missing samples can not be corrected automatically.\n");
        dftEmpty(&dft); return(2);
      }
    }
 
    /* Check for the peak */
    if(dft.voiNr==1) { // just one TAC inside DFT
      checked++;
      ret=dftVerifyPeak(&dft, 0, verbose-2, tmp);
      if(ret==0) { // All is well
        if(verbose>0) printf("  ok\n");
        dftEmpty(&dft); continue;
      } else if(ret==1) {
        fprintf(stderr, "Error: program can not verify the peak.\n");
        if(verbose>0) fprintf(stderr, "Error: %s\n", tmp);
        dftEmpty(&dft); return(1);
      } else if(ret<0) { // Warning
        /* Different messages depending on whether it is to fixed or not */
        if(doFix==0) { // no correction
          if(verbose<=0) fprintf(stdout, "%s:\n", dftfile);
          printf("  TAC may have minor problems\n");
        } else { // correction will be tried
          if(verbose>0) printf("  TAC may have minor problems\n");
        }
        errorCount++;
      } else { // Problem was found
        /* Different messages depending on whether it is to fixed or not */
        if(doFix==0) { // no correction
          if(verbose<=0) fprintf(stdout, "%s:\n", dftfile);
          printf("  Error: TAC must be corrected before using as input\n");
        } else { // correction will be tried
          if(verbose>0)
            printf("  TAC must be corrected before using as input\n");
        }
        errorCount++;
      }
      if(doFix) {
        if(verbose>1) printf("  trying to fix...\n");
        ret=dftFixPeak(&dft, verbose-2);
        if(ret==0) { // correction succeeded
          if(verbose>0) printf("  corrected.\n");
          fixed++; errorCount--;
        } else { // could not be corrected
          if(verbose<=0) fprintf(stdout, "%s:\n", dftfile);
          printf("  Error: peak can not be corrected automatically.\n");
          printf("  Correct the input TAC manually!\n");
        }
      }
    } else { // several TACs inside DFT
 
      int localErrorCount=0;
      for(ri=0; ri<dft.voiNr; ri++) {
        if(verbose>0) fprintf(stdout, "  TAC %s:\n", dft.voi[ri].name);
        checked++;
        ret=dftVerifyPeak(&dft, ri, verbose-2, tmp);
        if(ret==0) { // All is well
          if(verbose>0) printf("    ok\n");
        } else if(ret==1) {
          fprintf(stderr, "Error: program can not verify the peak.\n");
          if(verbose>0) fprintf(stderr, "Error: %s\n", tmp);
          dftEmpty(&dft); return(1);
        } else if(ret<0) { // Warning
          if(verbose>0) printf("    TAC may have minor problems\n");
          localErrorCount++;
        } else { // Problem was found
          if(verbose>0)
            printf("    TAC must be corrected before using as input\n");
          localErrorCount++;
        }
      } // next TAC
      if(localErrorCount>0 && doFix) {
        if(verbose>1) printf("  trying to fix...\n");
        ret=dftFixPeak(&dft, verbose-2);
        if(ret<=0) {
          if(verbose>0) printf("  corrected.\n");
          fixed++; localErrorCount=0;
        } else {
          printf("  peaks can not be corrected automatically.\n");
        }
      }
      errorCount+=localErrorCount;
    }
 
    /* Write edited file */
    if(fixed>0) {
      if(verbose>2) printf("  writing corrected data...\n");
      ret=dftWrite(&dft, dftfile);
      if(ret!=0) {
        fprintf(stderr, "Error in writing '%s': %s\n", dftfile, dfterrmsg);
        dftEmpty(&dft);
        return(11);
      }
      if(verbose>1) printf("  corrected TAC written in file.\n");
    }
 
    /* Prepare to the next file */
    dftEmpty(&dft);
  } /* next file */
 
  if(verbose>0 && checked>1) {
    if(errorCount==0)
      printf("All files appear as suitable to be used as model input.\n");
    else
      printf("%d possibly severe problem(s) encountered.\n", errorCount);
  }
  return(errorCount);
}
/*****************************************************************************/
 
/*****************************************************************************/
 
int dftFixPeak(
  DFT *dft,
  int verbose
) {
  int ri, ret, mini, maxi, n, warn=0, fixnr=0;
  double minx, maxx, miny, maxy, dif, d, slope, ic;
  char tmp[256];
 
  if(verbose>0) printf("dftFixPeak(dft, %d)\n", verbose);
  /* Check the arguments */
  if(dft==NULL || dft->frameNr<1 || dft->voiNr<1) return 1;
 
  /* If less than 3 samples, then this can not be suitable as input TAC */
  if(dft->frameNr<3) return 2;
 
  /* Check if any correction is needed for any TAC */
  if(dftVerifyPeak(dft, -1, 0, NULL)==0) {
    if(verbose>1) printf("TAC does not need to be corrected\n");
    return 0;
  }
 
  /* Check if problem was something else than missing initial rise */
  if(dft->x[0]<=0.001) {
    if(verbose>1) printf("TAC already has zero sample\n");
    return 2;
  }
 
  /* Allocate memory for new sample points to be added */
  double zx, zy[dft->voiNr];
  zx=0.0; if(dft->timetype==DFT_TIME_STARTEND) zx=0.5*dft->x1[0];
 
  /* Go through all TACs */
  for(ri=0; ri<dft->voiNr; ri++) {
    zy[ri]=0.0;
    if(verbose>1 && dft->voiNr>1)
      printf("checking region %d: %s\n", 1+ri, dft->voi[ri].name);
 
    /* Check if any correction is needed for this TAC */
    if(dft->voiNr>1 && dftVerifyPeak(dft, ri, 0, NULL)==0) continue; // next TAC
 
    /* Get the TAC min and max */
    ret=dftMinMaxTAC(dft, ri, &minx, &maxx, &miny, &maxy, NULL, NULL,
                     &mini, &maxi);
    if(ret!=0) {
      if(verbose>0) printf("Error %d in dftMinMaxTAC()\n", ret);
      return 1;
    }
 
    /* If peak is not much higher than lowest value, that may indicate
       a problem */
    if(maxy<1.5*miny) {
      if(verbose>0) printf("TAC does not have a clear peak.\n");
      return 2;
    }
    if(maxy<5.0*miny) {
      if(verbose>1) printf("TAC does not have a clear peak.\n");
      warn++;
    }
 
    /* If the first sample is the peak:
       then add zero sample unless sample time is not too late */
    if(maxi==0) {
      if(verbose>1) printf("peak at the first sample.\n");
      if(dft->timetype==DFT_TIME_STARTEND) {
        dif=dft->x1[maxi]/(dft->x2[maxi]-dft->x1[maxi]);
      } else {
        dif=dft->x[maxi]/(dft->x[maxi+1]-dft->x[maxi]);
      }
      if(verbose>2) printf("dif := %g\n", dif);
      if(dif>5.0) {
        if(verbose>1)
          printf("peak at the first sample, which is too late.\n");
        return 2;
      } else if(dif>1.0 && maxy<=10.*miny) {
        if(verbose>1)
          printf("peak at the first sample and bad peak/tail ratio.\n");
        return 2;
      } 
      if(verbose>1) printf("zero sample added.\n");
      // default zeroes are good for this TAC
      fixnr++; continue;
    }
 
    /* If first sample is closer to peak and relatively high, then
       it is not safe to guess the start sample */
    if(dft->x[0]>0.75*maxx) {
      if(dft->voi[ri].y[0]>0.50*maxy) {
        if(verbose>1) printf("The first sample is relatively late and high.\n");
        return 2;
      }
    }
 
    /* Try to estimate zero point from ascending part */
    n=maxi+1; if(n>6) n/=2; else if(n>3) n--;
    ret=highest_slope(dft->x, dft->voi[ri].y, maxi+1, n,
                      &slope, &ic, &d, NULL);
    if(ret!=0 || d>=maxx) {
      if(verbose>0) printf("ascending part of TAC not available.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      return 2;
    }
    if(verbose>2) {
      printf("based on ascending part:\n");
      printf("  slope: %g\n", slope);
      printf("  ic: %g\n", ic);
      printf("  new x intercept: %g\n", d);
    }
    if(dft->timetype==DFT_TIME_STARTEND) {
      /* If start and end times, then set the y value */
      zy[ri]=slope*zx+ic;
      if(zy[ri]>0.5*dft->voi[ri].y[0]) zy[ri]=0.5*dft->voi[ri].y[0];
      else if(zy[ri]<0.0) zy[ri]=0.0;
    } else if(dft->voiNr>1) {
      /* If more than one TAC, then only set the y value */
      zy[ri]=slope*zx+ic;
      if(zy[ri]>0.5*dft->voi[ri].y[0]) zy[ri]=0.5*dft->voi[ri].y[0];
      else if(zy[ri]<0.0) zy[ri]=0.0;
    } else {
      /* Otherwise change the sample time */
      if(d<0.0) {
        if(ic<0.5*dft->voi[ri].y[0]) d=0.5*dft->x[0]; else d=0.0;
      }
      if(d>=dft->x[0]) d=0.5*dft->x[0];
      zx=d;
    }
    fixnr++;
 
  } // next TAC
 
  /* If TAC(s) were fixed, only then add the zero sample */
  if(fixnr>0) {
    ret=dftFillInitialGap(dft); if(ret!=0) return 1;
    dft->x[0]=zx;
    for(ri=0; ri<dft->voiNr; ri++) {
      dft->voi[ri].y[0]=zy[ri];
      if(verbose>1) printf("zero sample added at (%g,%g).\n", zx, zy[ri]);
    }
  }
 
  /* Check TACs are fine after corrections */
  ret=dftVerifyPeak(dft, -1, 0, tmp);
  if(verbose>0) {
    if(verbose>2) printf("dftVerifyPeak()=%d\n", ret);
    if(verbose>1) fprintf(stderr, "Note: %s\n", tmp);
    if(ret>1) printf("correction was not successful\n");
  }
  return(ret);
}
/*****************************************************************************/
 
/*****************************************************************************/