Header file for libtpcmodext. More...

#include "tpcclibConfig.h"
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpcsvg.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"

Functions
int	bf_srtm (double t, double cr, int n, int bfNr, double t3min, double t3max, DFT *bf)

int	bfRadiowater (DFT input, DFT tissue, DFT bf, int bfNr, double k2min, double k2max, char status, int verbose)

int	bfIrr2TCM (DFT input, DFT tissue, DFT bf, int bfNr, double thetamin, double thetamax, char status, int verbose)

int	dftReadinput (DFT input, DFT tissue, char filename, int filetype, double ti1, double ti2, int verifypeak, char *status, int verbose)

int	dftReadReference (DFT tissue, char filename, int filetype, int ref_index, char *status, int verbose)

int	dftReadModelingData (char tissuefile, char inputfile1, char inputfile2, char inputfile3, double fitdur, int fitframeNr, DFT tis, DFT inp, FILE loginfo, int verbose, char status)

int	dftRobustMinMaxTAC (DFT dft, int tacindex, double minx, double maxx, double miny, double maxy, int mini, int maxi, int mins, int *maxs, int verbose)

int	dftVerifyPeak (DFT dft, int index, int verbose, char status)

int	imgReadModelingData (char petfile, char siffile, char inputfile1, char inputfile2, char inputfile3, double fitdur, int fitframeNr, IMG img, DFT inp, DFT iinp, int verifypeak, FILE loginfo, int verbose, char status)

int	extrapolate_monoexp (DFT dft, double fittime, int min_nr, int max_nr, double mintime, double extr_to, DFT ext, FILE loginfo, char status)

int	dftAutointerpolate (DFT dft, DFT dft2, double endtime, int verbose)

int	dftDoubleFrames (DFT dft, DFT dft2)

int	dftDivideFrames (DFT dft, int voi_index, int add_nr, DFT dft2, int verbose)

int	dft_end_line (DFT dft, double fittime, int min_nr, int max_nr, double mintime, int check_impr, FIT fit, FILE loginfo, char status)

int	dft_ln (DFT dft1, DFT dft2)

int	fittime_from_dft (DFT dft, double startTime, double endTime, int first, int *last, int verbose)

int	fittime_from_img (IMG img, double fittime, int verbose)

int	check_times_dft_vs_img (IMG img, DFT dft, int verbose)

int	check_times_dft_vs_dft (DFT dft1, DFT dft2, int verbose)

int	copy_times_from_img_to_dft (IMG img, DFT dft, int verbose)

int	getActualSamplenr (DFT *dft, int ri)

double	dftEndtime (DFT *dft)

double	imgEndtime (IMG *img)

int	dftMatchTimeunits (DFT dft1, DFT dft2, int *tunit2, int verbose)

int	img_patlak (DFT input, IMG dyn_img, int start, int end, linefit_range fit_range, float thrs, IMG ki_img, IMG ic_img, IMG nr_img, char status, int verbose)

int	img_logan (DFT input, IMG dyn_img, int start, int end, linefit_range fit_range, float thrs, double k2, IMG vt_img, IMG ic_img, IMG nr_img, char status, int verbose)

int	img_k1_using_ki (DFT input, IMG dyn_img, int frame_nr, IMG ki_img, IMG k1_img, IMG k2k3_img, char status, int verbose)

int	dftInterpolateCheckStart (DFT input, DFT output, char *status, int verbose)

int	dftInterpolateCheckEnd (DFT input, DFT output, char *status, int verbose)

int	dftInterpolate (DFT input, DFT tissue, DFT output, char status, int verbose)

int	dftInterpolateInto (DFT inp, DFT tis, char *status, int verbose)

int	dftTimeIntegral (DFT dft, double t1, double t2, DFT idft, int calc_mode, char *status, int verbose)

int	dftDerivative (DFT dft, DFT deriv, char *status)

int	dftDerivative_old (DFT dft, DFT deriv, char *status)

int	dftInterpolateForIMG (DFT input, IMG img, int frame_nr, DFT output, double ti1, double ti2, int verbose, char status)

int	imgTimeIntegral (IMG img, float t1, float t2, IMG iimg, int calc_mode, char *status, int verbose)

int	dftAllocateWithIMG (DFT dft, int tacNr, IMG img)

int	sif2dft (SIF sif, DFT dft)

int	sifAllocateWithIMG (SIF sif, IMG img, int doCounts, int verbose)

int	mrl_between_tacs (double y1, double y2, int n)

int	noiseSD4Simulation (double y, double t1, double dt, double hl, double a, double sd, char status, int verbose)

int	noiseSD4SimulationFromDFT (DFT dft, int index, double pc, double sd, char *status, int verbose)

int	plot_svg (DFT dft, RES res, int first, int last, char main_title, char x_title, char y_title, int color_scale, char fname, int verbose)

int	plotdata (DFT dft, RES res, int first, int last, char mtitle, char xtitle, char ytitle, char fname)

int	plotdata_as_dft (DFT dft, char fname)

int	plot_fit_svg (DFT dft1, DFT dft2, char main_title, char fname, int verbose)

int	plot_fitrange_svg (DFT dft1, DFT dft2, char main_title, double x1, double x2, double y1, double y2, char fname, int verbose)

int	cunit_check_dft_vs_img (DFT dft, IMG img, char *errmsg, int verbose)

int	dftWeightByFreq (DFT *dft)

int	imgSetWeights (IMG *img, int wmet, int verbose)

int	dftWSampleNr (DFT *tac)

int	clusterTACs (IMG dimg, IMG cimg, int nr, DFT *tac, int verbose)

Detailed Description

Header file for libtpcmodext.

Author: Vesa Oikonen

Definition in file libtpcmodext.h.

Function Documentation

◆ bf_srtm()

int bf_srtm	(	double *	t,
		double *	cr,
		int	n,
		int	bfNr,
		double	t3min,
		double	t3max,
		DFT *	bf )

Calculates set of basis functions for SRTM.

Returns: Returns 0 if successful, otherwise non-zero.

See also: bfRadiowater, bfIrr2TCM

Parameters

t	PET frame mid times.
cr	Non-decay corrected Cr(t).
n	Nr of PET frames.
bfNr	Nr of basis functions to calculate.
t3min	theta3 min.
t3max	theta3 max.
bf	data for basis functions is allocated and filled here.

Definition at line 16 of file bf_model.c.

  {
  int bi, fi, ret;
  double a, b, c;
 
  // printf("\n%s(*t, *cr, %d, %d, %g, %g, *bf)\n", __func__, n, bfNr, t3min, t3max);
 
  /* Check the parameters */
  if(t==NULL || cr==NULL || n<2 || bfNr<1 || t3min<1.0E-10 || t3min>=t3max) return(1);
  if(bf==NULL || bf->voiNr>0) return(1);
  
  /* Allocate meory for basis functions */
  ret=dftSetmem(bf, n, bfNr); if(ret) return(2);
 
  /* Copy and set information fields */
  bf->voiNr=bfNr; bf->frameNr=n;
  bf->_type=DFT_FORMAT_STANDARD;
  for(bi=0; bi<bf->voiNr; bi++) {
    snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
    strcpy(bf->voi[bi].hemisphere, ".");
    strcpy(bf->voi[bi].place, ".");
    strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
  }
  for(fi=0; fi<bf->frameNr; fi++) bf->x[fi]=t[fi];
 
  /* Compute theta3 values to size fields */
  a=log10(t3min); b=log10(t3max); c=(b-a)/(double)(bfNr-1);
  for(bi=0; bi<bf->voiNr; bi++) {
    bf->voi[bi].size=pow(10.0, (double)bi*c+a);
  }
  
  /* Calculate the functions */
  for(bi=0; bi<bf->voiNr; bi++) {
    a=bf->voi[bi].size;
    ret=simC1(t, cr, n, 1.0, a, bf->voi[bi].y);
    if(ret) return(4);
  }
 
  return(0);
}

◆ bfIrr2TCM()

int bfIrr2TCM	(	DFT *	input,
		DFT *	tissue,
		DFT *	bf,
		int	bfNr,
		double	thetamin,
		double	thetamax,
		char *	status,
		int	verbose )

Calculates set of basis functions for irreversible 2TCM.

Returns: Returns 0 if successful, otherwise non-zero.

See also: bf_srtm, bfRadiowater

Parameters

input	Arterial PTAC (not modified).
tissue	PET TTAC (not modified, just to get frame times).
bf	Place for basis functions (initiated DFT struct, allocated and filled here).
bfNr	Nr of basis functions to calculate.
thetamin	Minimum of theta=k2+k3 (sec-1 or min-1, corresponding to TAC time units).
thetamax	Maximum of theta=k2+k3 (sec-1 or min-1, corresponding to TAC time units).
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 217 of file bf_model.c.

  {
  int bi, fi, ret;
  double a, b, c;
 
  if(verbose>0)
    printf("\n%s(*inp, *tis, *bf, %d, %g, %g, status, %d)\n", __func__,
           bfNr, thetamin, thetamax, verbose);
 
  /* Check the parameters */
  if(input==NULL || tissue==NULL || bf==NULL) {
    if(status!=NULL) strcpy(status, "program error");
    else if(verbose>0) fprintf(stderr, "invalid function parameters\n");
    return(1);
  }
  if(input->frameNr<3 || input->voiNr<1) {
    if(status!=NULL) strcpy(status, "no input data");
    else if(verbose>0) fprintf(stderr, "invalid input data\n");
    return(2);
  }
  if(tissue->frameNr<1) {
    if(status!=NULL) strcpy(status, "no pet data");
    else if(verbose>0) fprintf(stderr, "invalid PET data\n");
    return(3);
  }
  if(input->timeunit!=tissue->timeunit) {
    if(status!=NULL) strcpy(status, "invalid time units");
    else if(verbose>0) fprintf(stderr, "invalid time units\n");
    return(4);
  }
  if(bfNr<2) {
    if(status!=NULL) strcpy(status, "invalid nr of basis functions");
    else if(verbose>0) fprintf(stderr, "invalid number of basis functions\n");
    return(5);
  }
  if(thetamin<0.0) thetamin=0.0;
  if(thetamin>=thetamax) {
    if(status!=NULL) strcpy(status, "invalid theta range");
    else if(verbose>0) fprintf(stderr, "invalid theta range\n");
    return(6);
  }
  if(verbose>1) {
    printf("input timerange: %g - %g\n", input->x[0], input->x[input->frameNr-1]);
    printf("tissue timerange: %g - %g\n", tissue->x[0], tissue->x[tissue->frameNr-1]);
  }
  
  /* Allocate memory for basis functions */
  if(verbose>1) printf("allocating memory for basis functions\n");
  ret=dftSetmem(bf, tissue->frameNr, bfNr);
  if(ret) {
    if(status!=NULL) strcpy(status, "out of memory");
    else if(verbose>0) fprintf(stderr, "out of memory\n");
    return(10);
  }
 
  /* Copy and set information fields */
  bf->voiNr=bfNr; bf->frameNr=tissue->frameNr;
  bf->_type=tissue->_type;
  dftCopymainhdr2(tissue, bf, 1);
  for(bi=0; bi<bf->voiNr; bi++) {
    snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
    strcpy(bf->voi[bi].hemisphere, ".");
    strcpy(bf->voi[bi].place, ".");
    strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
  }
  for(fi=0; fi<bf->frameNr; fi++) {
    bf->x[fi]=tissue->x[fi];
    bf->x1[fi]=tissue->x1[fi];
    bf->x2[fi]=tissue->x2[fi];
  }
  
  /* Compute the range of theta values to size fields */
  if(verbose>1) printf("computing theta values\n");
  a=thetamin; b=thetamax; c=(b-a)/(double)(bfNr-1);
  if(verbose>20) printf("a=%g b=%g, c=%g\n", a, b, c);
  for(bi=0; bi<bf->voiNr; bi++) bf->voi[bi].size=(double)bi*c+a;
  if(verbose>2) {
    printf("final BF theta range: %g - %g\n", bf->voi[0].size, bf->voi[bf->voiNr-1].size);
    printf("theta step size: %g\n", c);
  }
  
  /* Allocate memory for simulated TAC */
  double *sim;
  sim=(double*)malloc(input->frameNr*sizeof(double));
  if(sim==NULL) {
    if(status!=NULL) strcpy(status, "out of memory");
    else if(verbose>0) fprintf(stderr, "out of memory\n");
    dftEmpty(bf); return(11);
  }
  
  /* Calculate the basis functions at input time points */
  if(verbose>1) printf("computing basis functions at input sample times\n");
  for(bi=0; bi<bf->voiNr; bi++) {
    a=bf->voi[bi].size;
    ret=simC1(input->x, input->voi[0].y, input->frameNr, 1.0, a, sim);
    if(ret) {
      if(status!=NULL) strcpy(status, "simulation problem");
      else if(verbose>0) fprintf(stderr, "simulation problem\n");
      free(sim); dftEmpty(bf); return(20);
    }
    if(verbose>100) {
      printf("\ntheta := %g\n", a);
      printf("simulated TAC:\n");
      for(fi=0; fi<input->frameNr; fi++)
        printf("  %12.6f  %12.3f\n", input->x[fi], sim[fi]);
    }
    /* interpolate to PET time frames */
    if(tissue->timetype==DFT_TIME_STARTEND)
      ret=interpolate4pet(input->x, sim, input->frameNr, tissue->x1, tissue->x2,
                          bf->voi[bi].y, NULL, NULL, bf->frameNr);
    else
      ret=interpolate(input->x, sim, input->frameNr, tissue->x,
                      bf->voi[bi].y, NULL, NULL, bf->frameNr);
    if(ret) {
      if(status!=NULL) strcpy(status, "simulation problem");
      else if(verbose>0) fprintf(stderr, "simulation problem\n");
      free(sim); dftEmpty(bf); return(20);
    }
 
  } // next basis function
 
  free(sim);
  if(verbose>1) printf("%s() done.\n\n", __func__);
  if(status!=NULL) strcpy(status, "ok");
  return(0);
}

◆ bfRadiowater()

int bfRadiowater	(	DFT *	input,
		DFT *	tissue,
		DFT *	bf,
		int	bfNr,
		double	k2min,
		double	k2max,
		char *	status,
		int	verbose )

Calculates set of basis functions for generic radiowater model.

Returns: Returns 0 if successful, otherwise non-zero.

See also: bf_srtm, bfIrr2TCM

Parameters

input	Arterial blood input TAC (not modified).
tissue	PET TACs (not modified, just to get frame times).
bf	Place for basis functions (initiated DFT struct, allocated and filled here).
bfNr	Nr of basis functions to calculate.
k2min	Minimum of k2 (sec-1 or min-1, corresponding to TAC time units).
k2max	Maximum of k2 (sec-1 or min-1, corresponding to TAC time units).
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 77 of file bf_model.c.

  {
  if(verbose>0)
    printf("\n%s(*inp, *tis, *bf, %d, %g, %g, status, %d)\n", __func__,
           bfNr, k2min, k2max, verbose);
 
  /* Check the parameters */
  if(input==NULL || tissue==NULL || bf==NULL) {
    if(status!=NULL) strcpy(status, "program error");
    return 1;
  }
  if(input->frameNr<3 || input->voiNr<1) {
    if(status!=NULL) strcpy(status, "no input data");
    return 2;
  }
  if(tissue->frameNr<1) {
    if(status!=NULL) strcpy(status, "no pet data");
    return 3;
  }
  if(input->timeunit!=tissue->timeunit) {
    if(status!=NULL) strcpy(status, "invalid time units");
    return 4;
  }
  if(bfNr<2) {
    if(status!=NULL) strcpy(status, "invalid nr of basis functions");
    return 5;
  }
  if(k2min<1.0E-10) k2min=1.0E-10; // range calculation does not work otherwise
  if(k2min>=k2max || k2min<0.0) {
    if(status!=NULL) strcpy(status, "invalid k2 range");
    return 6;
  }
  if(verbose>1) {
    printf("input timerange: %g - %g\n", input->x[0], input->x[input->frameNr-1]);
    printf("tissue timerange: %g - %g\n", tissue->x[0], tissue->x[tissue->frameNr-1]);
  }
  
  /* Allocate memory for basis functions */
  if(verbose>1) printf("allocating memory for basis functions\n");
  if(dftSetmem(bf, tissue->frameNr, bfNr)) {
    if(status!=NULL) strcpy(status, "out of memory");
    return 10;
  }
 
  /* Copy and set information fields */
  bf->voiNr=bfNr; bf->frameNr=tissue->frameNr;
  bf->_type=tissue->_type;
  dftCopymainhdr2(tissue, bf, 1);
  for(int bi=0; bi<bf->voiNr; bi++) {
    snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
    strcpy(bf->voi[bi].hemisphere, ".");
    strcpy(bf->voi[bi].place, ".");
    strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
  }
  for(int fi=0; fi<bf->frameNr; fi++) {
    bf->x[fi]=tissue->x[fi];
    bf->x1[fi]=tissue->x1[fi];
    bf->x2[fi]=tissue->x2[fi];
  }
  
  /* Compute the range of k2 values to size fields */
  if(verbose>1) printf("computing k2 values\n");
  double a, b, c;
  a=log10(k2min); b=log10(k2max); c=(b-a)/(double)(bfNr-1);
  if(verbose>20) printf("a=%g b=%g, c=%g\n", a, b, c);
  for(int bi=0; bi<bf->voiNr; bi++) {
    bf->voi[bi].size=pow(10.0, (double)bi*c+a);
  }
  if(verbose>2) {
    printf("final BF k2 range: %g - %g\n", bf->voi[0].size, bf->voi[bf->voiNr-1].size);
  }
  
  /* Allocate memory for simulated TAC */
  double *sim;
  sim=(double*)malloc(input->frameNr*sizeof(double));
  if(sim==NULL) {
    if(status!=NULL) strcpy(status, "out of memory");
    dftEmpty(bf); return 11;
  }
  
  /* Calculate the basis functions at input time points */
  if(verbose>1) printf("computing basis functions at input sample times\n");
  for(int bi=0; bi<bf->voiNr; bi++) {
    a=bf->voi[bi].size;
    if(simC1(input->x, input->voi[0].y, input->frameNr, 1.0, a, sim)) {
      if(status!=NULL) strcpy(status, "simulation problem");
      free(sim); dftEmpty(bf);
      return(20);
    }
    if(verbose>100) {
      printf("\nk2 := %g\n", a);
      printf("simulated TAC:\n");
      for(int fi=0; fi<input->frameNr; fi++)
        printf("  %12.6f  %12.3f\n", input->x[fi], sim[fi]);
    }
    /* interpolate to PET time frames */
    int ret=0;
    if(tissue->timetype==DFT_TIME_STARTEND)
      ret=interpolate4pet(input->x, sim, input->frameNr, tissue->x1, tissue->x2,
                          bf->voi[bi].y, NULL, NULL, bf->frameNr);
    else
      ret=interpolate(input->x, sim, input->frameNr, tissue->x,
                      bf->voi[bi].y, NULL, NULL, bf->frameNr);
    if(ret) {
      if(status!=NULL) strcpy(status, "simulation problem");
      free(sim); dftEmpty(bf);
      return(20);
    }
 
  } // next basis function
 
  free(sim);
  if(verbose>1) printf("%s() done.\n\n", __func__);
  if(status!=NULL) strcpy(status, "ok");
  return(0);
}

◆ check_times_dft_vs_dft()

int check_times_dft_vs_dft	(	DFT *	dft1,
		DFT *	dft2,
		int	verbose )

Check whether sample times are the same (or very close to) in two DFT structs. Data sets are not edited. If DFT structs contain different sample number, then only common nr of samples are compared.

See also: check_times_dft_vs_img, copy_times_from_img_to_dft

Returns: Returns 0 in case of no match, 1 if times do match, and <1 in case of an error.

Parameters

dft1	Pointer to first DFT data; times can be both seconds or minutes, if unit is correctly set.
dft2	Pointer to second DFT data; times can be both seconds or minutes, if unit is correctly set; sample nr may be different than in dft1.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 169 of file fittime.c.

  {
  int fi, n, smaller_frameNr=0;
  double f, ts1, ts2;
  double accepted_timedif=2.2; // sec
 
  if(verbose>0) printf("%s(*img, *dft)\n", __func__);
  if(dft1==NULL || dft2==NULL) return -1;
  if(dft1->frameNr<=0 || dft2->frameNr<=0) return 0;
  
  /* Which has less samples? */
  if(dft1->frameNr<dft2->frameNr) smaller_frameNr=dft1->frameNr;
  else smaller_frameNr=dft2->frameNr;
 
  /* Convert times to sec if necessary (and possible) */
  ts1=ts2=1.0;
  if(dft1->timeunit!=TUNIT_UNKNOWN && dft1->timeunit!=TUNIT_UNKNOWN) {
    if(dft1->timeunit==TUNIT_MIN) ts1=60.0;
    if(dft2->timeunit==TUNIT_MIN) ts2=60.0;
  }
  if(verbose>1) {
    printf("dft1->timetype := %d\n", dft1->timetype);
    printf("dft2->timetype := %d\n", dft2->timetype);
    if(verbose>2) {
      printf("time range 1 := %g - %g %s\n", dft1->x[0],
             dft1->x[dft1->frameNr-1], petTunit(dft1->timeunit));
      printf("time range 2 := %g - %g %s\n", dft2->x[0],
             dft2->x[dft2->frameNr-1], petTunit(dft2->timeunit));
    }
  }
 
  /* With short study the accepted time difference must be shorter */
  f=0.01*dft1->x2[dft1->frameNr-1]*ts1;
  if(accepted_timedif>f) accepted_timedif=f;
  if(verbose>1) printf("accepted_timedif := %g [s]\n", accepted_timedif);
 
  /* Compare sample times frame-by-frame */
  for(fi=0, n=0; fi<smaller_frameNr; fi++) {
    if(dft1->timetype==DFT_TIME_MIDDLE) { // check frame mid times
      f=fabs(dft1->x[fi]*ts1-dft2->x[fi]*ts2);
      if(verbose>10) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>12) printf("  %g vs %g\n", ts1*dft1->x[fi], ts2*dft2->x[fi]);
      if(f>accepted_timedif) n++;
      continue;
    }
    if(dft1->timetype==DFT_TIME_START || dft1->timetype==DFT_TIME_STARTEND) {
      f=fabs(dft1->x1[fi]*ts1-dft2->x1[fi]*ts2);
      if(verbose>10) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>12) printf("  %g vs %g\n", ts1*dft1->x1[fi], ts2*dft2->x1[fi]);
      if(f>accepted_timedif) {n++; continue;}
    }
    if(dft1->timetype==DFT_TIME_END || dft1->timetype==DFT_TIME_STARTEND) {
      f=fabs(dft1->x2[fi]*ts1-dft2->x2[fi]*ts2);
      if(verbose>10) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>12) printf("  %g vs %g\n", ts1*dft1->x2[fi], ts2*dft2->x2[fi]);
      if(f>accepted_timedif) n++;
    }
  }
  if(verbose>2) printf("nr of different frame times := %d\n", n);
 
  if(n==0) return 1; else return 0; 
} 

Referenced by dftInterpolate(), dftInterpolateForIMG(), and dftInterpolateInto().

◆ check_times_dft_vs_img()

int check_times_dft_vs_img	(	IMG *	img,
		DFT *	dft,
		int	verbose )

Check whether DFT sample times are the same (or very close to) as the frame times in IMG. This would suggest that DFT data originates from the same or similar PET scan. Specified data sets are not edited. If frame nr is different, then only the common frames are compared.

See also: check_times_dft_vs_dft, copy_times_from_img_to_dft, dftEndtime, imgEndtime

Returns: Returns 0 in case of no match, 1 if times do match, and <1 in case of an error.

Parameters

img	Pointer to IMG data; times must be in sec as usual.
dft	Pointer to DFT data; times can be both seconds or minutes, if unit is correctly set; sample nr may be different than in IMG.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 109 of file fittime.c.

  {
  int smaller_dimt=0;
  double f, ts;
  double accepted_timedif=2.2; // sec
 
  if(verbose>0) printf("%s(*img, *dft)\n", __func__);
  if(img==NULL || dft==NULL) return -1;
  if(img->dimt<=0 || dft->frameNr<=0) return 0;
  
  /* Get the smaller frame nr */
  if(img->dimt<dft->frameNr) smaller_dimt=img->dimt; else smaller_dimt=dft->frameNr;
 
  /* With short study the accepted time difference must be shorter */
  f=0.01*img->end[img->dimt-1]; if(accepted_timedif>f) accepted_timedif=f;
  if(verbose>1) printf("accepted_timedif := %g [s]\n", accepted_timedif);
 
  /* Convert DFT times to sec if necessary */
  if(dft->timeunit==TUNIT_MIN) ts=60.0; else ts=1.0;
 
  /* Compare sample times frame-by-frame */
  int n=0;
  for(int fi=0; fi<smaller_dimt; fi++) {
    if(dft->timetype==DFT_TIME_MIDDLE) { // check frame mid times
      f=fabs(img->mid[fi]-dft->x[fi]*ts);
      if(verbose>10) printf("timedif[%d] := %g\n", fi, f);
      if(f>accepted_timedif) n++;
      continue;
    }
    if(dft->timetype==DFT_TIME_START || dft->timetype==DFT_TIME_STARTEND) {
      f=fabs(img->start[fi]-dft->x1[fi]*ts);
      if(verbose>10) printf("timedif[%d] := %g\n", fi, f);
      if(f>accepted_timedif) {n++; continue;}
    }
    if(dft->timetype==DFT_TIME_END || dft->timetype==DFT_TIME_STARTEND) {
      f=fabs(img->end[fi]-dft->x2[fi]*ts);
      if(verbose>10) printf("timedif[%d] := %g\n", fi, f);
      if(f>accepted_timedif) n++;
    }
  }
  if(verbose>2) printf("nr of different frame times := %d\n", n);
 
  if(n==0) return 1; else return 0; 
} 

Referenced by img_logan(), and img_patlak().

◆ clusterTACs()

int clusterTACs	(	IMG *	dimg,
		IMG *	cimg,
		int	nr,
		DFT *	tac,
		int	verbose )

Allocates memory and calculates the values for average TACs for clusters.

Returns: Returns 0 if ok, and >0 in case of an error.

Parameters

dimg	Dynamic image
cimg	Cluster image
nr	Highest cluster ID
tac	Pointer to initiated but empty DFT data
verbose	Verbose level; if zero, then only warnings are printed into stderr

Definition at line 15 of file cluster_tac.c.

  {
  if(verbose>0) printf("clusterTACs(dimg, cimg, %d, tac, %d)\n", nr, verbose);
 
  /* Check the arguments */
  if(dimg==NULL || cimg==NULL || nr<1 || tac==NULL) return(1);
  if(dimg->dimt<1 || cimg->dimt<1) return(1);
  if(cimg->dimx!=dimg->dimx || cimg->dimy!=dimg->dimy || cimg->dimz!=dimg->dimz)
    return(2);
 
  /* Allocate memory for the TACs */
  dftEmpty(tac);
  if(dftSetmem(tac, dimg->dimt, nr+1)) return(3);
  
  /* Set TAC info */
  tac->voiNr=0; tac->frameNr=dimg->dimt; tac->_type=1;
  for(int fi=0; fi<tac->frameNr; fi++) {
    tac->x1[fi]=dimg->start[fi];
    tac->x2[fi]=dimg->end[fi];
    tac->x[fi] =dimg->mid[fi];
  }
  tac->timetype=DFT_TIME_STARTEND;
  tac->timeunit=TUNIT_SEC;
  strcpy(tac->unit, imgUnit(dimg->unit));
 
  /* Calculate one cluster at a time */
  float y[dimg->dimt];
  int clusterID;
  for(clusterID=1; clusterID<=nr; clusterID++) {
    char buf[128]; snprintf(buf, 128, "%06d", clusterID);
    char *p=buf+strlen(buf)-6;
    snprintf(tac->voi[clusterID-1].voiname, MAX_REGIONSUBNAME_LEN+1, "%s", p);
    int n=imgsegmClusterMean(dimg, cimg, clusterID, y, verbose);
    if(verbose>1) printf("  clusterID%d -> %d pixels\n", clusterID, n);
    if(n<0) return(5); else if(n==0) return(6);
    for(int fi=0; fi<tac->frameNr; fi++) tac->voi[clusterID-1].y[fi]=(double)y[fi];
    tac->voi[clusterID-1].size=(double)n*dimg->sizex*dimg->sizey*dimg->sizez;
    tac->voiNr++;
  }
  /* and once more for cluster 0, i.e. the thresholded pixels;  */
  /* note that it is possible that there is no cluster 0 at all */
  clusterID=0;
  sprintf(tac->voi[tac->voiNr].voiname, "%06d", clusterID);
  strcpy(tac->voi[tac->voiNr].name, tac->voi[tac->voiNr].voiname);
  int n=imgsegmClusterMean(dimg, cimg, clusterID, y, verbose);
  if(n<0) return(7);
  if(n>0) {
    for(int fi=0; fi<tac->frameNr; fi++) tac->voi[tac->voiNr].y[fi]=(double)y[fi];
    tac->voi[tac->voiNr].size=(double)n*dimg->sizex*dimg->sizey*dimg->sizez;
    tac->voiNr++;
  }
 
  return(0);
}

◆ copy_times_from_img_to_dft()

int copy_times_from_img_to_dft	(	IMG *	img,
		DFT *	dft,
		int	verbose )

Copies frame times (especially start and end times, but also mid times) from an IMG data into DFT data, and sets DFT 'header' to indicate that frame start and end times are present.

See also: check_times_dft_vs_img, dftEndtime, imgEndtime, dftMatchTimeunits

Returns: Returns 0 if successful and <>0 in case of an error.

Parameters

img	Pointer to IMG data; times must be in sec as usual.
dft	Pointer to DFT data; times can be both seconds or minutes, if unit is correctly set; sample nr may be smaller than in IMG.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 246 of file fittime.c.

  {
  int times_changed=0;
 
  if(verbose>0) printf("%s(*img, *dft)\n", __func__);
  if(img==NULL || dft==NULL) return 1;
  if(img->dimt<=0 || dft->frameNr<=0) return 0;
  
  /* If more DFT samples, then that is an error */
  if(img->dimt<dft->frameNr) return 2;
 
  /* Convert DFT times to sec if necessary */
  if(dft->timeunit==TUNIT_MIN) {dftMin2sec(dft); times_changed=1;}
 
  /* Copy the frame times frame-by-frame */
  for(int fi=0; fi<dft->frameNr; fi++) {
    dft->x1[fi]=img->start[fi];
    dft->x2[fi]=img->end[fi];
    dft->x[fi]=img->mid[fi];
  }
  dft->timetype=DFT_TIME_STARTEND;
 
  /* Convert DFT times back to min if necessary */
  if(times_changed!=0) dftSec2min(dft);
 
  return 0; 
} 

Referenced by dftInterpolateForIMG(), img_logan(), and img_patlak().

◆ cunit_check_dft_vs_img()

int cunit_check_dft_vs_img	(	DFT *	dft,
		IMG *	img,
		char *	errmsg,
		int	verbose )

Check that calibration units in IMG (PET image) and DFT (input TAC) are the same, and if not, then try to convert DFT calibration unit to IMG unit. If input unit is unknown, then assume it is the same as the PET unit.

Returns: Returns 0 if successful, >0 in case of error, and <0 in case of a warning or error message to user is suggested.

Parameters

dft	Pointer to DFT.
img	Pointer to IMG.
errmsg	Char pointer to string (at least of length 128) where possible error description or warning is copied; set to NULL if not necessary.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 18 of file units_check.c.

  {
  int iunit, punit;
 
  if(verbose>0) printf("calibration_unit_check_dft_vs_img()\n");
  if(errmsg!=NULL) sprintf(errmsg, "program error");
  if(dft==NULL || img==NULL) return 1;
 
  iunit=petCunitId(dft->unit); // identify input file unit
  punit=img->unit;
 
  if(iunit==CUNIT_UNKNOWN) { // Input file unit is unknown
    // If PET unit is not known either, give a warning
    if(punit==CUNIT_UNKNOWN) {
      if(errmsg!=NULL) sprintf(errmsg, "unknown concentration units");
      return -1;
    } else { // Set to PET unit, and give a warning
      if(errmsg!=NULL)
        sprintf(errmsg, "unknown input concentration unit, now set to PET unit");
      strcpy(dft->unit, imgUnit(img->unit));
      return -2;
    }
  }
 
  // Input unit is known; if PET unit is not, then give a warning
  if(punit==CUNIT_UNKNOWN) {
    if(errmsg!=NULL) sprintf(errmsg, "unknown concentration units in PET data");
    return -3;
  }
 
  // Both units are known, so convert input data if necessary/possible
  if(iunit==CUNIT_KBQ_PER_ML) { // input is in units kBq/ml
    // If PET unit is the same, then everything is fine
    if(punit==CUNIT_KBQ_PER_ML) {
      if(errmsg!=NULL)
        sprintf(errmsg, "input and PET data have the same concentration units.\n");
      return 0;
    } else if(punit==CUNIT_BQ_PER_ML) { // image is in Bq/ml, convert input
      dftUnitConversion(dft, CUNIT_BQ_PER_ML);
      if(errmsg!=NULL)
        sprintf(errmsg, "input units converted to %s\n", dft->unit);
      return 0;
    } else { // image is in some other units, just give a warning for now
      if(errmsg!=NULL)
        sprintf(errmsg, "different concentration units in input and PET data");
      return -4;
    }
  } else if(iunit==CUNIT_BQ_PER_ML) { // input is in units Bq/ml
    // If PET unit is the same, then everything is fine
    if(punit==CUNIT_BQ_PER_ML) {
      if(errmsg!=NULL)
        sprintf(errmsg, "input and PET data have the same concentration units.\n");
      return 0;
    } else if(punit==CUNIT_KBQ_PER_ML) { // image is in kBq/ml, convert input
      dftUnitConversion(dft, CUNIT_KBQ_PER_ML);
      if(errmsg!=NULL)
        sprintf(errmsg, "input units converted to %s\n", dft->unit);
      return 0;
    } else { // image is in some other units, just give a warning for now
      if(errmsg!=NULL)
        sprintf(errmsg, "different concentration units in input and PET data");
      return -4;
    }
  } else { // input unit is known, but not kBq/ml or Bq/ml
    if(errmsg!=NULL)
      sprintf(errmsg, "check the concentration units in input and PET data");
    return -5;
  }
 
  return 0;
}

Referenced by imgReadModelingData().

◆ dft_end_line()

int dft_end_line	(	DFT *	dft,
		double *	fittime,
		int *	min_nr,
		int	max_nr,
		double	mintime,
		int	check_impr,
		FIT *	fit,
		FILE *	loginfo,
		char *	status )

Fits line to the end-part of TACs. By default, the included data points are determined based on maximum adjusted R^2 from at least three points; regression to the larger number of points is used in case difference in adjusted R^2 values is not larger than 0.0001.

Returns: Returns 0 when successful, otherwise <>0.

Parameters

dft	Pointer to original TAC data
fittime	By default, the search for the best line fit is started from the last sample towards the first sample; set fit time to -1 to use this default. However, if the end phase is unreliable or very noisy, you may want to set fittime to include only certain time range from the beginning. Function will write here the fittime that was actually used.
min_nr	The minimum number of samples used in searching the best fit; at least 2, but 3 is recommended. If data is very noisy, then this number may need to be increased. Function will write here the nr of samples that was actually used. This can be used as an alternative to mintime or in addition to it.
max_nr	The maximum number of samples used in searching the best fit; must be higher than min_nr, or set to -1 to not to limit the number.
mintime	Minimum time range used in searching the best fit. If data is very noisy, then this may need to be set, otherwise setting mintime to -1 will use the default. This can be used as an alternative to min_nr or in addition to it.
check_impr	Linear fitting can be applied to all data subsets in the fit time range (check_impr=0), or fitting is stopped when increasing n does not improve the adjusted R^2 (check_impr=1); the latter mode is for compatibility with WinNonlin.
fit	Pointer to data for fitted parameters. Struct must be initiated. Any existing data is deleted. Additionally, adjusted R^2 is written as 3rd (non-documented) parameter.
loginfo	Give file pointer (for example stdout) where log information is printed; NULL if not needed.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 514 of file extrapolate.c.

  {
  int fi, ret, n, to, from, to_min, from_min, orig_max_nr;
  int fitNr, snr;
  double *cx, *cy;
  double adj_r2, slope, slope_sd, ic, ic_sd, r, f, ySD_min;
  double adj_r2_max, adj_r2_prev=-10.0, ic_min, slope_min;
 
 
  /* Check the input */
  if(status!=NULL) sprintf(status, "program error");
  if(dft==NULL || fit==NULL) return -1;
  if(*min_nr<2) *min_nr=3;  // Setting erroneous value to a recommended value
  if(max_nr>-1 && max_nr<*min_nr) return -2;
  orig_max_nr=max_nr;
 
  /* Set the last sample to be fitted */
  fitNr=dft->frameNr;
  if(*fittime>0.0) {
    while(dft->x[fitNr-1]>*fittime && fitNr>0) fitNr--;
    if(loginfo!=NULL) fprintf(loginfo, "fitNr := %d\nTime range := %g - %g\n", 
      fitNr, dft->x[0], dft->x[fitNr-1]);
  }
  *fittime=dft->x[fitNr-1];
  /* Check that there are at least 3 samples */
  if(fitNr<3) { /* min_nr is checked later */
    if(status!=NULL) sprintf(status, "too few samples for linear fit");
    return(2);
  }
  
  /* If mintime was set, then get min_nr based on that */
  if(mintime>0.0) {
    fi=0; while(dft->x[fi]<(*fittime-mintime) && fi<fitNr) fi++;
    if(--fi<=0) {
      if(status!=NULL) sprintf(status, "required minimum fit range too large");
      return(2);
    }
    *min_nr=((fitNr-fi)>*min_nr?fitNr-fi:*min_nr);
  }
  if(loginfo!=NULL) fprintf(loginfo, "final_min_nr := %d\n", *min_nr);
 
  /* Allocate data for fits */
  if((ret=fit_allocate_with_dft(fit, dft))!=0) {
    if(loginfo!=NULL) fprintf(loginfo, "Error %d: cannot allocate memory for fits.\n", ret);
    if(status!=NULL) sprintf(status, "cannot allocate memory for fits");
    return 4;
  }
  /* and for x,y data to be fitted */  
  cx=(double*)malloc(2*fitNr*sizeof(double)); if(cx==NULL) {
    if(status!=NULL) sprintf(status, "out of memory");
    return(6);
  }
  cy=cx+fitNr;
 
  /* Fit each TAC */
  if(loginfo!=NULL) fprintf(loginfo, "linear fitting\n");
  for(int ri=0; ri<dft->voiNr; ri++) {
    max_nr=orig_max_nr;
    /* Print TAC name, if more than one was found */
    if(dft->voiNr>1 && loginfo!=NULL)
      fprintf(loginfo, "%s :\n", dft->voi[ri].name);
    /* Set header */
    fit->voi[ri].parNr=2;
    fit->voi[ri].type=101;
    /* Copy appropriate TAC data */
    for(fi=n=0; fi<fitNr; fi++)
      if(dft->x[fi]>0.0 && !isnan(dft->voi[ri].y[fi])) {
        cx[n]=dft->x[fi]; cy[n++]=dft->voi[ri].y[fi];}
    if(n<*min_nr) {
      if(status!=NULL) sprintf(status, "check the datafile (%d<%d)", n, *min_nr);
      free(cx); return(7);
    }
    if(max_nr<=0) max_nr=n;
    /* Search the plot range that gives the max adjusted R^2 */
    from_min=to_min=-1; adj_r2_max=-9.99E+99; ic_min=slope_min=0.0; ySD_min=0.0;
    for(from=n-*min_nr, to=n-1; from>=n-max_nr; from--) {
      snr=(to-from)+1;
      /* Calculation of linear regression using pearson() */
      ret=pearson(
        cx+from, cy+from, snr, &slope, &slope_sd, &ic, &ic_sd, &r, &f
      );
      if(ret==0) {
        adj_r2= 1.0 - ((1.0-r*r)*(double)(snr-1)) / (double)(snr-2); 
        if(adj_r2<0.0 && loginfo!=NULL)
          fprintf(loginfo, "  r=%g; snr=%d\n", r, snr);
      } else {
        adj_r2=-9.99E+99;
      }
      if(adj_r2>adj_r2_max-0.0001) {
        adj_r2_max=adj_r2; from_min=from; to_min=to;
        ic_min=ic; slope_min=slope; ySD_min=f;
      }
      if(loginfo!=NULL)
        fprintf(loginfo, "  adj_r2=%g from=%d (%g)\n",
          adj_r2, from, *(cx+from) );
      /* check for improvement, if required */
      if(check_impr!=0 && adj_r2_prev>-1.0 && adj_r2>0.0 && adj_r2<adj_r2_prev)
        break;      
 
      adj_r2_prev=adj_r2;
    }
    if(from_min<0) {
      if(status!=NULL) sprintf(status, "check the datafile");
      free(cx); return(7);
    }
    if(loginfo!=NULL) fprintf(loginfo, "  adj_r2_max=%g.\n", adj_r2_max );
    /* Set fit line parameters */
    fit->voi[ri].p[0]=ic_min;
    fit->voi[ri].p[1]=slope_min;
    fit->voi[ri].p[2]=adj_r2_max;
    fit->voi[ri].wss=ySD_min;
    fit->voi[ri].start=cx[from_min]; fit->voi[ri].end=cx[to_min];
    fit->voi[ri].dataNr=(to_min-from_min)+1;
  } // next curve
 
 
  free(cx);
  if(status!=NULL) sprintf(status, "ok");
  return 0;
}

◆ dft_ln()

int dft_ln	(	DFT *	dft1,
		DFT *	dft2 )

Natural logarithm (ln) transformation for TAC concentrations.

Returns: Returns 0 when successful, otherwise <>0.

Parameters

dft1	Pointer to original TAC data.
dft2	Pointer to allocated memory for ln transformed TAC data; enter NULL, if original data is to be overwritten by ln transformed values.

Definition at line 677 of file extrapolate.c.

  {
  DFT *out;
 
  if(dft1==NULL || dft1->voiNr<1 || dft1->frameNr<1) return 1;
  if(dft2!=NULL) out=dft2; else out=dft1;
 
  int ok_nr=0;
  for(int ri=0; ri<dft1->voiNr; ri++) for(int fi=0; fi<dft1->frameNr; fi++) {
    if(!isnan(dft1->voi[ri].y[fi]) && dft1->voi[ri].y[fi]>0.0) {
      out->voi[ri].y[fi]=log(dft1->voi[ri].y[fi]); ok_nr++;
    } else {
      out->voi[ri].y[fi]=nan("");
    }
  }
  if(ok_nr>0) return 0; else return 2;
}

◆ dftAllocateWithIMG()

int dftAllocateWithIMG	(	DFT *	dft,
		int	tacNr,
		IMG *	img )

Allocate memory for DFT based on information in IMG.

See also: sifAllocateWithIMG

Returns: Returns 0 if successful, otherwise <>0.

Parameters

dft	Pointer to initiated DFT struct which will be allocated here.
tacNr	Nr of TACs to be allocated in DFT; if zero, then TACs for all IMG pixels is allocated, but it may lead to out of memory error.
img	Pointer to IMG struct from where necessary information is read.

Definition at line 296 of file misc_model.c.

  {
  int ri, fi, ret;
 
  // Check the input data
  if(dft==NULL || img==NULL) return 1;
  if(img->status!=IMG_STATUS_OCCUPIED) return 2;
  if(img->dimt<1) return 3;
 
  // Get tacNr from image dimensions if necessary
  if(tacNr<1) {
    tacNr=img->dimz*img->dimx*img->dimy;
    if(tacNr<1) return 4;
  }
 
  // Allocate memory
  ret=dftSetmem(dft, img->dimt, tacNr);
  if(ret) return(100+ret);
  dft->voiNr=tacNr;
  dft->frameNr=img->dimt;
 
  // Set header contents
  dft->timetype=DFT_TIME_STARTEND;
  for(fi=0; fi<dft->frameNr; fi++) {
    dft->x[fi]=img->mid[fi];
    dft->x1[fi]=img->start[fi];
    dft->x2[fi]=img->end[fi];
  }
  dft->isweight=0;
  strncpy(dft->unit, imgUnit(img->unit), MAX_UNITS_LEN);
  dft->timeunit=TUNIT_SEC;
  dft->_type=DFT_FORMAT_STANDARD;
  for(ri=0; ri<dft->voiNr; ri++) {
    snprintf(dft->voi[ri].voiname, 6, "%06d", ri+1);
    strcpy(dft->voi[ri].name, dft->voi[ri].voiname);
  }
  strlcpy(dft->studynr, img->studyNr, MAX_STUDYNR_LEN);
 
  return 0;
}

Referenced by imgMaskPixelTACs().

◆ dftAutointerpolate()

int dftAutointerpolate	(	DFT *	dft,
		DFT *	dft2,
		double	endtime,
		int	verbose )

Interpolates TACs to automatically determined sample times with smaller intervals in the beginning. Only data in y arrays are interpolated; data in y2 and y3 are not used.

Returns: Function returns 0 when succesful, else a value >= 1.

Parameters

dft	Data to be interpolated is read from this structure.
dft2	Interpolated data is written in this struct; must be initiated; any previous content is deleted.
endtime	The length of interpolated/extrapolated data.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 243 of file extrapolate.c.

  {
  int i, newnr, maxNr=10000;
  double t, dt;
 
 
  if(verbose>0) printf("dftAutointerpolate(dft1, dft2, %g)\n", endtime);
  /* Check the arguments */
  if(dft->frameNr<1 || dft->voiNr<1) return 1;
  if(endtime<1.0 || endtime>1.0E12) return 2;
  if(dft->timetype!=DFT_TIME_STARTEND && dft->timetype!=DFT_TIME_MIDDLE)
    return 10;
 
  /* Calculate the number of interpolated data points */
  t=0.0; dt=0.02; newnr=1;
  while(t+dt<endtime && newnr<maxNr-1) {
    t+=dt; dt*=1.05; newnr++;
  }
  dt=endtime-t; t=endtime; newnr++;
  if(verbose>1) printf("newnr := %d\n", newnr);
 
  /* Allocate memory for interpolated data */
  dftEmpty(dft2); if(dftSetmem(dft2, newnr, dft->voiNr)) return 3;
  /* Copy header info */
  (void)dftCopymainhdr(dft, dft2); dft2->voiNr=dft->voiNr;
  for(i=0; i<dft->voiNr; i++) if(dftCopyvoihdr(dft, i, dft2, i)) return 4;
 
  /* Set times */
  if(dft->timetype==DFT_TIME_STARTEND) {
 
    dft2->timetype=DFT_TIME_STARTEND;
    t=0.0; dt=0.02; i=0;
    if(verbose>1) printf("%05d: %12.5f  %10.5f\n", i, t, dt);
    dft2->x1[i]=t; dft2->x2[i]=t+dt; dft2->x[i]=0.5*(dft2->x1[i]+dft2->x2[i]);
    i++;
    while((t+2.5*dt)<endtime && newnr<maxNr-2) {
      t+=dt; dt*=1.05;
      if(verbose>1) printf("%05d: %12.5f  %10.5f\n", i, t, dt);
      dft2->x1[i]=t; dft2->x2[i]=t+dt; dft2->x[i]=0.5*(dft2->x1[i]+dft2->x2[i]);
      i++;
    }
    t+=dt; dt=endtime-t;
    if(verbose>1) printf("%05d: %12.5f  %10.5f\n", i, t, dt);
    dft2->x1[i]=t; dft2->x2[i]=t+dt; dft2->x[i]=0.5*(dft2->x1[i]+dft2->x2[i]);
    i++; dft2->frameNr=i;
 
    /* Interpolate */
    for(i=0; i<dft->voiNr; i++)
      if(interpolate4pet(dft->x, dft->voi[i].y, dft->frameNr,
         dft2->x1, dft2->x2, dft2->voi[i].y, NULL, NULL, dft2->frameNr)) {
        dftEmpty(dft2); return 5;
      }
 
  } else if(dft->timetype==DFT_TIME_MIDDLE) {
 
    dft2->timetype=DFT_TIME_MIDDLE;
    t=0.0; dt=0.02; i=0;
    if(verbose>1) printf("%05d: %12.5f  %10.5f\n", i, t, dt);
    dft2->x1[i]=t; dft2->x2[i]=t+dt; dft2->x[i]=0.5*(dft2->x1[i]+dft2->x2[i]);
    i++;
    while((t+2.5*dt)<endtime && newnr<maxNr-2) {
      t+=dt; dt*=1.05;
      if(verbose>1) printf("%05d: %12.5f  %10.5f\n", i, t, dt);
      dft2->x1[i]=t; dft2->x2[i]=t+dt; dft2->x[i]=0.5*(dft2->x1[i]+dft2->x2[i]);
      i++;
    }
    t+=dt; dt=endtime-t;
    if(verbose>1) printf("%05d: %12.5f  %10.5f\n", i, t, dt);
    dft2->x1[i]=t; dft2->x2[i]=t+2.0*dt;
    dft2->x[i]=0.5*(dft2->x1[i]+dft2->x2[i]);
    i++; dft2->frameNr=i;
 
    /* Interpolate */
    for(i=0; i<dft->voiNr; i++)
      if(interpolate4pet(dft->x, dft->voi[i].y, dft->frameNr,
         dft2->x1, dft2->x2, dft2->voi[i].y, NULL, NULL, dft2->frameNr)) {
        dftEmpty(dft2); return 5;
      }
 
  }
 
  return 0;
}

◆ dftDerivative()

int dftDerivative	(	DFT *	dft,
		DFT *	deriv,
		char *	status )

Calculate simple derivatives from regional PET TACs. This must not be used for any quantitative purpose.

Returns: Returns 0 if successul.

Parameters

dft	DFT struct containing the regional tissue TACs
deriv	Allocated DFT struct for derivatives
status	Pointer to allocated string where the status or error message is written; set to NULL if not needed

Definition at line 635 of file dftint.c.

  {
  int ri, fi;
  double fdur;
 
  /* check input */
  if(status!=NULL) sprintf(status, "invalid input for dftDerivative()");
  if(dft==NULL || dft->frameNr<1 || dft->voiNr<1) return 1;
  if(deriv==NULL || deriv->frameNr<dft->frameNr || deriv->voiNr<dft->voiNr)
    return 2;
  if(dft->timetype!=DFT_TIME_MIDDLE && dft->timetype!=DFT_TIME_STARTEND) {
    if(status!=NULL)
      sprintf(status, "frame start and end times or mid times are required");
    return 3;
  }
 
  /* calculate frame mid times if necessary */
  if(dft->timetype==DFT_TIME_STARTEND)
    for(fi=0; fi<dft->frameNr; fi++)
      dft->x[fi]=0.5*(dft->x1[fi]+dft->x2[fi]);
 
  /* calculate derivative */
  if(dft->x[0]<=1.0E-20) for(ri=0; ri<dft->voiNr; ri++)
    deriv->voi[ri].y[0]=0.0;
  else for(ri=0; ri<dft->voiNr; ri++)
    deriv->voi[ri].y[0]=dft->voi[ri].y[0]/dft->x[0];
  for(fi=1; fi<dft->frameNr; fi++) {
    fdur=dft->x[fi]-dft->x[fi-1];
    if(fdur<=1.0E-20) for(ri=0; ri<dft->voiNr; ri++)
      deriv->voi[ri].y[fi]=0.0;
    else for(ri=0; ri<dft->voiNr; ri++)
      deriv->voi[ri].y[fi]=(dft->voi[ri].y[fi]-dft->voi[ri].y[fi-1])/fdur;
  }  
  for(ri=0; ri<dft->voiNr; ri++) for(fi=0; fi<dft->frameNr-1; fi++) {
    deriv->voi[ri].y[fi]+=deriv->voi[ri].y[fi+1];
    deriv->voi[ri].y[fi]*=0.5;
  }
 
  return 0;
}

◆ dftDerivative_old()

int dftDerivative_old	(	DFT *	dft,
		DFT *	deriv,
		char *	status )

Calculate simple derivatives from regional PET TACs. Old version. Requires that frame start and end times are known.

Returns: Returns 0 if successul.

Parameters

dft	DFT struct containing the regional tissue TACs
deriv	Allocated DFT struct for derivatives
status	Pointer to allocated string where the status or error message is written; set to NULL if not needed

Definition at line 591 of file dftint.c.

  {
  int ri, fi;
  double fdur;
 
  /* check input */
  if(status!=NULL) sprintf(status, "invalid input for dftDerivative()");
  if(dft==NULL || dft->frameNr<1 || dft->voiNr<1) return 1;
  if(deriv==NULL || deriv->frameNr<dft->frameNr || deriv->voiNr<dft->voiNr)
    return 2;
  if(dft->timetype!=3) {
    if(status!=NULL) sprintf(status, "frame start and end times are required");
    return 3;
  }
 
  /* calculate derivative */
  for(fi=0; fi<dft->frameNr; fi++) {
    fdur=dft->x2[fi]-dft->x1[fi];
    if(fdur<=1.0E-10) { 
      for(ri=0; ri<dft->voiNr; ri++) deriv->voi[ri].y[fi]=0.0;
      continue;
    }
    for(ri=0; ri<dft->voiNr; ri++) {
      deriv->voi[ri].y[fi]=dft->voi[ri].y[fi];
      if(fi>0) deriv->voi[ri].y[fi]-=dft->voi[ri].y[fi-1];
      deriv->voi[ri].y[fi]/=fdur;
    }
  }
  return 0;
}

◆ dftDivideFrames()

int dftDivideFrames	(	DFT *	dft,
		int	voi_index,
		int	add_nr,
		DFT *	dft2,
		int	verbose )

Interpolates TACs to automatically determined sample times with smaller intervals in the beginning. Only data in y arrays are interpolated; data in y2 and y3 are not used.

Returns: Function returns 0 when succesful, else a value >= 1.

Parameters

dft	Data to be interpolated into more time frames is read from this struct.
voi_index	Region index [0..voiNr-1] that is interpolated; <0 means all VOIs.
add_nr	Nr of extra time frames that are created from each original frame; valid numers are 1-100; 1 doubles the frame number.
dft2	Interpolated data is written in this struct; must be initiated, may be allocated but do not need to be; any previous content is deleted.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 406 of file extrapolate.c.

  {
  int ret, fi, fj, i, new_frameNr=0, new_voiNr=0;
  double fdur;
 
  if(verbose>0) {
    printf("dftDivideFrames(*dft, %d, %d, *dft2)\n", voi_index, add_nr);
    fflush(stdout);
  }
  /* Check the arguments */
  if(dft==NULL || dft2==NULL) return 1;
  if(dft->frameNr<1 || dft->voiNr<1) return 2;
  if(add_nr<1 || add_nr>100) return 3;
  if(voi_index>=dft->voiNr) return 4;
 
  /* Calculate how many frames and TACs will be needed */
  if(dft->timetype==DFT_TIME_STARTEND) new_frameNr=(add_nr+1)*dft->frameNr;
  else new_frameNr=(add_nr+1)*dft->frameNr-1;
  if(voi_index<0) new_voiNr=dft->voiNr; else new_voiNr=1;
  if(verbose>1) {
    printf("new_frameNr := %d\n", new_frameNr);
    printf("new_voiNr := %d\n", new_voiNr);
    fflush(stdout);
  }
 
  /* Allocate memory for interpolated data if necessary */
  if(new_frameNr>dft2->frameNr || new_voiNr>dft2->_voidataNr) {
    if(verbose>1) {
      printf("deleting old data and allocating new\n"); fflush(stdout);}
    dftEmpty(dft2);
    if(dftSetmem(dft2, new_frameNr, dft->voiNr)) return 11;
  }
 
  /* Copy header info */
  ret=dftCopymainhdr(dft, dft2); if(ret!=0) return 12;
  dft2->voiNr=new_voiNr; dft2->frameNr=new_frameNr;
  if(voi_index>=0) ret=dftCopyvoihdr(dft, voi_index, dft2, 0);
  else {
    for(int ri=0; ri<dft->voiNr; ri++) {
      ret=dftCopyvoihdr(dft, ri, dft2, ri); if(ret!=0) break;}
  }
  if(ret!=0) return 13;
 
  /* Set new frame times and interpolate */
  ret=0;
  if(dft->timetype==DFT_TIME_STARTEND) {
    for(fi=0, fj=0; fi<dft->frameNr; fi++) {
      fdur=(dft->x2[fi]-dft->x1[fi])/(double)(add_nr+1);
      for(i=0; i<add_nr+1; i++) {
        fj=fi*(add_nr+1)+i;
        dft2->x1[fj]=dft->x1[fi]+fdur*(double)i;
        dft2->x2[fj]=dft2->x1[fj]+fdur;
        dft2->x[fj]=0.5*(dft2->x1[fj]+dft2->x2[fj]);
        if(voi_index>=0)
          dft2->voi[0].y[fj]=dft->voi[voi_index].y[fi];
        else
          for(int ri=0; ri<dft->voiNr; ri++) dft2->voi[ri].y[fj]=dft->voi[ri].y[fi];
      }
    }
    dft2->frameNr=fj+1;
  } else {
    dft2->x[0]=dft->x[0];
    for(fi=1, fj=0; fi<dft->frameNr; fi++) {
      fdur=(dft->x[fi]-dft->x[fi-1])/(double)(add_nr+1);
      for(i=0; i<add_nr+1; i++) {
        fj=(fi-1)*(add_nr+1) + i + 1;
        dft2->x[fj]=dft->x[fi-1]+fdur*(double)(i+1);
      }
    }
    dft2->frameNr=fj+1;
    if(voi_index>=0) {
      ret=interpolate(dft->x, dft->voi[voi_index].y, dft->frameNr,
        dft2->x, dft2->voi[0].y, NULL, NULL, dft2->frameNr);
    } else {
      for(int ri=0; ri<dft->voiNr; ri++) {
        ret=interpolate(dft->x, dft->voi[ri].y, dft->frameNr,
                        dft2->x, dft2->voi[ri].y, NULL, NULL, dft2->frameNr);
        if(ret) break;
      }
    }
  }
  if(ret!=0) return 21;
 
  return 0;
}

◆ dftDoubleFrames()

int dftDoubleFrames	(	DFT *	dft,
		DFT *	dft2 )

Doubles the TAC sample number by making each sample/frame into two by linear interpolation. Only data in y arrays are interpolated; data in y2 and y3 are not used.

Returns: Function returns 0 when succesful, else a value >= 1.

Parameters

dft	Data to be interpolated is read from this structure.
dft2	Interpolated data is written in this struct; must be initiated; any previous content is deleted.

Definition at line 343 of file extrapolate.c.

  {
  int ri, fi, fj, ret;
  double f;
 
 
  /* Check the arguments */
  if(dft==NULL || dft2==NULL) return 1;
  if(dft->frameNr<1 || dft->voiNr<1) return 2;
  if(dft->timetype!=DFT_TIME_STARTEND && dft->x[0]<0.0) return 3;
 
  /* Allocate memory for interpolated data */
  dftEmpty(dft2); if(dftSetmem(dft2, 2*dft->frameNr, dft->voiNr)) return 11;
  /* Copy header info */
  (void)dftCopymainhdr(dft, dft2);
  dft2->voiNr=dft->voiNr; dft2->frameNr=2*dft->frameNr;
  for(ri=0; ri<dft->voiNr; ri++) if(dftCopyvoihdr(dft, ri, dft2, ri)) return 12;
 
  /* Set new time frames and interpolate */
  if(dft->timetype==DFT_TIME_STARTEND) {
    for(fi=fj=0; fi<dft->frameNr; fi++, fj+=2) {
      f=0.5*(dft->x1[fi]+dft->x2[fi]);
      dft2->x1[fj]=dft->x1[fi]; dft2->x2[fj]=f;
      dft2->x[fj]=0.5*(dft2->x1[fj]+dft2->x2[fj]);
      dft2->x1[fj+1]=f; dft2->x2[fj+1]=dft->x2[fi];
      dft2->x[fj+1]=0.5*(dft2->x1[fj+1]+dft2->x2[fj+1]);
      for(ri=0; ri<dft->voiNr; ri++)
        dft2->voi[ri].y[fj]=dft2->voi[ri].y[fj+1]=dft->voi[ri].y[fi];
    }
    ret=0;
  } else {
    for(fi=fj=0; fi<dft->frameNr; fi++) {
      if(dft->x[fi]<=0.0) {
        /* just copy, no doubles */
        dft2->x[fj++]=dft->x[fi]; continue;
      }
      if(fi==0) f=0.5*dft->x[fi];
      else f=0.5*(dft->x[fi-1]+dft->x[fi]);
      dft2->x[fj++]=f; dft2->x[fj++]=dft->x[fi];
    }
    dft2->frameNr=fj;
    for(ri=0, ret=0; ri<dft->voiNr; ri++) {
      ret=interpolate(dft->x, dft->voi[ri].y, dft->frameNr,
        dft2->x, dft2->voi[ri].y, NULL, NULL, dft2->frameNr);
      if(ret) break;
    }
  }
  if(ret) return 20+ret;
  return 0;
}

◆ dftEndtime()

double dftEndtime ( DFT * dft )

Get TAC end time. Sample times are assumed to be sorted to increasing order.

See also: imgEndtime, fittime_from_dft, getActualSamplenr

Returns: Returns the TAC end time, not converting the time units.

Parameters

dft	Pointer to DFT TAC structure.

Definition at line 315 of file fittime.c.

  {
  if(dft==NULL || dft->frameNr<1) return(0.0);
  for(int fi=dft->frameNr-1; fi>=0; fi--) {
    if(dft->timetype==DFT_TIME_MIDDLE) {
      if(!isnan(dft->x[fi])) return(dft->x[fi]);
    } else if(dft->timetype==DFT_TIME_STARTEND) {
      if(!isnan(dft->x2[fi])) return(dft->x2[fi]);
    } else if(dft->timetype==DFT_TIME_START) {
      if(!isnan(dft->x[fi])) return(dft->x[fi]);
    } else if(dft->timetype==DFT_TIME_END) {
      if(!isnan(dft->x[fi])) return(dft->x[fi]);
    } else
      return(0.0);
  }
  return(0.0);
}

Referenced by imgReadModelingData().

◆ dftInterpolate()

int dftInterpolate	(	DFT *	input,
		DFT *	tissue,
		DFT *	output,
		char *	status,
		int	verbose )

Interpolate (and integrate) TAC data to sample times that are given with another TAC data. PET frame lengths are taken into account if available in tissue DFT struct.

Input frame lengths are taken into account if the framing is same as with tissue data.

See also: dftInterpolateInto, dftTimeIntegral

Returns: Returns 0 if successful, and <>0 in case of an error.

Parameters

input	Data which is interpolated; make sure that time unit is the same as in tissue data; time range is checked to cover the tissue data
tissue	Data to which (sample times) the interpolation is done
output	Pointer to initiated DFT into which interpolated values and integrals will be written
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 162 of file dftint.c.

  {
  int fi, ri, ret;
 
  if(verbose>0) printf("dftInterpolate()");
 
  // Check the input
  if(input==NULL || tissue==NULL || output==NULL) {
    if(status!=NULL) sprintf(status, "program error");
    return 1;
  }
 
  // If input and tissue data have the same frame times already, then
  // copy frame times and timetype into input
  if(tissue->timetype==DFT_TIME_STARTEND && 
     check_times_dft_vs_dft(tissue, input, verbose)==1 &&
     input->frameNr<=tissue->frameNr)
  {
    for(fi=0; fi<input->frameNr; fi++) {
      input->x[fi]=tissue->x[fi];
      input->x1[fi]=tissue->x1[fi]; input->x2[fi]=tissue->x2[fi];
    }
    input->timetype=tissue->timetype;
  }
 
  // Check that there is no need for excess extrapolation
  ret=dftInterpolateCheckEnd(input, tissue, status, verbose);
  if(ret<0) return 2;
  ret=dftInterpolateCheckStart(input, tissue, status, verbose);
  if(ret<0) return 2;
 
  // Delete any previous output data
  dftEmpty(output);
 
  // Allocate memory for interpolated data
  if(dftSetmem(output, tissue->frameNr, input->voiNr)) {
    if(status!=NULL) strcpy(status, "memory allocation error");
    return 3;
  }
  output->voiNr=input->voiNr; output->frameNr=tissue->frameNr;
 
  // Copy header information
  dftCopymainhdr(input, output);
  for(ri=0; ri<input->voiNr; ri++) dftCopyvoihdr(input, ri, output, ri);
 
  // Copy frame information
  output->isweight=tissue->isweight;
  for(fi=0; fi<tissue->frameNr; fi++) {
    output->x[fi]=tissue->x[fi];
    output->x1[fi]=tissue->x1[fi]; output->x2[fi]=tissue->x2[fi];
    output->w[fi]=tissue->w[fi];
  }
  output->timetype=tissue->timetype;
 
  // Check if input and tissue data do have the same frame times already
  if(check_times_dft_vs_dft(tissue, input, verbose)==1 &&
     input->frameNr>=tissue->frameNr)
  {
    // copy the values directly and integrate in place
    for(ri=0, ret=0; ri<output->voiNr && ret==0; ri++) {
      for(fi=0; fi<tissue->frameNr; fi++)
        output->voi[ri].y[fi]=input->voi[ri].y[fi];
      if(output->timetype==3)
        ret=petintegral(output->x1, output->x2, output->voi[ri].y,
               output->frameNr, output->voi[ri].y2, output->voi[ri].y3);
      else
        ret=interpolate(output->x, output->voi[ri].y, output->frameNr,
               output->x, output->voi[ri].y, output->voi[ri].y2,
               output->voi[ri].y3, output->frameNr);
    }
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot interpolate (%d)", ret);
      dftEmpty(output); return 5;
    }
    return 0; // that's it then
  }
 
  // Interpolate and integrate input data to tissue sample times,
  // taking into account tissue frame lengths if available
  for(ri=0, ret=0; ri<output->voiNr && ret==0; ri++) {
    if(output->timetype==DFT_TIME_STARTEND)
      ret=interpolate4pet(input->x, input->voi[ri].y, input->frameNr,
                output->x1, output->x2, output->voi[ri].y, output->voi[ri].y2,
                output->voi[ri].y3, output->frameNr);
    else
      ret=interpolate(input->x, input->voi[ri].y, input->frameNr,
                output->x, output->voi[ri].y, output->voi[ri].y2,
                output->voi[ri].y3, output->frameNr);
  }
  if(ret) {
    if(status!=NULL) sprintf(status, "cannot interpolate (%d)", ret);
    dftEmpty(output); return 6;
  }
 
  return 0;
}

Referenced by dftReadinput().

◆ dftInterpolateCheckEnd()

int dftInterpolateCheckEnd	(	DFT *	input,
		DFT *	output,
		char *	status,
		int	verbose )

Verify that data to-be-interpolated does not need too much extrapolation in the end, and that samples are not too sparse.

Time units in DFTs can be different, if specified.

See also: dftInterpolateCheckStart, dftInterpolateInto, dftInterpolate

Returns: Returns 0 if data is fine, 1 if extrapolation can be done, but there may be too few samples, and -1 if extrapolation in the end is impossible.

Parameters

input	Data that will be verified for reliable interpolation
output	Data containing sample times for interpolated data
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 82 of file dftint.c.

  {
  int fi, n, itunit;
  double t1, t2;
 
  if(verbose>0) printf("dftInterpolateCheckEnd()\n");
 
  // Check the input
  if(status!=NULL) sprintf(status, "program error");
  if(input==NULL || output==NULL) return -1;
  if(input->frameNr<1 || output->frameNr<1) return -1;
  if(input->frameNr==1 && output->frameNr>1) {
    if(status!=NULL) sprintf(status, "too short data for interpolation");
    return -1;
  }
  if(status!=NULL) sprintf(status, "ok");
 
  /* Try to make sure that time units are the same */
  dftMatchTimeunits(output, input, &itunit, verbose);
 
  /* Check that input data extends almost to the end of output range */
  if(input->timetype==DFT_TIME_STARTEND && output->timetype==DFT_TIME_STARTEND) {
    t1=input->x2[input->frameNr-1]; t2=output->x2[output->frameNr-1]; 
  } else {
    t1=input->x[input->frameNr-1];  t2=output->x[output->frameNr-1];
  }
  if(t1<0.95*t2) {
    if(status!=NULL) strcpy(status, "too short data for interpolation");
    if(verbose>1) printf("t1 := %g\nt2 := %g\n", t1, t2);
    dftTimeunitConversion(input, itunit); // back to original time units
    return -1;
  }
 
  /* Check that input sample frequency is not too low in the end */
  if(output->frameNr>3) {
    if(output->timetype==DFT_TIME_STARTEND) {
      t1=output->x1[output->frameNr-3]; 
      t2=output->x2[output->frameNr-1];
    } else { 
      t1=output->x[output->frameNr-3]; 
      t2=output->x[output->frameNr-1];
    }
    for(fi=n=0; fi<input->frameNr; fi++)
      if(input->x[fi]>=t1 && input->x[fi]<=t2) n++;
    if(n<1 || (n<2 && t2>input->x[input->frameNr-1])) {
      if(status!=NULL) strcpy(status, "too sparse sampling for interpolation");
      if(verbose>1) printf("n=%d t1=%g t2=%g\n", n, t1, t2);
      dftTimeunitConversion(input, itunit); // back to original time units
      return -1; // Error
    }
    if(n<2 || (n<3 && t2>input->x[input->frameNr-1])) {
      if(status!=NULL) strcpy(status, "too sparse sampling for interpolation");
      if(verbose>1) printf("n=%d t1=%g t2=%g\n", n, t1, t2);
      dftTimeunitConversion(input, itunit); // back to original time units
      return 1; // Warning
    }
  } 
  dftTimeunitConversion(input, itunit); // back to original time units
  return 0;
}

Referenced by dftInterpolate(), dftInterpolateForIMG(), and dftInterpolateInto().

◆ dftInterpolateCheckStart()

int dftInterpolateCheckStart	(	DFT *	input,
		DFT *	output,
		char *	status,
		int	verbose )

Verify that data to-be-interpolated does not need too much extrapolation in the beginning.

See also: dftInterpolateCheckEnd, dftInterpolate, dftInterpolateInto, dftTimeIntegral

Returns: Returns 0 if data is fine, 1 if it starts late but extrapolation can be done reliably, and -1 if extrapolation in the beginning would be too risky.

Parameters

input	Data that will be verified for reliable interpolation
output	Data containing sample times for interpolated data
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 17 of file dftint.c.

  {
  int ri, itunit;
  double t1, t2, f;
 
  if(verbose>0) printf("dftInterpolateCheckStart()\n");
 
  // Check the input
  if(status!=NULL) sprintf(status, "program error");
  if(input==NULL || output==NULL) return -1;
  if(input->frameNr<1 || output->frameNr<1) return -1;
  if(input->frameNr==1 && output->frameNr>1) {
    if(status!=NULL) sprintf(status, "too short data for interpolation");
    return -1;
  }
  if(status!=NULL) sprintf(status, "ok");
 
  /* Try to make sure that time units are the same */
  dftMatchTimeunits(output, input, &itunit, verbose);
 
  /* Check the start */
  if(input->timetype==DFT_TIME_STARTEND && output->timetype==DFT_TIME_STARTEND) {
    t1=input->x1[0]; t2=output->x1[0];
  } else {
    t1=input->x[0];  t2=output->x[0];
  }
  if(0.95*t1>t2) {
    if(verbose>1) printf("t1 := %g\nt2 := %g\n", t1, t2);
    /* Check that first value is relatively low, so that there will not be any
       error when the initial part is assumed to be a straight line from
       (0,0) to the first sample point */
    f=0.25*dft_kBqMax(input);
    for(ri=0; ri<input->voiNr; ri++) if(input->voi[ri].y[0] > f) {
      if(status!=NULL) strcpy(status, "data starts too late");
      dftTimeunitConversion(input, itunit); // back to original time units
      return -1;
    }
    if(status!=NULL) strcpy(status, "data starts late");
    dftTimeunitConversion(input, itunit); // back to original time units
    return 1;
  }
  dftTimeunitConversion(input, itunit); // back to original time units
  return 0;
}

Referenced by dftInterpolate(), dftInterpolateForIMG(), and dftInterpolateInto().

◆ dftInterpolateForIMG()

int dftInterpolateForIMG	(	DFT *	input,
		IMG *	img,
		int	frame_nr,
		DFT *	output,
		double *	ti1,
		double *	ti2,
		int	verbose,
		char *	status )

Interpolate (and integrate) TAC data to sample times that are given with IMG data. Input frame lengths are taken into account if they are available in DFT or if the framing is the same as in IMG data.

Returns: Returns 0 if successful, and <>0 in case of an error.

Parameters

input	Data which is interpolated.
img	Data to which (sample times) the interpolation is done.
frame_nr	Number of IMG frames that are needed; can be set to 0 if all frames can be included.
output	Pointer to initiated DFT into which interpolated values and integrals will be written at input sample times and units.
ti1	First time of input data before interpolation (in input time units); use this to check that required time range was measured; NULL if not needed.
ti2	Last time of input data before interpolation (in input time units); use this to check that required time range was measured; NULL if not needed.
verbose	Verbose level; set to zero to not to print any comments.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 17 of file misc_model.c.

  {
  /* Check input */
  if(verbose>0) {
    printf("%s(*inp, *img, %d, *out, *ti1, *ti2, %d, status)\n", __func__, frame_nr, verbose);
    fflush(stdout);
  }
  if(input==NULL || img==NULL || output==NULL) {
    if(status!=NULL) strcpy(status, "program error");
    return 1;
  }
  if(input->frameNr<1 || input->voiNr<1 || img->dimt<1) {
    if(status!=NULL) strcpy(status, "no pet data");
    return 2;
  }
  if(input->timeunit!=TUNIT_MIN && input->timeunit!=TUNIT_SEC) {
    if(status!=NULL) strcpy(status, "unknown time units");
    return 3;
  }
  if(frame_nr<1 || frame_nr>img->dimt) frame_nr=img->dimt;
  /* Set ti1 and ti2 */
  if(ti1!=NULL) {
    if(input->timetype==DFT_TIME_STARTEND)
      *ti1=input->x1[0]; else *ti1=input->x[0];
  }
  if(ti2!=NULL) {
    if(input->timetype==DFT_TIME_STARTEND) *ti2=input->x2[input->frameNr-1];
    else *ti2=input->x[input->frameNr-1];
  }
 
  /* Delete any previous data */
  dftEmpty(output);
 
  /* Allocate memory for interpolated data */
  if(verbose>10) printf("allocating memory for interpolated data\n");
  if(dftAllocateWithHeader(output, frame_nr, input->voiNr, input)) {
    if(status!=NULL) strcpy(status, "memory allocation error");
    return 11;
  }
  output->voiNr=input->voiNr; output->frameNr=frame_nr;
 
  /* Set output times */
  if(copy_times_from_img_to_dft(img, output, verbose)!=0) {
    if(status!=NULL) strcpy(status, "frame time error");
    dftEmpty(output); return 12;
  }
  if(verbose>10) {
    printf("time range := %g - %g %s\n", output->x[0],
           output->x[output->frameNr-1], petTunit(output->timeunit));
    printf("  timetype := %d\n", output->timetype);
  }
 
  // Check if input and output data do have the same frame times already
  if(check_times_dft_vs_dft(input, output, verbose)==1 && input->frameNr>=output->frameNr) {
    // copy the values directly and integrate in place
    if(verbose>10) printf("frame times are assumed to be the same\n");
    int ret=0;
    for(int ri=0; ri<output->voiNr && ret==0; ri++) {
      for(int fi=0; fi<output->frameNr; fi++)
        output->voi[ri].y[fi]=input->voi[ri].y[fi];
      ret=petintegral(output->x1, output->x2, output->voi[ri].y,
             output->frameNr, output->voi[ri].y2, output->voi[ri].y3);
    }
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot interpolate (%d)", ret);
      dftEmpty(output); return 15;
    }
    return 0; // that's it then
  }
  if(verbose>10) printf("frame times are not the same\n");
 
  // Check that there is no need for extrapolation in the start
  if(dftInterpolateCheckStart(input, output, status, verbose)<0) {
    dftEmpty(output); return 16;
  }
  // Check that there is no need for excess extrapolation in the end either
  if(dftInterpolateCheckEnd(input, output, status, verbose)<0) { 
    dftEmpty(output); return 17;
  }
 
  // Interpolate and integrate input data to tissue sample times,
  // taking into account tissue frame lengths
  {
    int ret=0;
    for(int ri=0; ri<output->voiNr && ret==0; ri++) {
      ret=interpolate4pet(input->x, input->voi[ri].y, input->frameNr,
                  output->x1, output->x2, output->voi[ri].y, output->voi[ri].y2,
                  output->voi[ri].y3, output->frameNr);
    }
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot interpolate (%d)", ret);
      dftEmpty(output); return 18;
    }
  }
 
  return 0;
}

Referenced by imgReadModelingData().

◆ dftInterpolateInto()

int dftInterpolateInto	(	DFT *	inp,
		DFT *	tis,
		char *	status,
		int	verbose )

Interpolate (and integrate) TAC data to sample times that are given with another TAC data. New TAC is written into existing TAC data. PET frame lengths are taken into account if available in tissue DFT struct. Input frame lengths are taken into account if the framing is same as with tissue data.

See also: dftTimeIntegral, dftInterpolate, dftInterpolateCheckStart, dftInterpolateCheckEnd

Returns: Returns 0 if successful, and <>0 in case of an error.

Parameters

inp	Data which is interpolated; make sure that time unit is the same as in tissue data; time range is checked to cover the tissue data
tis	Data into which the interpolation is done
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 281 of file dftint.c.

  {
  int fi, ri, ret;
 
  if(verbose>0) printf("dftInterpolateInto()\n");
 
  // Check the input
  if(inp==NULL || tis==NULL) {
    if(status!=NULL) sprintf(status, "program error");
    return 1;
  }
  ret=dft_nr_of_NA(inp); if(ret>0) {
    if(status!=NULL) sprintf(status, "missing sample(s)");
    return 2;
  }
 
  // Check that there is no need for excess extrapolation
  ret=dftInterpolateCheckEnd(inp, tis, status, verbose);
  if(ret<0) return 3;
  ret=dftInterpolateCheckStart(inp, tis, status, verbose);
  if(ret<0) return 4;
 
  // Allocate memory for interpolated data
  if(dftAddmem(tis, inp->voiNr)) {
    if(status!=NULL) strcpy(status, "memory allocation error");
    return 5;
  }
 
  // Copy header information
  for(ri=0; ri<inp->voiNr; ri++)
    dftCopyvoihdr(inp, ri, tis, tis->voiNr+ri);
 
  // Check if input and tissue data do have the same frame times already
  if(check_times_dft_vs_dft(inp, tis, verbose)==1 &&
     inp->frameNr>=tis->frameNr)
  {
    // copy the values directly and integrate in place
    for(ri=0, ret=0; ri<inp->voiNr && ret==0; ri++) {
      for(fi=0; fi<tis->frameNr; fi++)
        tis->voi[tis->voiNr+ri].y[fi]=inp->voi[ri].y[fi];
      if(tis->timetype==DFT_TIME_STARTEND)
        ret=petintegral(tis->x1, tis->x2, tis->voi[tis->voiNr+ri].y,
               tis->frameNr, tis->voi[tis->voiNr+ri].y2, 
               tis->voi[tis->voiNr+ri].y3);
      else
        ret=interpolate(tis->x, tis->voi[tis->voiNr+ri].y, tis->frameNr,
               tis->x, tis->voi[tis->voiNr+ri].y, 
               tis->voi[tis->voiNr+ri].y2,
               tis->voi[tis->voiNr+ri].y3, tis->frameNr);
    }
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot interpolate (%d)", ret);
      dftEmpty(tis); return 7;
    }
    tis->voiNr+=inp->voiNr;
    return 0; // that's it then
  }
 
  // Interpolate and integrate input data to tissue sample times,
  // taking into account tissue frame lengths if available
  for(ri=0, ret=0; ri<inp->voiNr && ret==0; ri++) {
    if(tis->timetype==DFT_TIME_STARTEND)
      ret=interpolate4pet(inp->x, inp->voi[ri].y, inp->frameNr,
                tis->x1, tis->x2, tis->voi[tis->voiNr+ri].y, 
                tis->voi[tis->voiNr+ri].y2,
                tis->voi[tis->voiNr+ri].y3, tis->frameNr);
    else
      ret=interpolate(inp->x, inp->voi[ri].y, inp->frameNr,
                tis->x, tis->voi[tis->voiNr+ri].y, 
                tis->voi[tis->voiNr+ri].y2,
                tis->voi[tis->voiNr+ri].y3, tis->frameNr);
  }
  if(ret) {
    if(status!=NULL) sprintf(status, "cannot interpolate (%d)", ret);
    return 9;
  }
  tis->voiNr+=inp->voiNr;
 
  return 0;
}

Referenced by dftReadModelingData(), dftReadReference(), and imgReadModelingData().

◆ dftMatchTimeunits()

int dftMatchTimeunits	(	DFT *	dft1,
		DFT *	dft2,
		int *	tunit2,
		int	verbose )

Make sure that time units in two DFT structs are the same, converting units when necessary, optionally saving original units so that units can be converted back to what they were using dftTimeunitConversion().

See also: dftTimeunitConversion, dftEndtime, fittime_from_dft

Returns: Returns 0 if the was no need for time unit conversion, 1 if time units in dft2 were converted, and <0 if units were not identified or in case of an error.

Parameters

dft1	Pointer to DFT struct 1.
dft2	Pointer to DFT struct 2; sample time units are changed to match the data in dft1.
tunit2	Pointer for original time unit in DFT 2; enter NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 358 of file fittime.c.

  {
  if(verbose>0) printf("%s()\n", __func__);
  if(dft1==NULL || dft2==NULL) return(-1);
  /* Save original units if required */
  if(tunit2!=NULL) *tunit2=dft2->timeunit;
  /* Check that time units are available */
  if(dft1->timeunit!=TUNIT_MIN && dft1->timeunit!=TUNIT_SEC) {
    if(verbose>0) printf("  unknown time units in dft1\n");
    return(-2);
  }
  if(dft2->timeunit!=TUNIT_MIN && dft2->timeunit!=TUNIT_SEC) {
    if(verbose>0) printf("  unknown time units in dft2\n");
    return(-3);
  }
  /* Check if they are the same */
  if(dft1->timeunit==dft2->timeunit) {
    if(verbose>1) printf("  time units are the same in dft1 and dft2.\n");
    return(0);
  }
  /* Conversion to dft2 */
  if(verbose>1) printf("  time units in dft2 converted from %s to %s\n",
    petTunit(dft2->timeunit), petTunit(dft1->timeunit));
  int ret=dftTimeunitConversion(dft2, dft1->timeunit);
  if(ret) {
    if(verbose>0) printf("  time unit conversion failed\n");
    return(-10-ret);
  }
  return(1);
}

Referenced by dftInterpolateCheckEnd(), and dftInterpolateCheckStart().

◆ dftReadinput()

int dftReadinput	(	DFT *	input,
		DFT *	tissue,
		char *	filename,
		int *	filetype,
		double *	ti1,
		double *	ti2,
		int	verifypeak,
		char *	status,
		int	verbose )

Read DFT format input TAC data to match the time frames in the specified tissue DFT file. Instead of input filename, a reference region name can be given: then all the matching tacs (based on roi names) are moved from the roi data to the input data, with the best match as first tac. Input data is interpolated (if necessary) and also integrated to y2[]. Input time and concentration units are tried to be converted to be the same as in tissue data.

Returns: Returns 0 if successful, and >0 in case of an error, and specifically 101 in case input TAC is not valid.

See also: dftReadReference, dftRead, dftReadModelingData, dftVerifyPeak, imgReadModelingData

Parameters

input	Input data, previous contents are cleared.
tissue	PET data containing frames and possible input regions.
filename	Input data filename, or region name in PET data.
filetype	Type of input, as found out here; 1 and 2 =tac, 3=fit, 5=region name; enter NULL, if not needed.
ti1	First time of input data before interpolation (in tissue time units); use this to check that required time range was measured; NULL if not needed.
ti2	Last time of input data before interpolation (in tissue time units); use this to check that required time range was measured; NULL if not needed.
verifypeak	Set to <>0 to verify that input TAC does not start too late and has decent peak to provide reliable AUC0-T.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 25 of file dftinput.c.

  {
  int ri, n, ret, ftype=0;
  DFT temp;
  FILE *fp;
 
  if(verbose>0) printf("dftReadinput(inp, tis, %s, type, ti1, ti2, %d, status, %d)\n",
                       filename, verifypeak, verbose);
  if(filetype!=NULL) *filetype=0;
 
  /* Check that PET data is ok */
  if(input==NULL || tissue==NULL || filename==NULL) {
    if(status!=NULL) strcpy(status, "program error");
    return 1;
  }
  if(tissue->frameNr<1 || tissue->voiNr<1) {
    if(status!=NULL) strcpy(status, "no pet data");
    return 2;
  }
 
  /* Delete any previous input data and initiate temp data */
  dftEmpty(input); dftInit(&temp);
 
  /* Try to figure out what is the 'filename' */
  ftype=0;
  /* Can we open it as a file? */
  fp=fopen(filename, "r");
  if(fp!=NULL) {
    /* File can be opened for read; close it for now */
    if(verbose>1) printf("  file can be opened for reading.\n");
    if(filetype!=NULL) *filetype=1;
    fclose(fp);
 
    /* Try to identify the file format */
    ftype=dftFormat(filename);
    if(ftype==DFT_FORMAT_UNKNOWN) {
      if(filetype!=NULL) *filetype=0;
      if(status!=NULL) strcpy(status, "unknown file format");
      return 3;
    } else if(ftype==DFT_FORMAT_FIT) {
      if(filetype!=NULL) *filetype=3;
      if(status!=NULL) strcpy(status, "cannot read fit file"); 
      return 3;
    }
    if(verbose>2) printf("  fileformat=%d\n", ftype);
 
    /* Try to read it */
    if((ret=dftRead(filename, &temp))!=0) {
      if(filetype!=NULL) *filetype=0;
      if(status!=NULL) sprintf(status, "cannot read file (%d)", ret);
      return(2);
    }
    /* Convert input time units to the same as in tissue data */
    (void)dftTimeunitConversion(&temp, tissue->timeunit);
    /* Check the tissue and plasma TAC concentration units */
    ret=dftUnitConversion(&temp, petCunitId(tissue->unit));
    if(ret!=0) {sprintf(status, "check the units of input and tissue data");}
    /* Tell user what was the original input time range */
    if(temp.timetype==DFT_TIME_STARTEND) {
      if(ti1!=NULL) *ti1=temp.x1[0];
      if(ti2!=NULL) *ti2=temp.x2[temp.frameNr-1];
    } else {
      if(ti1!=NULL) *ti1=temp.x[0];
      if(ti2!=NULL) *ti2=temp.x[temp.frameNr-1];
    }
    /* Verify the peak if requested */
    if(verifypeak!=0) {
      ret=dftVerifyPeak(&temp, 0, verbose-2, status);
      //if(ret!=0) sprintf(status, "input TAC should start at time zero");
      if(ret>0) {dftEmpty(&temp); return 101;}
    }
    /* Interpolate and integrate data to pet times */
    ret=dftInterpolate(&temp, tissue, input, status, verbose);
    dftEmpty(&temp);
    if(ret!=0) return 4;
 
  } else {
 
    /* it's not a file, at least an accessible file, but is it a region name? */
    if(filetype!=NULL) *filetype=5;
    /* Select ROIs that match the specified input name */
    n=dftSelectRegions(tissue, filename, 1);
    if(n<=0) {
      if(status!=NULL) sprintf(status, "cannot find region");
      return 7;
    }
    if(n==tissue->voiNr) {
      if(status!=NULL) sprintf(status, "all regions do match");
      return 8;
    }
    /* one or more regions was found ; move them to input data */
    ret=dftdup(tissue, input); if(ret==0) {
      ri=0; ret=0; while(ri<input->voiNr && ret==0) {
        if(input->voi[ri].sw==0) ret=dftDelete(input, ri); else ri++;
      }
      if(ret==0) {
        ri=0; ret=0; while(ri<tissue->voiNr && ret==0) {
          if(tissue->voi[ri].sw!=0) ret=dftDelete(tissue, ri); else ri++;
        }
      }
    }
    if(ret!=0) {
      if(status!=NULL) sprintf(status, "cannot separate input regions");
      dftEmpty(input); return 9;
    }
    /* Try to select the best reference ROI */
    ri=dftSelectBestReference(input); if(ri<0) {
      if(status!=NULL) sprintf(status, "cannot separate input regions");
      dftEmpty(input); return 10;
    }
    /* And move it to the first place */
    if(ri>0) {dftMovevoi(input, ri, 0); ri=0;}
    if(verbose>1) printf("selected ref region := %s\n", input->voi[0].name);
    /* Verify the peak if requested */
    if(verifypeak!=0) {
      ret=dftVerifyPeak(input, 0, verbose-2, status);
      //if(ret!=0) sprintf(status, "input TAC should start at time zero");
      if(ret>0) {dftEmpty(input); return 101;}
    }
 
    /* Calculate integrals */
    for(ri=0; ri<input->voiNr; ri++) {
      if(input->timetype==DFT_TIME_STARTEND) {
        ret=petintegral(input->x1, input->x2, input->voi[ri].y, input->frameNr,
              input->voi[ri].y2, input->voi[ri].y3);
        if(ti1!=NULL) *ti1=input->x1[0];
        if(ti2!=NULL) *ti2=input->x2[input->frameNr-1];
      } else {
        ret=interpolate(input->x, input->voi[ri].y, input->frameNr,
               input->x, NULL, input->voi[ri].y2, input->voi[ri].y3, input->frameNr);
        if(ti1!=NULL) *ti1=input->x[0];
        if(ti2!=NULL) *ti2=input->x[input->frameNr-1];
      }
      if(ret) {
        if(status!=NULL) sprintf(status, "cannot integrate input");
        dftEmpty(input); return(11);
      }
    }
  }
 
  return(0);
}

◆ dftReadModelingData()

int dftReadModelingData	(	char *	tissuefile,
		char *	inputfile1,
		char *	inputfile2,
		char *	inputfile3,
		double *	fitdur,
		int *	fitframeNr,
		DFT *	tis,
		DFT *	inp,
		FILE *	loginfo,
		int	verbose,
		char *	status )

Read tissue and input data for modeling. Time units are converted to min and input calibration units to units of tissue data. Input data is NOT interpolated to tissue times. If input data extends much further than fit duration, the last part is removed to save computation time in simulations. Input data is verified not to end too early, and not to start too late.

Returns: Returns 0 when successful, otherwise a non-zero value.

See also: dftReadReference, dftRead, dftReadinput, imgReadModelingData

Parameters

tissuefile	Tissue data filename; one time sample is sufficient here.
inputfile1	1st input data filename.
inputfile2	2nd input data filename (or NULL if only not needed).
inputfile3	3rd input data filename (or NULL if only not needed).
fitdur	Fit duration (in minutes); shortened if longer than tissue data.
fitframeNr	Nr of time frames (samples) in tissue data that are inside fitdur will be written here.
tis	Pointer to initiated DFT into which tissue data will be written.
inp	Pointer to initiated DFT into which input data (plasma and/or blood) will be written.
loginfo	Give file pointer (for example stdout) where log information is printed; NULL if not needed; warnings will be printed in stderr anyway.
verbose	Verbose level; if zero, then only warnings are printed into stderr.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 342 of file dftinput.c.

  {
  int ret, fi, n, first, last, ii, input_nr=0;
  DFT tmpdft;
  double f, starttime, endtime;
  FILE *wfp;
  char *fname;
 
 
  if(loginfo!=NULL && verbose>0) {
    fprintf(loginfo, "dftReadModelingData(\n");
    fprintf(loginfo, "  %s,\n", tissuefile);
    fprintf(loginfo, "  %s,\n", inputfile1);
    fprintf(loginfo, "  %s,\n", inputfile2);
    fprintf(loginfo, "  %s,\n", inputfile3);
    fprintf(loginfo, "  %g,\n", *fitdur);
    fprintf(loginfo, "  *fitframeNr, *tis, *inp, *loginfo, %d, *status\n", verbose);
    fprintf(loginfo, ")\n");
  }
  /* Check the input */
  if(status!=NULL) sprintf(status, "program error");
  if(tis==NULL || inp==NULL) return -1;
  if(tissuefile==NULL || strlen(tissuefile)<1) return -2;
  ret=0;
  if(inputfile1==NULL || strlen(inputfile1)<1) return -3; else input_nr++;
  if(inputfile2!=NULL && strlen(inputfile2)>0) input_nr++;
  if(inputfile3!=NULL && strlen(inputfile3)>0) {
    if(input_nr<2) return -4;
    input_nr++;
  }
  if(status!=NULL) sprintf(status, "arguments validated");
  /* Warnings will be printed into stderr */
  wfp=stderr;
 
  /* Delete any previous data and initiate temp data */
  dftEmpty(inp); dftEmpty(tis); dftInit(&tmpdft);
 
  /*
   *  Read tissue data
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading tissue data in %s\n", tissuefile);
  if(dftRead(tissuefile, tis)) {
    if(status!=NULL) sprintf(status, "cannot read '%s': %s", tissuefile, dfterrmsg);
    return(2);
  }
  /* Check for NaN's */
  ret=dft_nr_of_NA(tis); if(ret>0) {
    if(status!=NULL) sprintf(status, "missing sample(s) in %s", tissuefile);
    dftEmpty(tis); return(2);
  }
  /* Sort the data by increasing sample times */
  dftSortByFrame(tis);
  /* Do not check frame number; static scan is fine here */
  /* Make sure that there is no overlap in frame times */
  if(tis->timetype==DFT_TIME_STARTEND) {
    if(loginfo!=NULL && verbose>0)
      fprintf(loginfo, "checking frame overlap in %s\n", tissuefile);
    ret=dftDeleteFrameOverlap(tis);
    if(ret) {
      if(status!=NULL) sprintf(status, "%s has overlapping frame times", tissuefile);
      dftEmpty(tis); return(2);
    }
  }
  if(tis->timeunit==TUNIT_UNKNOWN) {
    fprintf(wfp, "Warning: tissue sample time units not known.\n");
  }
 
  /*
   *  Read first input data
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading input data in %s\n", inputfile1);
  if(dftRead(inputfile1, inp)) {
    if(status!=NULL) sprintf(status, "cannot read '%s': %s", inputfile1, dfterrmsg);
    dftEmpty(tis); return(3);
  }
  /* Check and correct the sample time unit */
  if(tis->timeunit==TUNIT_UNKNOWN) tis->timeunit=inp->timeunit;
  else if(inp->timeunit==TUNIT_UNKNOWN) inp->timeunit=tis->timeunit;
  if(inp->timeunit==TUNIT_UNKNOWN) {
    fprintf(wfp, "Warning: input sample time units not known.\n");
  }
  if(tis->timeunit!=inp->timeunit &&
     dftTimeunitConversion(inp, tis->timeunit)!=0) {
    if(status!=NULL) sprintf(status, "tissue and plasma do have different time units");
    dftEmpty(tis); dftEmpty(inp);
    return(3);
  }
  if(inp->voiNr>1) {
    fprintf(wfp, "Warning: using only first TAC in %s\n", inputfile1);
    inp->voiNr=1;
  }
  if(inp->frameNr<4) {
    if(status!=NULL) sprintf(status, "%s contains too few samples", inputfile1);
    dftEmpty(tis); dftEmpty(inp); return(3);
  }
  /* Sort the data by increasing sample times */
  dftSortByFrame(inp);
 
  /*
   *  Read following input files, if required
   */
  for(ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=inputfile2;
    else fname=inputfile3;
 
    /* Make room for one more curve in input data */
    ret=dftAddmem(inp, 1);
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot allocate more memory");
      dftEmpty(tis); dftEmpty(inp); return(4);
    }
 
    /* Read blood data */
    dftInit(&tmpdft);
    if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading input data in %s\n", fname);
    if(dftRead(fname, &tmpdft)) {
      if(status!=NULL) sprintf(status, "cannot read '%s': %s", fname, dfterrmsg);
      dftEmpty(tis); dftEmpty(inp); return(4);
    }
    if(tmpdft.frameNr<4) {
      if(status!=NULL) sprintf(status, "%s contains too few samples", fname);
      dftEmpty(tis); dftEmpty(inp); dftEmpty(&tmpdft); return(4);
    }
 
    /* Check and correct the sample time unit */
    if(tis->timeunit==TUNIT_UNKNOWN) tis->timeunit=tmpdft.timeunit;
    else if(tmpdft.timeunit==TUNIT_UNKNOWN) tmpdft.timeunit=tis->timeunit;
    if(tmpdft.timeunit==TUNIT_UNKNOWN) {
      fprintf(wfp, "Warning: blood sample time units not known.\n");
    }
    if(inp->timeunit!=tmpdft.timeunit &&
       dftTimeunitConversion(&tmpdft, inp->timeunit)!=0) {
      if(status!=NULL) sprintf(status, "two input data are in different time units");
      dftEmpty(tis); dftEmpty(inp); dftEmpty(&tmpdft);
      return(4);
    }
    /* Sort the data by increasing sample times */
    dftSortByFrame(&tmpdft);
 
    /* Copy to input data */
    if(tmpdft.voiNr>1) {
      fprintf(wfp, "Warning: using only first TAC in %s\n", fname);
      tmpdft.voiNr=1;
    }
    if(loginfo!=NULL && verbose>1)
      fprintf(loginfo, "interpolating %d samples into %d samples.\n",
        tmpdft.frameNr, inp->frameNr);
    ret=dftInterpolateInto(&tmpdft, inp, status, verbose);
    if(ret) {
      if(loginfo!=NULL && verbose>0)
        fprintf(loginfo, "dftInterpolateInto() := %d\n", ret);
      dftEmpty(tis); dftEmpty(inp); dftEmpty(&tmpdft); return(4);
    }
    /* Remove the originally read blood data */
    dftEmpty(&tmpdft);
 
  } // next input file
 
  /* Set time unit to min */
  if(loginfo!=NULL && verbose>1)fprintf(loginfo, "setting time units to min.\n");
  ret=dftTimeunitConversion(tis, TUNIT_MIN); if(ret)
    fprintf(wfp, "Warning: check that regional data times are in minutes.\n");
  ret=dftTimeunitConversion(inp, TUNIT_MIN); if(ret)
    fprintf(wfp, "Warning: check that input data times are in minutes.\n");
 
  /* 
   *  Check the input data
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "checking input data\n");
  if(inp->frameNr<4) {
    if(status!=NULL) sprintf(status, "%s contains too few samples", inputfile1);
    dftEmpty(tis); dftEmpty(inp); return(4);
  }
  /* Check for NaN's */
  ret=dft_nr_of_NA(inp); if(ret>0) {
    if(status!=NULL) sprintf(status, "missing sample(s) in data");
    dftEmpty(tis); dftEmpty(inp); return(4);
  }
  /* Check that input and tissue time ranges are about the same */
  if(tis->x[tis->frameNr-1]>10.0*inp->x[inp->frameNr-1] ||
     tis->x[tis->frameNr-1]<0.10*inp->x[inp->frameNr-1]) {
    fprintf(wfp, "Warning: you might need to check the sample time units.\n");
  }
  /* Check and give a warning, if the first value is the highest */
  for(fi=1, n=0, f=inp->voi[0].y[0]; fi<inp->frameNr; fi++)
    if(inp->voi[0].y[fi]>=f) {f=inp->voi[0].y[fi]; n=fi;}
  if(n<2) fprintf(wfp, "Warning: check the first input sample values.\n");
 
  /* Check the tissue and blood TAC concentration units */
  ret=dftUnitConversion(inp, petCunitId(tis->unit));
  if(ret!=0) {fprintf(wfp, "Note: check the units of input and tissue data.\n");}
 
  /*
   *  Check and set fit time length
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "checking and setting fit time length\n");
  /* Set fit duration */
  starttime=0; endtime=*fitdur;
  *fitframeNr=fittime_from_dft(tis, &starttime, &endtime, &first, &last, verbose);
  if(loginfo!=NULL && verbose>1) {
    fprintf(loginfo, "tis.frameNr := %d\n", tis->frameNr);
    fprintf(loginfo, "starttime := %g\n", starttime);
    fprintf(loginfo, "endtime := %g\n", endtime);
    fprintf(loginfo, "first := %d\n", first);
    fprintf(loginfo, "last := %d\n", last);
    fprintf(loginfo, "fitframeNr := %d\n", *fitframeNr);
  }
  *fitdur=endtime;
  /* Check that input data does not end much_before fitdur */
  if(inp->timetype==DFT_TIME_STARTEND) f=inp->x2[inp->frameNr-1]; 
  else f=inp->x[inp->frameNr-1];
  if(*fitdur>1.2*f) {
    if(status!=NULL) sprintf(status, "input TAC is too short");
    dftEmpty(inp); dftEmpty(tis); return(5);
  }
  /* Cut off too many input samples to make calculation faster */
  f=*fitdur+1.0; // One minute extra to allow time delay correction
  if(loginfo!=NULL && verbose>0) printf("Input TAC cutoff at %g min\n", f); 
  for(fi=0; fi<inp->frameNr; fi++) if(inp->x[fi]>f) break;
  if(fi<inp->frameNr) fi++; 
  inp->frameNr=fi;
  if(inp->frameNr<4) {
    if(status!=NULL) sprintf(status, "too few samples in specified fit duration.\n");
    dftEmpty(inp); dftEmpty(tis); return(5);
  }
  if(loginfo!=NULL && verbose>1) {
    fprintf(loginfo, "dft.frameNr := %d\ninp.frameNr := %d\nfitdur := %g\n",
    tis->frameNr, inp->frameNr, *fitdur);
    fprintf(loginfo, "fitframeNr := %d\n", *fitframeNr);
  }
 
  if(status!=NULL) sprintf(status, "ok");
  return 0;
}

◆ dftReadReference()

int dftReadReference	(	DFT *	tissue,
		char *	filename,
		int *	filetype,
		int *	ref_index,
		char *	status,
		int	verbose )

Read DFT format reference region TAC data and add it into DFT already containing other tissue TACs, if reference region TAC(s) are be given in a separate file. Alternatively the reference region name can be given, which will then be selected from existing tissue TACs. Reference TACs are marked in DFT with sw=2 (best name match) or sw=1; if reference region file contains several TACs then the one which contains name 'mean' or 'avg' or has shortest total name length is given value sw=2. When necessary, reference data is interpolated and units converted to match the existing tissue DFT. Reference TAC is also integrated into y2[].

Returns: Returns the number of reference TACs, and <=0 in case of an error.

See also: dftReadinput, dftRead, dftReadModelingData, imgReadModelingData

Parameters

tissue	PET data containing existing tissue TACs, possibly also ref regions.
filename	Reference TAC filename, or region name in previous data.
filetype	Type of input, as found out here; 1 and 2 =tac, 3=fit, 5=region name; enter NULL, if not needed.
ref_index	Index of the best reference region; enter NULL if not needed.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 204 of file dftinput.c.

  {
  int ri, n, ret;
  DFT temp;
  FILE *fp;
  int ftype=0;
 
  if(verbose>0) printf("dftReadReference(tis, %s, type, i, status, %d)\n", filename, verbose);
 
  /* Check that input is ok */
  if(tissue==NULL || filename==NULL || strlen(filename)<1) {
    if(status!=NULL) strcpy(status, "program error");
    return -1;
  }
  if(tissue->frameNr<1 || tissue->voiNr<1) {
    if(status!=NULL) strcpy(status, "no pet data");
    return -2;
  }
 
  /* Check if we can identify the reference as an unsupported file */
  ret=dftFormat(filename);
  if(ret==DFT_FORMAT_FIT) {
    if(filetype!=NULL) *filetype=3;
    if(status!=NULL) strcpy(status, "cannot read fit file");
    return -3;
  }
 
 
  /* Can we open it as a file? */
  fp=fopen(filename, "r");
  if(fp!=NULL) {
    fclose(fp);
    /* Try to read the reference as a TAC file */
    dftInit(&temp);
    if((ret=dftRead(filename, &temp))!=0) {
      if(status!=NULL) sprintf(status, "cannot read file (%d)", ret);
      return -4;
    }
    if(filetype!=NULL) *filetype=1;
 
    /* Convert ref time units to the same as in tissue data */
    (void)dftTimeunitConversion(&temp, tissue->timeunit);
    /* Check the tissue and ref TAC concentration units */
    if((ret=dftUnitConversion(&temp, petCunitId(tissue->unit)))!=0) {
      sprintf(status, "check the units of reference and tissue data");
    }
    /* Interpolate and integrate data to pet times */
    /* this also verifies that ref data does not need extensive extrapolation */
    ret=dftInterpolateInto(&temp, tissue, status, verbose);
    if(ret!=0) {dftEmpty(&temp); return -5;}
    /* Set switches to define which are reference regions and which are not */
    for(ri=0; ri<tissue->voiNr-temp.voiNr; ri++) tissue->voi[ri].sw=0;
    for(; ri<tissue->voiNr; ri++) tissue->voi[ri].sw=1;
    /* Find the best reference region */
    n=temp.voiNr; dftEmpty(&temp);
    if(n==1) {
      ri=tissue->voiNr-n;
    } else {
      if((ri=dftSelectBestReference(tissue))<0) {
        sprintf(status, "cannot select the best reference region");
        dftEmpty(&temp); return -6;
      }
    }
    tissue->voi[ri].sw=2; if(ref_index!=NULL) *ref_index=ri;
    if(verbose>1) printf("selected ref region := %s\n", tissue->voi[ri].name);
    if(status!=NULL) sprintf(status, "%d reference curve(s) read", n);
 
    return n;
 
  } // reference file was read and processed 
 
 
  /* So it's not a file, at least an accessible file, but is it region name? */
  ftype=5; if(filetype!=NULL) *filetype=ftype;
  /* Select ROIs that match the specified input name */
  n=dftSelectRegions(tissue, filename, 1);
  if(verbose>1) printf("nr of ref regions := %d/%d\n", n, tissue->voiNr);
  if(n<=0) {
    if(status!=NULL) sprintf(status, "cannot find region");
    return -7;
  }
  if(n==tissue->voiNr && tissue->voiNr>1) {
    if(status!=NULL) sprintf(status, "all regions do match");
    return -8;
  }
 
  /* Try to select the best reference ROI */
  ri=dftSelectBestReference(tissue); if(ri<0) {
    if(status!=NULL) sprintf(status, "cannot select the best reference region");
    return -9;
  }
  tissue->voi[ri].sw=2; if(ref_index!=NULL) *ref_index=ri;
  if(verbose>1) printf("selected ref region := %s\n", tissue->voi[ri].name);
 
  /* Calculate integrals */
  for(ri=0; ri<tissue->voiNr; ri++) if(tissue->voi[ri].sw>0) {
    if(tissue->timetype==DFT_TIME_STARTEND) {
      ret=petintegrate(tissue->x1, tissue->x2, tissue->voi[ri].y,
              tissue->frameNr, tissue->voi[ri].y2, NULL);
    } else {
      ret=integrate(tissue->x, tissue->voi[ri].y, tissue->frameNr, tissue->voi[ri].y2);
    }
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot integrate input");
      return(-11);
    }
  }
 
  if(status!=NULL) sprintf(status, "%d reference curve(s) read", n);
  return(n);
}

◆ dftRobustMinMaxTAC()

int dftRobustMinMaxTAC	(	DFT *	dft,
		int	tacindex,
		double *	minx,
		double *	maxx,
		double *	miny,
		double *	maxy,
		int *	mini,
		int *	maxi,
		int *	mins,
		int *	maxs,
		int	verbose )

Robust search of the min and max values of DFT TAC data. Data may contain NaN's, and individual outliers are not taken as min or max.

See also: dftMinMaxTAC, dftMinMax

Returns: Returns 0 if successful.

Parameters

dft	Pointer to the DFT TAC data to search.
tacindex	Index of the only TAC which is searched for min and max; <0 if all.
minx	Pointer to X at TAC min; set to NULL if not needed.
maxx	Pointer to X at TAC max; set to NULL if not needed.
miny	Pointer to min Y; set to NULL if not needed.
maxy	Pointer to max Y; set to NULL if not needed.
mini	Index of min TAC; set to NULL if not needed.
maxi	Index of max TAC; set to NULL if not needed.
mins	Index of min sample; set to NULL if not needed.
maxs	Index of max sample; set to NULL if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 608 of file dftinput.c.

  {
  int ri, fi, i1, i2, s1, s2;
  double x, x1, x2, y1, y2;
  int n, runnh, maxrunnh, runns, maxrunns;
  double ym;
  int run1h, run2h, maxrun1h, maxrun2h, run1s, run2s, maxrun1s, maxrun2s;
 
  if(verbose>0)
    printf("dftRobustMinMaxTAC(dft, %d, minx, maxx, miny, maxy, mini, maxi, mins, maxs, %d)\n",
           tacindex, verbose);
 
  if(dft==NULL) return(1);
  if(tacindex>=dft->voiNr) return(2);
  if(dft->voiNr<1 || dft->frameNr<1) return(3);
 
  /* One TAC at a time */
  double list[dft->frameNr];
  x1=x2=y1=y2=nan(""); i1=i2=s1=s2=0;
  for(ri=0; ri<dft->voiNr; ri++) {
    if(tacindex>=0 && ri!=tacindex) continue;
    //fprintf(stderr, "Region: %s\n", dft->voi[ri].name);
    /* Calculate median of y values */
    for(fi=n=0; fi<dft->frameNr; fi++) {
      if(isnan(dft->voi[ri].y[fi])) continue;
      if(dft->timetype==DFT_TIME_STARTEND) {
        if(isnan(dft->x1[fi]) || isnan(dft->x2[fi])) continue;}
      else {if(isnan(dft->x[fi])) continue;}
      list[n++]=dft->voi[ri].y[fi];
    }
    if(n<10) {
      maxrunnh=maxrunns=dft->frameNr;
      maxrun1h=maxrun1s=0;
      maxrun2h=maxrun2s=dft->frameNr-1;
    } else {
      ym=dmedian(list, n);
      //fprintf(stderr, "median := %g\n", ym);
      /* Search the longest run of values which are larger/smaller than the median */
      runnh=maxrunnh=run1h=run2h=maxrun1h=maxrun2h=0;
      runns=maxrunns=run1s=run2s=maxrun1s=maxrun2s=0;
      for(fi=0; fi<dft->frameNr; fi++) {
        if(isnan(dft->voi[ri].y[fi])) continue;
        if(dft->timetype==DFT_TIME_STARTEND) {
          if(isnan(dft->x1[fi]) || isnan(dft->x2[fi])) continue;}
        else {if(isnan(dft->x[fi])) continue;}
        /* max */
        if(dft->voi[ri].y[fi]>ym) {
          runnh++; if(runnh==1) {run1h=run2h=fi;} else {run2h=fi;}
        } else {
          if(runnh>maxrunnh) {maxrunnh=runnh; maxrun1h=run1h; maxrun2h=run2h;}
          runnh=0;
        }
        /* min */
        if(dft->voi[ri].y[fi]<ym) {
          runns++; if(runns==1) {run1s=run2s=fi;} else {run2s=fi;}
        } else {
          if(runns>maxrunns) {maxrunns=runns; maxrun1s=run1s; maxrun2s=run2s;}
          runns=0;
        }
      }
      /* Get range also from the end part */
      if(runnh>maxrunnh) {maxrunnh=runnh; maxrun1h=run1h; maxrun2h=run2h;}
      if(runns>maxrunns) {maxrunns=runns; maxrun1s=run1s; maxrun2s=run2s;}
      /* If longest range includes just one sample, then take the whole range */
      if(maxrunnh<2) {maxrunnh=dft->frameNr; maxrun1h=0; maxrun2h=dft->frameNr-1;}
      if(maxrunns<2) {maxrunns=dft->frameNr; maxrun1s=0; maxrun2s=dft->frameNr-1;}
    }
    if(verbose>12) {
      fprintf(stderr, "longest run for max: %g - %g\n", dft->x[maxrun1h], dft->x[maxrun2h]);
      fprintf(stderr, "longest run for min: %g - %g\n", dft->x[maxrun1s], dft->x[maxrun2s]);
    }
    /* Inside the range, search for max */
    for(fi=maxrun1h; fi<=maxrun2h; fi++) {
      if(isnan(dft->voi[ri].y[fi])) continue;
      if(dft->timetype==DFT_TIME_STARTEND) {
        if(isnan(dft->x1[fi]) || isnan(dft->x2[fi])) continue;
        x=0.5*(dft->x1[fi]+dft->x2[fi]);
      } else {
        if(isnan(dft->x[fi])) continue;
        x=dft->x[fi];
      }
      if(isnan(y2) || y2<dft->voi[ri].y[fi]) {y2=dft->voi[ri].y[fi]; i2=ri; x2=x; s2=fi;}
    }
    /* Inside the range, search for min */
    for(fi=maxrun1s; fi<=maxrun2s; fi++) {
      if(isnan(dft->voi[ri].y[fi])) continue;
      if(dft->timetype==DFT_TIME_STARTEND) {
        if(isnan(dft->x1[fi]) || isnan(dft->x2[fi])) continue;
        x=0.5*(dft->x1[fi]+dft->x2[fi]);
      } else {
        if(isnan(dft->x[fi])) continue;
        x=dft->x[fi];
      }
      if(isnan(y1) || y1>dft->voi[ri].y[fi]) {y1=dft->voi[ri].y[fi]; i1=ri; x1=x; s1=fi;}
    }
 
  } /* next TAC */
 
  if(minx!=NULL) {if(isnan(x1)) return(11); else *minx=x1;}
  if(maxx!=NULL) {if(isnan(x2)) return(12); else *maxx=x2;}
  if(miny!=NULL) {if(isnan(y1)) return(13); else *miny=y1;}
  if(maxy!=NULL) {if(isnan(y2)) return(14); else *maxy=y2;}
  if(mini!=NULL) {if(isnan(y1)) return(13); else *mini=i1;}
  if(maxi!=NULL) {if(isnan(y2)) return(14); else *maxi=i2;}
  if(mins!=NULL) {if(isnan(y1)) return(13); else *mins=s1;}
  if(maxs!=NULL) {if(isnan(y2)) return(14); else *maxs=s2;}
  return(0);
}

◆ dftTimeIntegral()

int dftTimeIntegral	(	DFT *	dft,
		double	t1,
		double	t2,
		DFT *	idft,
		int	calc_mode,
		char *	status,
		int	verbose )

Integration of regional TAC data from time1 to time2, i.e. AUC(t1,t2).

You may want to test the time range before applying this routine, because this function accepts relatively large extrapolation large.

See also: dftInterpolateInto, dftInterpolate, dftInterpolateCheckStart, dftInterpolateCheckEnd

Returns: Returns 0 when call was successful, and >0 in case of an error.

Parameters

dft	Regional TAC data in DFT struct. Number of samples must be at least one. If only one sample, then the integration time range must match with the possible frame start and times. Frames do not have to be continuous in time.
t1	Time where to start integration (same unit as in TACs)
t2	Time to stop integration (same unit as in TACs); must be higher than t1, except that t1=t2 is acceptable when calc_mode=average
idft	Pointer to initiated but empty AUC DFT data (output)
calc_mode	Calculate integral or average: 0=integral, 1=average
status	Pointer to a string (allocated for at least 128 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 382 of file dftint.c.

  {
  int ri, fi, ret, f1, f2;
  double fdur, aucroi, t[2], auc[2];
  double accept_tdif=1.0; // in seconds
 
  if(verbose>0)
    printf("dftTimeIntegral(dft, %g, %g, idft, %d, status, %d)\n", t1, t2, calc_mode, verbose);
 
  /* Check the arguments */
  if(status!=NULL) sprintf(status, "program error");
  if(dft==NULL || idft==NULL || t1<0.0 || t2<0.0) return STATUS_FAULT;
  fdur=t2-t1; if(fdur<0.0) return STATUS_FAULT;
  if(fdur==0.0 && (calc_mode!=1 || dft->frameNr>1)) return STATUS_FAULT;
  if(dft->frameNr<1 || dft->voiNr<1) return STATUS_FAULT;
 
  /* Set time unit based things */
  if(dft->timeunit==TUNIT_MIN) accept_tdif/=60.0;
 
  /* Create the DFT struct for AUC values */
  ret=dftAllocateWithHeader(idft, 1, dft->voiNr, dft);
  if(ret!=0) {
    if(status!=NULL) sprintf(status, "cannot setup AUC data (%d)", ret);
    return STATUS_FAULT;
  }
  /* 'frame' start and end time are taken from integration time */
  idft->timetype=DFT_TIME_STARTEND; 
  if(calc_mode==1 && fdur==0.0) idft->timetype=DFT_TIME_MIDDLE;
  idft->_type=DFT_FORMAT_STANDARD;
  /* Initially integral is zero */
  for(ri=0; ri<idft->voiNr; ri++) idft->voi[ri].y[0]=0.0;
  //dftPrint(dft);
 
  /* Different paths if frame start and end times are taken into account or not */
  if(dft->timetype==DFT_TIME_STARTEND) {
 
    /* Check that time range matches with pet frame */
    if(dft->frameNr==1) { // static data (one frame only)
      /* check that specified times match with frame times, but
         accept 1 s differences */
      if(fabs(dft->x2[0]-t2)>accept_tdif || fabs(dft->x1[0]-t1)>accept_tdif) {
        if(status!=NULL) {
           strcpy(status, "for static data the integration time range must");
           strcat(status, " be exactly as long as the scan");
        }
        return STATUS_FAULT;
      }
      /* Set the times, in case there was a small difference */
      dft->x2[0]=t2; dft->x1[0]=t1;
    } else { // Dynamic data (more than one frame)
      /* First frame must start before 1/3 of the integration time */
      /* Last frame must end after 2/3 of integration time */
      if(dft->x1[0]>(0.66*t1+0.34*t2) ||
         dft->x2[dft->frameNr-1]<(0.34*t1+0.66*t2)) {
        if(status!=NULL)
          strcpy(status, "integration time range oversteps data range");
        return STATUS_FAULT;
      }
    }
 
    /* Get the first and last frame index that resides inside integration time*/
    if(dft->frameNr==1) {
      f1=f2=0;
    } else {
      for(fi=0, f1=f2=-1; fi<dft->frameNr; fi++) {
        if(f1<0) {if(dft->x1[fi]>=t1 && dft->x2[fi]<=t2) f1=f2=fi;}
        if(f1>=0) {if(t2>=dft->x2[fi]) f2=fi;}
      }
    }
    if(verbose>1) printf("f1=%d f2=%d\n", f1, f2);
    if(f1>=0 && f2>=0) {
      /* Integrate over the frames that are included in time range as a whole */
      fi=f1;
      for(ri=0; ri<dft->voiNr; ri++) {
        idft->voi[ri].y[0]+=(dft->x2[fi]-dft->x1[fi])*dft->voi[ri].y[fi];
      }
      for(fi=f1+1; fi<=f2; fi++) {
        for(ri=0; ri<dft->voiNr; ri++) {
          idft->voi[ri].y[0]+=(dft->x2[fi]-dft->x1[fi])*dft->voi[ri].y[fi];
        }
        /* Check whether frames are contiguous */
        if(dft->x1[fi]==dft->x2[fi-1]) continue;
        /* When not, calculate the area of an imaginary frame */
        double x, a;
        x=(dft->x1[fi]+dft->x2[fi-1])/2.0;
        for(ri=0; ri<dft->voiNr; ri++) {
          a=(dft->x1[fi]-dft->x2[fi-1])*
            (dft->voi[ri].y[fi]-(dft->voi[ri].y[fi]-dft->voi[ri].y[fi-1])
             *(dft->x2[fi]+dft->x1[fi]-2.0*x)
             /(dft->x2[fi]+dft->x1[fi]-dft->x2[fi-1]-dft->x1[fi-1]));
          idft->voi[ri].y[0]+=a;
        }
      }
 
      /* If necessary, add the partial integrals */
      if(dft->x1[f1]>t1) {
        t[0]=t1; t[1]=dft->x1[f1];
        if(verbose>5) printf("t[0]=%g t[1]=%g\n", t[0], t[1]); 
        for(ri=0; ri<dft->voiNr; ri++) {
          ret=interpolate(dft->x, dft->voi[ri].y, dft->frameNr, t, NULL, auc, NULL, 2);
          if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
          idft->voi[ri].y[0]+=aucroi;
        }
      }
      if(t2>dft->x2[f2]) {
        t[0]=dft->x2[f2]; t[1]=t2;
        if(verbose>5) printf("t[0]=%g t[1]=%g\n", t[0], t[1]);
        for(ri=0; ri<dft->voiNr; ri++) {
          ret=interpolate(dft->x, dft->voi[ri].y, dft->frameNr, t, NULL, auc, NULL, 2);
          if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
          idft->voi[ri].y[0]+=aucroi;
        }
      }
 
    } else { // no full frames inside integration range
 
      t[0]=t1; t[1]=t2; 
      for(ri=0; ri<dft->voiNr; ri++) {
        ret=interpolate(dft->x, dft->voi[ri].y, dft->frameNr, t, NULL, auc, NULL, 2);
        if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
        idft->voi[ri].y[0]+=aucroi;
      }
 
    }
 
    /* Set output image time frame */
    idft->x2[0]=t2; idft->x1[0]=t1; idft->x[0]=0.5*(t1+t2);
 
  } else if(dft->timetype==DFT_TIME_MIDDLE) { // do not consider frame lengths
 
    /* If average calculation was required, and sample times match the required
       time range, then directly take the values; otherwise, calculate AUC */
    if(calc_mode==1 && dft->x[0]==t1 && dft->x[0]==t2 && dft->frameNr==1) {
      for(ri=0; ri<dft->voiNr; ri++)
        idft->voi[ri].y[0]=dft->voi[ri].y[0];
    } else {
 
      /* First sample time must be before 1/3 of the integration time */
      /* Last sample time must end after 2/3 of integration time */
      if(dft->x[0]>(0.66*t1+0.34*t2) ||
         dft->x[dft->frameNr-1]<(0.34*t1+0.66*t2)) {
        if(status!=NULL) sprintf(status, "integration time range oversteps data range");
        return STATUS_FAULT;
      }
 
      t[0]=t1; t[1]=t2; 
      if(verbose>5) printf("t[0]=%g t[1]=%g\n", t[0], t[1]); 
      for(ri=0; ri<dft->voiNr; ri++) {
        ret=interpolate(dft->x, dft->voi[ri].y, dft->frameNr, t, NULL, auc, NULL, 2);
        if(ret) idft->voi[ri].y[0]=0.0; else idft->voi[ri].y[0]=auc[1]-auc[0];
      }
    }
    /* Set output image time frame */
    idft->x2[0]=t2; idft->x1[0]=t1; idft->x[0]=0.5*(t1+t2);
 
  } else {
    if(status!=NULL)
      sprintf(status, "frame mid times or start and end times required");
    return STATUS_FAULT;
  }
 
  /* If required, then calculate average by dividing integral with time */
    /* Set units accordingly */
  if(calc_mode!=0) { // average
    if(fdur>0.0)
      for(ri=fi=0; ri<idft->voiNr; ri++) 
        idft->voi[ri].y[fi]/=fdur;
    if(status!=NULL) sprintf(status, "average TAC [%g,%g] calculated", t1, t2);
  } else { // integral
    int unit;
    unit=petCunitId(idft->unit);
    strcpy(idft->unit, dftUnit(CUNIT_UNKNOWN));
    if(unit==CUNIT_KBQ_PER_ML) {
      if(idft->timeunit==TUNIT_MIN)
        strcpy(idft->unit, dftUnit(CUNIT_MIN_KBQ_PER_ML));
      else if(idft->timeunit==TUNIT_SEC)
        strcpy(idft->unit, dftUnit(CUNIT_SEC_KBQ_PER_ML));
    }
    if(status!=NULL) sprintf(status, "TAC integral [%g,%g] calculated", t1, t2);
  }
  
  return(0);
}

◆ dftVerifyPeak()

int dftVerifyPeak	(	DFT *	dft,
		int	index,
		int	verbose,
		char *	status )

Verify that specified (input) TAC contains peak and that it does not start too late to get reliable estimate of AUC.

Returns: Returns 0 in case data is suitable, < 0 if there may be a problem, 1 in case of an error, and 2 if peak seems to be missed.

See also: dftReadModelingData, imgReadModelingData

Parameters

dft	Pointer to TAC data which is verified.
index	Index of TAC to be verified; <0 in case all are (separately) verified.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 747 of file dftinput.c.

  {
  int ri, fi, ret, mini, maxi, starti, warn=0;
  double minx, maxx, miny, maxy, startx, dif, lowesty;
 
  if(verbose>0) printf("dftVerifyPeak(dft, %d, %d)\n", index, verbose);
  /* Check the arguments */
  if(status!=NULL) strcpy(status, "program error");
  if(dft==NULL || dft->frameNr<1 || dft->voiNr<1) return 1;
  if(index>=dft->voiNr) return 1;
 
  /* If less than 3 samples, then this can not be suitable as input TAC */
  if(dft->frameNr<3) {
    if(status!=NULL) strcpy(status, "too few samples");
    return 2;
  }
 
  /* Make sure data is sorted by increasing sample time */
  dftSortByFrame(dft);
 
  /* Go through all TACs that are requested */
  for(ri=0; ri<dft->voiNr; ri++) {
    if(index>=0 && index!=ri) continue;
    if(verbose>1) printf("checking region %d: %s\n", 1+ri, dft->voi[ri].name);
 
    /* 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);
      if(status!=NULL) strcpy(status, "invalid TAC");
      return 1;
    }
 
    /* Check that there are no big negative values */
    if(maxy<=0.0) {
      if(verbose>0) printf("TAC has no positive values.\n");
      if(status!=NULL) strcpy(status, "no positive TAC values");
      return 2;
    }
    if(miny<0.0) {
      if((-miny)>0.40*maxy) {
        if(verbose>0) printf("TAC has high negative value(s).\n");
        if(status!=NULL) strcpy(status, "too high negative TAC values");
        return 2;
      }
      if((-miny)>0.02*maxy) {
        if(verbose>1) printf("TAC has negative value(s).\n");
        warn++;
      }
    }
 
    /* Get the first sample time, considering missing values */
    startx=1.0E+10; starti=dft->frameNr-1;
    for(fi=0; fi<dft->frameNr; fi++) {
      if(isnan(dft->voi[ri].y[fi])) continue;
      if(dft->timetype==DFT_TIME_STARTEND) {
        if(isnan(dft->x1[fi])) continue;
        if(isnan(dft->x2[fi])) continue;
        startx=dft->x1[fi]; starti=fi; break;
      } else {
        if(isnan(dft->x[fi])) continue;
        startx=dft->x[fi]; starti=fi; break;
      }
    }
    if(verbose>2) printf("first time sample at %g\n", startx);
 
    /* If the first sample is the peak, then check that sample time is not too late */
    if(maxi==starti) {
      if(verbose>2) printf("Peak at the first sample.\n");
      if(status!=NULL) strcpy(status, "input TAC should start at time zero");
      /* Calculate peak time to initial sample frequency ratio */
      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]);
        for(fi=maxi+1; fi<dft->frameNr; fi++)
          if(!isnan(dft->x[fi]) && dft->x[fi]>dft->x[maxi]) break;
        if(fi==dft->frameNr) dif=1.0E+10;
        else dif=dft->x[maxi]/(dft->x[fi]-dft->x[maxi]);
      }
      if(dif>0.3) {
        if(verbose>0) printf("Peak at the first sample which is not at zero.\n");
        if(verbose>1) printf("dif := %g\n", dif);
      }
      if(dif>5.0) {
        return 2;
      } else if(dif>1.0) {
        /* Not bad, if peak is still much higher than tale */
        if(maxy>20.*miny) warn++; else return 2;
        if(verbose>1) printf("good peak/tail -ratio\n");
      } else if(dif>0.3) warn++;
    }
 
    /* Search the lowest value before the peak */
    lowesty=dft->voi[ri].y[starti];
    for(fi=starti+1; fi<maxi; fi++)
      if(dft->voi[ri].y[fi]<lowesty) lowesty=dft->voi[ri].y[fi];
    if(verbose>2) printf("lowest value before peak: %g\n", lowesty);
 
    /* If first sample is closer to peak and relatively high, then that is or may be a problem */
    if(maxi>starti && startx>0.001 && startx>0.75*maxx) {
      if(dft->voi[ri].y[starti]>0.66*maxy && lowesty>0.05*maxy && mini>maxi) {
        if(verbose>0) printf("The first sample is relatively late and high.\n");
        if(status!=NULL) strcpy(status, "input TAC should start at time zero");
        if(verbose>2) {
          printf("startx=%g\n", startx);
          printf("starty=%g\n", dft->voi[ri].y[starti]);
          printf("maxx=%g\n", maxx);
          printf("maxy=%g\n", maxy);
        }
        return 2;
      }
    }
    if(maxi>starti && startx>0.001 && startx>0.5*maxx) {
      if(dft->voi[ri].y[starti]>0.5*maxy && lowesty>0.05*maxy && mini>maxi) {
        if(verbose>1) printf("The first sample is relatively late and high.\n");
        warn++;
      }
    }
    if(verbose>5) {
      printf("startx=%g\n", startx);
      printf("starty=%g\n", dft->voi[ri].y[starti]);
      printf("maxx=%g\n", maxx);
      printf("maxy=%g\n", maxy);
    }
 
    /* 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");
      if(status!=NULL) strcpy(status, "input TAC peak missing");
      return 2;
    }
    if(maxy<5.0*miny) {
      if(verbose>1) printf("TAC does not have a clear peak.\n");
      warn++;
    }
 
  } // next TAC
 
  if(verbose>0 && warn>0) printf("%d warning(s)\n", warn);
  if(warn>0) {
    if(status!=NULL) strcpy(status, "input TAC is not optimal");
    return -1;
  }
  if(status!=NULL) strcpy(status, "input TAC ok");
  return 0;
}

Referenced by dftFixPeak(), dftReadinput(), and imgReadModelingData().

◆ dftWeightByFreq()

int dftWeightByFreq ( DFT * dft )

Add weights based on sample frequency or frame length.

Existing weights are overwritten.

Returns: Returns 0 when successful, otherwise <>0.

See also: dftRead, dftSortByFrame, sifWeight, sifWeightByFrames, dftDecayCorrection, imgSetWeights, dftWSampleNr

Parameters

dft	Samples/frames must be sorted by sample time, but duplicate samples are allowed

Definition at line 20 of file weight_model.c.

  {
  int fi, fi1, fi2;
  double f, t, t1, t2, sumw=0.0;
 
  if(dft==NULL || dft->frameNr<1) return 1;
  if(dft->frameNr==1) {
    dft->w[0]=1.0;
    dft->isweight=1;
    return 0;
  }
 
  if(dft->timetype==DFT_TIME_STARTEND) { /* weights by frame lengths */
 
    for(fi=0; fi<dft->frameNr; fi++) {
      dft->w[fi]=dft->x2[fi]-dft->x1[fi]; 
      sumw+=dft->w[fi];
    }
    /* scale weights so that sum of weights equals sample number */
    sumw/=(double)dft->frameNr;
    for(fi=0; fi<dft->frameNr; fi++) dft->w[fi]/=sumw;
 
  } else if(dft->timetype==DFT_TIME_MIDDLE) { /* weights by sample distance */
 
    for(fi=0; fi<dft->frameNr; fi++) {
      t=t1=t2=dft->x[fi]; //printf("fi=%d t=%g\n", fi, t);
      /* Find the closest sample time before this one */
      for(fi1=fi; fi1>=0; fi1--) {t1=dft->x[fi1]; if(t1<t) break;} 
      /* Find the closest sample time after this one */ 
      for(fi2=fi; fi2<dft->frameNr; fi2++) {t2=dft->x[fi2]; if(t2>t) break;} 
      //printf("  t1=%g t2=%g\n", t1, t2);
      /* Mean sample distance */
      f=0.0;
      if(t1<t) f+=t-t1; else f+=t2-t;
      if(t2>t) f+=t2-t; else f+=t-t1;
      f*=0.5; if(f<=0.0) f=1.0;
      /* Set initial weight */
      //printf("    f=%g\n", f);
      dft->w[fi]=f; sumw+=f;
    }
    /* Scale weights */
    sumw/=(double)dft->frameNr;
    for(fi=0; fi<dft->frameNr; fi++) dft->w[fi]/=sumw;
 
  } else
    return 2;
 
  dft->isweight=1;
  return(0);
}

◆ dftWSampleNr()

int dftWSampleNr ( DFT * tac )

Get the number of samples in DFT that have weight > 0. Missing (NaN) sample values are included as long as weight is not missing. If weights are not set, then nr of all samples is returned.

Returns: The number of samples with weight above zero.

See also: dftWeightByFreq, aicSS, parFreeNr

Parameters

tac	Pointer to the DFT struct.

Definition at line 147 of file weight_model.c.

  {
  if(tac==NULL || tac->voiNr<1 || tac->frameNr<1) return(0);
  if(tac->isweight==0) return(tac->frameNr);
  int n=0;
  for(int i=0; i<tac->frameNr; i++) if(tac->w[i]>0.0) n++;
  return(n);
}

◆ extrapolate_monoexp()

int extrapolate_monoexp	(	DFT *	dft,
		double *	fittime,
		int *	min_nr,
		int	max_nr,
		double	mintime,
		double	extr_to,
		DFT *	ext,
		FILE *	loginfo,
		char *	status )

Extrapolation of exponentially decreasing tail of PET radiotracer plasma curves. This is accomplished by fitting line to the end-part of the plot of the natural logarithm of tracer concentration against time. By default, the included data points are determined based on maximum adjusted R^2 from at least three points; fit to the larger number of points is used in case difference in adjusted R^2 values is not larger than 0.0001. This function is used and tested with program extrapol.

Returns: Returns 0 when successful, otherwise <>0.

Parameters

dft	Pointer to original TAC data.
fittime	By default, the search for the best line fit is started from the last sample towards the first sample; set fit time to -1 to use this default. However, if the end phase is unreliable or very noisy, you may want to set fittime to include only certain time range from the beginning. Function will write here the fittime that was actually used.
min_nr	The minimum number of samples used in searching the best fit; at least 2, but 3 is recommended. If data is very noisy, then this number may need to be increased. Function will write here the nr of samples that was actually used. This can be used as an alternative to mintime or in addition to it.
max_nr	The maximum number of samples used in searching the best fit; must be higher than min_nr, or set to -1 to not to limit the number.
mintime	Minimum time range used in searching the best fit. If data is very noisy, then this may need to be set, otherwise setting mintime to -1 will use the default. This can be used as an alternative to min_nr or in addition to it.
extr_to	Last extrapolated sample time in same time units than in the data.
ext	Pointer to data for extrapolated TACs. Struct must be initiated. Any existing data is deleted.
loginfo	Give file pointer (for example stdout) where log information is printed; NULL if not needed.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 22 of file extrapolate.c.

  {
  int fi, ret, n, to, from, to_min, from_min, orig_max_nr;
  int fitNr, snr;
  double kel, c0=0.0, *cx, *cy;
  double slope, slope_sd, ic, ic_sd, r, f, adj_r2, adj_r2_max;
  double extr_sampl;
 
 
  /* Check the input */
  if(status!=NULL) sprintf(status, "program error");
  if(dft==NULL || ext==NULL) return -1;
  if(*min_nr<2) *min_nr=3;  // Setting erroneous value to a recommended value
  if(max_nr>-1 && max_nr<*min_nr) return -2;
  orig_max_nr=max_nr;
 
  /* Set the last sample to be fitted */
  fitNr=dft->frameNr;
  if(*fittime>0.0) {
    while(dft->x[fitNr-1]>*fittime && fitNr>0) fitNr--;
    if(loginfo!=NULL) fprintf(loginfo, "fitNr := %d\nTime range := %g - %g\n", 
      fitNr, dft->x[0], dft->x[fitNr-1]);
  }
  *fittime=dft->x[fitNr-1];
  /* Check that there are at least 3 samples */
  if(fitNr<3) { /* min_nr is checked later */
    if(status!=NULL) sprintf(status, "too few samples for extrapolation");
    return(2);
  }
  
  /* If mintime was set, then get min_nr based on that */
  if(mintime>0.0) {
    fi=0; while(dft->x[fi]<(*fittime-mintime) && fi<fitNr) fi++;
    if(--fi<=0) {
      if(status!=NULL) sprintf(status, "required minimum fit range too large");
      return(2);
    }
    *min_nr=((fitNr-fi)>*min_nr?fitNr-fi:*min_nr);
  }
  if(loginfo!=NULL) fprintf(loginfo, "final_min_nr := %d\n", *min_nr);
 
 
  /*
   *  Initiate data for extrapolated data
   */
  /* Determine the sampling time for extrapolated range */
  extr_sampl=0.5*(dft->x[dft->frameNr-1]-dft->x[dft->frameNr-3]);
  if(loginfo!=NULL) fprintf(loginfo, "extr_sampl=%g\n", extr_sampl);
  if(extr_sampl<1.0E-8) {
    if(status!=NULL) sprintf(status, "check sample times");
    return(2);
  }
  /* Determine the nr of extrapolated samples */
  f=extr_to-dft->x[dft->frameNr-1];
  if(f<=0.0) {
    if(status!=NULL) sprintf(status, "no extrapolation is needed");
    return(2);
  }
  n=ceil(f/extr_sampl);
  if(loginfo!=NULL) fprintf(loginfo, "  extr_sampl=%g n=%d\n", extr_sampl, n);
  /* Allocate memory */
  dftEmpty(ext);
  ret=dftSetmem(ext, dft->frameNr+n, dft->voiNr);
  if(ret) {
    if(status!=NULL) sprintf(status, "error in memory allocation.\n");
    return(4);
  }
  ext->frameNr=dft->frameNr+n; ext->voiNr=dft->voiNr;
  /* Set sample times */
  for(fi=0; fi<dft->frameNr; fi++) {
    ext->x[fi]=dft->x[fi]; ext->x1[fi]=dft->x1[fi]; ext->x2[fi]=dft->x2[fi];}
  if(dft->timetype==DFT_TIME_MIDDLE) {
    for(; fi<ext->frameNr; fi++) {
      ext->x[fi]=ext->x[fi-1]+extr_sampl;
      ext->x1[fi]=ext->x2[fi-1]; ext->x2[fi]=ext->x[fi]+0.5*extr_sampl;
    }
  } else {
    for(; fi<ext->frameNr; fi++) {
      ext->x1[fi]=ext->x2[fi-1]; ext->x2[fi]=ext->x1[fi]+extr_sampl;
      ext->x[fi]=0.5*(ext->x1[fi]+ext->x2[fi]);
    }
  }
  /* Copy "header" information */
  dftCopymainhdr(dft, ext);
  for(int ri=0; ri<dft->voiNr; ri++) dftCopyvoihdr(dft, ri, ext, ri);
  /* Copy existing values */
  for(int ri=0; ri<dft->voiNr; ri++)
    for(int fi=0; fi<dft->frameNr; fi++)
      ext->voi[ri].y[fi]=dft->voi[ri].y[fi];
 
 
  /*
   *  Make ln transformation for TACs
   */
  if(loginfo!=NULL) fprintf(loginfo, "ln transformation\n");
  for(int ri=0; ri<dft->voiNr; ri++)
    for(int fi=0; fi<dft->frameNr; fi++) {
      if(!isnan(dft->voi[ri].y[fi]) && dft->voi[ri].y[fi]>0.0)
        dft->voi[ri].y2[fi]=log(dft->voi[ri].y[fi]);
      else
        dft->voi[ri].y2[fi]=nan("");
    }
 
 
  /*
   *  Compute best linear fit to the end of ln transformed TACs
   */
  if(loginfo!=NULL) fprintf(loginfo, "linear fitting\n");
  cx=(double*)malloc(2*fitNr*sizeof(double)); if(cx==NULL) {
    if(status!=NULL) sprintf(status, "out of memory");
    return(6);
  }
  cy=cx+fitNr;
  for(int ri=0; ri<dft->voiNr; ri++) {
    kel=0.0; max_nr=orig_max_nr;
    /* Print TAC name, if more than one was found */
    if(dft->voiNr>1 && loginfo!=NULL)
      fprintf(loginfo, "%s :\n", dft->voi[ri].name);
 
    /* Copy appropriate TAC data */
    for(fi=n=0; fi<fitNr; fi++)
      if(dft->x[fi]>0.0 && !isnan(dft->voi[ri].y2[fi])) {
        cx[n]=dft->x[fi]; cy[n++]=dft->voi[ri].y2[fi];}
    if(n<*min_nr) {
      if(status!=NULL) sprintf(status, "check the datafile (%d<%d)", n, *min_nr);
      free(cx); return(7);
    }
    if(max_nr<=0) max_nr=n;
    /* Search the plot range that gives the max adjusted R^2 */
    from_min=to_min=-1; adj_r2_max=-9.99E+99;
    for(from=n-max_nr, to=n-1; from<1+n-*min_nr; from++) {
      snr=(to-from)+1;
      /* Calculation of linear regression using pearson() */
      ret=pearson(
        cx+from, cy+from, snr, &slope, &slope_sd, &ic, &ic_sd, &r, &f
      );
      if(ret==0) {
        adj_r2= 1.0 - ((1.0-r*r)*(double)(snr-1)) / (double)(snr-2); 
      } else {
        adj_r2=-9.99E+99;
      }
      //if(cv<cv_min) {
      if(adj_r2>adj_r2_max+0.0001) {
        adj_r2_max=adj_r2; from_min=from; to_min=to; kel=-slope; c0=exp(ic);
      }
      if(loginfo!=NULL)
        fprintf(loginfo, "  adj_r2=%g from=%d (%g)\n", adj_r2, from, *(cx+from) );
    }
    if(from_min<0) {
      if(status!=NULL) sprintf(status, "check the datafile");
      free(cx); return(7);
    }
    /* Check the sign of slope */
    if(kel>=0.0) {
      /* negative slope i.e. positive kel is ok */
      if(loginfo!=NULL)
         fprintf(loginfo, "  k(el)=%g adj_r2=%g C(0)=%g; %d data points.\n",
           kel, adj_r2_max, c0, (to_min-from_min)+1 );
      /* Calculate the extrapolated values */
      for(fi=fitNr; fi<ext->frameNr; fi++) {
        ext->voi[ri].y[fi]=c0*exp(-kel*ext->x[fi]);
      }
    } else {
      /* If slope is positive, then extrapolate with line f(x)=b+0*t */
      for(fi=fitNr-1, f=0.0, n=0; fi>=0 && n<3; fi--)
        if(!isnan(dft->voi[ri].y[fi])) {f+=dft->voi[ri].y[fi]; n++;}
      if(n>0) f/=(double)n;
      if(status!=NULL)
        sprintf(status, "curve end is not descending; extrapolation with horizontal line determined as avg of %d samples", n);
      /* Calculate the extrapolated values */
      for(fi=fitNr; fi<ext->frameNr; fi++) {
        ext->voi[ri].y[fi]=f;
      }
    }
  } /* next curve */
  free(cx);
 
  if(status!=NULL) sprintf(status, "ok");
  return 0;
}

◆ fittime_from_dft()

int fittime_from_dft	(	DFT *	dft,
		double *	startTime,
		double *	endTime,
		int *	first,
		int *	last,
		int	verbose )

Reset user-defined fit time range to comply with DFT data.

See also: fittime_from_img, getActualSamplenr, dftEndtime, dftReadModelingData

Returns: Returns the number of samples included in the fit range, or <0 in case of an error.

Parameters

dft	Pointer to DFT containing (regional tissue) data; times can be in minutes or seconds, as long as units are defined.
startTime	Pointer containing originally the requested fit start time (min). This is changed to contain the time of the first included frame. Unit must be minutes. Initially, set to <0 to start from the beginning of the data.
endTime	Pointer containing originally the requested fit end time (min). This is changed to contain the time of the last included frame. Unit must be minutes. Initially, set to <0 or to a very large value to reach to the end of data.
first	Function writes the index of the first included sample (frame) here.
last	Function writes the index of the last included sample (frame) here.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 16 of file fittime.c.

  {
  int fi, dataNr;
  
  if(verbose>0) printf("%s(*dft, %g, %g, first, last)\n", __func__, *startTime, *endTime);
  if(dft==NULL || startTime==NULL || endTime==NULL) return(-1);
  if(!isfinite(*startTime)) *startTime=-1.0;
  if(!isfinite(*endTime)) *endTime=-1.0;
  *first=*last=0;
  if(dft->frameNr<=0) {*startTime=*endTime=0.0; return 0;}
  /* Change start and end times to seconds if necessary */
  if(dft->timeunit==TUNIT_SEC) {*startTime*=60.; *endTime*=60.;}
  /* Check that data range is not outside required range */
  if(dft->x[dft->frameNr-1]<*startTime || dft->x[0]>*endTime) {
    *startTime=*endTime=0.0; return 0;}
  /* Get first and last data point */
  for(fi=0, *first=dft->frameNr; fi<dft->frameNr; fi++)
    if(dft->x[fi]>=*startTime) {*first=fi; break;}
  for(*last=fi=*first; fi<dft->frameNr; fi++)
    if(dft->x[fi]<=*endTime) *last=fi; else break;
  if(*first>=dft->frameNr) {*startTime=*endTime=0.0; return 0;}
  /* Correct fit range to frame start and end times */
  *startTime=(dft->timetype==DFT_TIME_STARTEND?dft->x1[*first]:dft->x[*first]);
  *endTime=(dft->timetype==DFT_TIME_STARTEND?dft->x2[*last]:dft->x[*last]);
  /* Calculate the number of data points in the fit range */
  dataNr=(*last)-(*first)+1;
  /* Change start and end times back to minutes if necessary */
  if(dft->timeunit==TUNIT_SEC) {*startTime/=60.; *endTime/=60.;}
  return dataNr;
}

Referenced by dftReadModelingData(), and imgReadModelingData().

◆ fittime_from_img()

int fittime_from_img	(	IMG *	img,
		double *	fittime,
		int	verbose )

Get the IMG frame end time of the last frame that is inside (mid time before) the specified maximum fittime.

See also: fittime_from_dft, imgEndtime, imgReadModelingData

Returns: Returns the number of IMG frames included in the fittime, or <0 in case of an error.

Parameters

img	Pointer to IMG
fittime	Pointer containing originally the fit time maximum, after this the last included IMG frame end time. Unit must be seconds. Initially, set to <0 or to a very large value to include all frames.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 74 of file fittime.c.

  {
  if(verbose>0) printf("%s(*img, %g)\n", __func__, *fittime);
  if(img==NULL || fittime==NULL) return(-1);
  if(!isfinite(*fittime)) *fittime=-1.0;
  if(img->dimt<=0) {*fittime=0.0; return 0;}
  if(*fittime<0) {*fittime=img->end[img->dimt-1]; return img->dimt;}
  int fitdimt=0;
  for(int fi=0; fi<img->dimt; fi++)
    if(img->mid[fi]>*fittime) break; else fitdimt++;
  *fittime=img->end[fitdimt-1];
  if(verbose>1) {
    printf("  fitdimt := %d\n", fitdimt);
    printf("  fittime := %g\n", *fittime);
  }
  return fitdimt;
}

Referenced by imgReadModelingData().

◆ getActualSamplenr()

int getActualSamplenr	(	DFT *	dft,
		int	ri )

Returns the actual TAC sample number, not including NaNs, samples with negative x, duplicate samples, or samples with zero weights (if data is weighted).

See also: dftEndtime, fittime_from_dft

Returns: Returns sample number.

Parameters

dft	Pointer to TAC data in DFT struct; must be sorted by increasing x.
ri	Region index [0..voiNr-1].

Definition at line 288 of file fittime.c.

  {
  int n, fi;
  double last_x=-1.0E+99;
 
  if(dft==NULL || dft->frameNr<1 || ri<0 || ri>=dft->voiNr) return 0;
  for(fi=n=0; fi<dft->frameNr; fi++) {
    if(isnan(dft->x[fi])) continue;
    if(isnan(dft->voi[ri].y[fi])) continue;
    if(dft->x[fi]<0.0) continue;
    if(dft->isweight!=0 && dft->w[fi]<=0.0) continue;
    if(dft->x[fi]==last_x) continue;
    n++; last_x=dft->x[fi];
  }
  return n;
}

◆ img_k1_using_ki()

int img_k1_using_ki	(	DFT *	input,
		IMG *	dyn_img,
		int	frame_nr,
		IMG *	ki_img,
		IMG *	k1_img,
		IMG *	k2k3_img,
		char *	status,
		int	verbose )

Computing pixel-by-pixel the K1 for irreversible PET tracers using previously determined Ki (K1*k3/(k2+k3) and bilinear regression.

Returns: Returns 0 if succesful, and >1 in case of an error.

Parameters

input	Pointer to the TAC data to be used as model input. Sample times in minutes. Curve is interpolated to PET frame times, if necessary
dyn_img	Pointer to dynamic PET image data. Image and input data must be in the same calibration units
frame_nr	Nr of frames that will be included in the fit [3-frame_nr]
ki_img	Pointer to previously calculated Ki image
k1_img	Pointer to initiated IMG structure where K1 values will be placed
k2k3_img	Pointer to initiated IMG structure where (k2+k3) values will be placed; enter NULL, if not needed
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 16 of file img_k1.c.

  {
  int zi, yi, xi, fi, ret=0;
  int nnls_n=2, nnls_m, m;
  DFT tac;
  clock_t fitStart, fitFinish;
 
 
  if(verbose>0) {
    printf("%s(input, dyn_img, %d, ki_img, k1_img, ", __func__, frame_nr);
    if(k2k3_img==NULL) printf("NULL, "); else printf("k2k3_img, ");
    if(status==NULL) printf("NULL)\n"); else printf("status)\n");
  }
  /* Initial check for the arguments */
  if(status!=NULL) sprintf(status, "invalid data");
  if(dyn_img->status!=IMG_STATUS_OCCUPIED || dyn_img->dimt<1) return(1);
  if(ki_img->status!=IMG_STATUS_OCCUPIED || ki_img->dimt<1) return(1);
  if(input==NULL || input->frameNr<1) return(1);
  if(frame_nr>dyn_img->dimt) return(1);
  if(k1_img==NULL) return(1);
  if(frame_nr<3) {
    if(status!=NULL) sprintf(status, "invalid fit time"); 
    return(1);
  }
  
  /* Check that input contains samples until at least 80% of line fit range */
  if(verbose>1) {
    printf("input_last_sample_time := %g\n", input->x[input->frameNr-1]);
    printf("k1_start_time := %g\n", dyn_img->start[0]/60.0);
    printf("k1_end_time := %g\n", dyn_img->end[frame_nr-1]/60.0);
  }
  if(input->x[input->frameNr-1] <
      (0.2*dyn_img->mid[0]+0.8*dyn_img->mid[frame_nr-1])/60.0)
  {
    sprintf(status, "too few input samples"); return(2);
  }
 
  fitStart=clock();
 
  /* Allocate memory for interpolated and integrated input and 
     tissue pixel TACs */
  dftInit(&tac); if(dftSetmem(&tac, frame_nr, 2)!=0) {
    sprintf(status, "out of memory"); return(3);
  }
  strcpy(tac.voi[0].voiname, "input"); strcpy(tac.voi[1].voiname, "tissue");
  tac.voiNr=2; tac.frameNr=frame_nr;
  for(fi=0; fi<tac.frameNr; fi++) {
    tac.x1[fi]=dyn_img->start[fi]/60.;
    tac.x2[fi]=dyn_img->end[fi]/60.;
    tac.x[fi]=dyn_img->mid[fi]/60.;
  }
 
  /* If input sample times are much different from PET frame times,
     as they usually are if arterial plasma is used as input, then use 
     interpolate4pet(), otherwise with image-derived input use petintegral()
  */
  ret=0;
  if(input->frameNr<dyn_img->dimt) ret=1;
  for(fi=0; fi<tac.frameNr; fi++) {
    if(verbose>8) printf("  %g  vs   %g\n", input->x1[fi], tac.x1[fi]);
    if(input->x1[fi]>tac.x1[fi]+0.034 || input->x1[fi]<tac.x1[fi]-0.034) {
      ret++; continue;}
    if(input->x2[fi]>tac.x2[fi]+0.034 || input->x2[fi]<tac.x2[fi]-0.034) ret++;
  }
  if(ret>0) {
    if(verbose>1) printf("using interpolate4pet() for input curve\n");
    ret=interpolate4pet(
      input->x, input->voi[0].y, input->frameNr,
      tac.x1, tac.x2, tac.voi[0].y, tac.voi[0].y2, tac.voi[0].y3, tac.frameNr
    );
  } else {
    if(verbose>1) printf("copying input curve and using petintegral()\n");
    /* because times are the same, the values can be copied directly */
    for(fi=0; fi<tac.frameNr; fi++) tac.voi[0].y[fi]=input->voi[0].y[fi];
    ret=petintegral(tac.x1, tac.x2, tac.voi[0].y, tac.frameNr,
      tac.voi[0].y2, tac.voi[0].y3);
  }
  if(ret) {
    dftEmpty(&tac); sprintf(status, "cannot interpolate input data");
    if(verbose>0) printf("  ret := %d\n", ret);
    return(4);
  }
  if(verbose>3) dftPrint(&tac);
 
  /*
   *  Allocate result images and fill the header info
   */
  /* K1 image */
  imgEmpty(k1_img);
  ret=imgAllocate(k1_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1);
  if(ret) {
    sprintf(status, "cannot allocate memory for K1 image");
    imgEmpty(k1_img); dftEmpty(&tac); return(5);
  }
  ret=imgCopyhdr(dyn_img, k1_img);
  if(ret) {
    sprintf(status, "cannot copy header info for result image");
    imgEmpty(k1_img); dftEmpty(&tac); return(5);
  }
  k1_img->unit=CUNIT_ML_PER_ML_PER_MIN;
  k1_img->decayCorrection=IMG_DC_NONCORRECTED; k1_img->isWeight=0;
  k1_img->start[0]=dyn_img->start[0]; k1_img->end[0]=dyn_img->end[frame_nr-1];
  if(verbose>9) imgInfo(k1_img);
 
  /* (k2+k3) image */
  if(k2k3_img!=NULL) {
    imgEmpty(k2k3_img);
    ret=imgAllocate(k2k3_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1);
    if(ret) {
      sprintf(status, "cannot allocate memory for (k2+k3) image");
      imgEmpty(k2k3_img); imgEmpty(k1_img); dftEmpty(&tac); return(6);
    }
    ret=imgCopyhdr(dyn_img, k2k3_img);
    if(ret) {
      sprintf(status, "cannot copy header info for result image");
      imgEmpty(k2k3_img); imgEmpty(k1_img); dftEmpty(&tac); return(6);
    }
    k2k3_img->unit=CUNIT_PER_MIN;
    k2k3_img->decayCorrection=IMG_DC_NONCORRECTED; k2k3_img->isWeight=0;
    k2k3_img->start[0]=dyn_img->start[0];
    k2k3_img->end[0]=dyn_img->end[frame_nr-1];
    if(verbose>9) imgInfo(k2k3_img);
  }
 
  /*
   *  Allocate memory required by NNLS; C99 !!!
   */
  if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
  nnls_m=frame_nr;
  double *nnls_a[nnls_n], nnls_mat[nnls_n][nnls_m], nnls_b[nnls_m],
          nnls_zz[nnls_m];
  double nnls_x[nnls_n], nnls_wp[nnls_n], nnls_rnorm;
  int nnls_index[nnls_n];
  nnls_a[0]=nnls_mat[0];
  nnls_a[1]=nnls_mat[1];
 
  /*
   *  Compute pixel-by-pixel
   */
  if(verbose>1) printf("computing K1 pixel-by-pixel\n");
  for(zi=0; zi<dyn_img->dimz; zi++) {
    for(yi=0; yi<dyn_img->dimy; yi++) {
      for(xi=0; xi<dyn_img->dimx; xi++) {
        /* Initiate pixel output values */
        k1_img->m[zi][yi][xi][0]=0.0;
        if(k2k3_img!=NULL) k2k3_img->m[zi][yi][xi][0]=0.0;
 
        /* Copy and integrate pixel curve */
        for(m=0; m<nnls_m; m++) tac.voi[1].y[m]=dyn_img->m[zi][yi][xi][m];
        ret=petintegral(tac.x1, tac.x2, tac.voi[1].y, tac.frameNr,
                        tac.voi[1].y2, NULL);
        if(ret) continue;
        /* if AUC at the end is <= 0, then do nothing */
        if(tac.voi[1].y2[nnls_m-1]<=0.0) continue;       
 
        /* Fill the NNLS data matrix */
        for(m=0; m<nnls_m; m++) {
          nnls_a[0][m]=tac.voi[0].y2[m];
          nnls_a[1][m]=ki_img->m[zi][yi][xi][0]*tac.voi[0].y3[m]-tac.voi[1].y2[m];
          nnls_b[m]=tac.voi[1].y[m];
        }
        
        /* NNLS */
        ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm,
                 nnls_wp, nnls_zz, nnls_index);
        if(ret>1) { /* no solution is possible */
          continue;
        } else if(ret==1) { /* max iteration count exceeded */ }
        k1_img->m[zi][yi][xi][0]=nnls_x[0];
        if(k2k3_img!=NULL) k2k3_img->m[zi][yi][xi][0]=nnls_x[1];
        
      } /* next column */
    } /* next row */
  } /* next plane */
  dftEmpty(&tac);
 
  fitFinish=clock(); //printf("CLOCKS_PER_SEC=%ld\n", CLOCKS_PER_SEC);
  //printf("%ld - %ld\n", fitFinish, fitStart);
  if(verbose>0 && fitFinish!=(clock_t)(-1)) printf("done in %g seconds.\n",
    (double)(fitFinish - fitStart) / (double)CLOCKS_PER_SEC );
  
  return(0);
}

◆ img_logan()

int img_logan	(	DFT *	input,
		IMG *	dyn_img,
		int	start,
		int	end,
		linefit_range	fit_range,
		float	thrs,
		double	k2,
		IMG *	vt_img,
		IMG *	ic_img,
		IMG *	nr_img,
		char *	status,
		int	verbose )

Computing pixel-by-pixel the graphical analysis for reversible PET tracers (Logan plot).

See also: img_patlak, llsqperp, logan_data, mtga_best_perp

Returns: Returns 0 if successful, and >0 in case of an error.

Parameters

input	Pointer to the TAC data to be used as model input. Sample times in minutes. Curve is interpolated to PET frame times, if necessary.
dyn_img	Pointer to dynamic PET image data. Image and input data must be in the same calibration units.
start	The range of frames where line is fitted, given as the frame start here and next the end index, i.e. [0..frame_nr-1].
end	The range of frames where line is fitted, given as the frame start above and here the end index, i.e. [0..frame_nr-1].
fit_range	Use the whole range or based on data leave out points from the beginning; PRESET or EXCLUDE_BEGIN.
thrs	Threshold as fraction of input AUC.
k2	Reference region k2; set to <=0 if not needed.
vt_img	Pointer to initiated IMG structure where Vt (or DVR) values will be placed.
ic_img	Pointer to initiated IMG structure where plot y axis intercept values times -1 will be placed; enter NULL, if not needed.
nr_img	Pointer to initiated IMG structure where the number of plot data points actually used in the fit is written; enter NULL, when not needed.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 228 of file img_mtga.c.

  {
  int zi, yi, xi, fi, nr, ret=0;
  DFT tac;
  double *plotData, *xaxis, *yaxis, slope, ic, f;
  double aucrat;
 
 
  if(verbose>0) {
    printf("img_logan(%s, dyn_img, %d, %d, %g, %g, ki_img, ", __func__, start, end, thrs, k2);
    if(ic_img==NULL) printf("NULL, "); else printf("ic_img, ");
    if(nr_img==NULL) printf("NULL, "); else printf("nr_img, ");
    if(status==NULL) printf("NULL)\n"); else printf("status)\n");
  }
  /* Initial check for the arguments */
  if(status!=NULL) sprintf(status, "invalid data");
  if(dyn_img->status!=IMG_STATUS_OCCUPIED || dyn_img->dimt<1) return(1);
  if(input==NULL || input->frameNr<1) return(2);
  nr=1+end-start; if(nr<2) return(3);
  if(end>dyn_img->dimt-1 || start<0) return(4);
  if(vt_img==NULL) return(5);
  /* Convert input time units to min */
  if(input->timeunit==TUNIT_SEC) dftTimeunitConversion(input, TUNIT_MIN);
 
  /* Check that input contains samples until at least 80% of line fit range */
  if(verbose>1) {
    printf("input_last_sample_time := %g\n", input->x[input->frameNr-1]);
    printf("logan_start_time := %g\n", dyn_img->start[start]/60.0);
    printf("logan_end_time := %g\n", dyn_img->end[end]/60.0);
  }
  if(input->x[input->frameNr-1] < (0.2*dyn_img->mid[start]+0.8*dyn_img->mid[end])/60.0) {
    sprintf(status, "too few input samples"); return(6);
  }
 
  /* Allocate memory for interpolated and integrated input and 
     tissue pixel TACs */
  dftInit(&tac); if(dftSetmem(&tac, nr, 2)!=0) {sprintf(status, "out of memory"); return(11);}
  strcpy(tac.voi[0].voiname, "input"); strcpy(tac.voi[1].voiname, "tissue");
  tac.voiNr=2; tac.frameNr=nr;
  for(int fi=0; fi<tac.frameNr; fi++) {
    tac.x1[fi]=dyn_img->start[start+fi]/60.;
    tac.x2[fi]=dyn_img->end[start+fi]/60.;
    tac.x[fi]=dyn_img->mid[start+fi]/60.;
  }
 
  /* If input sample times are much different from PET frame times,
     as they usually are if arterial plasma is used as input, then use 
     interpolate4pet(), otherwise with image-derived input use petintegral()
  */
  int equal_times;
  equal_times=check_times_dft_vs_img(dyn_img, input, verbose-1);
  if(equal_times==1) {
    if(verbose>1) printf("copying input curve and using petintegral()\n");
    /* Copy frame start and end times from image to input data */
    ret=copy_times_from_img_to_dft(dyn_img, input, verbose-1);
    if(ret) {
      sprintf(status, "cannot set input sample times");
      dftEmpty(&tac); return(12);
    }
    /* integral must be calculated from the zero time */
    ret=petintegral(input->x1, input->x2, input->voi[0].y, input->frameNr,
      input->voi[0].y2, input->voi[0].y3);
    /* because times are the same, the values can be copied directly */
    if(ret==0) for(fi=0; fi<tac.frameNr; fi++) {
      tac.voi[0].y[fi]=input->voi[0].y[start+fi];
      tac.voi[0].y2[fi]=input->voi[0].y2[start+fi];
      tac.voi[0].y3[fi]=input->voi[0].y3[start+fi];
    }
  } else {
    if(verbose>1) printf("using interpolate4pet() for input curve\n");
    ret=interpolate4pet(
      input->x, input->voi[0].y, input->frameNr,
      tac.x1, tac.x2, tac.voi[0].y, tac.voi[0].y2, tac.voi[0].y3, tac.frameNr
    );
  }
  if(ret) {
    dftEmpty(&tac);
    sprintf(status, "cannot interpolate input data"); return(12);
  }
  if(verbose>3) dftPrint(&tac);
 
  /*
   *  Allocate result images and fill the header info
   */
  /* Vt image */
  imgEmpty(vt_img);
  ret=imgAllocateWithHeader(vt_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1, dyn_img);
  if(ret) {
    sprintf(status, "cannot setup memory for Vt image");
    imgEmpty(vt_img); dftEmpty(&tac); return(21);
  }
  vt_img->unit=CUNIT_ML_PER_ML;
  vt_img->decayCorrection=IMG_DC_NONCORRECTED; vt_img->isWeight=0;
  vt_img->start[0]=dyn_img->start[start]; vt_img->end[0]=dyn_img->end[end];
  if(verbose>9) imgInfo(vt_img);
  /* Ic image */
  if(ic_img!=NULL) {
    imgEmpty(ic_img);
    ret=imgAllocateWithHeader(ic_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1, dyn_img);
    if(ret) {
      sprintf(status, "cannot setup memory for Ic image");
      imgEmpty(ic_img); imgEmpty(vt_img); dftEmpty(&tac); return(23);
    }
    ic_img->unit=CUNIT_UNITLESS; /* mL/mL */
    ic_img->decayCorrection=IMG_DC_NONCORRECTED; ic_img->isWeight=0;
    ic_img->start[0]=dyn_img->start[start]; ic_img->end[0]=dyn_img->end[end];
    if(verbose>9) imgInfo(ic_img);
  }
  /* Nr image */
  if(nr_img!=NULL) {
    imgEmpty(nr_img);
    ret=imgAllocateWithHeader(nr_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1, dyn_img);
    if(ret) {
      sprintf(status, "cannot setup memory for nr image");
      imgEmpty(nr_img); imgEmpty(ic_img); imgEmpty(vt_img); dftEmpty(&tac);
      return(25);
    }
    nr_img->unit=CUNIT_UNITLESS;
    nr_img->decayCorrection=IMG_DC_NONCORRECTED; nr_img->isWeight=0;
    nr_img->start[0]=dyn_img->start[start]; nr_img->end[0]=dyn_img->end[end];
    if(verbose>9) imgInfo(ic_img);
  }
 
  /* Allocate memory for graphical analysis plot data */
  plotData=malloc(2*nr*sizeof(double));
  if(plotData==NULL) {
    sprintf(status, "cannot allocate memory for plots");
    imgEmpty(ic_img); imgEmpty(vt_img); dftEmpty(&tac); 
    return(25);
  }
  xaxis=plotData; yaxis=plotData+nr;
  float pxlauc[dyn_img->dimt];
 
  /* Calculate threshold */
  thrs*=tac.voi[0].y2[tac.frameNr-1];
  if(verbose>1) printf("  threshold-AUC := %g\n", thrs);
 
  /*
   *  Compute pixel-by-pixel
   */
  if(verbose>1) printf("computing MTGA pixel-by-pixel\n");
  int best_nr;
  for(zi=0; zi<dyn_img->dimz; zi++) {
    for(yi=0; yi<dyn_img->dimy; yi++) {
      for(xi=0; xi<dyn_img->dimx; xi++) {
        /* Initiate pixel output values */
        vt_img->m[zi][yi][xi][0]=0.0;
        if(ic_img!=NULL) ic_img->m[zi][yi][xi][0]=0.0;
        if(nr_img!=NULL) nr_img->m[zi][yi][xi][0]=0.0;
        /* Copy TTAC values */
        for(fi=0; fi<tac.frameNr; fi++) tac.voi[1].y[fi]=dyn_img->m[zi][yi][xi][start+fi];
        /* Compute TTAC AUC(0-t) and check for threshold */
        ret=fpetintegral(dyn_img->start, dyn_img->end, dyn_img->m[zi][yi][xi], 
                         dyn_img->dimt, pxlauc, NULL);
        if(ret) continue;
        if((pxlauc[dyn_img->dimt-1]/60.0) < thrs) continue;
        for(fi=0; fi<tac.frameNr; fi++)
          tac.voi[1].y2[fi]=pxlauc[start+fi]/60.0; // conc*sec -> conc*min
        /* Calculate Logan plot data */
        int pn=logan_data(nr, tac.voi[0].y, tac.voi[0].y2, tac.voi[1].y, tac.voi[1].y2, k2, xaxis, yaxis);
        /* Line fit */
        if(fit_range==PRESET || pn<MTGA_BEST_MIN_NR) {
          ret=llsqperp(xaxis, yaxis, pn, &slope, &ic, &f);
          best_nr=pn;
        } else {
          ret=mtga_best_perp(xaxis, yaxis, pn, &slope, &ic, &f, &best_nr);
        }
        if(ret!=0) continue; // line fit failed
        /* Use 10xAUCratio as upper limit to prevent image where only 
           noise-induced hot spots can be seen */
        aucrat=tac.voi[1].y2[tac.frameNr-1]/tac.voi[0].y2[tac.frameNr-1];
        if(slope>10.0*aucrat) {
          if(verbose>50) printf("%g > 10 x %g\n", slope, aucrat);
          slope=10.0*aucrat;
        }
        /* copy result to pixels */
        vt_img->m[zi][yi][xi][0]=slope;
        if(ic_img!=NULL) ic_img->m[zi][yi][xi][0]=-ic;
        if(nr_img!=NULL) nr_img->m[zi][yi][xi][0]=best_nr;
      } /* next column */
    } /* next row */
  } /* next plane */
 
  free(plotData); dftEmpty(&tac);
  return(0);
}

◆ img_patlak()

int img_patlak	(	DFT *	input,
		IMG *	dyn_img,
		int	start,
		int	end,
		linefit_range	fit_range,
		float	thrs,
		IMG *	ki_img,
		IMG *	ic_img,
		IMG *	nr_img,
		char *	status,
		int	verbose )

Computing pixel-by-pixel the graphical analysis for irreversible PET tracers (Gjedde-Patlak plot).

See also: img_logan, llsqperp, patlak_data, mtga_best_perp

Returns: Returns 0 if successful, and >0 in case of an error.

Parameters

input	Pointer to the TAC data to be used as model input. Sample times in minutes. Curve is interpolated to PET frame times, if necessary.
dyn_img	Pointer to dynamic PET image data. Image and input data must be in the same calibration units.
start	The range of frames where line is fitted, given as the frame start here and next the end index, i.e. [0..frame_nr-1].
end	The range of frames where line is fitted, given as the frame start above and here the end index, i.e. [0..frame_nr-1].
fit_range	Use the whole range or based on data leave out points from the beginning; PRESET or EXCLUDE_BEGIN.
thrs	Threshold as fraction of input AUC.
ki_img	Pointer to initiated IMG structure where Ki values will be placed.
ic_img	Pointer to initiated IMG structure where plot y axis intercept values will be placed; enter NULL, if not needed.
nr_img	Pointer to initiated IMG structure where the number of plot data points actually used in the fit is written; enter NULL, when not needed.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 17 of file img_mtga.c.

  {
  int zi, yi, xi, fi, nr, ret=0;
  DFT tac;
  double *plotData, *xaxis, *yaxis, slope, ic, f;
 
 
  if(verbose>0) {
    printf("%s(input, dyn_img, %d, %d, range, %g, ki_img, ", __func__, start, end, thrs);
    if(ic_img==NULL) printf("NULL, "); else printf("ic_img, ");
    if(nr_img==NULL) printf("NULL, "); else printf("nr_img, ");
    if(status==NULL) printf("NULL)\n"); else printf("status)\n");
  }
  /* Initial check for the arguments */
  if(status!=NULL) sprintf(status, "invalid data");
  if(dyn_img->status!=IMG_STATUS_OCCUPIED || dyn_img->dimt<1) return(1);
  if(input==NULL || input->frameNr<1) return(2);
  nr=1+end-start; if(nr<2) return(3);
  if(end>dyn_img->dimt-1 || start<0) return(4);
  if(ki_img==NULL) return(5);
  /* Convert input time units to min */
  if(input->timeunit==TUNIT_SEC) dftTimeunitConversion(input, TUNIT_MIN);
 
  /* Check that input contains samples until at least 80% of line fit range */
  if(verbose>1) {
    printf("input_last_sample_time := %g\n", input->x[input->frameNr-1]);
    printf("patlak_start_time := %g\n", dyn_img->start[start]/60.0);
    printf("patlak_end_time := %g\n", dyn_img->end[end]/60.0);
  }
  if(input->x[input->frameNr-1] < (0.2*dyn_img->mid[start]+0.8*dyn_img->mid[end])/60.0) {
    sprintf(status, "too few input samples"); return(6);
  }
 
  /* Allocate memory for interpolated and integrated input and 
     tissue pixel TACs */
  dftInit(&tac); if(dftSetmem(&tac, nr, 2)!=0) {sprintf(status, "out of memory"); return(11);}
  strcpy(tac.voi[0].voiname, "input"); strcpy(tac.voi[1].voiname, "tissue");
  tac.voiNr=2; tac.frameNr=nr;
  for(fi=0; fi<tac.frameNr; fi++) {
    tac.x1[fi]=dyn_img->start[start+fi]/60.;
    tac.x2[fi]=dyn_img->end[start+fi]/60.;
    tac.x[fi]=dyn_img->mid[start+fi]/60.;
  }
 
  /* If input sample times are much different from PET frame times,
     as they usually are if arterial plasma is used as input, then use 
     interpolate4pet(), otherwise with image-derived input use petintegral()
  */
  int equal_times;
  equal_times=check_times_dft_vs_img(dyn_img, input, verbose-1);
  if(equal_times==1) {
    if(verbose>1) printf("copying input curve and using petintegral()\n");
    /* Copy frame start and end times from image to input data */
    ret=copy_times_from_img_to_dft(dyn_img, input, verbose-1);
    if(ret) {
      sprintf(status, "cannot set input sample times");
      dftEmpty(&tac); return(12);
    }
    /* integral must be calculated from the zero time */
    ret=petintegral(input->x1, input->x2, input->voi[0].y, input->frameNr,
      input->voi[0].y2, input->voi[0].y3);
    /* because times are the same, the values can be copied directly */
    if(ret==0) for(fi=0; fi<tac.frameNr; fi++) {
      tac.voi[0].y[fi]=input->voi[0].y[start+fi];
      tac.voi[0].y2[fi]=input->voi[0].y2[start+fi];
      tac.voi[0].y3[fi]=input->voi[0].y3[start+fi];
    }
  } else {
    if(verbose>1) printf("using interpolate4pet() for input curve\n");
    ret=interpolate4pet(
      input->x, input->voi[0].y, input->frameNr,
      tac.x1, tac.x2, tac.voi[0].y, tac.voi[0].y2, tac.voi[0].y3, tac.frameNr
    );
  }
  if(ret) {
    dftEmpty(&tac);
    sprintf(status, "cannot interpolate input data"); return(12);
  }
  if(verbose>3) dftPrint(&tac);
 
  /*
   *  Allocate result images and fill the header info
   */
  /* Ki image */
  imgEmpty(ki_img);
  ret=imgAllocateWithHeader(ki_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1, dyn_img);
  if(ret) {
    sprintf(status, "cannot setup memory for Ki image");
    imgEmpty(ki_img); dftEmpty(&tac); return(21);
  }
  ki_img->unit=CUNIT_ML_PER_ML_PER_MIN; /* mL/(mL*min) */
  ki_img->decayCorrection=IMG_DC_NONCORRECTED; ki_img->isWeight=0;
  ki_img->start[0]=dyn_img->start[start]; ki_img->end[0]=dyn_img->end[end];
  if(verbose>9) imgInfo(ki_img);
  /* Ic image */
  if(ic_img!=NULL) {
    imgEmpty(ic_img);
    ret=imgAllocateWithHeader(ic_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1, dyn_img);
    if(ret) {
      sprintf(status, "cannot setup memory for Ic image");
      imgEmpty(ic_img); imgEmpty(ki_img); dftEmpty(&tac); return(23);
    }
    ic_img->unit=CUNIT_ML_PER_ML; /* mL/mL */
    ic_img->decayCorrection=IMG_DC_NONCORRECTED; ic_img->isWeight=0;
    ic_img->start[0]=dyn_img->start[start]; ic_img->end[0]=dyn_img->end[end];
    if(verbose>9) imgInfo(ic_img);
  }
  /* Nr image */
  if(nr_img!=NULL) {
    imgEmpty(nr_img);
    ret=imgAllocateWithHeader(nr_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1, dyn_img);
    if(ret) {
      sprintf(status, "cannot setup memory for nr image");
      imgEmpty(nr_img); imgEmpty(ic_img); imgEmpty(ki_img); dftEmpty(&tac);
      return(25);
    }
    nr_img->unit=CUNIT_UNITLESS;
    nr_img->decayCorrection=IMG_DC_NONCORRECTED; nr_img->isWeight=0;
    nr_img->start[0]=dyn_img->start[start]; nr_img->end[0]=dyn_img->end[end];
    if(verbose>9) imgInfo(ic_img);
  }
 
  /* Allocate memory for graphical analysis plot data */
  plotData=malloc(2*nr*sizeof(double));
  if(plotData==NULL) {
    sprintf(status, "cannot allocate memory for plots");
    imgEmpty(ic_img); imgEmpty(ki_img); dftEmpty(&tac); return(25);
  }
  xaxis=plotData; yaxis=plotData+nr;
 
  /* Calculate threshold */
  thrs*=tac.voi[0].y2[tac.frameNr-1];
  if(verbose>1) printf("  threshold-AUC := %g\n", thrs);
 
  /*
   *  Compute pixel-by-pixel
   */
  if(verbose>1) printf("computing MTGA pixel-by-pixel\n");
  float pxlauc[dyn_img->dimt];
  int best_nr;
  for(zi=0; zi<dyn_img->dimz; zi++) {
    for(yi=0; yi<dyn_img->dimy; yi++) {
      for(xi=0; xi<dyn_img->dimx; xi++) {
        /* Initiate pixel output values */
        ki_img->m[zi][yi][xi][0]=0.0;
        if(ic_img!=NULL) ic_img->m[zi][yi][xi][0]=0.0;
        if(nr_img!=NULL) nr_img->m[zi][yi][xi][0]=0;
        /* Check for threshold */
        ret=fpetintegral(dyn_img->start, dyn_img->end, dyn_img->m[zi][yi][xi], 
                         dyn_img->dimt, pxlauc, NULL);
        if(ret) continue;
        if((pxlauc[dyn_img->dimt-1]/60.0) < thrs) continue;
        /* Calculate Patlak plot data */
        for(fi=0; fi<tac.frameNr; fi++) tac.voi[1].y[fi]=dyn_img->m[zi][yi][xi][start+fi];
        int pn=patlak_data(nr, tac.voi[0].y, tac.voi[0].y2, tac.voi[1].y, xaxis, yaxis);
        /* Line fit */
        if(fit_range==PRESET || pn<MTGA_BEST_MIN_NR) {
          ret=llsqperp(xaxis, yaxis, pn, &slope, &ic, &f);
          best_nr=pn;
        } else {
          ret=mtga_best_perp(xaxis, yaxis, pn, &slope, &ic, &f, &best_nr);
        }
        if(ret==0) {
          ki_img->m[zi][yi][xi][0]=slope;
          if(ic_img!=NULL) ic_img->m[zi][yi][xi][0]=ic;
          if(nr_img!=NULL) nr_img->m[zi][yi][xi][0]=best_nr;
        }
      } /* next column */
    } /* next row */
  } /* next plane */
 
  free(plotData); dftEmpty(&tac);
  return(0);
}

◆ imgEndtime()

double imgEndtime ( IMG * img )

Get IMG end time. Frame times are assumed to be sorted to increasing order.

See also: dftEndtime, fittime_from_img

Returns: Returns the last frame end time. By default, times are in sec.

Parameters

img	Pointer to IMG structure.

Definition at line 339 of file fittime.c.

  {
  if(img==NULL || img->dimt<1) return(0.0);
  for(int fi=img->dimt-1; fi>=0; fi--)
    if(!isnan(img->end[fi])) return(img->end[fi]);
  return(0.0);
}

Referenced by imgReadModelingData().

◆ imgReadModelingData()

int imgReadModelingData	(	char *	petfile,
		char *	siffile,
		char *	inputfile1,
		char *	inputfile2,
		char *	inputfile3,
		double *	fitdur,
		int *	fitframeNr,
		IMG *	img,
		DFT *	inp,
		DFT *	iinp,
		int	verifypeak,
		FILE *	loginfo,
		int	verbose,
		char *	status )

Read PET image and input TAC for modelling.

Input calibration units are converted to units of tissue data. Input data is optionally interpolated to tissue times. If input data extends much further than fit duration, the last part is removed to save computation time in simulations. Input data is optionally verified not to end too early, and not to start too late.

Returns: Returns 0 when successful, otherwise a non-zero value.

See also: imgRead, dftRead, dftReadinput, dftReadModelingData, dftEndtime, imgEndtime

Parameters

petfile	PET dynamic image filename; one time frame is sufficient here.
siffile	Optional SIF. SIF is automatically read with NIfTI and Analyze formats, if found in the same folder and files are consistently named. SIF filename needs to be specified here only if that is not the case, or if other image formats do not contain information on frame times or tracer isotope, or if weights from SIF are needed. Enter NULL if not needed.
inputfile1	1st input data filename.
inputfile2	2nd input data filename (or NULL if only not needed).
inputfile3	3rd input data filename (or NULL if only not needed).
fitdur	Fit duration (in minutes); shortened if longer than PET data; enter <=0 to use all data in image; input is verified to not be much shorter than fitdur.
fitframeNr	Nr of time frames (samples) in PET data that are inside fitdur will be written here; pointer MUST be provided.
img	Pointer to initiated IMG structure into which PET data will be written.
inp	Pointer to initiated DFT into which input data (plasma, blood, and/or reference tissue TAC) will be written without interpolation; enter NULL if not needed. Times will be in seconds, like in image.
iinp	Pointer to initiated DFT into which input data (plasma, blood, and/or reference tissue TAC) will be written, interpolated to PET frames; enter NULL if not needed. Times will be in seconds, like in image.
verifypeak	Set to <>0 to verify that input TAC does not start too late and has decent peak to provide reliable AUC0-T.
loginfo	Give file pointer (for example stdout) where log information is printed; NULL if not needed; warnings will be printed in stderr anyway.
verbose	Verbose level; if zero, then only warnings are printed into stderr.
status	Pointer to a string (allocated for at least 256 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 24 of file imginput.c.

  {
  int ret, fi, n, first, last, ii, input_nr=0;
  DFT tmpdft, dft;
  SIF sif;
  double f, starttime, endtime;
  FILE *wfp;
  char *fname, tmp[512];
 
 
  if(loginfo!=NULL && verbose>0) {
    fprintf(loginfo, "imgReadModelingData(\n");
    fprintf(loginfo, "  %s,\n", petfile);
    fprintf(loginfo, "  %s,\n", inputfile1);
    fprintf(loginfo, "  %s,\n", inputfile2);
    fprintf(loginfo, "  %s,\n", inputfile3);
    fprintf(loginfo, "  %g,\n", *fitdur);
    fprintf(loginfo, "  *fitframeNr, *tis, *inp, *iinp *loginfo, %d, *status\n", verbose);
    fprintf(loginfo, ")\n");
  }
  /* Check the function input */
  if(status!=NULL) sprintf(status, "program error");
  if(img==NULL || (inp==NULL && iinp==NULL)) return -1;
  if(petfile==NULL || strlen(petfile)<1) return -2;
  ret=0;
  if(inputfile1==NULL || strlen(inputfile1)<1) return -3; else input_nr++;
  if(inputfile2!=NULL && strlen(inputfile2)>0) input_nr++;
  if(inputfile3!=NULL && strlen(inputfile3)>0) {
    if(input_nr<2) return -4;
    input_nr++;
  }
  if(fitframeNr==NULL || fitdur==NULL) return -5;
  if(status!=NULL) sprintf(status, "arguments validated");
  /* Warnings will be printed into stderr */
  wfp=stderr;
 
  /* Check that input file(s) exist, because image is read first, and it might take a while */
  if(loginfo!=NULL && verbose>1) fprintf(loginfo, "checking access to input files\n");
  if(access(inputfile1, 0) == -1) {
    sprintf(tmp, "cannot read '%s'", inputfile1);
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    return(1);
  }
  for(ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=inputfile2; else fname=inputfile3;
    if(access(fname, 0) == -1) {
      sprintf(tmp, "cannot read '%s'", fname);
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      return(1);
    }
  }
 
  /* Delete any previous data and initiate temp data */
  dftEmpty(inp); dftEmpty(iinp); imgEmpty(img);
 
 
  /*
   *  Read PET image
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading image %s\n", petfile);
  imgInit(img);
  ret=imgRead(petfile, img); fflush(stderr); fflush(stdout);
  if(ret) {
    sprintf(tmp, "cannot read '%s': %s", petfile, img->statmsg);
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    if(loginfo!=NULL) fprintf(loginfo, "imgRead(%s, img) := %d\n", petfile, ret);
    return(2);
  }
  /* Check if PET data is raw or image */
  if(img->type!=IMG_TYPE_IMAGE) {
    sprintf(tmp, "%s is not an image.\n", petfile);
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); return(3);
  }
  if(loginfo!=NULL && verbose>2)
    fprintf(loginfo, "image contains %d frames and %d planes.\n", img->dimt, img->dimz);
 
  /* Read SIF, if specified */
  if(siffile!=NULL && siffile[0]) {
    if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading SIF %s\n", siffile);
    sifInit(&sif);
    ret=sifRead(siffile, &sif);
    if(ret) {
      sprintf(tmp, "cannot read '%s': %s", siffile, siferrmsg);
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); return(4);
    }
    /* Get isotope from SIF, if available */
    {
      double hl;
      hl=hlFromIsotope(sif.isotope_name);
      if(hl>0.0) {
        img->isotopeHalflife=60.0*hl;
        if(loginfo!=NULL && verbose>0) fprintf(loginfo, "isotope code read from %s\n", siffile);
      }
    }
    /* Get frame times from SIF if not available in image */
    if(!imgExistentTimes(img) && img->dimt<=sif.frameNr) {
      for(fi=0; fi<img->dimt; fi++) {
        img->start[fi]=sif.x1[fi]; img->end[fi]=sif.x2[fi];
        img->mid[fi]=0.5*(sif.x1[fi]+sif.x2[fi]);
      }
      if(loginfo!=NULL && verbose>0) fprintf(loginfo, "image frame times read from %s\n", siffile);
    }
    sifEmpty(&sif);
  }
 
  /* Check that image frame times are available */
  if(loginfo!=NULL && verbose>1) fprintf(loginfo, "checking image contents\n");
  if(!imgExistentTimes(img)) {
    strcpy(tmp, "image frame times not available");
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); return(5);
  }
  /* Make sure that there is no overlap in image frames */ 
  if(loginfo!=NULL && verbose>1) fprintf(loginfo, "checking frame overlap in %s\n", petfile);
  ret=imgDeleteFrameOverlap(img);
  if(ret) {
    strcpy(tmp, "image has overlapping frame times");
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); return(5);
  }
  /* Check fit end time with the PET image */
  if(!(*fitdur>0.0)) *fitdur=1.0E+99;
  else if(*fitdur<1.0E+10) *fitdur*=60.0; // fit time must be in sec
  *fitframeNr=fittime_from_img(img, fitdur, verbose-2);
  if(*fitframeNr<1 || *fitdur<=3.5) {
    strcpy(tmp, "image has no data in fit time range");
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); return(5);
  }
  *fitdur/=60.0; // back to minutes
  if(loginfo!=NULL && verbose>3) {
    fprintf(loginfo, "fittimeFinal := %g min\n", *fitdur);
    fprintf(loginfo, "fitframeNr := %d\n", *fitframeNr);
  }
 
 
  /*
   *  Read first input data
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading input data in %s\n", inputfile1);
  dftInit(&dft);
  if(dftRead(inputfile1, &dft)) {
    sprintf(tmp, "cannot read '%s': %s", inputfile1, dfterrmsg);
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); return(11);
  }
  /* Check TAC nr */
  if(dft.voiNr>1) {
    if(verbose>0) {
      strcpy(tmp, "only first TAC is used as input");
      if(loginfo!=NULL) fprintf(loginfo, "Warning: %s.\n", tmp);
      else fprintf(wfp, "Warning: %s.\n", tmp);
    }
    dft.voiNr=1;
  }
  /* check if file contains NAs (missing values) */
  if(dft_nr_of_NA(&dft)>0) {
    strcpy(tmp, "missing values in input file");
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s.\n", tmp);
    else fprintf(wfp, "Error: %s.\n", tmp);
    imgEmpty(img); dftEmpty(&dft); return(12);
  }
  /* Check the concentration units */
  ret=cunit_check_dft_vs_img(&dft, img, tmp, verbose-2);
  if(ret==0) {
    if(loginfo!=NULL && verbose>3) fprintf(loginfo, "%s\n", tmp);
  } else if(ret<0) {
    fprintf(wfp, "Warning: %s\n", tmp);
  } else {
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s.\n", tmp);
    else fprintf(wfp, "Error: %s.\n", tmp);
    imgEmpty(img); dftEmpty(&dft); return(13);
  }
  if(verbose>3)
    printf("input time range := %g - %g %s\n",
      dft.x[0], dft.x[dft.frameNr-1], petTunit(dft.timeunit));
  /* Set time unit to sec like in IMG */
  if(dft.timeunit==TUNIT_UNKNOWN) {
    if(dftEndtime(&dft)<0.2*imgEndtime(img)) dft.timeunit=TUNIT_MIN;
    else dft.timeunit=TUNIT_SEC;
    fprintf(wfp, "Warning: assuming that times are in %s in %s\n",
            petTunit(dft.timeunit), inputfile1);
  }
  /* Set time unit to sec like in IMG */
  ret=dftTimeunitConversion(&dft, TUNIT_SEC);
  if(verbose>1)
    printf("input time range := %g - %g %s\n",
           dft.x[0], dft.x[dft.frameNr-1], petTunit(dft.timeunit));
  
  /* Verify the peak if requested; only for the first input TAC, because
     for example metabolite TAC may not show bolus shape */
  if(verifypeak!=0) {
    if(loginfo!=NULL && verbose>1) fprintf(loginfo, "verifying input peak\n");
    ret=dftVerifyPeak(&dft, 0, verbose-5, tmp);
    if(ret>0) {
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s.\n", tmp);
      else fprintf(wfp, "Error: %s.\n", tmp);
      imgEmpty(img); dftEmpty(&dft); return(14);
    }
  }
  if(verbose>5)
    printf("input time range := %g - %g %s\n",
           dft.x[0], dft.x[dft.frameNr-1], petTunit(dft.timeunit));
 
 
  /*
   *  Read following input files, if required
   */
  dftInit(&tmpdft);
  for(ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=inputfile2; else fname=inputfile3;
 
    /* Make room for one more curve in input data */
    ret=dftAddmem(&dft, 1);
    if(ret) {
      strcpy(tmp, "cannot allocate more memory");
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); return(15);
    }
    /* Read input TAC */
    if(loginfo!=NULL && verbose>0) fprintf(loginfo, "reading input data in %s\n", fname);
    if(dftRead(fname, &tmpdft)) {
      sprintf(tmp, "cannot read '%s': %s", fname, dfterrmsg);
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); return(16);
    }
    /* Check TAC nr */
    if(tmpdft.voiNr>1) {
      if(verbose>0) {
        strcpy(tmp, "only first TAC is used as input");
        if(loginfo!=NULL) fprintf(loginfo, "Warning: %s.\n", tmp);
        else fprintf(wfp, "Warning: %s.\n", tmp);
      }
      dft.voiNr=1;
    }
    /* check if file contains NAs (missing values) */
    if(dft_nr_of_NA(&tmpdft)>0) {
      strcpy(tmp, "missing values in input file");
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); dftEmpty(&tmpdft); return(17);
    }
    /* Check and correct the sample time unit */
    if(tmpdft.timeunit==TUNIT_UNKNOWN) {
      if(dftEndtime(&tmpdft)<0.2*imgEndtime(img)) tmpdft.timeunit=TUNIT_MIN;
      else tmpdft.timeunit=TUNIT_SEC;
      fprintf(wfp, "Warning: assuming that times are in %s in %s\n",
              petTunit(tmpdft.timeunit), fname);
    }
    dftTimeunitConversion(&tmpdft, dft.timeunit);
    /* Check the concentration units */
    ret=cunit_check_dft_vs_img(&tmpdft, img, tmp, verbose-2);
    if(ret==0) {
      if(loginfo!=NULL && verbose>3) fprintf(loginfo, "%s\n", tmp);
    } else if(ret<0) {
      fprintf(wfp, "Warning: %s\n", tmp);
    } else {
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); dftEmpty(&tmpdft); return(18);
    }
    /* Copy to input data */
    if(loginfo!=NULL && verbose>1)
      fprintf(loginfo, "interpolating %d samples into %d samples.\n", tmpdft.frameNr, dft.frameNr);
    ret=dftInterpolateInto(&tmpdft, &dft, tmp, verbose);
    if(ret) {
      if(loginfo!=NULL) fprintf(loginfo, "dftInterpolateInto() := %d\n", ret);
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); dftEmpty(&tmpdft); return(19);
    }
    /* Remove the originally read input data */
    dftEmpty(&tmpdft);
 
  } // next input file
  if(verbose>10)
    printf("input time range := %g - %g %s\n",
           dft.x[0], dft.x[dft.frameNr-1], petTunit(dft.timeunit));
 
 
  /*
   *  Check and set input data range vs fit time length
   */
  if(loginfo!=NULL && verbose>0) fprintf(loginfo, "checking and setting input sample time range\n");
  /* Check that input data does not end _much_ before fitdur; 
     dftInterpolateForIMG() below tests this again. */
  starttime=0; endtime=*fitdur; // start and end times always in min
  n=fittime_from_dft(&dft, &starttime, &endtime, &first, &last, verbose-1);
  if(loginfo!=NULL && verbose>2) {
    fprintf(loginfo, "starttime := %g min\n", starttime);
    fprintf(loginfo, "endtime := %g min\n", endtime);
    fprintf(loginfo, "sample_nr := %d\n", n);
  }
  if(*fitdur>1.5*endtime && (*fitdur-endtime)>0.15) {
    strcpy(tmp, "input TAC is too short");
    if(status!=NULL) strcpy(status, tmp);
    else if(verbose>0) fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); dftEmpty(&dft); return(21);
  }
#if(0)
  /* Do NOT test this here! It would prevent using this function in
     reading steady-state studies etc */ 
  if(n<4) {
    strcpy(tmp, "too few input samples in specified fit duration");
    if(status!=NULL) strcpy(status, tmp);
    else if(loginfo!=NULL) fprintf(loginfo, "Error: %s\n", tmp);
    else fprintf(wfp, "Error: %s\n", tmp);
    imgEmpty(img); dftEmpty(&dft); return(22);
  }
#endif
  /* Cut off too many input samples to make calculation faster */
  if(verbose>10) {
    printf("input time range := %g - %g %s\n",
           dft.x[0], dft.x[dft.frameNr-1], petTunit(dft.timeunit));
    printf("fitdur := %g min\n", *fitdur);
  }
  f=*fitdur*60.0+60.0; // add 60 s to allow time delay correction
  if(dftEndtime(&dft)>f) {
    if(loginfo!=NULL && verbose>0) fprintf(loginfo, "Input TAC cutoff at %g sec\n", f); 
    for(fi=0; fi<dft.frameNr; fi++) if(dft.x[fi]>f) break;
    if(fi<dft.frameNr) fi++; 
    dft.frameNr=fi;
  }
  if(verbose>10) printf("input time range := %g - %g %s\n",
                        dft.x[0], dft.x[dft.frameNr-1], petTunit(dft.timeunit));
 
  /* Copy input data to given pointer, if required */
  if(inp!=NULL) {
    if(loginfo!=NULL && verbose>0) fprintf(loginfo, "copying input TAC\n");
    ret=dftdup(&dft, inp);
    if(ret!=0) {
      strcpy(tmp, "cannot copy TAC contents");
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s.\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); return(31);
    }
  }
    
  /* Interpolate input data to given pointer, if required */
  if(iinp!=NULL) {
    if(loginfo!=NULL && verbose>0) fprintf(loginfo, "interpolating input TAC\n");
    ret=dftInterpolateForIMG(&dft, img, *fitframeNr, iinp, NULL, NULL, verbose, tmp);
    if(ret!=0) {
      if(status!=NULL) strcpy(status, tmp);
      else if(loginfo!=NULL) fprintf(loginfo, "Error: %s.\n", tmp);
      else fprintf(wfp, "Error: %s\n", tmp);
      imgEmpty(img); dftEmpty(&dft); return(33);
    }
  }
 
  dftEmpty(&dft);
  if(status!=NULL) sprintf(status, "ok");
  return(0);
}

◆ imgSetWeights()

int imgSetWeights	(	IMG *	img,
		int	wmet,
		int	verbose )

Add weights based on average voxel values and frame lengths, just frame lengths, or set all weights to 1 (no weighting).

Existing weights are overwritten. Decay correction is not considered in calculation of weights.

Returns: Returns 0 when successful, otherwise <>0.

See also: dftWeightByFreq, sifWeight, sifWeightByFrames

Parameters

img	Pointer to IMG struct; image data is used to calculate relational frame weights, which are also written into IMG struct.
wmet	Weighting method: 0=based on mean voxel values times frame lengths; 1=based on frame lengths, 2=set weights to 1 (no weighting).
verbose	Verbose level; if zero, then only warnings are printed into stderr

Definition at line 85 of file weight_model.c.

  {
  if(verbose>0) printf("imgSetWeights(*img, %d, ...)\n", wmet);
  /* Check the function input */
  if(img==NULL) return(1);
  if(img->status!=IMG_STATUS_OCCUPIED) return(2);
  if(img->dimt<1 || img->dimx<1 || img->dimy<1 || img->dimz<1) return(3);
  if(wmet<0 || wmet>2) return(4);
  img->isWeight=0;
 
  /* If one time frame, then weight is 1 */
  if(img->dimt==1) {
    img->weight[0]=1.0; 
    if(wmet==0 || wmet==1) img->isWeight=1; 
    return(0);
  }
  
  int fi, ret;
  double f, fm;
 
  if(wmet==2) { /* If weights are to be set to 1 (removed) */
    for(fi=0; fi<img->dimt; fi++) img->weight[fi]=1.0;
    return(0);
  }
 
  if(wmet==0) { /* weight based on frame mean values and durations */
    SIF sif;
    sifInit(&sif);
    ret=sifAllocateWithIMG(&sif, img, 1, verbose); if(ret) return(100+ret);
    sifModerateTrues(&sif, 100.0);
    sifWeight(&sif, 0.0); if(verbose>2) sifPrint(&sif);
    for(fi=0; fi<img->dimt; fi++) img->weight[fi]=sif.weights[fi];
    sifEmpty(&sif);
  } else if(wmet==1) { /* weight based on frame durations */
    for(fi=0, fm=0.0; fi<img->dimt; fi++) {
      f=img->end[fi]-img->start[fi]; fm+=f;
      img->weight[fi]=f;
    }
    /* scale weights so that sum of weights equals sample number */
    if(fm>0.0) f=(double)img->dimt/fm; else f=1.0;
    for(fi=0; fi<img->dimt; fi++) img->weight[fi]*=f;
  }
  img->isWeight=1;
  return(0);
}

◆ imgTimeIntegral()

int imgTimeIntegral	(	IMG *	img,
		float	t1,
		float	t2,
		IMG *	iimg,
		int	calc_mode,
		char *	status,
		int	verbose )

Integration of dynamic image from time1 to time2, storing integrals in iimg, which is allocated here. Frames do not have to be continuous in time. Time unit in integral is sec. Raw data (sinogram) must be divided by frame durations before calling this. You may want to test the time range before applying this routine, because this function accepts relatively large extrapolation.

See also: imgFrameIntegral, dftInterpolateForIMG, imgExistentTimes

Returns: Returns errstatus, which is STATUS_OK (0) when call was successful, and >0 in case of an error.

Parameters

img	IMG data; preferably dynamic, if static, then the specified time range must match with the frame time.
t1	Time where to start integration (sec).
t2	Time to stop integration (sec); must be higher than t1.
iimg	Pointer to initiated but empty AUC IMG data.
calc_mode	Calculate integral or average: 0=integral, 1=average.
status	Pointer to a string (allocated for at least 128 chars) where error message or other execution status will be written; enter NULL, if not needed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 147 of file misc_model.c.

  {
  int ret;
  double accept_tdif=1.0; // in seconds
 
  if(verbose>0) {
    printf("%s(img, %g, %g, iimg, %d, status, %d)\n", __func__, t1, t2, calc_mode, verbose);
    fflush(stdout);
  }
 
  /* Check the arguments */
  if(status!=NULL) sprintf(status, "program error");
  if(img==NULL || iimg==NULL || t1<0.0 || t2<0.0) return STATUS_FAULT;
  if(img->status!=IMG_STATUS_OCCUPIED) return STATUS_FAULT;
  float fdur=t2-t1; if(fdur<=0.0) return STATUS_FAULT;
  if(img->dimt<1) return STATUS_FAULT;
 
  /* Check that time range matches with image frame */
  if(img->dimt==1) { // static image (one frame only)
    /* check that specified times match with image frame times, but accept 1 s differences */
    if(fabs(img->end[0]-t2)>accept_tdif || fabs(img->start[0]-t1)>accept_tdif) {
      if(status!=NULL) sprintf(status, 
        "for static image the integration time range must be exactly as long as the scan");
      return STATUS_FAULT;
    }
    /* Set the times, in case there was a small difference */
    img->end[0]=t2; img->start[0]=t1; img->mid[0]=0.5*(t1+t2);
  } else { // Dynamic image (more than one frame)
    /* First PET frame must start before 1/3 of the integration time */
    /* Last PET frame must end after 2/3 of integration time */
    if(img->start[0]>(0.66*t1+0.34*t2) || img->end[img->dimt-1]<(0.34*t1+0.66*t2)) {
      if(status!=NULL) sprintf(status, "integration time range oversteps data range");
      return STATUS_FAULT;
    }
  }
  if(verbose>10) printf("t1=%g t2=%g fdur=%g\n", t1, t2, fdur);
 
  /* Get the first and last frame index that resides inside integration time */
  int f1, f2;
  if(img->dimt==1) {
    f1=f2=0;
  } else {
    f1=f2=-1;
    for(int fi=0; fi<img->dimt; fi++) {
      if(f1<0) {if(img->start[fi]>=t1 && img->end[fi]<=t2) f1=f2=fi;}
      if(f1>=0) {if(t2>=img->end[fi]) f2=fi;}
    }
  }
  if(verbose>10) printf("f1=%d f2=%d\n", f1, f2);
 
  float aucroi, t[2], auc[2];
  if(f1>=0 && f2>=0) {
 
    /* Integrate over the frames that are included in time range as a whole */
    ret=imgFrameIntegral(img, f1, f2, iimg, verbose-1);
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot integrate (%d)", ret);
      return STATUS_FAULT;
    }
 
    /* If necessary, add the partial integrals */
    if(img->start[f1]>t1) {
      t[0]=t1; t[1]=img->start[f1];
      if(verbose>20) printf("t[0]=%g t[1]=%g\n", t[0], t[1]); 
      for(int zi=0; zi<img->dimz; zi++)
        for(int yi=0; yi<img->dimy; yi++)
          for(int xi=0; xi<img->dimx; xi++) {
            ret=finterpolate(img->mid, img->m[zi][yi][xi], img->dimt, t, NULL, auc, NULL, 2);
            if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
            iimg->m[zi][yi][xi][0]+=aucroi;
          }
    }
    if(t2>img->end[f2]) {
      t[0]=img->end[f2]; t[1]=t2;
      if(verbose>20) printf("t[0]=%g t[1]=%g\n", t[0], t[1]);
      for(int zi=0; zi<img->dimz; zi++)
        for(int yi=0; yi<img->dimy; yi++)
          for(int xi=0; xi<img->dimx; xi++) {
            ret=finterpolate(img->mid, img->m[zi][yi][xi], img->dimt, t, NULL, auc, NULL, 2);
            if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
            iimg->m[zi][yi][xi][0]+=aucroi;
          }
    }
 
  } else { // no full frames inside integration range
 
    /* Create the AUC image */
    ret=imgAllocateWithHeader(iimg, img->dimz, img->dimy, img->dimx, 1, img);
    if(ret!=0) {
      if(status!=NULL) sprintf(status, "cannot setup integral image");
      return STATUS_FAULT;
    }
 
    t[0]=t1; t[1]=t2;
    for(int zi=0; zi<img->dimz; zi++)
      for(int yi=0; yi<img->dimy; yi++)
        for(int xi=0; xi<img->dimx; xi++) {
          ret=finterpolate(img->mid, img->m[zi][yi][xi], img->dimt, t, NULL, auc, NULL, 2);
          if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
          iimg->m[zi][yi][xi][0]+=aucroi;
        }
 
  }
 
  /* Set output image time frame */
  iimg->end[0]=t2; iimg->start[0]=t1; iimg->mid[0]=0.5*(t1+t2);
 
  /* If required, then calculate average by dividing integral with time. */
  /* Set units accordingly */
  if(calc_mode!=0) { // average
    ret=imgArithmConst(iimg, fdur, ':', 1.0E+10, verbose-1);
    if(ret!=0) {
      if(status!=NULL) sprintf(status, "cannot divide integral image");
      return STATUS_FAULT;
    }
    iimg->unit=img->unit;
    if(status!=NULL) sprintf(status, "average image [%g,%g] calculated",t1,t2);
  } else { // integral
    if(img->unit==CUNIT_KBQ_PER_ML) iimg->unit=CUNIT_SEC_KBQ_PER_ML;
    else iimg->unit=CUNIT_UNKNOWN;
    if(status!=NULL) sprintf(status, "integral image [%g,%g] calculated",t1,t2);
  }
 
  imgSetStatus(iimg, STATUS_OK);
  return(0);
}

◆ mrl_between_tacs()

int mrl_between_tacs	(	double *	y1,
		double *	y2,
		int	n )

Return the maximum run length between given n length arrays of data

Parameters

y1	Array of data1; may contain NaNs
y2	Array of data2; may contain NaNs
n	Nr of samples in array 1 and 2

Definition at line 13 of file mrl.c.

  {
  int i, mrl=0, rl=0;
  char last_sign=0, sign;
 
  if(n<1 || y1==NULL || y2==NULL) return(0);
  for(i=0; i<n; i++) {
    if(isnan(y1[i]) || isnan(y2[i])) continue;
    if(y1[i]>y2[i]) sign=1; else if(y1[i]<y2[i]) sign=-1; else sign=0;
    if(sign!=last_sign) {
      rl=0; last_sign=sign;
    } else {
      if(sign!=0) {rl++; if(rl>mrl) mrl=rl;}
    }
  }
  return(mrl);
}

◆ noiseSD4Simulation()

int noiseSD4Simulation	(	double	y,
		double	t1,
		double	dt,
		double	hl,
		double	a,
		double *	sd,
		char *	status,
		int	verbose )

Calculate SD for PET radioactivity concentration data to be used to simulate noise.

Note that SD is dependent on the time units.

Reference: Varga & Szabo. J Cereb Blood Flow Metab 2002;22(2):240-244.

Returns: Returns 0 when successful, otherwise <>0.

Parameters

y	Sample radioactivity concentration (decay corrected to zero time)
t1	Radioactivity measurement (frame) start time (in same units as the halflife)
dt	Radioactivity measurement (frame) duration (in same units as the halflife)
hl	Isotope halflife (in same units as the sample time); enter 0 if not to be considered
a	Proportionality factor. Note that it is inside square root.
sd	Pointer in which SD is written
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 21 of file noise.c.

  {
  double d, lambda, f;
 
  if(verbose>0) printf("noiseSD4Simulation(%g, %g, %g, %g, %g, *sd, ...)\n",
                       y, t1, dt, hl, a);
  if(status!=NULL) strcpy(status, "invalid data");
  if(sd==NULL) return 1;
  if(t1<0.0) return 2;
  if(dt<=0.0) return 3;
 
  /* Decay factor (<=1) */
  if(hl<=0.0) {
    d=1.0; lambda=0.0;
  } else {
    if(status!=NULL) strcpy(status, "invalid half-life");
    lambda=hl2lambda(hl); if(lambda<0.0) return 4;
    d=hlLambda2factor(-lambda, t1, dt); if(d<0.0) return 5;
  }
  /* SD */
  f=y*d*dt;
  if(f<=0.0) *sd=0.0;
  else *sd=y*sqrt(a/f);
 
  if(status!=NULL) strcpy(status, "ok");
  return(0);
}

Referenced by noiseSD4SimulationFromDFT().

◆ noiseSD4SimulationFromDFT()

int noiseSD4SimulationFromDFT	(	DFT *	dft,
		int	index,
		double	pc,
		double *	sd,
		char *	status,
		int	verbose )

Calculate SD for noise simulation from TAC data.

Sample times will be converted to minutes if necessary.

Reference: Varga & Szabo. J Cereb Blood Flow Metab 2002;22(2):240-244.

Returns: Returns 0 when successful, otherwise <>0.

Parameters

dft	Pointer to TAC data in DFT struct, based on which the SD for each frame is calculated. Contents are not changed. Struct must contain correct values for the isotope, time unit, and frame times (start and end). Status of decay correction is used, and if not known, then data is assumed to be decay corrected.
index	TAC index [0..voiNr-1] which is used to calculate the SD in case DFT contains more than one TAC. Non-effective if DFT contains only one TAC. Set to <0 if mean of all TACs is to be used.
pc	Proportionality factor
sd	Pointer to allocated array in which SD is written
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 79 of file noise.c.

  {
  int ret, fi;
  double hl=0.0, t1, deltat;
  double *y;
  DFT mean;
 
  if(verbose>0) printf("noiseSD4SimulationFromDFT(DFT, %d, %g, sd[], ...)\n",
                       index, pc);
  if(status!=NULL) strcpy(status, "invalid data");
  if(dft==NULL || sd==NULL) return 1;
  if(dft->voiNr<1 || dft->frameNr<1) return 2;
  if(dft->voiNr>1 && index>=dft->voiNr) return 3;
 
  /* Check that valid frame times are available */
  if(status!=NULL) strcpy(status, "invalid frame times");
  for(fi=0; fi<dft->frameNr; fi++)
    if(dft->x2[fi]<=dft->x1[fi] || dft->x1[fi]<0.0) return 4;
  if(status!=NULL) strcpy(status, "missing time unit");
  if(dft->timeunit!=TUNIT_SEC && dft->timeunit!=TUNIT_MIN) return 5;
 
  /* Check if isotope is available, if it is needed */
  if(dft->decayCorrected==DFT_DECAY_CORRECTED || 
     dft->decayCorrected==DFT_DECAY_UNKNOWN)
  {
    hl=hlFromIsotope(dft->isotope);
    if(status!=NULL) strcpy(status, "missing isotope halflife");
    if(hl<=0.0) return 6;
  }
  if(verbose>1) printf("halflife := %g\n", hl);
 
  /* Set pointer to activity concentration data */
  dftInit(&mean);
  if(dft->voiNr==1) {
    y=dft->voi[0].y;
  } else if(index>=0) {
    y=dft->voi[index].y;
  } else {
    /* Compute TAC mean */
    ret=dftMeanTAC(dft, &mean); if(ret!=0) {
      if(status!=NULL) strcpy(status, "cannot calculate mean TAC");
      return 8;
    }
    /* and use that */
    y=mean.voi[0].y;
  }
 
  /* Compute SD for each sample time */
  if(pc<=0.0) pc=1.0;
  for(fi=0; fi<dft->frameNr; fi++) {
    t1=dft->x1[fi]; deltat=dft->x2[fi]-dft->x1[fi];
    if(dft->timeunit==TUNIT_SEC) {t1/=60.0; deltat/=60.0;}
    ret=noiseSD4Simulation(y[fi], t1, deltat, hl, pc, sd+fi, status, verbose);
    if(ret!=0) break;
  }
  dftEmpty(&mean);
  if(ret!=0) {
    if(status!=NULL)
      sprintf(status, "cannot calculate SD for noise simulation (%d)", ret);
    return(100+ret);
  }
 
  if(status!=NULL) strcpy(status, "ok");
  return 0;
}

◆ plot_fit_svg()

int plot_fit_svg	(	DFT *	dft1,
		DFT *	dft2,
		char *	main_title,
		char *	fname,
		int	verbose )

Writes plots of original and fitted TACs in SVG 1.1 format. Data must not contain NaNs.

Returns: Returns 0 if successful, otherwise nonzero.

See also: plot_svg, plot_fitrange_svg

Parameters

dft1	Measured data points.
dft2	Fitted data points. Times can be different but unit must be the same.
main_title	String for plot main title, or "".
fname	SVG filename; existing file is backed up.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 216 of file plotfit.c.

  {
  int ret, n, ri, si, ei;
  char x_title[64], y_title[64], tac_id[32], tac_title[64];
  double minx, maxx, miny, maxy, tx1, tx2, ty1, ty2, f;
  struct svg_viewports viewports; svg_init_viewports(&viewports);
  int max_color_nr, color_nr;
  int max_symbol_nr, symbol_nr;
  SVG_LEGENDS legends; svg_init_legends(&legends);
  FILE *fp_svg=NULL;
 
 
  if(verbose>0) {
    printf("%s(dft1, dft2, mt, fn, %d)\n", __func__, verbose);
  }
 
  /* Check data */
  if(dft1==NULL || dft1->voiNr<1) return(1);
  if(dft2==NULL) return(1);
  if(dft2->voiNr!=dft1->voiNr) return(1);
 
  int is_label=0; if(dft1->voiNr>1) is_label=1;
 
  /* Check if file exists; backup, if necessary */
  ret=backupExistingFile(fname, NULL, NULL); if(ret) return(2);
 
  /* Determine the plot min and max values */
  minx=maxx=miny=maxy=0;
  ret=dftMinMax(dft1, &tx1, &tx2, &ty1, &ty2); if(ret) return(3);
  minx=tx1; maxx=tx2; miny=ty1; maxy=ty2;
  ret=dftMinMax(dft2, &tx1, &tx2, &ty1, &ty2); if(ret) return(3);
  if(minx>tx1) minx=tx1;
  if(maxx<tx2) maxx=tx2;
  if(miny>ty1) miny=ty1;
  if(maxy<ty2) maxy=ty2;
  if(verbose>1) printf("minx:=%g\nmaxx:=%g\nminy:=%g\nmaxy:=%g\n", minx, maxx, miny, maxy);
  if(miny>0.0) {f=maxy-miny; miny-=0.01*f;}
 
  /* Calculate the axis ticks */
  viewports.label_area_viewport.is=is_label; // needed for x axis ticks
  viewports.x.fixed_min=0; viewports.y.fixed_min=0;
  viewports.x.min=minx; viewports.x.max=maxx;
  viewports.y.min=miny; viewports.y.max=maxy;
  ret=svg_calculate_axes(&viewports, verbose-3); if(ret) return(4);
 
  /* Set x and y axis titles based on activity and time units */
  if(verbose>2) printf("set x title\n");
  strcpy(x_title, "");
  if(dft1->timeunit==DFTTIME_SEC || dft1->timeunit==DFTTIME_MIN) 
    sprintf(x_title, "Time (%s)", dftTimeunit(dft1->timeunit));
  else if(dft1->timeunit!=DFTTIME_UNKNOWN)
    strcpy(x_title, dftTimeunit(dft1->timeunit));
  if(verbose>2) printf("set y title\n");
  strcpy(y_title, "");
  if(dftUnitId(dft1->unit)!=DFTUNIT_UNKNOWN)
    strcpy(y_title, dftUnit(dftUnitId(dft1->unit)));
 
  /* Set the plot window and window area sizes */
  if(verbose>2) printf("set window sizes\n");
  ret=svg_define_viewports(0, 0, strlen(main_title), strlen(y_title),
    strlen(x_title), is_label, &viewports, verbose-3);
  if(ret) return(5);
 
  /* Initiate graphics file */
  fp_svg=svg_initiate(fname, 0, 0, &viewports, NULL, verbose-3);
  if(fp_svg==NULL) return(6);
 
  /* Put the graph titles into their own viewports */
  ret=svg_create_main_title(fp_svg, main_title, "", &viewports, NULL,verbose-3);
  if(ret) return(7);
  ret=svg_create_yaxis_title(fp_svg, y_title, &viewports, NULL, verbose-3);
  if(ret) return(8);
  ret=svg_create_xaxis_title(fp_svg, x_title, &viewports, NULL, verbose-3);
  if(ret) return(9);
 
  /* Put the plot into its own viewport */
  ret=svg_start_plot_viewport(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(10);
 
  /*  Start coordinate area viewport */
  ret=svg_start_coordinate_viewport(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(11);
 
  /* Write plot axes */
  ret=svg_write_axes(fp_svg, &viewports, NULL, verbose-3);
  if(ret) {
    if(verbose>0) printf("svg_write_axes() := %d\n", ret);
    return(12);
  }
 
  /*
   *  Draw the plots
   */
  max_color_nr=0; while(svgColorName(max_color_nr)!=NULL) max_color_nr++;
  if(max_color_nr<1) max_color_nr=1; // remainder works only if 2nd operator>0
  if(verbose>3) printf("max_color_nr := %d\n", max_color_nr);
  max_symbol_nr=0; while(svgSymbolName(max_symbol_nr)!=NULL) max_symbol_nr++;
  if(max_symbol_nr<1) max_symbol_nr=1; // remainder works only if 2nd operator>0
  if(verbose>3) printf("max_symbol_nr := %d\n", max_symbol_nr);
  if(dft1->voiNr==1) color_nr=0; else color_nr=1;
  symbol_nr=0;
  for(ri=0, n=0; ri<dft1->voiNr; ri++) {
    sprintf(tac_id, "plot_%d", n);
/*
    if(strlen(dft1->studynr)>0 && strcmp(dft1->studynr, ".")!=0)
      sprintf(tac_title, "%s: %s", dft1->studynr, dft1->voi[ri].name);
    else strcpy(tac_title, dft1->voi[ri].name);
*/
    rnameRmDots(dft1->voi[ri].name, tac_title);
    /* Draw the fitted line */
    for(si=0; si<dft2->frameNr; si++) if(!isnan(dft2->voi[ri].y[si])) break;
    for(ei=dft2->frameNr-1; ei>=si; ei--) if(!isnan(dft2->voi[ri].y[ei])) break;
    if((ei-si)>0) {
      ret=svg_write_tac(fp_svg, &viewports, 1, tac_id, tac_title,
            dft2->x+si, dft2->voi[ri].y+si, 1+ei-si,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLFILLED,
            NULL, verbose-3);
      if(ret) {svg_legend_empty(&legends); return(21);}
    }
    /* Draw the measured points */
    ret=svg_write_tac(fp_svg, &viewports, 2, tac_id, tac_title,
            dft1->x, dft1->voi[ri].y, dft1->frameNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLFILLED,
            NULL, verbose-3);
    if(ret) {svg_legend_empty(&legends); return(22);}
    /* Set legend too, if requested */
    if(is_label!=0) {
      svg_legend_add(&legends, 0, symbol_nr%max_symbol_nr, SYMBOLFILLED,
                     color_nr%max_color_nr, tac_title);
    }
    /* Prepare for the next plot */
    color_nr++; n++;
    if(color_nr==max_color_nr) {symbol_nr++; color_nr=0;}
    if(symbol_nr==max_symbol_nr) symbol_nr=0;
  }
 
  /* Close the coordinate viewport */
  ret=svg_end_coordinate_viewport(fp_svg, NULL, verbose-3);
  if(ret) {svg_legend_empty(&legends); return(91);}
 
  /* Write the axis ticks */
  if(svg_write_xticks(fp_svg, &viewports, NULL, verbose-3)!=0) {
    svg_legend_empty(&legends); return(92);}
  if(svg_write_yticks(fp_svg, &viewports, NULL, verbose-3)!=0) {
    svg_legend_empty(&legends); return(93);}
 
  /* Close the plot viewport */
  ret=svg_end_plot_viewport(fp_svg, NULL, verbose-3);
  if(ret) {svg_legend_empty(&legends); return(94);}
 
  /* Make the plot legends into their own viewport */
  if(viewports.label_area_viewport.is!=0) {
    if(verbose>2) printf("creating plot legends\n");
    ret=svg_create_legends(fp_svg, &viewports, &legends, NULL, verbose-3);
    if(ret) {svg_legend_empty(&legends); return(95);}
  }
  svg_legend_empty(&legends);
 
  /* Close the SVG file */
  ret=svg_close(fp_svg, NULL, verbose-3); if(ret) return(101);
 
  return(0);
}

◆ plot_fitrange_svg()

int plot_fitrange_svg	(	DFT *	dft1,
		DFT *	dft2,
		char *	main_title,
		double	x1,
		double	x2,
		double	y1,
		double	y2,
		char *	fname,
		int	verbose )

Writes specified range of plots of original and fitted TACs in SVG 1.1 format.

Returns: Returns 0 if successful, otherwise nonzero.

See also: plot_fit_svg, plot_svg

Parameters

dft1	Measured data points
dft2	Fitted data points. Times can be different but unit must be the same.
main_title	String for plot main title, or NULL
x1	Start time; NaN if determined from data
x2	End time; NaN if determined from data
y1	Minimum y value; NaN if determined from data
y2	Maximum y value; NaN if determined from data
fname	SVG filename; existing file is backed up
verbose	Verbose level; set to zero to not to print any comments

Definition at line 16 of file plotfit.c.

  {
  int ret, n, ri, si, ei;
  char x_title[64], y_title[64], tac_id[32], tac_title[64];
  double minx, maxx, miny, maxy, tx1, tx2, ty1, ty2, f;
  struct svg_viewports viewports; svg_init_viewports(&viewports);
  int max_color_nr, color_nr;
  int max_symbol_nr, symbol_nr;
  SVG_LEGENDS legends; svg_init_legends(&legends);
  FILE *fp_svg=NULL;
 
  if(verbose>0) {
    printf("%s(dft1, dft2, mt, x1, x2, y1, y2, fn, %d)\n", __func__, verbose);
  }
 
  /* Check data */
  if(dft1==NULL || dft1->voiNr<1) return(1);
  if(dft2==NULL) return(1);
  if(dft2->voiNr!=dft1->voiNr) return(1);
 
  /* Check if file exists; backup, if necessary */
  ret=backupExistingFile(fname, NULL, NULL); if(ret) return(2);
 
  int is_label=0; if(dft1->voiNr>1) is_label=1;
 
  /* Determine the plot min and max x values */
  ret=dftMinMax(dft1, &tx1, &tx2, NULL, NULL); if(ret) return(3);
  minx=tx1; maxx=tx2;
  ret=dftMinMax(dft2, &tx1, &tx2, NULL, NULL); if(ret) return(3);
  if(minx>tx1) minx=tx1; 
  if(maxx<tx2) maxx=tx2;
  if(minx>0.0) {
    f=maxx-minx; minx-=0.05*f; 
    if(minx<0.0) minx=0.0;
  }
  if(!isnan(x1)) minx=x1; 
  if(!isnan(x2)) maxx=x2;
 
  /* Determine the plot min and max y values */
  ret=dftMaxY(dft1, minx, maxx, &ty1, &ty2); if(ret) return(3);
  miny=ty1; maxy=ty2;
  ret=dftMaxY(dft2, minx, maxx, &ty1, &ty2); if(ret) return(3);
  if(miny>ty1) miny=ty1; 
  if(maxy<ty2) maxy=ty2;
  if(miny>0.0) {
    f=maxy-miny; miny-=0.05*f; 
    if(miny<0.0) miny=0.0;
  }
  if(!isnan(y1)) miny=y1; 
  if(!isnan(y2)) maxy=y2;
 
  if(verbose>1) printf("minx:=%g\nmaxx:=%g\nminy:=%g\nmaxy:=%g\n", minx, maxx, miny, maxy);
 
  /* Calculate the axis ticks */
  viewports.label_area_viewport.is=is_label; // needed for x axis ticks
  viewports.x.min=minx; viewports.x.max=maxx;
  viewports.y.min=miny; viewports.y.max=maxy;
  if(isnan(x1) || isnan(x2)) viewports.x.fixed_min=0; 
  if(isnan(y1) || isnan(y2)) viewports.y.fixed_min=0; 
  else viewports.y.fixed_min=1; 
  ret=svg_calculate_axes(&viewports, verbose-3); if(ret) return(4);
 
  /* Set x and y axis titles based on activity and time units */
  strcpy(x_title, "");
  if(dft1->timeunit==DFTTIME_SEC || dft1->timeunit==DFTTIME_MIN) 
    sprintf(x_title, "Time (%s)", dftTimeunit(dft1->timeunit));
  else if(dft1->timeunit!=DFTTIME_UNKNOWN)
    strcpy(x_title, dftTimeunit(dft1->timeunit));
  strcpy(y_title, "");
  if(dftUnitId(dft1->unit)!=DFTUNIT_UNKNOWN)
    strcpy(y_title, dftUnit(dftUnitId(dft1->unit)));
 
  /* Set the plot window and window area sizes */
  ret=svg_define_viewports(0, 0, strlen(main_title), strlen(y_title),
    strlen(x_title), is_label, &viewports, verbose-3);
  if(ret) return(5);
 
  /* Initiate graphics file */
  fp_svg=svg_initiate(fname, 0, 0, &viewports, NULL, verbose-3);
  if(fp_svg==NULL) return(6);
 
  /* Put the graph titles into their own viewports */
  ret=svg_create_main_title(fp_svg, main_title, "", &viewports, NULL,verbose-3);
  if(ret) return(7);
  ret=svg_create_yaxis_title(fp_svg, y_title, &viewports, NULL, verbose-3);
  if(ret) return(8);
  ret=svg_create_xaxis_title(fp_svg, x_title, &viewports, NULL, verbose-3);
  if(ret) return(9);
 
  /* Put the plot into its own viewport */
  ret=svg_start_plot_viewport(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(10);
 
  /*  Start coordinate area viewport */
  ret=svg_start_coordinate_viewport(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(11);
 
  /* Write plot axes */
  ret=svg_write_axes(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(12);
 
  /*
   *  Draw the plots
   */
  max_color_nr=0; while(svgColorName(max_color_nr)!=NULL) max_color_nr++;
  if(max_color_nr==0) max_color_nr=1; // remainder works only if 2nd operator>0
  if(verbose>3) printf("max_color_nr := %d\n", max_color_nr);
  max_symbol_nr=0; while(svgSymbolName(max_symbol_nr)!=NULL) max_symbol_nr++;
  if(max_symbol_nr==0) max_symbol_nr=1; // remainder works only if 2nd operator>0
  if(verbose>3) printf("max_symbol_nr := %d\n", max_symbol_nr);
  if(dft1->voiNr==1) color_nr=0; else color_nr=1;
  symbol_nr=0;
  for(ri=0, n=0; ri<dft1->voiNr; ri++) {
    sprintf(tac_id, "plot_%d", n);
/*
    if(strlen(dft1->studynr)>0 && strcmp(dft1->studynr, ".")!=0)
      sprintf(tac_title, "%s: %s", dft1->studynr, dft1->voi[ri].name);
    else strcpy(tac_title, dft1->voi[ri].name);
*/
    rnameRmDots(dft1->voi[ri].name, tac_title);
    /* Draw the fitted line */
    for(si=0; si<dft2->frameNr; si++) if(!isnan(dft2->voi[ri].y[si])) break;
    for(ei=dft2->frameNr-1; ei>=si; ei--) if(!isnan(dft2->voi[ri].y[ei])) break;
    if((ei-si)>0) {
      ret=svg_write_tac(fp_svg, &viewports, 1, tac_id, tac_title,
            dft2->x+si, dft2->voi[ri].y+si, 1+ei-si,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLFILLED,
            NULL, verbose-3);
      if(ret) {svg_legend_empty(&legends); return(21);}
    }
    /* Draw the measured points */
    ret=svg_write_tac(fp_svg, &viewports, 2, tac_id, tac_title,
            dft1->x, dft1->voi[ri].y, dft1->frameNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLFILLED,
            NULL, verbose-3);
    if(ret) {svg_legend_empty(&legends); return(22);}
    /* Set legend too, if requested */
    if(is_label!=0) {
      svg_legend_add(&legends, 0, symbol_nr%max_symbol_nr, SYMBOLFILLED,
                     color_nr%max_color_nr, tac_title);
    }
    /* Prepare for the next plot */
    color_nr++; n++;
    if(color_nr==max_color_nr) {symbol_nr++; color_nr=0;}
    if(symbol_nr==max_symbol_nr) symbol_nr=0;
  }
 
  /* Close the coordinate viewport */
  ret=svg_end_coordinate_viewport(fp_svg, NULL, verbose-3);
  if(ret) {svg_legend_empty(&legends); return(91);}
 
  /* Write the axis ticks */
  if(svg_write_xticks(fp_svg, &viewports, NULL, verbose-3)!=0) {
    svg_legend_empty(&legends); return(92);}
  if(svg_write_yticks(fp_svg, &viewports, NULL, verbose-3)!=0) {
    svg_legend_empty(&legends); return(93);}
 
  /* Close the plot viewport */
  ret=svg_end_plot_viewport(fp_svg, NULL, verbose-3);
  if(ret) {svg_legend_empty(&legends); return(94);}
 
  /* Make the plot legends into their own viewport */
  if(viewports.label_area_viewport.is!=0) {
    if(verbose>2) printf("creating plot legends\n");
    ret=svg_create_legends(fp_svg, &viewports, &legends, NULL, verbose-3);
    if(ret) {svg_legend_empty(&legends); return(95);}
  }
  svg_legend_empty(&legends);
 
  /* Close the SVG file */
  ret=svg_close(fp_svg, NULL, verbose-3); if(ret) return(101);
 
  return(0);
}

◆ plot_svg()

int plot_svg	(	DFT *	dft,
		RES *	res,
		int	first,
		int	last,
		char *	main_title,
		char *	x_title,
		char *	y_title,
		int	color_scale,
		char *	fname,
		int	verbose )

Writes graphical analysis plots in SVG 1.1 format. Assumes that line slope and ic are in res->parameter[0] and [1].

Returns: Returns 0 if successful, otherwise nonzero.

See also: plot_fit_svg, plot_fitrange_svg

Parameters

dft	Plot points: X in y2, Y in y3.
res	Results containing parameters of line.
first	First sample (starting from 0) used in linear fit.
last	last sample (starting from 0) used in linear fit.
main_title	String for plot main title, or NULL.
x_title	String for X axis title, or NULL.
y_title	String for Y axis title, or NULL.
color_scale	Colour-scale: 0=colour, 2=black-and-white.
fname	File name for SVG; existing file is renamed as *.bak.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 17 of file plotdata.c.

  {
  int n, ri, fi, ret;
  char tac_id[32], tac_title[64];
  double maxPlotX=0, maxy;
  double px[2], py[2];
  struct svg_viewports viewports; svg_init_viewports(&viewports);
  int max_color_nr, color_nr;
  int max_symbol_nr, symbol_nr;
  SVG_LEGENDS legends; svg_init_legends(&legends);
  FILE *fp_svg=NULL;
 
  if(verbose>0) {
    printf("%s(dft, res, %d, %d, mt, xt, yt, fn, cs, %d)\n", __func__, first, last, verbose);
  }
 
  /* Check data */
  if(dft==NULL || dft->voiNr<1) return(1);
  if(res==NULL || res->voiNr!=dft->voiNr) return(1);
  if(first>last) return(1);
 
  int is_label=0; if(dft->voiNr>1) is_label=1;
  if(color_scale!=0 && color_scale!=2) color_scale=0; 
 
  /* Check if file exists; backup, if necessary */
  backupExistingFile(fname, NULL, NULL);
 
  /* Search the largest plot x-value to be used as line end point */
  for(ri=0; ri<dft->voiNr; ri++) for(fi=0; fi<dft->frameNr; fi++)
    if(dft->voi[ri].y2[fi]>maxPlotX) maxPlotX=dft->voi[ri].y2[fi];
  /* Get maxy */
  maxy=0.0;
  for(ri=0; ri<dft->voiNr; ri++) for(fi=0; fi<dft->frameNr; fi++)
    if(dft->voi[ri].y3[fi]>maxy) maxy=dft->voi[ri].y3[fi];
 
  /* Calculate the axis ticks */
  viewports.label_area_viewport.is=is_label; // needed for x axis ticks
  viewports.x.fixed_min=0; viewports.y.fixed_min=0;
  viewports.x.min=0.0; viewports.x.max=maxPlotX;
  viewports.y.min=0.0; viewports.y.max=maxy;
  ret=svg_calculate_axes(&viewports, verbose-3);
  if(ret) return(2);
 
  /* Set the plot window and window area sizes */
  ret=svg_define_viewports(0, 0, strlen(main_title), strlen(y_title),
                           strlen(x_title), is_label, &viewports, verbose-3);
  if(ret) return(3);
 
  /* Initiate graphics file */
  fp_svg=svg_initiate(fname, 0, 0, &viewports, NULL, verbose-3);
  if(fp_svg==NULL) return(4);
 
  /* Put the graph titles into their own viewports */
  ret=svg_create_main_title(fp_svg, main_title, "", &viewports, NULL,verbose-3);
  if(ret) return(5);
  ret=svg_create_yaxis_title(fp_svg, y_title, &viewports, NULL, verbose-3);
  if(ret) return(6);
  ret=svg_create_xaxis_title(fp_svg, x_title, &viewports, NULL, verbose-3);
  if(ret) return(7);
 
  /* Put the plot into its own viewport */
  ret=svg_start_plot_viewport(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(8);
 
  /*  Start coordinate area viewport */
  ret=svg_start_coordinate_viewport(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(9);
 
  /* Write plot axes */
  ret=svg_write_axes(fp_svg, &viewports, NULL, verbose-3);
  if(ret) return(10);
 
  /*
   *  Draw the plots
   */
  int symbol_fill=SYMBOLFILLED;
  if(color_scale==2) { // black-and-white
    symbol_fill=SYMBOLOPEN;
    max_color_nr=1;
  } else { // colour
    max_color_nr=0; while(svgColorName(max_color_nr)!=NULL) max_color_nr++;
  }
  if(max_color_nr==0) max_color_nr=1; // remainder works only if 2nd operator>0
  if(verbose>3) printf("max_color_nr := %d\n", max_color_nr);
  max_symbol_nr=0; while(svgSymbolName(max_symbol_nr)!=NULL) max_symbol_nr++;
  if(max_symbol_nr==0) max_symbol_nr=1; // remainder works only if 2nd operator>0
  if(verbose>3) printf("max_symbol_nr := %d\n", max_symbol_nr);
  if(dft->voiNr==1) color_nr=0; else color_nr=1;
  symbol_nr=0;
  for(ri=0, n=0; ri<dft->voiNr; ri++) {
    sprintf(tac_id, "plot_%d", n);
    //printf("ri=%d color_nr=%d symbol_nr=%d\n", ri, color_nr, symbol_nr);
/*
    if(strlen(dft->studynr)>0 && strcmp(dft->studynr, ".")!=0)
      sprintf(tac_title, "%s: %s", dft->studynr, dft->voi[ri].name);
    else strcpy(tac_title, dft->voi[ri].name);
*/
    rnameRmDots(dft->voi[ri].name, tac_title);
    /*printf("tac_title := %s ; color := %s ; symbol := %d\n",
      tac_title, svgcolor[color_nr%max_color_nr], symbol_nr%(max_symbol_nr+1));*/
    /* Draw symbols */
    if(dft->frameNr<150) {
      /* plot symbols only if less than 150 samples */
      ret=svg_write_tac(fp_svg, &viewports, 2, tac_id, tac_title,
            dft->voi[ri].y2, dft->voi[ri].y3, dft->frameNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, symbol_fill,
            NULL, verbose-3);
    } else {
      /* plot samples as line if 150 samples or more */
      ret=svg_write_tac(fp_svg, &viewports, 1, tac_id, tac_title,
            dft->voi[ri].y2, dft->voi[ri].y3, dft->frameNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, symbol_fill,
            NULL, verbose-3);
    }
    if(ret) {svg_legend_empty(&legends); return(21);}
    /* Draw the line */
    px[0]=0.0; py[0]=res->voi[ri].parameter[1];
    px[1]=maxPlotX;
    py[1]=maxPlotX*res->voi[ri].parameter[0]+res->voi[ri].parameter[1];
    sprintf(tac_id, "line_%d", n);
    ret=svg_write_tac(fp_svg, &viewports, 1, tac_id, tac_title,
            px, py, 2,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, symbol_fill,
            NULL, verbose-3);
    if(ret) {svg_legend_empty(&legends); return(22);}
    /* Set legend too, if requested */
    if(is_label!=0) {
      svg_legend_add(&legends, 0, symbol_nr%max_symbol_nr, symbol_fill, color_nr%max_color_nr, tac_title);
    }
 
    /* Prepare for the next plot */
    if(color_scale==0) {
      color_nr++;
      if(color_nr==max_color_nr) {symbol_nr++; color_nr=0;}
      if(symbol_nr==max_symbol_nr) symbol_nr=0;
    }
    if(color_scale==2) {
      symbol_nr++;
      if(symbol_nr==max_symbol_nr) {
        if(symbol_fill==SYMBOLOPEN) symbol_fill=SYMBOLFILLED; else symbol_fill=SYMBOLOPEN;
        symbol_nr=0;
      }
    }
    n++;
  }
 
  /* Close the coordinate viewport */
  ret=svg_end_coordinate_viewport(fp_svg, NULL, verbose-3);
  if(ret) {svg_legend_empty(&legends); return(91);}
 
  /* Write the axis ticks */
  if(svg_write_xticks(fp_svg, &viewports, NULL, verbose-3)!=0) {
    svg_legend_empty(&legends); return(92);}
  if(svg_write_yticks(fp_svg, &viewports, NULL, verbose-3)!=0) {
    svg_legend_empty(&legends); return(93);}
 
  /* Close the plot viewport */
  ret=svg_end_plot_viewport(fp_svg, NULL, verbose-3);
  if(ret) {svg_legend_empty(&legends); return(94);}
 
  /* Make the plot legends into their own viewport */
  if(viewports.label_area_viewport.is!=0) {
    if(verbose>2) printf("creating plot legends\n");
    ret=svg_create_legends(fp_svg, &viewports, &legends, NULL, verbose-3);
    if(ret) {svg_legend_empty(&legends); return(95);}
  }
  svg_legend_empty(&legends);
 
  /* Close the SVG file */
  ret=svg_close(fp_svg, NULL, verbose-3); if(ret) return(101);
  return(0);
}

◆ plotdata()

int plotdata	(	DFT *	dft,
		RES *	res,
		int	first,
		int	last,
		char *	mtitle,
		char *	xtitle,
		char *	ytitle,
		char *	fname )

Write plot and line fit data in XHTML 1.1. Strict table format. Assumes that line slope and ic are in res->parameter[0] and [1].

Returns: Returns 0 if successful, otherwise non-zero.

See also: plotdata_as_dft

Parameters

dft	Plot points: X in y2, Y in y3.
res	Results containing parameters of line.
first	First sample (starting from 0) used in linear fit.
last	last sample (starting from 0) used in linear fit.
mtitle	String for plot main title, or NULL.
xtitle	String for X axis title, or NULL.
ytitle	String for Y axis title, or NULL.
fname	Filename for plot data; existing file is renamed as *%. If extension is .dft, plot data (excluding lines) is written in DFT format with x values as separate columns before corresponding y values.

Definition at line 217 of file plotdata.c.

  {
  int n, ri, row, fi;
  char tmp[FILENAME_MAX], *cptr=NULL;
  FILE *fp;
  double maxPlotX=0, maxRegX=0, maxFitX=0, f;
 
 
  /* Check data */
  if(dft==NULL || dft->voiNr<1) return(1);
  if(res==NULL || res->voiNr!=dft->voiNr) return(1);
  if(first>last) return(1);
 
  /* Get filename extension to determine output type */
  cptr=strrchr(fname, '.');
  if(cptr!=NULL && strcasecmp(cptr, ".DFT")==0) {
    /* Write in DFT if required */
    return(plotdata_as_dft(dft, fname));
  }
 
  /* Check if file exists; backup, if necessary */
  backupExistingFile(fname, NULL, NULL);
 
  /* Search the largest plot x-value to be used as line end point */
  for(ri=0; ri<dft->voiNr; ri++) for(fi=0; fi<dft->frameNr; fi++)
    if(dft->voi[ri].y2[fi]>maxPlotX) maxPlotX=dft->voi[ri].y2[fi];
 
  /* Open output file */
  if((fp = fopen(fname, "w")) == NULL) return(3);
 
  /* Write XHTML doctype and head */
//  n=fprintf(fp, "<!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML 1.1//EN\" \"ht
//tp://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd\">\n");
  n=fprintf(fp, "<!DOCTYPE html PUBLIC \"-//W3C//DTD XHTML 1.1//EN\"");
  if(n<10) return(4);
  n=fprintf(fp, " \"https://www.w3.org/TR/xhtml11/DTD/xhtml11.dtd\">\n");
  if(n<10) return(4);
  n=fprintf(fp, "<html xmlns=\"https://www.w3.org/1999/xhtml\" xml:lang=\"en\">\n\n");
  if(n<20) return(4);
  /* Write XHTML header */
  n=fprintf(fp, "<head>\n"); if(n<6) return(4);
  fprintf(fp, "  <title>Graphical analysis plot</title>\n");
//  fprintf(fp, "  <meta http-equiv=\"content-type\" cont
//ent=\"text/html; charset=iso-8859-1\" />\n");
  fprintf(fp, "  <meta http-equiv=\"content-type\" content=\"text/html;");
  fprintf(fp, " charset=iso-8859-1\" />\n");
 
  fprintf(fp, "  <meta http-equiv=\"content-language\" content=\"en-gb\" />\n");
  fprintf(fp, "  <meta name=\"ProgId\" content=\"Excel.Sheet\" />\n");
/*
  fprintf(fp, "  <link rel=\"icon\" href=\"https://www.turkupetcentre.
net/favicon.ico\" type=\"image/x-icon\" />\n");
  fprintf(fp, "  <link rel=\"shortcut icon\" href=\"https://www.turkupetcentre.
net/favicon.ico\" type=\"image/x-icon\" />\n");
*/
  fprintf(fp, "  <link rel=\"icon\" href=\"https://www.turkupetcentre.net/");
  fprintf(fp, "favicon.ico\" type=\"image/x-icon\" />\n");
  fprintf(fp, "  <link rel=\"shortcut icon\" href=\"https://www.turkupet");
  fprintf(fp, "centre.net/favicon.ico\" type=\"image/x-icon\" />\n");
 
  fprintf(fp, "  <style type=\"text/css\">\n");
  fprintf(fp, "    thead {background-color:#999999; color:black;}\n");
/*  fprintf(fp, "    table {text-align:left; width:100%%; 
border-collapse:collapse; empty-cells:show;}\n");*/
  fprintf(fp, "    table {text-align:left; width:100%%;");
  fprintf(fp, " border-collapse:collapse; empty-cells:show;}\n");
 
  fprintf(fp, "    td {border:1px solid black;}\n");
  fprintf(fp, "    <!--table\n");
  fprintf(fp, "     {mso-displayed-decimal-separator:\"\\.\";\n");
  fprintf(fp, "     mso-displayed-thousand-separator:\" \";}\n");
  fprintf(fp, "    -->\n");
  fprintf(fp, "  </style>\n");
  fprintf(fp, "</head>\n");
 
  /* Start writing the body of the HTML file */
  fprintf(fp, "\n<body>\n");
 
  /* Start the div for tables */
  fprintf(fp, "\n<div id=\"tables\">\n");
 
  /* Write information on the graphical analysis */
  fprintf(fp, "<table>\n");
  fprintf(fp, "<tbody>\n");
  fprintf(fp, "<tr><th>Main title</th><th>%s</th></tr>\n", mtitle);
  fprintf(fp, "<tr><th>X title</th><th>%s</th></tr>\n", xtitle);
  fprintf(fp, "<tr><th>Y title</th><th>%s</th></tr>\n", ytitle);
  if(ctime_r_int(&res->time, tmp))
    fprintf(fp, "<tr><th>Date</th><th>%s</th></tr>\n", tmp);
  fprintf(fp, "</tbody>\n");
  fprintf(fp, "</table>\n");
 
  /* Write the plots, each to their own table */
  for(ri=0; ri<dft->voiNr; ri++) {
    /* Search the largest regional plot x-value to be used as line end points */
    for(fi=0, maxRegX=maxFitX=0; fi<dft->frameNr; fi++) {
      if(dft->voi[ri].y2[fi]>maxRegX) maxRegX=dft->voi[ri].y2[fi];
      if(fi>=first && fi<=last && dft->voi[ri].y2[fi]>maxFitX) maxFitX=dft->voi[ri].y2[fi];
    }
    /* Begin a new table */
    fprintf(fp, "<table>\n");
    /* Write the title row */
    fprintf(fp, "<thead>\n");
    fprintf(fp, "<tr><th>%s %s %s</th>", dft->voi[ri].voiname,
            dft->voi[ri].hemisphere, dft->voi[ri].place);
    fprintf(fp, "<th>symbol open</th><th>symbol filled</th><th>text</th>");
    fprintf(fp, "<th>X</th><th>line</th>");
    fprintf(fp, "</tr>\n");
    fprintf(fp, "</thead>\n");
    /* Write the plot rows */
    fprintf(fp, "<tbody>\n");
    row=0;
    for(fi=0; fi<(dft->frameNr>2?dft->frameNr:2); fi++)
      if(!isnan(dft->voi[ri].y2[fi]) && !isnan(dft->voi[ri].y3[fi]))
    {
        fprintf(fp, "<tr>");
        if(fi<dft->frameNr) {
          fprintf(fp, "<th>%g</th>", dft->voi[ri].y2[fi]); /* x-axis value */
          fprintf(fp, "<th>%g</th>", dft->voi[ri].y3[fi]); /* y-axis value */
        } else {
          fprintf(fp, "<th> </th>");
          fprintf(fp, "<th> </th>");
        }
      /* If included in the fit, y-axis value again */
      if(fi>=first && fi<=last)
        fprintf(fp, "<th>%g</th>", dft->voi[ri].y3[fi]);
      else fprintf(fp, "<th></th>");
      if(fi<dft->frameNr)
        fprintf(fp, "<th>%g</th>", dft->x[fi]); /* Frame time as text */
      else
        fprintf(fp, "<th> </th>");
      /* Line points */
      if(row==0) { /* line start */
        fprintf(fp, "<th>0</th>"); /* x-axis value */
        fprintf(fp, "<th>%g</th>", res->voi[ri].parameter[1]); /* y-axis value */
      } else if(row==1) { /* line point at the end of fitted range */
        fprintf(fp, "<th>%g</th>", maxFitX); /* x-axis value */
        f=maxFitX*res->voi[ri].parameter[0]+res->voi[ri].parameter[1];
        fprintf(fp, "<th>%g</th>", f); /* y-axis value */
      } else if(row==2) { /* line "mid" point at the end of regional data */
        fprintf(fp, "<th>%g</th>", maxRegX); /* x-axis value */
        f=maxRegX*res->voi[ri].parameter[0]+res->voi[ri].parameter[1];
        fprintf(fp, "<th>%g</th>", f); /* y-axis value */
      } else if(row==3) { /* line end point at the end of all plots */
        fprintf(fp, "<th>%g</th>", maxPlotX); /* x-axis value */
        f=maxPlotX*res->voi[ri].parameter[0]+res->voi[ri].parameter[1];
        fprintf(fp, "<th>%g</th>", f); /* y-axis value */
      }
      fprintf(fp, "</tr>\n");
      row++;
    }
    fprintf(fp, "</tbody>\n");
 
    /* End the data table */
    fprintf(fp, "</table>\n");
 
  } /* next region plot */
 
  /* End the div for tables */
  fprintf(fp, "</div>\n");
 
  /* Stop writing the body of the HTML file, and end the file */
  n=fprintf(fp, "</body></html>\n");
  if(n==0) return(4);
 
  /* Close file */
  fclose(fp);
 
  return(0);
}

◆ plotdata_as_dft()

int plotdata_as_dft	(	DFT *	dft,
		char *	fname )

Write plot data in DFT format with x values as separate columns before corresponding y values.

Returns: Returns 0 if successful, otherwise non-zero.

See also: plotdata

Parameters

dft	Plot points: X in y2, Y in y3.
fname	File name for plot data.

Definition at line 412 of file plotdata.c.

  {
  int ri, rj, fi, ret;
  DFT plot;
 
  /* Check input data */
  if(dft==NULL || dft->voiNr<1) return(1);
  /* Create the plot data */
  dftInit(&plot);
  ret=dftSetmem(&plot, dft->frameNr, 2*dft->voiNr); if(ret) return(ret);
  ret=dftCopymainhdr(dft, &plot); if(ret) {dftEmpty(&plot); return(ret);}
  for(ri=rj=0; ri<dft->voiNr; ri++) {
    /* x */
    strcpy(plot.voi[rj].voiname, "X");
    strcpy(plot.voi[rj].name, plot.voi[rj].voiname);
    for(fi=0; fi<dft->frameNr; fi++) plot.voi[rj].y[fi]=dft->voi[ri].y2[fi];
    rj++;
    /* y */
    dftCopyvoihdr(dft, ri, &plot, rj);
    for(fi=0; fi<dft->frameNr; fi++) plot.voi[rj].y[fi]=dft->voi[ri].y3[fi];
    rj++;
  }
  for(fi=0; fi<dft->frameNr; fi++) {
    plot.x[fi]=dft->x[fi]; plot.x1[fi]=dft->x1[fi]; plot.x2[fi]=dft->x2[fi];
  }
  plot.voiNr=2*dft->voiNr; plot.frameNr=dft->frameNr;
  /* Save plot data */
  strcpy(plot.comments, "");
  ret=dftWrite(&plot, fname);
  dftEmpty(&plot);
  return(ret);
}

Referenced by plotdata().

◆ sif2dft()

int sif2dft	(	SIF *	sif,
		DFT *	dft )

Convert SIF data to DFT data.

See also: sifAllocateWithIMG, dftAllocateWithIMG

Returns: Returns 0 if successful, otherwise >0.

Parameters

sif	Pointer to SIF struct containing data to be converted.
dft	Pointer to initiated but empty DFT struct.

Definition at line 350 of file misc_model.c.

  {
  int fi, ri, tacNr, ret;
 
  /* Check input */
  if(dft==NULL || sif==NULL) return 1;
  if(sif->frameNr<1) return 2;
  tacNr=sif->colNr-2; if(tacNr<=0) tacNr=1;
 
  /* Allocate memory */
  ret=dftSetmem(dft, sif->frameNr, tacNr);
  if(ret) return(100+ret);
  dft->voiNr=tacNr;
  dft->frameNr=sif->frameNr;
 
  // Set header contents
  dft->_type=DFT_FORMAT_STANDARD;
  dft->timetype=DFT_TIME_STARTEND;
  for(fi=0; fi<dft->frameNr; fi++) {
    dft->x[fi]=0.5*(sif->x1[fi]+sif->x2[fi]);
    dft->x1[fi]=sif->x1[fi];
    dft->x2[fi]=sif->x2[fi];
  }
  dft->timeunit=TUNIT_SEC;
  dft->isweight=0;
  strcpy(dft->unit, imgUnit(CUNIT_COUNTS));
  for(ri=0; ri<dft->voiNr; ri++) {
    /* Set ROI name */
    if(ri==0) sprintf(dft->voi[ri].voiname, "Prompt");
    else if(ri==1) sprintf(dft->voi[ri].voiname, "Random");
    else if(ri==2) sprintf(dft->voi[ri].voiname, "True");
    else if(ri==3) sprintf(dft->voi[ri].voiname, "Weight");
    else snprintf(dft->voi[ri].voiname, 7, "%06d", ri+1);
    strcpy(dft->voi[ri].name, dft->voi[ri].voiname);
    /* Copy TAC values */
    if(ri==0 && ri<sif->colNr-2) {
      for(fi=0; fi<dft->frameNr; fi++) dft->voi[ri].y[fi]=sif->prompts[fi];
    } else if(ri==1 && ri<sif->colNr-2) {
      for(fi=0; fi<dft->frameNr; fi++) dft->voi[ri].y[fi]=sif->randoms[fi];
    } else if(ri==2 && ri<sif->colNr-2) {
      for(fi=0; fi<dft->frameNr; fi++) dft->voi[ri].y[fi]=sif->trues[fi];
    } else if(ri==3 && ri<sif->colNr-2) {
      for(fi=0; fi<dft->frameNr; fi++) dft->voi[ri].y[fi]=sif->weights[fi];
    } else {
      for(fi=0; fi<dft->frameNr; fi++) dft->voi[ri].y[fi]=0.0;
    }
  }
  strcpy(dft->studynr, sif->studynr);
  strcpy(dft->isotope, hlCorrectIsotopeCode(sif->isotope_name));
 
  return 0;
}

◆ sifAllocateWithIMG()

int sifAllocateWithIMG	(	SIF *	sif,
		IMG *	img,
		int	doCounts,
		int	verbose )

Allocate memory for SIF based on information in IMG.

Frame times are set, 'counts' are optionally set to mean of voxel values times frame lengths, scaled to 0-1E7. Weights are not calculated.

See also: sif2dft, dftAllocateWithIMG

Returns: Returns 0 if successful, otherwise <>0.

Parameters

sif	Pointer to initiated SIF struct which will be allocated here.
img	Pointer to IMG struct from where necessary information is read.
doCounts	Set (1) or do not set (0) counts, based on sum of all image voxel values and frame lengths.
verbose	Verbose level; if zero, then only warnings are printed into stderr.

Definition at line 418 of file misc_model.c.

  {
  int fi, ret;
  float *cs;
  double f, mf;
 
  if(verbose>0) {printf("%s(*sif, *img, %d, ...)\n", __func__, doCounts); fflush(stdout);}
  /* Check the input data */
  if(sif==NULL || img==NULL) return 1;
  if(img->status!=IMG_STATUS_OCCUPIED) return 2;
  if(img->dimt<1) return 3;
  if(doCounts<0 || doCounts>1) return 4;
  
  /* Delete any previous SIF contents */
  sifEmpty(sif);
  
  /* Allocate memory for SIF data */
  ret=sifSetmem(sif, img->dimt); if(ret!=0) return(10+ret);
 
  /* Set SIF information */
  sif->version=1;
  sif->colNr=4;
  strcpy(sif->isotope_name, imgIsotope(img));
  strcpy(sif->studynr, img->studyNr);
  sif->scantime=img->scanStart;
  for(fi=0; fi<img->dimt; fi++) {sif->x1[fi]=img->start[fi]; sif->x2[fi]=img->end[fi];}
  
  if(doCounts==0) return 0; // that's it then
  
  /* Calculate average curve of all pixels */
  if(verbose>1) printf("calculate image average curve.\n");
  cs=(float*)malloc(img->dimt*sizeof(float));
  if(cs==NULL) {sifEmpty(sif); return(20);}
  ret=imgAverageTAC(img, cs); if(ret!=0) {sifEmpty(sif); return(20+ret);}
  /* Multiply average curve with frame durations, and get the max */
  for(fi=0, mf=0.0; fi<sif->frameNr; fi++) {
    f=sif->x2[fi]-sif->x1[fi]; if(f<=0.1) f=0.1;
    cs[fi]*=f; if(cs[fi]>mf) mf=cs[fi];
  }
  /* Put scaled counts into SIF */
  if(mf>0.0) f=1.0E+007/mf; else f=1.0;
  for(fi=0, mf=0.0; fi<sif->frameNr; fi++) {
    sif->prompts[fi]=sif->trues[fi]=cs[fi]*f;
    sif->randoms[fi]=0.0;
  }
  free(cs);
  if(verbose>2) sifPrint(sif);
  return 0;
}

Referenced by imgSetWeights().

Functions

Detailed Description

Function Documentation

◆ bf_srtm()

◆ bfIrr2TCM()

◆ bfRadiowater()

◆ check_times_dft_vs_dft()

◆ check_times_dft_vs_img()

◆ clusterTACs()

◆ copy_times_from_img_to_dft()

◆ cunit_check_dft_vs_img()

◆ dft_end_line()

◆ dft_ln()

◆ dftAllocateWithIMG()

◆ dftAutointerpolate()

◆ dftDerivative()

◆ dftDerivative_old()

◆ dftDivideFrames()

◆ dftDoubleFrames()

◆ dftEndtime()

◆ dftInterpolate()

◆ dftInterpolateCheckEnd()

◆ dftInterpolateCheckStart()

◆ dftInterpolateForIMG()

◆ dftInterpolateInto()

◆ dftMatchTimeunits()

◆ dftReadinput()

◆ dftReadModelingData()

◆ dftReadReference()

◆ dftRobustMinMaxTAC()

◆ dftTimeIntegral()

◆ dftVerifyPeak()

◆ dftWeightByFreq()

◆ dftWSampleNr()

◆ extrapolate_monoexp()

◆ fittime_from_dft()

◆ fittime_from_img()

◆ getActualSamplenr()

◆ img_k1_using_ki()

◆ img_logan()

◆ img_patlak()

◆ imgEndtime()

◆ imgReadModelingData()

◆ imgSetWeights()

◆ imgTimeIntegral()

◆ mrl_between_tacs()

◆ noiseSD4Simulation()

◆ noiseSD4SimulationFromDFT()

◆ plot_fit_svg()

◆ plot_fitrange_svg()

◆ plot_svg()

◆ plotdata()

◆ plotdata_as_dft()

◆ sif2dft()

◆ sifAllocateWithIMG()