inpstart.c File Reference

Verifies that peak is not missed in an input TAC file, and optionally adds a guess of first zero sample. More...

Go to the source code of this file.

Functions
int	dftFixPeak (DFT *dft, int verbose)

Detailed Description

Verifies that peak is not missed in an input TAC file, and optionally adds a guess of first zero sample.

Copyright

(c) Turku PET Centre

Author

Vesa Oikonen

Definition in file inpstart.c.

Function Documentation

dftFixPeak()

int dftFixPeak	(	DFT *	dft,
		int	verbose )

Tries to fix the situation when specified (input) TAC does not seem to contain peak or it starts too late to get reliable estimate of AUC.

Returns

Returns 0 in case data could be fixed, < 0 if there still may be a problem, 1 in case of an error, and 2 if TAC can not be corrected automatically.

Parameters

dft	Pointer to TAC data which is fixed and verified.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 266 of file inpstart.c.

  {
  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);
}