tpctacmod.h File Reference

Header file for libtpctacmod. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpcli.h"
#include "tpctac.h"
#include "tpcmodels.h"
#include "tpcpar.h"
#include "tpcfunc.h"

Go to the source code of this file.

Data Structures
struct	DELAYCMFITDATA

Functions
int	tacIntegrate (TAC *inp, TAC *out, TPCSTATUS *status)
	Integrate TACs from one TAC structure into a new TAC structure.
int	tacIntegrateFE (TAC *inp, TAC *out, TAC *out2, TPCSTATUS *status)
	Integrate TACs from one TAC structure into a new TAC structure. Integrals are calculated at frame end times.
int	tacInterpolate (TAC *inp, TAC *xinp, TAC *tac, TAC *itac, TAC *iitac, TPCSTATUS *status)
	Interpolate and/or integrate TACs from one TAC structure into a new TAC structure, using sample times from another TAC structure which is not changed.
int	tacInterpolateInto (TAC *inp, TAC *tac, TAC *itac, TAC *iitac, TPCSTATUS *status)
	Add TACs from one TAC structure into another TAC structure, interpolating the input TACs and allocating space if necessary.
double	tacAUC (TAC *tac, int ti, double t1, double t2, TPCSTATUS *status)
	Calculates TAC AUC from t1 to t2.
int	tacInput2sim (TAC *itac, TAC *ttac, TAC *stac, TPCSTATUS *status)
	Modify input TAC based on tissue TAC, for use in simulations.
int	tacVb (TAC *ttac, const int i, TAC *btac, double Vb, const int simVb, const int petVolume, TPCSTATUS *status)
	Correct TTACs for vascular blood, or simulate its effect.
int	tacInterpolateToEqualLengthFrames (TAC *inp, double minfdur, double maxfdur, TAC *tac, TPCSTATUS *status)
int	tacFittime (TAC *d, double *startTime, double *endTime, int *first, int *last, TPCSTATUS *status)
int	tacReadModelingData (const char *tissuefile, const char *inputfile1, const char *inputfile2, const char *inputfile3, double *fitdur, int cutInput, int *fitSampleNr, TAC *tis, TAC *inp, TPCSTATUS *status)
	Read tissue and input data for modelling.
int	tacReadReference (TAC *tis, const char *reference, TAC *ref, int *refIndex, TPCSTATUS *status)
	Read reference tissue TAC.
int	tacReadModelingInput (const char *inputfile1, const char *inputfile2, const char *inputfile3, TAC *inp, TPCSTATUS *status)
	Read arterial input data for modelling.
int	tacPlotFitSVG (TAC *tac1, TAC *tac2, const char *main_title, const double x1, const double x2, const double y1, const double y2, const char *fname, TPCSTATUS *status)
int	tacPlotLineSVG (TAC *tac, const char *main_title, const double x1, const double x2, const double y1, const double y2, const char *fname, TPCSTATUS *status)
int	mtgaPlotSVG (MTAC *mtac, const char *main_title, const char *fname, TPCSTATUS *status)
int	tacPlotHistogramSVG (TAC *d, const char *main_title, const double x1, const double x2, const double y1, const double y2, const char *fname, TPCSTATUS *status)
int	tacAllocateWithPAR (TAC *tac, PAR *par, int sampleNr, TPCSTATUS *status)
	Allocate TAC based on data in PAR.
int	parAllocateWithTAC (PAR *par, TAC *tac, int parNr, TPCSTATUS *status)
	Allocate PAR based on data in TAC.
int	tacToPAR (TAC *tac, PAR *par, TPCSTATUS *status)
	Copy the contents of TAC struct into PAR struct.
int	tacDelay (TAC *tac, double dt, int ti, TPCSTATUS *status)
	Move TAC y values (concentrations) in time, keeping sample times (x values) intact.
void	initDelayCMFitData (DELAYCMFITDATA *d)
	Before first use, initiate the data structure for delay time estimation.
void	freeDelayCMFitData (DELAYCMFITDATA *d)
	After last use, free memory in the data structure for delay time estimation.
int	tacDelayCMFit (TAC *btac, TAC *ttac, int ci, double dtmin, double dtmax, double fitend, double dtstep, double *dt, int mode, int model, DELAYCMFITDATA *tdata, TPCSTATUS *status)
	Fit time delay between PET tissue and plasma or blood curve.
int	mfCreateTAC (PAR *par, double endx, double dx, TAC *tac, TPCSTATUS *status)
	Make TAC(s) based on mathematical functions in PAR format.

Detailed Description

Header file for libtpctacmod.

libtpctacmod contains functions needed in reading and processing PET TAC data for modelling and fitting.

Author

Vesa Oikonen

Definition in file tpctacmod.h.

Function Documentation

freeDelayCMFitData()

void freeDelayCMFitData ( DELAYCMFITDATA * d )

extern

After last use, free memory in the data structure for delay time estimation.

See also

tacDelayCMFit, initDelayCMFitData

Parameters

d	Pointer to the data structure to be initialized.

Definition at line 150 of file delay.c.

  {
  if(d==NULL) return;
  d->iNr=d->tNr=d->moveNr=0;
  tacFree(&d->dtac); tacFree(&d->ditac); tacFree(&d->diitac); tacFree(&d->ddtac); tacFree(&d->dddtac);
  d->mode=0; d->model=0;
  d->llsq_n=d->llsq_m=0;
  free(d->llsq_a); free(d->llsq_mat);
  free(d->llsq_b); free(d->llsq_x); free(d->llsq_wp); free(d->llsq_zz); 
  free(d->llsq_i);
  return;
}

Referenced by tacDelayCMFit().

initDelayCMFitData()

void initDelayCMFitData ( DELAYCMFITDATA * d )

extern

Before first use, initiate the data structure for delay time estimation.

See also

tacDelayCMFit, freeDelayCMFitData

Parameters

d	Pointer to the data structure to be initialized.

Definition at line 131 of file delay.c.

  {
  if(d==NULL) return;
  d->iNr=d->tNr=d->moveNr=0;
  tacInit(&d->dtac); tacInit(&d->ditac); tacInit(&d->diitac); tacInit(&d->ddtac); tacInit(&d->dddtac);
  d->mode=0; d->model=0;
  d->llsq_n=d->llsq_m=0;
  d->llsq_mat=NULL;
  d->llsq_a=NULL;
  d->llsq_b=d->llsq_x=d->llsq_wp=d->llsq_zz=NULL; 
  d->llsq_i=NULL;
  return;
}

Referenced by tacDelayCMFit().

mfCreateTAC()

int mfCreateTAC ( PAR * par, double endx, double dx, TAC * tac, TPCSTATUS * status )

extern

Make TAC(s) based on mathematical functions in PAR format.

See also

mfEvalY, parExampleTTACs

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

par	Pointer to PAR structure containing functions for TACs.
endx	End time (x) of created TAC data. TAC will always start at zero.
dx	Step size (delta x).
tac	Pointer to initiated TAC structure. Any previous contents are deleted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 26 of file mftac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>2) {
    printf("%s(par, %g, %g, tac, status)\n", __func__, endx, dx);
    fflush(stdout);
  }
 
  /* Check provided data */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
  if(par==NULL || tac==NULL || !(endx>0.0) || !(dx>0.0)) return(status->error);
  if(par->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(status->error);
  }
 
  /* Make sure that TAC is empty */
  tacFree(tac);
 
  /* Allocate memory for TAC(s) */
  int xNr=1+endx/dx;  
  int ret=tacAllocate(tac, xNr, par->tacNr); 
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(status->error);
  }
  tac->sampleNr=xNr;
  tac->tacNr=par->tacNr;
 
  /* Fill x (times) */
  tac->isframe=0;
  for(int i=0; i<xNr; i++) tac->x[i]=(double)i*dx;
 
  /* Calculate curves */
  for(int ci=0; ci<par->tacNr; ci++) {
    strcpy(tac->c[ci].name, par->r[ci].name);
    if(verbose>6)
      printf("'%s' model=%u code=%s\n", par->r[ci].name, par->r[ci].model, modelCode(par->r[ci].model));
    if(mfEvalY(modelCode(par->r[ci].model), modelParNr(par->r[ci].model), par->r[ci].p, 
               tac->sampleNr, tac->x, tac->c[ci].y, verbose-3))
    {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
      tacFree(tac); return(status->error);
    }
  }
 
  /* Try to find units */
  int i;
  i=iftFindKey(&par->h, "unit", 0); if(i<0) i=iftFindKey(&par->h, "calibration_unit", 0);
  if(i>=0) tac->cunit=unitIdentify(par->h.item[i].value);
  i=iftFindKey(&par->h, "timeunit", 0); if(i<0) i=iftFindKey(&par->h, "time_unit", 0);
  if(i>=0) tac->tunit=unitIdentify(par->h.item[i].value);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

mtgaPlotSVG()

int mtgaPlotSVG ( MTAC * mtac, const char * main_title, const char * fname, TPCSTATUS * status )

extern

MTGA plots in SVG format.

Returns

enum tpcerror (TPCERROR_OK when successful).

See also

tacPlotLineSVG, tacPlotFitSVG

Author

Vesa Oikonen

Parameters

mtac	Pointer to MTAC structure, which contains plot data of tissue regions, each in its own TAC structure. TAC must contain the plot x axis values, and three y columns: points included in regression, points not included, and points representing the fitted line.
main_title	String for plot main title, or "".
fname	SVG file name.
status	Pointer to status data; enter NULL if not needed.

Definition at line 489 of file tacfitplot.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(MTAC, '%s', '%s')\n", __func__, main_title, fname); fflush(stdout);}
  /* Check data */
  if(mtac==NULL || mtac->nr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(fname==NULL || strnlen(fname, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return TPCERROR_INVALID_FILENAME;
  }
  if(status!=NULL && status->forgiving==0) {
    // any strict tests in here
  }
  if(verbose>1) {printf("%d plot(s)\n", mtac->nr); fflush(stdout);}
 
  /* Do we show plot legends? */
  int is_label=0; if(mtac->nr>1) is_label=1;
 
  /* Get the plot x and y range */
  double minx=nan(""), maxx=nan(""), miny=nan(""), maxy=nan(""), maxly=nan("");
  for(int i=0; i<mtac->nr; i++) {
    if(mtac->tac[i].tacNr<3) continue;
    double x1, x2, y1, y2;
    if(tacSampleXRange(&mtac->tac[i], &x1, &x2)!=0) continue;
    if(isfinite(x1) && !(x1>minx)) minx=x1;
    if(isfinite(x2) && !(x2<maxx)) maxx=x2;
    if(tacYRange(mtac->tac+i, -1, &y1, &y2, NULL, NULL, NULL, NULL)!=0) continue;
    if(isfinite(y1) && !(y1>miny)) miny=y1;
    if(isfinite(y2) && !(y2<maxy)) maxy=y2;
    /* Get the y maximum from fitted line alone */
    if(tacYRange(mtac->tac+i, 2, &y1, &y2, NULL, NULL, NULL, NULL)!=0) continue;
    if(isfinite(y2) && !(y2<maxly)) maxly=y2;
  }
  if(verbose>1) {
    printf("x-range := %g - %g\n", minx, maxx);
    printf("y-range := %g - %g (%g for the line)\n", miny, maxy, maxly);
  }
  if(!(maxx>minx) || !(maxy>miny)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  /* Limit the y range to maximal line value plus 1/40 of y range */
  {
    double yrange=maxy-miny;
    double ylim=maxly+0.025*yrange;
    if(maxy>ylim) maxy=ylim;
    if(verbose>2) printf("  yrange := %g\n  ylim := %g\n  maxy := %g\n", yrange, ylim, maxy);
  }
 
  /* Start setting viewports */
  struct svg_viewports viewports; svg_init_viewports(&viewports);
  viewports.x.min=minx; viewports.x.max=maxx;
  viewports.y.min=miny; viewports.y.max=maxy;
  viewports.label_area_viewport.is=is_label; // needed for x axis ticks
 
  /* Calculate the axis ticks */
  if(svg_calculate_axes(&viewports, verbose-3)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Set x and y axis titles based on activity and time units */
  char x_title[64], y_title[64];
  strcpy(x_title, ""); strcpy(y_title, "");
 
  /* Set the plot window and window area sizes; ticks should be set before this */
  if(svg_define_viewports(0, 0, strlen(main_title), strlen(y_title),
    strlen(x_title), is_label, &viewports, verbose-3)) 
  {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Initiate graphics file */
  FILE *fp=svg_initiate(fname, 0, 0, &viewports, NULL, verbose-3);
  if(fp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
 
  /* Put the graph titles into their own viewports */
  if(svg_create_main_title(fp, main_title, "", &viewports, NULL,verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  if(svg_create_yaxis_title(fp, y_title, &viewports, NULL, verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
  if(svg_create_xaxis_title(fp, x_title, &viewports, NULL, verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Put the plot into its own viewport */
  if(svg_start_plot_viewport(fp, &viewports, NULL, verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /*  Start coordinate area viewport */
  if(svg_start_coordinate_viewport(fp, &viewports, NULL, verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Write plot axes */
  if(svg_write_axes(fp, &viewports, NULL, verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
 
  /*
   *  Draw the plots
   */
  SVG_LEGENDS legends; svg_init_legends(&legends);
 
  int max_color_nr=0; while(svgColorName(max_color_nr)!=NULL) max_color_nr++;
  if(verbose>3) {printf("max_color_nr := %d\n", max_color_nr); fflush(stdout);}
  int color_nr; if(mtac->nr==1) color_nr=0; else color_nr=1;
  int max_symbol_nr=0; while(svgSymbolName(max_symbol_nr)!=NULL) max_symbol_nr++;
  if(verbose>3) printf("max_symbol_nr := %d\n", max_symbol_nr);
  int n=0, symbol_nr=0;
  for(int ri=0; ri<mtac->nr; ri++) {
    if(mtac->tac[ri].tacNr<3) continue;
    char tac_id[32], tac_title[64];
    sprintf(tac_id, "plot_%d", n);
    strcpy(tac_title, mtac->tac[ri].c[0].name);
 
 
    /* Draw the measured points */
    int ret=svg_write_tac(fp, &viewports, 2, tac_id, tac_title,
            mtac->tac[ri].x, mtac->tac[ri].c[1].y, mtac->tac[ri].sampleNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLOPEN,
            NULL, verbose-3);
    if(!ret) ret=svg_write_tac(fp, &viewports, 2, tac_id, tac_title,
            mtac->tac[ri].x, mtac->tac[ri].c[0].y, mtac->tac[ri].sampleNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLFILLED,
            NULL, verbose-3);
    /* Draw the regression line */
    if(!ret) ret=svg_write_tac(fp, &viewports, 1, tac_id, tac_title,
            mtac->tac[ri].x, mtac->tac[ri].c[2].y, mtac->tac[ri].sampleNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLOPEN,
            NULL, verbose-3);
    if(ret) {
      svg_legend_empty(&legends); fclose(fp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
 
 
    /* Set legend too, if requested */
    if(is_label!=0) {
      unsigned int ci, si; // remainder works only if 2nd operator > 0
      if(max_color_nr<1) ci=0; else ci=color_nr % max_color_nr;
      if(max_symbol_nr<1) si=0; else si=symbol_nr % max_symbol_nr;
      svg_legend_add(&legends, 0, si, SYMBOLFILLED, ci, tac_title);
    }
    /* Prepare for the next plot */
    color_nr++; n++;
    if(color_nr==max_color_nr) color_nr=0;
  }
 
 
 
  /* Close the coordinate viewport */
  if(svg_end_coordinate_viewport(fp, NULL, verbose-3)) {
    fclose(fp); svg_legend_empty(&legends);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Write the axis ticks */
  if(svg_write_xticks(fp, &viewports, NULL, verbose-3)!=0 ||
     svg_write_yticks(fp, &viewports, NULL, verbose-3)!=0) 
  {
    fclose(fp); svg_legend_empty(&legends);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Close the plot viewport */
  if(svg_end_plot_viewport(fp, NULL, verbose-3)) {
    fclose(fp); svg_legend_empty(&legends);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  /* Make the plot legends into their own viewport */
  if(viewports.label_area_viewport.is!=0) {
    if(verbose>2) {printf("creating plot legends\n"); fflush(stdout);}
    int ret=svg_create_legends(fp, &viewports, &legends, NULL, verbose-3);
    if(ret) {
      fclose(fp); svg_legend_empty(&legends);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
      return TPCERROR_CANNOT_WRITE;
    }
  }
  svg_legend_empty(&legends);
 
  /* Close the SVG file */
  if(svg_close(fp, NULL, verbose-3)) {
    fclose(fp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    return TPCERROR_CANNOT_WRITE;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

parAllocateWithTAC()

int parAllocateWithTAC ( PAR * par, TAC * tac, int parNr, TPCSTATUS * status )

extern

Allocate PAR based on data in TAC.

Region number, names, etc are copied from given TAC.

See also

parAllocate, parAllocateMore, parInit, tacToPar

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

par	Pointer to initiated PAR struct; any old contents are deleted. PAR tacNr and parNr will be set, although contents will be empty.

Precondition

Struct must be initiated with parInit.

Parameters

tac	Pointer to TAC struct which contains at least the TAC number and names.
parNr	Nr of parameters to allocate.
status	Pointer to status data; enter NULL if not needed.

Definition at line 90 of file partac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(par, tac, %d)\n", __func__, parNr);
 
  /* Check the input */
  if(tac==NULL || par==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(parNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  int i, ret;
 
  /* Allocate memory */
  if(verbose>2) printf("allocate memory for PAR\n");
  ret=parAllocate(par, parNr, tac->tacNr);
  if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return ret;}
  par->tacNr=tac->tacNr;
  par->parNr=parNr;
 
  /* Copy TAC names */
  if(verbose>3) printf("copy TAC names\n");
  for(int i=0; i<par->tacNr; i++) strcpy(par->r[i].name, tac->c[i].name);
 
  /* Copy header information */
  if(verbose>3) printf("copy header\n");
 
  i=iftFindKey(&tac->h, "studynr", 0);
  if(i<0) i=iftFindKey(&tac->h, "study_number", 0);
  if(i>=0) iftPut(&par->h, tac->h.item[i].key, tac->h.item[i].value, 0, NULL);
 
  i=iftFindKey(&tac->h, "injection_time", 0);
  if(i<0) i=iftFindKey(&tac->h, "injection time", 0);
  if(i>=0) iftPut(&par->h, tac->h.item[i].key, tac->h.item[i].value, 0, NULL);
 
  i=iftFindKey(&tac->h, "isotope", 0);
  if(i>=0) iftPut(&par->h, tac->h.item[i].key, tac->h.item[i].value, 0, NULL);
 
  i=iftFindKey(&tac->h, "radiopharmaceutical", 0);
  if(i>=0) iftPut(&par->h, tac->h.item[i].key, tac->h.item[i].value, 0, NULL);
 
  i=iftFindKey(&tac->h, "scan_start_time", 0);
  if(i>=0) iftPut(&par->h, tac->h.item[i].key, tac->h.item[i].value, 0, NULL);
 
  if(tacIsWeighted(tac)) iftPut(&par->h, "weighting", "yes", 0, NULL);
  else iftPut(&par->h, "weighting", "no", 0, NULL);
 
  if(tac->tunit!=UNIT_UNKNOWN)
    iftPut(&par->h, "timeunit", unitName(tac->tunit), 0, NULL);
  if(tac->cunit!=UNIT_UNKNOWN)
    iftPut(&par->h, "unit", unitName(tac->cunit), 0, NULL);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by tacToPAR().

tacAllocateWithPAR()

int tacAllocateWithPAR ( TAC * tac, PAR * par, int sampleNr, TPCSTATUS * status )

extern

Allocate TAC based on data in PAR.

TAC region number, names, and units are copied from given PAR.

See also

tacAllocate, tacAllocateMore, tacDuplicate, tacAllocateWithIMG

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

tac	Pointer to initiated TAC structure; any old contents are deleted. TAC sampleNr and tacNr will be set, although contents will be empty.

Precondition

Structure must be initiated with tacInit.

Parameters

par	Pointer to PAR structure which contains at least the TAC number and names.
sampleNr	Nr of samples (time frames) to allocate.
status	Pointer to status data; enter NULL if not needed.

Definition at line 28 of file partac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(tac, par, %d)\n", __func__, sampleNr);
 
  /* Check the input */
  if(tac==NULL || par==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(sampleNr<1 || par->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  int ret, i;
 
  /* Allocate */
  ret=tacAllocate(tac, sampleNr, par->tacNr);
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  if(ret!=TPCERROR_OK) return(ret);
  tac->sampleNr=sampleNr; tac->tacNr=par->tacNr;
 
  /* Copy TAC names */
  for(int i=0; i<par->tacNr; i++) strcpy(tac->c[i].name, par->r[i].name);
 
  /* Try to find the data units */
  i=iftFindKey(&par->h, "unit", 0);
  if(i<0) i=iftFindKey(&par->h, "calibration_unit", 0);
  if(i>0) tac->cunit=unitIdentify(par->h.item[i].value);
  i=iftFindKey(&par->h, "timeunit", 0);
  if(i<0) i=iftFindKey(&par->h, "time_unit", 0);
  if(i>0) tac->tunit=unitIdentify(par->h.item[i].value);
  /* Try to find study number */
  i=iftFindKey(&par->h, "studynr", 0);
  if(i<0) i=iftFindKey(&par->h, "study", 0);
  if(i>=0) {
    iftPut(&tac->h, par->h.item[i].key, par->h.item[i].value, 1, NULL);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacAUC()

double tacAUC ( TAC * tac, int ti, double t1, double t2, TPCSTATUS * status )

extern

Calculates TAC AUC from t1 to t2.

See also

tacInterpolate, liInterpolate, liIntegrate

Returns

AUC, or NaN in case of an error.

Parameters

tac	Pointer to TAC structure.

Precondition

Data must be sorted by increasing x. No gaps or overlap between frames allowed.

See also

tacSortByTime, tacSetXContiguous

Parameters

ti	Index [0..tacNr-1] of the TAC to calculate AUC from.
t1	AUC start time.
t2	AUC end time.
status	Pointer to status data; enter NULL if not needed.

Definition at line 486 of file litac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(tac, %d, %g, %g)\n", __func__, ti, t1, t2);
  /* Check the function input */
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(nan(""));
  }
  if(ti<0 || ti>=tac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(nan(""));
  }
  if(t1>=t2) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return(nan(""));
  }
 
  /* Get and check the x range */
  double xmin=nan(""), xmax=nan("");
  if(tacXRange(tac, &xmin, &xmax)!=0 || t1>=xmax || t2<=xmin) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return(nan(""));
  }
 
  /* Make temporary x and yi array */
  double tx[2]; tx[0]=t1; tx[1]=t2;
  double tyi[2];
 
  /* If TAC contains just frame mid times, then simple dot-to-dot integration */
  if(tac->isframe==0) {
    if(liInterpolate(tac->x, tac->c[ti].y, tac->sampleNr, tx, NULL, tyi, NULL, 2, 3, 1, verbose-10))
    {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(nan(""));
    }
    /* That's it then */
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(tyi[1]-tyi[0]);
  }
 
  /* We have frame start and end times */
  /* Convert the data into stepwise dot-to-dot data */
  TAC stac; tacInit(&stac);
  {
    int ret=tacFramesToSteps(tac, &stac, NULL);
    if(ret!=TPCERROR_OK) {
      tacFree(&stac);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(nan(""));
    }
  }
  /* Then simple dot-to-dot integration */
  if(liInterpolate(stac.x, stac.c[ti].y, stac.sampleNr, tx, NULL, tyi, NULL, 2, 3, 1, verbose-10))
  {
    tacFree(&stac);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(nan(""));
  }
  tacFree(&stac);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(tyi[1]-tyi[0]);
}

tacDelay()

int tacDelay ( TAC * tac, double dt, int ti, TPCSTATUS * status )

extern

Move TAC y values (concentrations) in time, keeping sample times (x values) intact.

If TAC has frame start and end times, this function tries to retain the AUC with stepwise frame interpretation. If you wish to prevent this, then set tac->isframe=0 before calling this function and after it set back tac->isframe=1.

See also

tacDelayCMFit, tacInterpolate, tacFramesToSteps, simDispersion, liInterpolate, tacVb

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

tac	Pointer to TAC structure.

Precondition

Data must be sorted by increasing x. No gaps or overlap between frames allowed.

See also

tacSortByTime, tacSetXContiguous

Postcondition

Data y values are modified, x values will stay the same.

Parameters

dt	Time (in TAC tunit's) to move the TAC forward (negative time moves TAC backward).
ti	Index [0..tacNr-1] of the TAC to move, or <0 to move all TACs.
status	Pointer to status data; enter NULL if not needed.

Definition at line 29 of file delay.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(tac, %g, %d)\n", __func__, dt, ti);
  /* Check the function input */
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(ti>=tac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(fabs(dt)<1.0E-12) {
    if(verbose>1) printf("  requested delay time is zero; nothing done.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
  /* Check that input does not have any missing values */
  if(tacXNaNs(tac)>0 || tacYNaNs(tac, ti)>0) {
    if(verbose>1) printf("  data contains NaN(s).\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
 
  /* Get and check the x range */
  double xmin=nan(""), xmax=nan("");
  if(tacXRange(tac, &xmin, &xmax)!=0 || dt>=xmax || dt<=-xmax) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  if(verbose>2) printf("  xrange: %g - %g\n", xmin, xmax);
 
 
  /* If TAC contains just frame mid times, then simple dot-to-dot interpolation */
  if(tac->isframe==0) {
    /* Make temporary x array */
    double tx[tac->sampleNr];
    for(int i=0; i<tac->sampleNr; i++) tx[i]=tac->x[i]+dt;
    /* One TAC at a time */
    for(int j=0; j<tac->tacNr; j++) if(ti<0 || ti==j) {
      /* Make temporary y array */
      double ty[tac->sampleNr];
      for(int i=0; i<tac->sampleNr; i++) ty[i]=tac->c[j].y[i];
      /* Interpolate temporary data to original times */
      if(liInterpolate(tx, ty, tac->sampleNr, tac->x, tac->c[j].y, NULL, NULL, tac->sampleNr, 0, 1, 
                       verbose-10)!=0)
      {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
        return TPCERROR_INVALID_VALUE;
      }
    }
    /* That's it then */
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
 
  /* We have frame start and end times */
  /* Convert the data into stepwise dot-to-dot data */
  TAC stac; tacInit(&stac);
  {
    int ret=tacFramesToSteps(tac, &stac, NULL);
    if(ret!=TPCERROR_OK) {
      tacFree(&stac);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
  }
  /* Change times of the stepwise data */
  for(int i=0; i<stac.sampleNr; i++) stac.x[i]+=dt;
  /* Interpolate stepwise data into original frame times, one TAC at a time */
  for(int j=0; j<tac->tacNr; j++) if(ti<0 || ti==j) {
    if(liInterpolateForPET(stac.x, stac.c[j].y, stac.sampleNr, tac->x1, tac->x2, tac->c[j].y, 
                           NULL, NULL, tac->sampleNr, 0, 1, verbose-10))
    {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
  }
  tacFree(&stac);
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacDelayCMFit()

int tacDelayCMFit ( TAC * btac, TAC * ttac, int ci, double dtmin, double dtmax, double fitend, double dtstep, double * dt, int mode, int model, DELAYCMFITDATA * tdata, TPCSTATUS * status )

extern

Fit time delay between PET tissue and plasma or blood curve.

Initial part of tissue curve is fitted with a range of input curves moved in time, using linearised one or two-tissue compartment model with blood volume term. The time that provides the best fit is selected. It is assumed that tissue heterogeneity, dispersion of the blood curve, metabolites, and variable plasma-to-blood ratio will only minimally impact the time delay estimate, since only the initial part of the data is fitted, especially if reversible two-tissue model (R2tCM) is used.

See also

tacDelay, simDispersion

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

btac	Pointer to input data (blood or plasma curve). Units of input and tissue TAC must match. All but first curve in TAC are ignored. Enter NULL in subsequent runs
ttac	Pointer to tissue data. Units of input and tissue TAC must match. All but first curve in TAC are ignored.
ci	Index [0..ttac->tacNr-1] of tissue tac to use.
dtmin	Minimum delay to be tested, for example -10 s, in the same units as in TACs. Not used in subsequent runs with tdata.
dtmax	Maximum delay to be tested, for example +50 s, in the same units as in TACs. Not used in subsequent runs with tdata.
fitend	Fitting uses data from 0 to end time, for example 120 s, in the same units as in TACs. In case of short PET scan, end time may need to be reduced so that input TAC with negative delay extends to the fit end. Not used in subsequent runs with tdata.
dtstep	The step length for moving the input curve, for example 1 s, in the same units as in TACs. Not used in subsequent runs with tdata.
dt	Pointer for the time delay estimate. Positive delay time means that tissue curve is delayed as compared to the input curve, and vice versa. Thus, input curve needs to be moved by the delay time to match the tissue curve.
mode	Mode for multilinear equations: 0=integrals and 2nd integrals, 1=derivates and 2nd derivates (experimental), 2=integrals and derivates (experimental).
model	Model: 0=R1TCM, 1=I2TCM, 2=R2TCM.
tdata	Pointer to temporary data structure, which quickens subsequent runs with other tissue TACs; enter NULL if not needed.

See also

initDelayCMFitData, freeDelayCMFitData

Parameters

status

Pointer to status data; enter NULL if not needed.

Definition at line 179 of file delay.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>2) {
    printf("%s(inp, tis, %d, %g, %g, %g, %g, %d, ...)\n", __func__, ci, dtmin, dtmax, fitend, dtstep, mode);
    fflush(stdout);
  }
 
 
  /* Check provided data */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
  if(ttac==NULL || (btac==NULL && tdata==NULL)) return(TPCERROR_FAIL);
  if(btac!=NULL && (btac->tacNr<1 || btac->sampleNr<4)) return(TPCERROR_FAIL);
  if(ci<0 || ci>=ttac->tacNr || ttac->sampleNr<4) return(TPCERROR_FAIL);
  if(model<0 || model>2) return(TPCERROR_FAIL);
  if(dt!=NULL) *dt=nan("");
 
  /* Initiate (possibly reused) data structure */
  DELAYCMFITDATA *b, ltdata; initDelayCMFitData(&ltdata);
  if(tdata!=NULL) {
    b=tdata;
    if(verbose>3) {
      if(b->moveNr==0) printf("  given temp data will be filled\n"); 
      else printf("  filled temp data will be used\n");
      fflush(stdout);
    }
  } else {
    b=&ltdata;
    if(verbose>4) {printf("  local temp data created\n"); fflush(stdout);}
  }
 
  /* Check and set temp data when run first time with the given TTAC */
  if(btac!=NULL || b->moveNr==0) {
    /* Delay range and steps*/
    if(!(dtmax>dtmin) || !(dtstep>0.001)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    b->dtmin=dtmin; b->dtmax=dtmax; b->dtstep=dtstep;
    /* How many delayed curves needed */
    b->moveNr=1+(b->dtmax-b->dtmin)/b->dtstep;
    if(verbose>3) printf("  moveNr := %d\n", b->moveNr);
    if(b->moveNr<5 || b->moveNr>5000) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    /* Get x range */
    if(tacXRange(ttac, &b->ttac_xmin, &b->ttac_xmax)!=0 || tacXRange(btac, &b->btac_xmin, &b->btac_xmax)!=0) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    if(verbose>3) printf("  TTAC range := %g - %g\n  BTAC range := %g - %g\n",
                         b->ttac_xmin, b->btac_xmax, b->ttac_xmin, b->btac_xmax);
    /* Check that time units do match */
    if(btac->tunit!=ttac->tunit) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return(TPCERROR_INVALID_X);
    }
    /* Check that input curve is long enough to fit duration */
    if(!(fitend>0.0)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    b->fitend=fitend;
    if(b->fitend>b->ttac_xmax) b->fitend=b->ttac_xmax;
    if(b->fitend>b->btac_xmax) b->fitend=b->btac_xmax;
    if(b->dtmin<0.0 && b->btac_xmax<(b->fitend+b->dtmin))
      b->fitend+=b->dtmin;
    if(verbose>3) printf("fitend=%g -> %g\n", fitend, b->fitend);
    b->iNr=0; for(unsigned int i=0; i<(unsigned int)btac->sampleNr; i++) if(btac->x[i]<=b->fitend) b->iNr++; else break;
    b->tNr=0; for(unsigned int i=0; i<(unsigned int)ttac->sampleNr; i++) if(ttac->x[i]<=b->fitend) b->tNr++; else break;
    //printf("iNr=%u tNr=%u\n", b->iNr, b->tNr);
    if(b->iNr<5 || b->tNr<5) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
 
    /* Further, when run the first time, make input TACs moved in time, interpolated and integrated 
       to tissue times */
    if(verbose>5) {printf("  preparing delayed input data\n"); fflush(stdout);}
    b->mode=mode;
    b->model=model;
    int ret;
    ret=tacAllocate(&b->dtac, b->tNr, b->moveNr+1); // additional place for tissue curve
    if(ret!=TPCERROR_OK) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(TPCERROR_OUT_OF_MEMORY);
    }
    b->dtac.sampleNr=b->tNr;
    b->dtac.tacNr=b->moveNr+1;  // tissue as the last curve
    b->dtac.isframe=ttac->isframe;
    for(int i=0; i<b->dtac.sampleNr; i++) {
      b->dtac.x1[i]=ttac->x1[i]; b->dtac.x[i]=ttac->x[i]; b->dtac.x2[i]=ttac->x2[i];}
    ret=tacDuplicate(&b->dtac, &b->ditac);
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->diitac); 
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->ddtac); 
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->dddtac); 
    if(ret!=TPCERROR_OK) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
    /* Integrate/derivate input data, change times, and interpolate to tissue data */
    double bd[btac->sampleNr], bdd[btac->sampleNr];
    double bi[btac->sampleNr], bii[btac->sampleNr];
    if(liDerivate3(btac->x, btac->c[0].y, btac->sampleNr, bd, bdd, 0)) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    if(btac->isframe==0) {
      if(liIntegrate(btac->x, btac->c[0].y, btac->sampleNr, bi, 3, 0) ||
         liIntegrate(btac->x, bi, btac->sampleNr, bii, 3, 0))
        ret=TPCERROR_INVALID_VALUE;
    } else {
      if(liIntegrateFE(btac->x1, btac->x2, btac->c[0].y, btac->sampleNr, bi, bii, 0))
        ret=TPCERROR_INVALID_VALUE;
    }
    if(ret) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
    double dx[btac->sampleNr], dx2[btac->sampleNr], *dxp;
    if(btac->isframe==0) dxp=dx; else dxp=dx2;
    for(int di=0; di<b->dtac.tacNr-1; di++) {
      for(int i=0; i<btac->sampleNr; i++) {
        dx[i]=btac->x[i]+b->dtmin+di*b->dtstep;
        dx2[i]=btac->x2[i]+b->dtmin+di*b->dtstep;
      }
      if(b->dtac.isframe==0) {
        ret=liInterpolate(dx, btac->c[0].y, btac->sampleNr, b->dtac.x, b->dtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
    
        if(!ret) ret=liInterpolate(dxp, bi, btac->sampleNr, b->dtac.x, b->ditac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dxp, bii, btac->sampleNr, b->dtac.x, b->diitac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dx, bd, btac->sampleNr, b->dtac.x, b->ddtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dx, bdd, btac->sampleNr, b->dtac.x, b->dddtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
      } else {
        ret=liInterpolateForPET(dx, btac->c[0].y, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->dtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dxp, bi, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->ditac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dxp, bii, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->diitac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dx, bd, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->ddtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dx, bdd, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->dddtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
      }
      if(ret) {
        freeDelayCMFitData(&ltdata);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
        return(TPCERROR_INVALID_XRANGE);
      }
    }
    if(verbose>6) {
      printf("integrated input curve 1:\n");
      for(int i=0; i<b->dtac.sampleNr; i++) 
        printf("%g  %g %g %g\n", b->dtac.x[i], b->dtac.c[0].y[i], b->ditac.c[0].y[i], b->diitac.c[0].y[i]);
      fflush(stdout);
    }
 
 
    /* Allocate memory for linear LSQ */
    if(verbose>5) {printf("  preparing data for LLSQ\n"); fflush(stdout);}
    b->llsq_n=3+b->model; 
    b->llsq_m=b->dtac.sampleNr;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    b->llsq_mat=(double*)malloc((b->llsq_n*b->llsq_m)*sizeof(double));
    if(b->llsq_mat==NULL) {freeDelayCMFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
    b->llsq_a=(double**)malloc(b->llsq_n*sizeof(double*));
    if(b->llsq_a==NULL) {freeDelayCMFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
    for(int ni=0; ni<b->llsq_n; ni++) b->llsq_a[ni]=b->llsq_mat+ni*b->llsq_m;
    b->llsq_b=(double*)malloc(b->llsq_m*sizeof(double));
    b->llsq_x=(double*)malloc(b->llsq_n*sizeof(double));
    b->llsq_wp=(double*)malloc(b->llsq_n*sizeof(double));
    b->llsq_zz=(double*)malloc(b->llsq_m*sizeof(double));
    b->llsq_i=(int*)malloc(b->llsq_n*sizeof(int));
    if(b->llsq_b==NULL || b->llsq_x==NULL || b->llsq_wp==NULL || b->llsq_zz==NULL || b->llsq_i==NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      freeDelayCMFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
  }
 
  /* Integrate tissue curve */
  {
    if(verbose>5) {printf("  preparing tissue data\n"); fflush(stdout);}
    int ti=b->dtac.tacNr-1; // index of tissue curve after the input curves
    for(int i=0; i<b->dtac.sampleNr; i++) b->dtac.c[ti].y[i]=ttac->c[ci].y[i];
 
    int ret=liDerivate(b->dtac.x, b->dtac.c[ti].y, b->dtac.sampleNr, b->ddtac.c[ti].y, b->dddtac.c[ti].y, 0);
    if(b->dtac.isframe==0) {
      if(!ret) ret=liIntegrate(b->dtac.x, b->dtac.c[ti].y, b->dtac.sampleNr, b->ditac.c[ti].y, 3, 0);
      if(!ret) ret=liIntegrate(b->dtac.x, b->ditac.c[ti].y, b->dtac.sampleNr, b->diitac.c[ti].y, 3, 0);
    } else {
      if(!ret) ret=liIntegrateFE(b->dtac.x1, b->dtac.x2, b->dtac.c[ti].y, b->dtac.sampleNr, 
                                 b->ditac.c[ti].y, b->diitac.c[ti].y, 0);
    }
    if(ret) {
      freeDelayCMFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
 
    if(verbose>10) {
      printf("integrated tissue curves:\n");
      for(int i=0; i<b->dtac.sampleNr; i++) 
        printf("%g  %g %g %g\n", b->dtac.x[i], b->dtac.c[ti].y[i], b->ditac.c[ti].y[i], b->diitac.c[ti].y[i]);
      fflush(stdout);
    }
  }
 
 
  /* Multilinear fitting with each delayed input */
  if(verbose>5) {printf("  LLSQ\n"); fflush(stdout);}
  double ss, ssmin=nan("");
  int mini=-1;
  for(int di=0; di<b->dtac.tacNr-1; di++) {
    /* Setup data matrix A and vector B */
    for(int mi=0; mi<b->llsq_m; mi++) {
      b->llsq_b[mi]=b->dtac.c[b->dtac.tacNr-1].y[mi]; // TTAC
      b->llsq_mat[mi]=b->dtac.c[di].y[mi]; // BTAC
      if(mode==1) {
        b->llsq_mat[mi+1*b->llsq_m]=b->ddtac.c[di].y[mi]; // BTAC derivate
        b->llsq_mat[mi+2*b->llsq_m]=-b->ddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->dddtac.c[di].y[mi]; // BTAC derivate
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->dddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
      } else if(mode==2) {
        b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
        b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->ddtac.c[di].y[mi]; // BTAC derivate
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->ddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
      } else { // mode 0
        b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
        b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->diitac.c[di].y[mi]; // BTAC 2nd integral
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->diitac.c[b->dtac.tacNr-1].y[mi]; // TTAC 2nd integral
      }
      b->llsq_mat[mi]=b->dtac.c[di].y[mi]; // BTAC
      b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
      b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
      if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->diitac.c[di].y[mi]; // BTAC 2nd integral
      if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->diitac.c[b->dtac.tacNr-1].y[mi]; // TTAC 2nd integral
    }
    if(verbose>3) {
      int failnr=0;
      for(int mi=0; mi<b->llsq_m; mi++) {
        if(!isfinite(b->llsq_b[mi])) failnr++;
        if(!isfinite(b->llsq_mat[mi])) failnr++;
        if(!isfinite(b->llsq_mat[mi+1*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+2*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+3*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+4*b->llsq_m])) failnr++;
      }
      if(failnr>0) printf("%d bad values\n", failnr);
    }
    if(verbose>9 && b->llsq_m<100) {
      printf("---- B x A ----\n");
      for(int mi=0; mi<b->llsq_m; mi++) {
        printf("  %g ", b->llsq_b[mi]);
        for(int ni=0; ni<b->llsq_n; ni++)
          printf(" %g", b->llsq_mat[mi+ni*b->llsq_m]);
         printf("\n");
      }
    }
 
    int ret=nnls(b->llsq_a, b->llsq_m, b->llsq_n, b->llsq_b, b->llsq_x, &ss, b->llsq_wp, b->llsq_zz, b->llsq_i);
    if(ret<2 && isfinite(ss)) {
      if(!isfinite(ssmin) || ssmin>ss) {ssmin=ss; mini=di;}
      if(verbose>2) printf("   %d %d %g %g\n", di, mini, ss, ssmin);
    } else {
      if(verbose>1) printf("  failed NNLS\n");
    }
  }
  if(mini<0) {
    if(verbose>1) fprintf(stderr, "Error: all fits failed.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
    freeDelayCMFitData(&ltdata); return(TPCERROR_BAD_FIT);
  }
  if(dt!=NULL && mini>=0) *dt=b->dtmin+mini*b->dtstep;
 
  freeDelayCMFitData(&ltdata); // the local one, not the one possibly given by caller 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacFittime()

int tacFittime ( TAC * d, double * startTime, double * endTime, int * first, int * last, TPCSTATUS * status )

extern

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

Returns

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

See also

tacReadModelingData, tacVerifyTimeOrder, tacSortByTime

Parameters

d	Pointer to TAC containing (regional tissue) data; times can be in minutes or seconds, as long as units are defined. TAC sampleNr is not changed in this function.
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; enter NULL if not needed.
last	Function writes the index of the last included sample (frame) here; enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed.

Definition at line 25 of file tacmodelinput.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(first!=NULL) *first=0; 
  if(last!=NULL) *last=0;
  if(d==NULL || startTime==NULL || endTime==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(-1);
  }
 
  if(d->sampleNr<=0) {
    *startTime=*endTime=0.0; 
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(0);
  }
  if(*endTime<0.0) *endTime=1.0E+100;
 
  /* Change start and end times to seconds if necessary */
  double cf;
  cf=unitConversionFactor(UNIT_MIN, d->tunit);
  if(!isnan(cf)) {*startTime*=cf; *endTime*=cf;}
  if(verbose>1) {
    printf("startTime := %g\n", *startTime);
    printf("endTime := %g\n", *endTime);
  }
 
#if(0)
  /* Check that data range is not outside required range */
  {
    double s, e; 
    if(tacXRange(d, &s, &e) || e<*startTime || s>*endTime) {
      *startTime=*endTime=0.0; 
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X_RANGE);
      return(0);
    }
  }
#endif
 
  /* Get first and last data point inside the range */
  int s, e;
  {
    double f;
    s=e=-1;
    for(int i=0; i<d->sampleNr; i++) {
      if(d->isframe) f=0.5*(d->x1[i]+d->x2[i]); else f=d->x[i];
      if(verbose>3) printf("  f := %g\n", f);
      if(s<0 && f>=*startTime) s=i;
      if(f<=*endTime) e=i; else break;
    }
    if(s<0 || e<0) {
      *startTime=*endTime=0.0; 
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(0);
    }
    if(first!=NULL) *first=s; 
    if(last!=NULL) *last=e;
  }
 
  /* Correct fit range to frame start and end times */
  *startTime=(d->isframe ? d->x1[s] : d->x[s]);
  *endTime=  (d->isframe ? d->x2[e] : d->x[e]);
  /* Change start and end times back to minutes if necessary */
  if(!isnan(cf)) {*startTime/=cf; *endTime/=cf;}
  /* Return the sample number in the fit range */
  return(1+e-s);
}

Referenced by tacReadModelingData().

tacInput2sim()

int tacInput2sim ( TAC * itac, TAC * ttac, TAC * stac, TPCSTATUS * status )

extern

Modify input TAC based on tissue TAC, for use in simulations.

Note

NOT FUNCTIONAL/TESTED YET!

See also

tacReadModelingInput, tacReadModelingData, tacInterpolateInto, tacInterpolateToEqualLengthFrames

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

itac	Pointer to original input TAC structure. Not modified. Data must not have missing values.
ttac	Pointer to tissue TAC structure. Not modified. Data must not have missing values.
stac	Pointer to the new input data, containing input data with new sample times.
status	Pointer to status data; obligatory.

Definition at line 29 of file lisim.c.

  {
  if(status==NULL) return TPCERROR_FAIL;
  int verbose=status->verbose;
  if(verbose>0) printf("%s(itac, ttac, ...)\n", __func__);
 
  /* Check the input */
  if(itac==NULL || ttac==NULL || stac==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(itac->sampleNr<1 || ttac->sampleNr<1 || itac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Delete any previous data */
  tacFree(stac);
 
  /* Make a working copy of input data */
  TAC inp; tacInit(&inp);
  status->error=tacDuplicate(itac, &inp); if(status->error!=TPCERROR_OK) return(status->error);
  /* Make sure that input data contains also frame mid times */
  if(inp.isframe) tacSetX(&inp, NULL);
  /* Convert time units to TTAC units */
  tacXUnitConvert(&inp, ttac->tunit, status);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
 
  /* Sort the input data by sample times */
  tacSortByTime(&inp, status); if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
  /* Add zero sample if necessary */
  tacAddZeroSample(&inp, status);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
 
  /* Make a combined list of sample times based on input and tissue TACs */
  unsigned int maxn=inp.sampleNr+ttac->sampleNr;
  if(inp.isframe) maxn+=2*inp.sampleNr;
  if(ttac->isframe) maxn+=ttac->sampleNr;
  double x[maxn];
 
  int n=0;
  for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x[i])) {
    int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x[i])<1.0E-20) break;
    if(j==n) x[n++]=inp.x[i];
  }
  if(inp.isframe) {
    for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x1[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x1[i])<1.0E-20) break;
      if(j==n) x[n++]=inp.x1[i];
    }
    for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x2[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x2[i])<1.0E-20) break;
      if(j==n) x[n++]=inp.x2[i];
    }
  }
  if(ttac->isframe) {
    for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x2[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x2[i])<1.0E-20) break;
      if(j==n) x[n++]=ttac->x2[i];
    }
    for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x2[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x2[i])<1.0E-20) break;
      if(j==n) x[n++]=ttac->x2[i];
    }
  } else {
    for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x[i])<1.0E-20) break;
      if(j==n) x[n++]=ttac->x[i];
    }
  }
  if(verbose>2) printf("n=%d\n", n);
  statSortDouble(x, n, 0);
 
  /* Put new sample times into output data */
  status->error=tacDuplicate(itac, stac);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
  status->error=tacAllocateMoreSamples(stac, n-stac->sampleNr);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
  stac->sampleNr=n; stac->tunit=inp.tunit; stac->isframe=0;
  for(int i=0; i<n; i++) stac->x[i]=x[i];
 
  /* Interpolate input data into the new sample times */
  for(int i=0; i<inp.tacNr; i++) {
    status->error=liInterpolate(inp.x, inp.c[i].y, inp.sampleNr,
                            stac->x, stac->c[i].y, NULL, NULL, stac->sampleNr, 4, 1, verbose-10);
    if(status->error!=TPCERROR_OK) break;
  }
  if(status->error!=TPCERROR_OK) {tacFree(&inp); tacFree(stac); return(status->error);}
 
  tacFree(&inp);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by bfm1TCM().

tacIntegrate()

int tacIntegrate ( TAC * inp, TAC * out, TPCSTATUS * status )

extern

Integrate TACs from one TAC structure into a new TAC structure.

PET frame lengths are taken into account, if available.

See also

liIntegrate, liIntegratePET, tacInterpolate, tacInterpolateInto

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

inp	Pointer to source TAC structure.

Precondition

Data must be sorted by increasing x.

Parameters

out	Pointer to target TAC structure into which the integrated TAC(s) are written. Structure must be initiated; any previous contents are deleted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 27 of file litac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || out==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  tacFree(out);
  int ret=TPCERROR_OK;
  ret=tacDuplicate(inp, out);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(inp->isframe==0) {
    for(int i=0; i<inp->tacNr; i++) {
      ret=liIntegrate(inp->x, inp->c[i].y, inp->sampleNr, out->c[i].y, 3, 0);
      if(ret!=TPCERROR_OK) break;
    }
  } else {
    for(int i=0; i<inp->tacNr; i++) {
      ret=liIntegratePET(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, NULL, 0);
      if(ret!=TPCERROR_OK) break;
    }
  }
  if(ret!=TPCERROR_OK) tacFree(out);
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  return(ret);
}

tacIntegrateFE()

int tacIntegrateFE ( TAC * inp, TAC * out, TAC * out2, TPCSTATUS * status )

extern

Integrate TACs from one TAC structure into a new TAC structure. Integrals are calculated at frame end times.

TAC must contain frame start and end times.

See also

tacIntegrate, liIntegrateFE, tacIsXContiguous, tacSetXContiguous, tacMultipleSamples

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

inp	Pointer to source TAC structure.

Precondition

Data must be sorted by increasing x. Frames must not overlap.

Parameters

out	Pointer to target TAC structure into which the integrated TAC(s) are written. Structure must be initiated; any previous contents are deleted.
out2	Pointer to target TAC structure into which the 2nd integrals are written. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed.

Definition at line 75 of file litac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || out==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(inp->isframe==0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_INVALID_X);
  }
 
  /* Make space */
  tacFree(out); if(out2!=NULL) tacFree(out2);
  int ret=TPCERROR_OK;
  ret=tacDuplicate(inp, out);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(out2!=NULL) {
    ret=tacDuplicate(inp, out2);
    if(ret!=TPCERROR_OK) {
      tacFree(out);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
  }
 
  for(int i=0; i<inp->tacNr; i++) {
    if(out2==NULL)
      ret=liIntegrateFE(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, NULL, 0);
    else
      ret=liIntegrateFE(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, out2->c[i].y, 0);
    if(ret!=TPCERROR_OK) break;
  }
  if(ret!=TPCERROR_OK) {tacFree(out); if(out2!=NULL) tacFree(out2);}
 
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  return(ret);
}

tacInterpolate()

int tacInterpolate ( TAC * inp, TAC * xinp, TAC * tac, TAC * itac, TAC * iitac, TPCSTATUS * status )

extern

Interpolate and/or integrate TACs from one TAC structure into a new TAC structure, using sample times from another TAC structure which is not changed.

PET frame lengths of output TACs are taken into account in interpolation, if available. Input frame lengths are taken into account if the framing is same as with output TACs, otherwise frame middle times are used.

See also

tacInterpolateInto, tacAUC, tacFramesToSteps, tacXCopy, liInterpolate

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

inp	Pointer to source TAC structure; make sure that x (time) and y (concentration) units are the same as in output TAC structure, because input time range is verified to cover the output data range. Detailed verification of suitability of TAC for interpolation must be done elsewhere.

Precondition

Data must be sorted by increasing x.

Parameters

xinp	Pointer to TAC structure which contains the x values (sample times) for the interpolation and integration. Data in this TAC structure is not modified.

Precondition

Data must be sorted by increasing x.

Parameters

tac	Pointer to target TAC structure into which the interpolated input TAC(s) are added. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
itac	Pointer to target TAC structure into which the integrated input TAC(s) are added. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
iitac	Pointer to target TAC structure into which 2nd integrals of input TAC(s) are added. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed.

Definition at line 141 of file litac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || xinp==NULL || xinp==tac || xinp==itac || xinp==iitac) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(tac==NULL && itac==NULL && iitac==NULL) { // there's nothing to do
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
  /* Check the function input */
  if(inp->sampleNr<1 || inp->tacNr<1 || xinp->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check that input does not have any missing values */
  if(tacNaNs(inp)>0 || tacXNaNs(xinp)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
 
 
  /* If input data have similar x values as requested already, then
     copy x1 and x2 into input if available */
  if(xinp->isframe && tacXMatch(inp, xinp, verbose-2) && inp->sampleNr<=xinp->sampleNr) {
    for(int i=0; i<inp->sampleNr; i++) {
      inp->x[i]=xinp->x[i]; inp->x1[i]=xinp->x1[i]; inp->x2[i]=xinp->x2[i];
    }
    inp->isframe=xinp->isframe;
  }
 
 
  /* Check that there is no need for excess extrapolation */
  {
    double start1, end1, start2, end2, range2;
    if(tacXRange(inp, &start1, &end1) || tacXRange(xinp, &start2, &end2)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return TPCERROR_INVALID_X;
    }
    range2=end2-start2;
    if(verbose>2) printf("time ranges: %g - %g vs %g - %g\n", start1, end1, start2, end2);
    if(start1>start2+0.2*range2 || end1<end2-0.3*range2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
  }
 
  /* Delete any previous contents in output structs */
  if(tac!=NULL) tacFree(tac);
  if(itac!=NULL) tacFree(itac);
  if(iitac!=NULL) tacFree(iitac);
 
  /* Allocate memory for output data */
  int ret=TPCERROR_OK;
  if(tac!=NULL) ret=tacAllocate(tac, xinp->sampleNr, inp->tacNr); 
  if(ret==TPCERROR_OK && itac!=NULL) 
    ret=tacAllocate(itac, xinp->sampleNr, inp->tacNr);
  if(ret==TPCERROR_OK && iitac!=NULL) 
    ret=tacAllocate(iitac, xinp->sampleNr, inp->tacNr);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  if(tac!=NULL) {tac->tacNr=inp->tacNr; tac->sampleNr=xinp->sampleNr;}
  if(itac!=NULL) {itac->tacNr=inp->tacNr; itac->sampleNr=xinp->sampleNr;}
  if(iitac!=NULL) {iitac->tacNr=inp->tacNr; iitac->sampleNr=xinp->sampleNr;}
 
 
  /* Copy header information */
  if(tac!=NULL) {
    tacCopyHdr(inp, tac);
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+i);
  }
  if(itac!=NULL) {
    tacCopyHdr(inp, itac);
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, itac->c+i);
  }
  if(iitac!=NULL) {
    tacCopyHdr(inp, iitac);
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, iitac->c+i);
  }
  /* Copy X (time) information */
  if(tac!=NULL) {
    tac->isframe=xinp->isframe; tacXCopy(xinp, tac, 0, tac->sampleNr-1);
  }
  if(itac!=NULL) {
    itac->isframe=xinp->isframe; tacXCopy(xinp, itac, 0, itac->sampleNr-1);
  }
  if(iitac!=NULL) {
    iitac->isframe=xinp->isframe; tacXCopy(xinp, iitac, 0, iitac->sampleNr-1);
  }
 
  /* Check if TACs have the same frame times already */
  if(tacXMatch(inp, tac, verbose-2) && inp->sampleNr>=xinp->sampleNr) {
    if(verbose>2) printf("  same frame times\n");
    int ret=0, ri, fi;
    double *i1, *i2;
    /* copy the values directly and integrate in place, if requested */
    for(ri=0; ri<inp->tacNr && ret==0; ri++) {
      if(tac!=NULL) 
        for(fi=0; fi<tac->sampleNr; fi++)
          tac->c[ri].y[fi]=inp->c[ri].y[fi];
      /* If integrals not requested then that's it for this TAC */
      if(itac==NULL && iitac==NULL) continue;
      /* otherwise set pointers for integrals and then integrate */
      if(itac!=NULL) i1=itac->c[ri].y; else i1=NULL;      
      if(iitac!=NULL) i2=iitac->c[ri].y; else i2=NULL;
      if(inp->isframe)
        ret=liIntegratePET(inp->x1, inp->x2, inp->c[ri].y,
                           xinp->sampleNr, i1, i2, verbose-2);
      else
        ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
                          xinp->x, NULL, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
    }
    if(ret) {
      if(tac!=NULL) tacFree(tac); 
      if(itac!=NULL) tacFree(itac);
      if(iitac!=NULL) tacFree(iitac);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK); // that's it then
  }
 
  /* Interpolate and/or integrate inp data to tac sample times,
     taking into account tac frame lengths if available */
  if(verbose>2) printf("  different frame times\n");
  ret=0;
  double *i0, *i1, *i2;
  for(int ri=0; ri<inp->tacNr && ret==0; ri++) {
    /* set pointers for output y arrays */
    if(tac!=NULL) i0=tac->c[ri].y; else i0=NULL;      
    if(itac!=NULL) i1=itac->c[ri].y; else i1=NULL;      
    if(iitac!=NULL) i2=iitac->c[ri].y; else i2=NULL;
    if(xinp->isframe)
      ret=liInterpolateForPET(inp->x, inp->c[ri].y, inp->sampleNr,
                xinp->x1, xinp->x2, i0, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
    else
      ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
                xinp->x, i0, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
  }
  if(ret) {
    if(tac!=NULL) tacFree(tac); 
    if(itac!=NULL) tacFree(itac);
    if(iitac!=NULL) tacFree(iitac);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacInterpolateInto()

int tacInterpolateInto ( TAC * inp, TAC * tac, TAC * itac, TAC * iitac, TPCSTATUS * status )

extern

Add TACs from one TAC structure into another TAC structure, interpolating the input TACs and allocating space if necessary.

PET frame lengths of output TACs are taken into account in interpolation, if available. Input frame lengths are taken into account if the framing is same as with output TACs, otherwise frame middle times are used.

See also

tacInterpolate, tacAUC, tacInterpolateToEqualLengthFrames, liInterpolateForPET

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

inp	Pointer to source TAC struct; make sure that x (time) and y (concentration) units are the same as in output TAC struct, because input time range is verified to cover the output data range. Detailed verification of suitability of TAC for interpolation must be done elsewhere.

Precondition

Data must be sorted by increasing x.

Parameters

tac	Pointer to target TAC struct into which the input TAC(s) are added. Struct must contain at least one TAC, and the x or x1 and x2 values.

Precondition

Data must be sorted by increasing x.

Parameters

itac	Pointer to target TAC struct into which integrated input TAC(s) are added; enter NULL if not needed.
iitac	Pointer to target TAC struct into which 2nd integrals of input TAC(s) are added; enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed.

Definition at line 330 of file litac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tac==NULL || inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  
  /* Check the function input */
  if(tac->sampleNr<1 || inp->sampleNr<1 || inp->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(itac!=NULL && itac->sampleNr!=tac->sampleNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
  if(iitac!=NULL && iitac->sampleNr!=tac->sampleNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
 
 
  /* Check that input does not have any missing values */
  /* Check that there are no missing sample times in output TACs */
  if(tacNaNs(inp)>0 || tacXNaNs(tac)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
  
  /* Check that there is no need for excess extrapolation */
  {
    double start1, end1, start2, end2, range2;
    if(tacXRange(inp, &start1, &end1) || tacXRange(tac, &start2, &end2)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return TPCERROR_INVALID_X;
    }
    range2=end2-start2;
    if(verbose>2) 
      printf("time ranges: %g - %g vs %g - %g\n", start1, end1, start2, end2);
    if(start1>(start2+0.2*range2) || end1<(end2-0.3*range2)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
  }
 
  /* Allocate memory for interpolated data */
  if(tacAllocateMore(tac, inp->tacNr)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  /* Copy header information */
  for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+tac->tacNr+i);
  /* Same for the integrals */
  if(itac!=NULL) {
    if(tacAllocateMore(itac, inp->tacNr)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      return TPCERROR_OUT_OF_MEMORY;
    }
    for(int i=0; i<inp->tacNr; i++) 
      tacCopyTacchdr(inp->c+i, itac->c+itac->tacNr+i);
  }
  if(iitac!=NULL) {
    if(tacAllocateMore(iitac, inp->tacNr)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      return TPCERROR_OUT_OF_MEMORY;
    }
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, iitac->c+iitac->tacNr+i);
  }
 
  /* Check if TACs have the same frame times already */
  if(tacXMatch(inp, tac, verbose-2) && inp->sampleNr>=tac->sampleNr) {
    if(verbose>2) printf("same frame times\n");
    int ret=0, ri, fi;
    double *i1, *i2;
    /* copy the values directly and integrate in place, if requested */
    for(ri=0; ri<inp->tacNr && ret==0; ri++) {
      for(fi=0; fi<tac->sampleNr; fi++)
        tac->c[tac->tacNr+ri].y[fi]=inp->c[ri].y[fi];
      /* If integrals not requested then that's it for this TAC */
      if(itac==NULL && iitac==NULL) continue;
      /* otherwise set pointers for integrals and then integrate */
      if(itac!=NULL) i1=itac->c[itac->tacNr+ri].y; else i1=NULL;      
      if(iitac!=NULL) i2=iitac->c[iitac->tacNr+ri].y; else i2=NULL;
      if(tac->isframe)
        ret=liIntegratePET(tac->x1, tac->x2, tac->c[tac->tacNr+ri].y,
                           tac->sampleNr, i1, i2, verbose-2);
      else
        ret=liInterpolate(tac->x, tac->c[tac->tacNr+ri].y, tac->sampleNr,
               tac->x, NULL, i1, i2, tac->sampleNr, 4, 1, verbose-2);
    }
    if(ret) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
    tac->tacNr+=inp->tacNr;
    if(itac!=NULL) itac->tacNr+=inp->tacNr;
    if(iitac!=NULL) iitac->tacNr+=inp->tacNr;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK); // that's it then
  }
 
  /* Interpolate and, if requested, integrate inp data to tac sample times,
     taking into account tac frame lengths if available */
  int ret=0;
  double *i1, *i2;
  for(int ri=0; ri<inp->tacNr && ret==0; ri++) {
    /* set pointers for integrals */
    if(itac!=NULL) i1=itac->c[itac->tacNr+ri].y; else i1=NULL;
    if(iitac!=NULL) i2=iitac->c[iitac->tacNr+ri].y; else i2=NULL;
    if(tac->isframe)
      ret=liInterpolateForPET(inp->x, inp->c[ri].y, inp->sampleNr,
                tac->x1, tac->x2, tac->c[tac->tacNr+ri].y, 
                i1, i2, tac->sampleNr, 4, 1, verbose-2);
    else
      ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
                tac->x, tac->c[tac->tacNr+ri].y, 
                i1, i2, tac->sampleNr, 4, 1, verbose-2);
  }
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  tac->tacNr+=inp->tacNr;
  if(itac!=NULL) itac->tacNr+=inp->tacNr;
  if(iitac!=NULL) iitac->tacNr+=inp->tacNr;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by tacReadModelingData(), and tacReadModelingInput().

tacInterpolateToEqualLengthFrames()

int tacInterpolateToEqualLengthFrames ( TAC * inp, double minfdur, double maxfdur, TAC * tac, TPCSTATUS * status )

extern

Interpolate TAC data into data with equal frame lengths, starting from zero.

By default, the shortest frame length or sampling interval in the input data is used.

Remarks

This is needed for TAC convolution.

See also

tacInterpolate, liInterpolateForPET, tacGetSampleInterval, simIsSteadyInterval

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

inp	Pointer to source TAC struct. Missing values are not allowed.

Precondition

Data must be sorted by increasing x.

Parameters

minfdur	Minimum frame length; if NaN or <=0 then not applied, otherwise frame lengths are not allowed to be lower.
maxfdur	Maximum frame length; if NaN or <=0 then not applied, otherwise frame lengths are not allowed to be higher.
tac	Pointer to target TAC struct into which the interpolated TACs are written.

Precondition

Struct must be initiated.

Parameters

status

Pointer to status data; enter NULL if not needed.

Definition at line 214 of file lisim.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(tac1, %g, %g, ...)\n", __func__, minfdur, maxfdur); fflush(stdout);}
  if(tac==NULL || inp==NULL || inp->sampleNr<1 || inp->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(minfdur>0.0 && maxfdur>0.0 && minfdur>maxfdur) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
  if(!tacIsX(inp)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
 
  /* Get the minimum sample interval (that is higher than 0) */
  double freq=nan("");
  int ret=tacGetSampleInterval(inp, 1.0E-08, &freq, NULL);
  //printf("dataset based freq: %g\n", freq);
  if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  if(minfdur>0.0 && freq<minfdur) freq=minfdur;
  if(maxfdur>0.0 && freq>maxfdur) freq=maxfdur;
  //printf("refined freq: %g\n", freq);
 
  /* Get the x range */
  double xmax=nan("");
  if(tacXRange(inp, NULL, &xmax)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  //printf("xmax: %g\n", xmax);
  if(freq>0.5*xmax) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
 
 
  /* Calculate the number of required interpolated samples */
  int nreq=ceil(xmax/freq); //printf("required_n := %d\n", nreq);
  /* Check that it is not totally stupid */
  if(nreq<1 || nreq>2E+06) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_BIG);
    return TPCERROR_TOO_BIG;
  }
 
  /* Setup the output TAC */
  if(tacAllocate(tac, nreq, inp->tacNr)!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  tac->tacNr=inp->tacNr; tac->sampleNr=nreq;
  /* Copy header information */
  tacCopyHdr(inp, tac);
  for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+i);
 
  /* Set X (sample time) information */
  tac->isframe=1;
  for(int i=0; i<nreq; i++) tac->x1[i]=(double)i*freq;
  for(int i=0; i<nreq; i++) tac->x2[i]=(double)(i+1)*freq;
  for(int i=0; i<nreq; i++) tac->x[i]=0.5*(tac->x1[i]+tac->x2[i]);
 
  /* Calculate mean value during each interval */
  if(verbose>2) {printf("  interpolating...\n"); fflush(stdout);}
  for(int i=0; i<inp->tacNr; i++) {
    if(liInterpolateForPET(inp->x, inp->c[i].y, inp->sampleNr, tac->x1, tac->x2, 
                           tac->c[i].y, NULL, NULL, tac->sampleNr, 2, 1, 0)!=0) {
      tacFree(tac);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
  }
  if(verbose>2) {printf("  ... done.\n"); fflush(stdout);}
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacPlotFitSVG()

int tacPlotFitSVG ( TAC * tac1, TAC * tac2, const char * main_title, const double x1, const double x2, const double y1, const double y2, const char * fname, TPCSTATUS * status )

extern

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

Returns

enum tpcerror (TPCERROR_OK when successful).

See also

tacPlotLineSVG, mtgaPlotSVG

Author

Vesa Oikonen

Parameters

tac1	Measured data points. Frame mid times must be set.
tac2	Fitted data points. Times can be different but unit must be the same as in tac1.
main_title	String for plot main title, or "".
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 file name.
status	Pointer to status data; enter NULL if not needed.

Definition at line 27 of file tacfitplot.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check data */
  if(tac1==NULL || tac2==NULL || tac1->tacNr<1 || tac2->tacNr<1 || 
     tac1->sampleNr<1 || tac2->sampleNr<1) 
  {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(fname==NULL || strnlen(fname, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return TPCERROR_INVALID_FILENAME;
  }
  if(tac1->tacNr!=tac2->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
    return TPCERROR_INCOMPATIBLE_DATA;
  }
  if(status!=NULL && status->forgiving==0) {
    // any strict tests in here
  }
 
  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=0, color_nr=0;
  int max_symbol_nr=0, symbol_nr=0;
  SVG_LEGENDS legends; svg_init_legends(&legends);
  FILE *fp_svg=NULL;
 
  int is_label=0; if(tac1->tacNr>1) is_label=1;
 
  /* Determine the plot min and max x values */
  if(tacSampleXRange(tac1, &tx1, &tx2)) {
    if(verbose>0) {printf("  tac1 failed\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  minx=tx1; maxx=tx2;
  if(tacSampleXRange(tac2, &tx1, &tx2)) {
    if(verbose>0) {printf("  tac2 failed\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  if(minx>tx1) minx=tx1; 
  if(maxx<tx2) maxx=tx2;
  if(minx>0.0) {double f=maxx-minx; minx-=0.05*f; if(minx<0.0) minx=0.0;}
  if(!isnan(x1)) minx=x1; 
  if(!isnan(x2)) maxx=x2;
  if(verbose>10) printf("  minx := %g\n  maxx:=%g\n", minx, maxx);
 
  /* Determine the plot min and max y values inside the x range */
  if(tacYRangeInXRange(tac1, -1, minx, maxx, &ty1, &ty2, NULL, NULL, NULL, NULL)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  miny=ty1; maxy=ty2;
  if(tacYRangeInXRange(tac2, -1, minx, maxx, &ty1, &ty2, NULL, NULL, NULL, NULL)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  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("  miny := %g\n  maxy:=%g\n", miny, maxy);
 
  /* Calculate the axis ticks */
  viewports.x.min=minx; viewports.x.max=maxx;
  viewports.y.min=miny; viewports.y.max=maxy;
  viewports.label_area_viewport.is=is_label; // needed for x axis ticks
  ret=svg_calculate_axes(&viewports, verbose-11); 
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Set x and y axis titles based on activity and time units */
  strcpy(x_title, "");
  if(unitIsTime(tac1->tunit)) 
    sprintf(x_title, "Time (%s)", unitName(tac1->tunit));
  else if(tac1->tunit!=UNIT_UNKNOWN)
    strcpy(x_title, unitName(tac1->tunit));
  strcpy(y_title, "");
  if(tac1->cunit!=UNIT_UNKNOWN) strcpy(y_title, unitName(tac1->cunit));
 
  /* 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-13);
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Initiate graphics file */
  fp_svg=svg_initiate(fname, 0, 0, &viewports, NULL, verbose-13);
  if(fp_svg==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Put the graph titles into their own viewports */
  ret=svg_create_main_title(fp_svg, main_title, "", &viewports, NULL,verbose-13);
  if(!ret) 
    ret=svg_create_yaxis_title(fp_svg, y_title, &viewports, NULL, verbose-13);
  if(!ret) 
    ret=svg_create_xaxis_title(fp_svg, x_title, &viewports, NULL, verbose-13);
 
  /* Put the plot into its own viewport */
  if(!ret) 
    ret=svg_start_plot_viewport(fp_svg, &viewports, NULL, verbose-13);
 
  /*  Start coordinate area viewport */
  if(!ret) 
    ret=svg_start_coordinate_viewport(fp_svg, &viewports, NULL, verbose-13);
 
  /* Write plot axes */
  if(!ret) ret=svg_write_axes(fp_svg, &viewports, NULL, verbose-13);
 
  if(ret) {
    fclose(fp_svg);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /*
   *  Draw the plots
   */
  max_color_nr=0; while(svgColorName(max_color_nr)!=NULL) max_color_nr++;
  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(verbose>3) printf("max_symbol_nr := %d\n", max_symbol_nr);
  if(tac1->tacNr==1) color_nr=0; else color_nr=1;
  symbol_nr=0;
  for(ri=0, n=0; ri<tac1->tacNr; ri++) {
    sprintf(tac_id, "plot_%d", n);
    strcpy(tac_title, tac1->c[ri].name);
    /* Draw the fitted line */
    for(si=0; si<tac2->sampleNr; si++) if(!isnan(tac2->c[ri].y[si])) break;
    for(ei=tac2->sampleNr-1; ei>=si; ei--) if(!isnan(tac2->c[ri].y[ei])) break;
    if((ei-si)>0) {
      ret=svg_write_tac(fp_svg, &viewports, 1, tac_id, tac_title,
            tac2->x+si, tac2->c[ri].y+si, 1+ei-si,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLOPEN,
            NULL, verbose-13);
      if(ret) {
        svg_legend_empty(&legends); fclose(fp_svg);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
        return TPCERROR_INVALID_VALUE;
      }
    }
    /* Draw the measured points */
    ret=svg_write_tac(fp_svg, &viewports, 2, tac_id, tac_title,
            tac1->x, tac1->c[ri].y, tac1->sampleNr,
            svgColorName(color_nr%max_color_nr), symbol_nr%max_symbol_nr, SYMBOLFILLED,
            NULL, verbose-13);
    if(ret) {
      svg_legend_empty(&legends); fclose(fp_svg);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
    /* Set legend too, if requested */
    if(is_label!=0) {
      int si=0, ci=0; // remainder works only if 2nd operator > 0
      if(max_symbol_nr>0) si=symbol_nr%max_symbol_nr;
      if(max_color_nr>0) ci=color_nr%max_color_nr;
 
      svg_legend_add(&legends, 0, si, SYMBOLFILLED, ci, 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-13);
 
  /* Write the axis ticks */
  if(!ret) {
    if(svg_write_xticks(fp_svg, &viewports, NULL, verbose-13)!=0) ret=1;}
  if(!ret) {
    if(svg_write_yticks(fp_svg, &viewports, NULL, verbose-13)!=0) ret=1;}
 
  /* Close the plot viewport */
  if(!ret) ret=svg_end_plot_viewport(fp_svg, NULL, verbose-13);
 
  if(ret) {
    fclose(fp_svg); svg_legend_empty(&legends);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Make the plot legends into their own viewport */
  ret=0;
  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-13);
  }
  svg_legend_empty(&legends);
  if(ret) {
    fclose(fp_svg);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Close the SVG file */
  ret=svg_close(fp_svg, NULL, verbose-13);
  if(ret) {
    fclose(fp_svg);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacPlotHistogramSVG()

int tacPlotHistogramSVG ( TAC * d, const char * main_title, const double x1, const double x2, const double y1, const double y2, const char * fname, TPCSTATUS * status )

extern

SVG plot of histogram from data given in TAC data structure.

Returns

enum tpcerror (TPCERROR_OK when successful).

See also

tacPlotLineSVG, tacPlotFitSVG, mtgaPlotSVG

Author

Vesa Oikonen

Parameters

d	Pointer to TAC structure. Only the first dataset is used.
main_title	String for plot main title, or "".
x1	Start x value; NaN if determined from data.
x2	End x value; 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 file name.
status	Pointer to status data; obligatory.

Definition at line 30 of file histplot.c.

  {
  if(status==NULL) return TPCERROR_FAIL;
  if(status->verbose>1)
    printf("%s(tac, '%s', %g, %g, %g, %g, '%s')\n", __func__, main_title, x1, x2, y1, y2, fname);
  if(d==NULL || d->sampleNr<1 || d->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA); return(status->error);}
  if(fname==NULL || strnlen(fname, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME); return(status->error);}
 
  /* Determine the plot min and max x and y values */
  double minx, maxx, miny, maxy, tx1, tx2, ty1, ty2;
  if(tacSampleXRange(d, &tx1, &tx2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE); return(status->error);}
  minx=tx1; maxx=tx2;
  if(minx>0.0) {double f=maxx-minx; minx-=0.05*f; if(minx<0.0) minx=0.0;}
  if(isfinite(x1)) minx=x1; 
  if(isfinite(x2)) maxx=x2;
  if(status->verbose>1) {printf("  minx := %g\n  maxx := %g\n", minx, maxx); fflush(stdout);}
  /* Determine the plot min and max y values inside the x range */
  if(tacYRangeInXRange(d, 0, minx, maxx, &ty1, &ty2, NULL, NULL, NULL, NULL)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE); return(status->error);}
//  if(ty1>0.0) ty1=0.0;
//  if(ty2<0.0) ty2=0.0;
  miny=ty1; maxy=ty2;
  if(miny>0.0) {double f=maxy-miny; miny-=0.05*f; if(miny<0.0) miny=0.0;}
  else if(maxy<0.0) maxy=0.0;
  if(isfinite(y1)) miny=y1; 
  if(isfinite(y2)) maxy=y2;
  if(status->verbose>1) {printf("  miny := %g\n  maxy := %g\n", miny, maxy); fflush(stdout);}
 
  /* Calculate the axis ticks */
  struct svg_viewports viewports; svg_init_viewports(&viewports);
  viewports.x.min=minx; viewports.x.max=maxx;
  viewports.y.min=miny; viewports.y.max=maxy;
  viewports.label_area_viewport.is=0; // no plot legends
  if(svg_calculate_axes(&viewports, status->verbose-1)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE); return(status->error);}
 
  /* Set the plot window and window area sizes */
  if(svg_define_viewports(0, 0, strlen(main_title), 0, 0, 0, &viewports, status->verbose-1)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE); return(status->error);}
 
  /* Initiate graphics file */
  FILE *fp=svg_initiate(fname, 0, 0, &viewports, NULL, status->verbose-1);
  if(fp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN); return(status->error);}
 
  /* Put the graph titles into their own viewports */
  if(svg_create_main_title(fp, main_title, "", &viewports, NULL, status->verbose-2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
  /* Put the plot into its own viewport */
  if(svg_start_plot_viewport(fp, &viewports, NULL, status->verbose-2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
  /*  Start coordinate area viewport */
  if(svg_start_coordinate_viewport(fp, &viewports, NULL, status->verbose-3)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
  /* Write plot axes */
  if(svg_write_axes(fp, &viewports, NULL, status->verbose-3)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
 
  /*
   *  Draw the data
   */
  char ilc[9], tmp[1024];
  if(SVG_INLINE) strcpy(ilc, "svg:"); else strcpy(ilc, "");
  /* Initiate the data object group */
  sprintf(tmp, "\n<%sg stroke=\"black\" stroke-width=\"25\" fill=\"black\">\n", ilc);
  if(svg_write(fp, tmp, NULL, status->verbose-5)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
  if(d->isframe==0) { // vertical lines
    /* Open the path */
    sprintf(tmp, "<%spath fill=\"none\" d=\"", ilc);
    svg_write(fp, tmp, NULL, status->verbose-5);
    for(int i=0; i<d->sampleNr; i++) {
      if(!isfinite(d->x[i]) || !isfinite(d->c[0].y[i])) continue;
      double nx1, ny1, nx2, ny2;
      nx1=nx2=viewports.x.origo+d->x[i]*viewports.x.scale;
      ny1=viewports.coordinate_area_viewport.h-(viewports.y.origo+viewports.y.scale*0.0);
      ny2=viewports.coordinate_area_viewport.h-(viewports.y.origo+viewports.y.scale*d->c[0].y[i]);
      /* Do not plot if x is outside viewport */
      if(nx1<0 && nx2<=0) continue;
      if(nx1>viewports.coordinate_area_viewport.w+1 || nx1>viewports.coordinate_area_viewport.w+1)
        continue;
      /* Do not plot if both y's are outside viewport */
      if((ny1<0 || ny2>viewports.coordinate_area_viewport.h+1) &&
         (ny2<0 || ny1>viewports.coordinate_area_viewport.h+1)) continue;
      /* Set y limit to the viewport border */
      if(ny1<0) ny1=0.0; else if(ny2<0) ny2=0.0;
      if(ny1>viewports.coordinate_area_viewport.h+1) ny1=viewports.coordinate_area_viewport.h+1;
      else if(ny2>viewports.coordinate_area_viewport.h+1) ny2=viewports.coordinate_area_viewport.h+1;
      /* Draw */
      sprintf(tmp, " M%.0f %.0f L%.0f %.0f", nx1, ny1, nx2, ny2);
      svg_write(fp, tmp, NULL, status->verbose-5);
    }
    /* Close the path */
    strcpy(tmp, "\" />\n");
    svg_write(fp, tmp, NULL, status->verbose-5);
  } else { // rectangles
    for(int i=0; i<d->sampleNr; i++) {
      if(!isfinite(d->x[i]) || !isfinite(d->c[0].y[i])) continue;
      double nx1, ny1, nx2, ny2;
      nx1=viewports.x.origo+d->x1[i]*viewports.x.scale;
      nx2=viewports.x.origo+d->x2[i]*viewports.x.scale;
      ny1=viewports.coordinate_area_viewport.h-(viewports.y.origo+viewports.y.scale*0.0);
      ny2=viewports.coordinate_area_viewport.h-(viewports.y.origo+viewports.y.scale*d->c[0].y[i]);
      /* Do not plot if x is outside viewport */
      if(nx1<0 && nx2<=0) continue;
      if(nx1>viewports.coordinate_area_viewport.w+1 || nx1>viewports.coordinate_area_viewport.w+1)
        continue;
      /* Do not plot if both y's are outside viewport */
      if((ny1<0 || ny2>viewports.coordinate_area_viewport.h+1) &&
         (ny2<0 || ny1>viewports.coordinate_area_viewport.h+1)) continue;
      /* Set y limit to the viewport border */
      if(ny1<0) ny1=0.0; else if(ny2<0) ny2=0.0;
      if(ny1>viewports.coordinate_area_viewport.h+1) ny1=viewports.coordinate_area_viewport.h+1;
      else if(ny2>viewports.coordinate_area_viewport.h+1) ny2=viewports.coordinate_area_viewport.h+1;
      if(ny1>ny2) {double f=ny1; ny1=ny2; ny2=f;}
      /* Draw */
      sprintf(tmp, " <rect x=\"%.0f\" y=\"%.0f\" width=\"%.0f\" height=\"%.0f\" />\n",
              nx1, ny1, nx2-nx1, fabs(ny2-ny1));
      svg_write(fp, tmp, NULL, status->verbose-5);
    }
  }
  /* Close the data object group */
  sprintf(tmp, "</%sg>\n", ilc);
  if(svg_write(fp, tmp, NULL, status->verbose-5)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
 
  /* Close the coordinate viewport */
  if(svg_end_coordinate_viewport(fp, NULL, status->verbose-1)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
 
  /* Write the axis ticks */
  if(svg_write_xticks(fp, &viewports, NULL, status->verbose-3)!=0 ||
     svg_write_yticks(fp, &viewports, NULL, status->verbose-3)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
 
  /* Close the plot viewport */
  if(svg_end_plot_viewport(fp, NULL, status->verbose-2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
 
  /* Close the SVG file */
  if(svg_close(fp, NULL, status->verbose-1)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(status->error);
  }
 
  if(status->verbose>0) fprintf(stdout, "  written file %s\n", fname);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacPlotLineSVG()

int tacPlotLineSVG ( TAC * tac, const char * main_title, const double x1, const double x2, const double y1, const double y2, const char * fname, TPCSTATUS * status )

extern

Writes specified range of line plots of TACs in SVG format.

Returns

enum tpcerror (TPCERROR_OK when successful).

See also

tacPlotFitSVG, mtgaPlotSVG

Author

Vesa Oikonen

Parameters

tac	TAC data points.
main_title	String for plot main title, or "".
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 file name.
status	Pointer to status data; enter NULL if not needed.

Definition at line 273 of file tacfitplot.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  /* Check data */
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(fname==NULL || strnlen(fname, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return TPCERROR_INVALID_FILENAME;
  }
  if(status!=NULL && status->forgiving==0) {
    // any strict tests in here
  }
 
  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);
  unsigned int max_color_nr=0, color_nr=0;
  unsigned int max_symbol_nr=0, symbol_nr=0;
  SVG_LEGENDS legends; svg_init_legends(&legends);
  FILE *fp_svg=NULL;
 
 
  int is_label=0; if(tac->tacNr>1) is_label=1;
 
  /* Determine the plot min and max x values */
  if(tacSampleXRange(tac, &tx1, &tx2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  minx=tx1; maxx=tx2;
  if(minx>0.0) {double f=maxx-minx; minx-=0.05*f; if(minx<0.0) minx=0.0;}
  if(!isnan(x1)) minx=x1; 
  if(!isnan(x2)) maxx=x2;
  if(verbose>1) {
    printf("  minx := %g\n  maxx := %g\n", minx, maxx); fflush(stdout);}
 
  /* Determine the plot min and max y values inside the x range */
  if(tacYRangeInXRange(tac, -1, minx, maxx, &ty1, &ty2, NULL, NULL, NULL, NULL)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  miny=ty1; 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("  miny := %g\n  maxy := %g\n", miny, maxy); fflush(stdout);}
 
  /* Calculate the axis ticks */
  viewports.x.min=minx; viewports.x.max=maxx;
  viewports.y.min=miny; viewports.y.max=maxy;
  viewports.label_area_viewport.is=is_label; // needed for x axis ticks
  ret=svg_calculate_axes(&viewports, verbose-3); 
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Set x and y axis titles based on activity and time units */
  strcpy(x_title, "");
  if(unitIsTime(tac->tunit)) 
    sprintf(x_title, "Time (%s)", unitName(tac->tunit));
  else if(tac->tunit!=UNIT_UNKNOWN)
    strcpy(x_title, unitName(tac->tunit));
  strcpy(y_title, "");
  if(tac->cunit!=UNIT_UNKNOWN) strcpy(y_title, unitName(tac->cunit));
 
  /* 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) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Initiate graphics file */
  fp_svg=svg_initiate(fname, 0, 0, &viewports, NULL, verbose-3);
  if(fp_svg==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Put the graph titles into their own viewports */
  ret=svg_create_main_title(fp_svg, main_title, "", &viewports, NULL,verbose-3);
  if(!ret) 
    ret=svg_create_yaxis_title(fp_svg, y_title, &viewports, NULL, verbose-3);
  if(!ret) 
    ret=svg_create_xaxis_title(fp_svg, x_title, &viewports, NULL, verbose-3);
 
  /* Put the plot into its own viewport */
  if(!ret) 
    ret=svg_start_plot_viewport(fp_svg, &viewports, NULL, verbose-3);
 
  /*  Start coordinate area viewport */
  if(!ret) 
    ret=svg_start_coordinate_viewport(fp_svg, &viewports, NULL, verbose-3);
 
  /* Write plot axes */
  if(!ret) ret=svg_write_axes(fp_svg, &viewports, NULL, verbose-3);
 
  if(ret) {
    fclose(fp_svg);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /*
   *  Draw the plots
   */
  max_color_nr=0; while(svgColorName(max_color_nr)!=NULL) max_color_nr++;
  if(verbose>3) {printf("max_color_nr := %d\n", max_color_nr); fflush(stdout);}
  if(tac->tacNr==1) color_nr=0; else color_nr=1;
  symbol_nr=0;
  for(ri=0, n=0; ri<tac->tacNr; ri++) {
    sprintf(tac_id, "plot_%d", n);
    strcpy(tac_title, tac->c[ri].name);
    /* Draw the line */
    for(si=0; si<tac->sampleNr; si++) if(!isnan(tac->c[ri].y[si])) break;
    for(ei=tac->sampleNr-1; ei>=si; ei--) if(!isnan(tac->c[ri].y[ei])) break;
    if((ei-si)>0) {
      unsigned int ci, si; // remainder works only if 2nd operator > 0
      if(max_color_nr<1) ci=0; else ci=color_nr % max_color_nr;
      if(max_symbol_nr<1) si=0; else si=symbol_nr % max_symbol_nr;
      ret=svg_write_tac(fp_svg, &viewports, 1, tac_id, tac_title,
            tac->x+si, tac->c[ri].y+si, 1+ei-si,
            svgColorName(ci), si, SYMBOLFILLED,
            NULL, verbose-3);
      if(ret) {
        svg_legend_empty(&legends); fclose(fp_svg);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
        return TPCERROR_INVALID_VALUE;
      }
    }
    /* Set legend too, if requested */
    if(is_label!=0) {
      unsigned int ci, si; // remainder works only if 2nd operator > 0
      if(max_color_nr<1) ci=0; else ci=color_nr % max_color_nr;
      if(max_symbol_nr<1) si=0; else si=symbol_nr % max_symbol_nr;
      svg_legend_add(&legends, 0, si, SYMBOLFILLED, ci, tac_title);
    }
    /* Prepare for the next plot */
    color_nr++; n++;
    if(color_nr==max_color_nr) color_nr=0;
  }
 
  /* Close the coordinate viewport */
  ret=svg_end_coordinate_viewport(fp_svg, NULL, verbose-3);
 
  /* Write the axis ticks */
  if(!ret) {
    if(svg_write_xticks(fp_svg, &viewports, NULL, verbose-3)!=0) ret=1;}
  if(!ret) {
    if(svg_write_yticks(fp_svg, &viewports, NULL, verbose-3)!=0) ret=1;}
 
  /* Close the plot viewport */
  if(!ret) ret=svg_end_plot_viewport(fp_svg, NULL, verbose-3);
 
  if(ret) {
    fclose(fp_svg); svg_legend_empty(&legends);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Make the plot legends into their own viewport */
  ret=0;
  if(viewports.label_area_viewport.is!=0) {
    if(verbose>2) {printf("creating plot legends\n"); fflush(stdout);}
    ret=svg_create_legends(fp_svg, &viewports, &legends, NULL, verbose-3);
  }
  svg_legend_empty(&legends);
  if(ret) {
    fclose(fp_svg);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  /* Close the SVG file */
  ret=svg_close(fp_svg, NULL, verbose-3);
  if(ret) {
    fclose(fp_svg);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacReadModelingData()

int tacReadModelingData ( const char * tissuefile, const char * inputfile1, const char * inputfile2, const char * inputfile3, double * fitdur, int cutInput, int * fitSampleNr, TAC * tis, TAC * inp, TPCSTATUS * status )

extern

Read tissue and input data for modelling.

Time units are converted to min and input calibration units to the units of tissue data. Input data is NOT interpolated to tissue times, but original sample times are kept. If input data extends much further than fit duration, the last part is removed to save computation time in simulations. Input data ranges or TAC shapes are NOT verified.

See also

tacReadReference, tacReadModelingInput, tacFittime, tacRead, tacInit, tacInterpolate

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

tissuefile	Tissue data file name; one time sample is sufficient here; required.
inputfile1	1st input data file name; required.
inputfile2	2nd input data file name (or NULL or empty string if not needed); required, if inputfile3 is given.
inputfile3	3rd input data file name (or NULL or empty string if not needed).
fitdur	Fit duration (in minutes); shortened if longer than tissue data; input data is cut (if requested) so that it will not be much longer. Tissue TACs are NOT cut to this time.
cutInput	Cut off too many input samples to make calculation faster by entering <>0 here, or 0 to keep all input samples (which may be needed for delay correction).
fitSampleNr	Nr of time frames (samples) in tissue data that are inside fitdur will be written here; enter NULL if not needed.
tis	Pointer to initiated TAC struct into which tissue data will be written; required.
inp	Pointer to initiated TAC struct into which input data (plasma and/or blood) TACs will be written; required.
status	Pointer to status data; enter NULL if not needed.

Definition at line 130 of file tacmodelinput.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tis==NULL || inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(tissuefile==NULL || inputfile1==NULL || strnlen(inputfile1, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check the function input */
  int input_nr=1;
  if(inputfile2!=NULL && strnlen(inputfile2, 1)>0) input_nr++;
  if(inputfile3!=NULL && strnlen(inputfile3, 1)>0) {
    if(input_nr<2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      return TPCERROR_NO_DATA;
    }
    input_nr++;
  }
  if(verbose>2) printf("input_nr := %d\n", input_nr);
 
  /* Delete any previous data */
  tacFree(tis); tacFree(inp);
  if(fitSampleNr!=NULL) *fitSampleNr=0;
 
  int ret;
  /*
   *  Read tissue data
   */
  if(verbose>1) printf("reading tissue data in %s\n", tissuefile);
  ret=tacRead(tis, tissuefile, status);
  if(ret!=TPCERROR_OK) return(ret);
 
  /* Do not check frame number; static scan is fine here */
  
  /* Check for NaN's */
  if(tacNaNs(tis)>0) {
    tacFree(tis);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
 
  /* Sort the tissue data by increasing sample times */
  ret=tacSortByTime(tis, status);
  if(ret!=TPCERROR_OK) {
    tacFree(tis);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return ret;
  }
  
  /* Make sure that there is no overlap in sample frame times; samples must
     be sorted before this */
  ret=tacCorrectFrameOverlap(tis, status);
  if(ret!=TPCERROR_OK) {tacFree(tis); return ret;}
 
  /*
   *  Read first input data
   */
  if(verbose>1) printf("reading input data 1 in %s\n", inputfile1);
  ret=tacRead(inp, inputfile1, status);
  if(ret!=TPCERROR_OK) {tacFree(tis); return ret;}
  /* Check and correct the sample time unit */
  if(tis->tunit==UNIT_UNKNOWN) tis->tunit=inp->tunit;
  else if(inp->tunit==UNIT_UNKNOWN) inp->tunit=tis->tunit;
  if(inp->tunit==UNIT_UNKNOWN && verbose>0) {
    fprintf(stderr, "Warning: input sample time units not known.\n");}
  if(tis->tunit!=inp->tunit && tacXUnitConvert(inp, tis->tunit, status)) {
    tacFree(tis); tacFree(inp); return ret;
  }
  /* Check TAC nr */
  if(inp->tacNr>1) {
    if(verbose>0)
      fprintf(stderr, "Warning: using only first TAC in %s\n", inputfile1);
    inp->tacNr=1;
  }
  /* Check sample nr */
  if(inp->sampleNr<4) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    tacFree(tis); tacFree(inp); return TPCERROR_TOO_FEW;
  }
  /* Sort the data by increasing sample times */
  tacSortByTime(inp, NULL);
  /* If inp contains isotope and tis does not, then copy it */
  if(tacGetIsotope(tis)!=ISOTOPE_UNKNOWN) {
    int isotope=tacGetIsotope(inp);
    if(isotope!=ISOTOPE_UNKNOWN) tacSetIsotope(tis, isotope);
  }
 
  /*
   *  Read following input files, if available
   */
  char *fname;
  TAC tmptac; tacInit(&tmptac);
  for(int ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=(char*)inputfile2; else fname=(char*)inputfile3;
    if(verbose>1) printf("reading input data %d in %s\n", ii, fname);
    ret=tacRead(&tmptac, fname, status);
    if(ret!=TPCERROR_OK) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;}
    /* Check TAC nr */
    if(tmptac.tacNr>1) {
      if(verbose>0)
        fprintf(stderr, "Warning: using only first TAC in %s\n", fname);
      tmptac.tacNr=1;
    }
    /* Check sample nr */
    if(tmptac.sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return TPCERROR_TOO_FEW;
    }
    /* Sort the data by increasing sample times */
    tacSortByTime(&tmptac, NULL);
 
    /* Check and correct the sample time unit */
    if(tis->tunit==UNIT_UNKNOWN) tis->tunit=tmptac.tunit;
    else if(tmptac.tunit==UNIT_UNKNOWN) tmptac.tunit=tis->tunit;
    if(inp->tunit!=tmptac.tunit && tacXUnitConvert(&tmptac, inp->tunit, status))
    {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
 
    /* Check and correct the sample concentration unit */
    if(inp->cunit==UNIT_UNKNOWN) inp->cunit=tmptac.cunit;
    else if(tmptac.cunit==UNIT_UNKNOWN) tmptac.cunit=inp->cunit;
    if(inp->cunit!=tmptac.cunit && tacYUnitConvert(&tmptac, inp->cunit, status))
    {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
 
    /* Copy to input data */
    if(tacInterpolateInto(&tmptac, inp, NULL, NULL, status)!=TPCERROR_OK) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
    
    tacFree(&tmptac);
  } // next input file
 
  /* Set time unit to min */
  ret=tacXUnitConvert(tis, UNIT_MIN, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check that regional data times are in minutes.\n");
  }
  ret=tacXUnitConvert(inp, UNIT_MIN, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check that input data times are in minutes.\n");
  }
  /* Check that input and tissue time ranges are about the same */
  double iend, tend;
  {
    ret=tacXRange(inp, NULL, &iend); if(ret==0) ret=tacXRange(tis, NULL, &tend);
    if(ret || iend<=0.0 || tend<=0.0) {
      tacFree(tis); tacFree(inp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
    if(tend>10.0*iend || tend<0.10*iend) {
      if(verbose>0) fprintf(stderr, "Warning: check the sample time units.\n");
    }
  }
  ret=tacYUnitConvert(inp, tis->cunit, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check the calibration units.\n");
  }
 
  /*
   *  Check and set fit time length
   */
  if(verbose>1) printf("checking and setting fit time length\n");
  /* Set fit duration */
  double starttime=0, endtime=*fitdur;
  int fnr=tacFittime(tis, &starttime, &endtime, NULL, NULL, status);
  if(verbose>3) {
    fprintf(stdout, "tis.sampleNr := %d\n", tis->sampleNr);
    fprintf(stdout, "starttime := %g\n", starttime);
    fprintf(stdout, "endtime := %g\n", endtime);
    //fprintf(stdout, "first := %d\n", first);
    //fprintf(stdout, "last := %d\n", last);
    fprintf(stdout, "fitSampleNr := %d\n", fnr);
  }
  *fitdur=endtime;
  if(fitSampleNr!=NULL) *fitSampleNr=fnr;
 
  /* Check that input data does not end much before fitdur */
  if(*fitdur>1.2*iend) {
    tacFree(tis); tacFree(inp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return TPCERROR_TOO_FEW;
  }
 
  /* Cut off too many input samples to make calculation faster */
  if(cutInput && iend>*fitdur) {
    if(verbose>1) printf("cut off too many input samples\n");
    int i; double f;
    for(i=0; i<inp->sampleNr; i++) {
      if(inp->isframe) f=0.5*(inp->x1[i]+inp->x2[i]); else f=inp->x[i];
      if(f>(*fitdur)) break;
    }
    if(i<inp->sampleNr) i++; 
    inp->sampleNr=i;
    if(inp->sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(tis); tacFree(inp); return TPCERROR_TOO_FEW;
    }
  }
  if(verbose>2) fprintf(stdout, "inp.sampleNr := %d\n", inp->sampleNr);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacReadModelingInput()

int tacReadModelingInput ( const char * inputfile1, const char * inputfile2, const char * inputfile3, TAC * inp, TPCSTATUS * status )

extern

Read arterial input data for modelling.

Time and calibration units to the units of the first data. Input data ranges or TAC shapes are NOT verified.

See also

tacReadModelingData, tacReadReference, tacFittime, tacInterpolate

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

inputfile1	1st input data file name; required.
inputfile2	2nd input data file name (or NULL or empty string if not needed); required, if inputfile3 is given.
inputfile3	3rd input data file name (or NULL or empty string if not needed).
inp	Pointer to initiated TAC data structure into which input data (plasma and/or blood) TACs will be written; required.

See also

tacInit, tacFree

Parameters

status

Pointer to status data; enter NULL if not needed.

Definition at line 560 of file tacmodelinput.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(inputfile1==NULL || strnlen(inputfile1, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }  
 
  /* Check the function input */
  int input_nr=1;
  if(inputfile2!=NULL && strnlen(inputfile2, 1)>0) input_nr++;
  if(inputfile3!=NULL && strnlen(inputfile3, 1)>0) {
    if(input_nr<2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      return TPCERROR_NO_DATA;
    }
    input_nr++;
  }
  if(verbose>2) {printf("input_nr := %d\n", input_nr); fflush(stdout);}
 
  /* Delete any previous data */
  tacFree(inp);
 
  int ret;
  /*
   *  Read first input data
   */
  if(verbose>1) {
    printf("reading input data 1 in %s\n", inputfile1); fflush(stdout);
  }
  ret=tacRead(inp, inputfile1, status);
  if(ret!=TPCERROR_OK) {return ret;}
  /* Check for NaN's */
  if(tacNaNs(inp)>0) {
    tacFree(inp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
  /* Sort the data by increasing sample times */
  ret=tacSortByTime(inp, status);
  if(ret!=TPCERROR_OK) {
    tacFree(inp); return ret;
  }
  /* Check TAC nr */
  if(inp->tacNr>1) {
    if(verbose>0)
      fprintf(stderr, "Warning: using only first TAC in %s\n", inputfile1);
    inp->tacNr=1;
  }
  /* Check sample nr */
  if(inp->sampleNr<4) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    tacFree(inp); return TPCERROR_TOO_FEW;
  }
 
  /*
   *  Read following input files, if available
   */
  char *fname;
  TAC tmptac; tacInit(&tmptac);
  for(int ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=(char*)inputfile2; else fname=(char*)inputfile3;
    if(verbose>1) {
      printf("reading input data %d in %s\n", ii, fname);
      fflush(stdout);
    }
    ret=tacRead(&tmptac, fname, status);
    if(ret!=TPCERROR_OK) {
      tacFree(inp); tacFree(&tmptac); return ret;}
    /* Check TAC nr */
    if(tmptac.tacNr>1) {
      if(verbose>0)
        fprintf(stderr, "Warning: using only first TAC in %s\n", fname);
      tmptac.tacNr=1;
    }
    /* Check sample nr */
    if(tmptac.sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(inp); tacFree(&tmptac); return TPCERROR_TOO_FEW;
    }
    /* Sort the data by increasing sample times */
    tacSortByTime(&tmptac, status);
    if(ret!=TPCERROR_OK) {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    /* Check and correct the sample time unit */
    if(inp->tunit==UNIT_UNKNOWN) inp->tunit=tmptac.tunit;
    else if(tmptac.tunit==UNIT_UNKNOWN) tmptac.tunit=inp->tunit;
    if(inp->tunit!=tmptac.tunit && tacXUnitConvert(&tmptac, inp->tunit, status))
    {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    /* Check and correct the sample concentration unit */
    if(inp->cunit==UNIT_UNKNOWN) inp->cunit=tmptac.cunit;
    else if(tmptac.cunit==UNIT_UNKNOWN) tmptac.cunit=inp->cunit;
    if(inp->cunit!=tmptac.cunit && tacYUnitConvert(&tmptac, inp->cunit, status))
    {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    /* Copy to input data */
    if(tacInterpolateInto(&tmptac, inp, NULL, NULL, status)!=TPCERROR_OK) {
      tacFree(inp); tacFree(&tmptac); return ret;
    }
    tacFree(&tmptac);
  } // next input file
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

tacReadReference()

int tacReadReference ( TAC * tis, const char * reference, TAC * ref, int * refIndex, TPCSTATUS * status )

extern

Read reference tissue TAC.

Reference tissue TAC is read either from a TAC file or from TAC structure containing all regional TACs based on its name. Reference tissue TACs are then either placed in a separate TAC structure, if provided by the user, or added/marked inside the existing TAC structure. Reference TAC(s) are marked with sw=2 (best or only 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 assumed to be the best guess of true reference region and marked with value sw=2. When necessary, reference data units are converted to match the existing tissue TAC. If reference TAC(s) are read from a file, then this function verifies that the sample times do match the existing TACs. Reference TAC shapes are NOT verified.

See also

tacReadModelingData, tacFittime, tacRead, tacInit, tacInterpolate, tacSelectedTACs, tacSelectBestReference, tacReadModelingInput.

Returns

The number of reference TACs, and 0 in case of an error.

Author

Vesa Oikonen

Parameters

tis	Pointer to TAC structure containing the regional tissue TACs, possibly also the reference TAC(s). If pointer to reference TAC structure is not given (below), then the reference TAC(s) are marked/added in this TAC structure; otherwise reference TAC(s) may be moved from this structure into the reference TAC structure.
reference	String containing either the name of reference TAC file, or reference region name or number in the TAC structure (above).
ref	Pointer to output TAC structure into which the reference tissue TACs are placed. Enter NULL, if reference TACs are to be placed/kept in the TAC structure with all tissue TACs. TACs in this structure will be marked with either sw=2 (best guess of true reference region) or sw=1 (less probable reference TACs). Any previous contents are deleted.

Precondition

Struct must be initiated before calling this function.

Parameters

refIndex	Index of the best reference region; enter NULL if not needed.
status	Pointer to status data; enter NULL if not needed.

Definition at line 390 of file tacmodelinput.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tis==NULL || reference==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return 0;
  }
  if(strnlen(reference, 1)<1 || tis->sampleNr<1 || !tacIsX(tis)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return 0;
  }
 
  /* First, try to read reference from file */
  TAC temp; tacInit(&temp);
  if(verbose>1) printf("trying to read file %s\n", reference);
  if(tacRead(&temp, reference, status)==TPCERROR_OK) {
    /* Check sample nr */
    if(temp.sampleNr<tis->sampleNr) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(&temp); return 0;
    }
    int n, ret=0, refi=0;
    /* Convert units */
    if(tacXUnitConvert(&temp, tis->tunit, status)) ret++;
    if(tacYUnitConvert(&temp, tis->cunit, status)) ret++;
    if(verbose>0 && ret>0)
      fprintf(stderr, "Warning: check the units of reference and tissue data.\n");
    /* Check that sample times do match */
    if(!tacXMatch(tis, &temp, verbose-1)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      tacFree(&temp); return 0;
    } else if(tis->isframe && !temp.isframe) {
      /* Copy frame start and end times, in case those had been lost */
      tacXCopy(tis, &temp, 0, tis->sampleNr-1);
      temp.isframe=1;
    }
    /* Select the best reference region */
    if(temp.tacNr==1) {
      temp.c[0].sw=2; refi=0; // just one
    } else {
      for(int i=0; i<temp.tacNr; i++) temp.c[i].sw=1;
      refi=tacSelectBestReference(&temp); if(refi<0) refi=0;
      temp.c[refi].sw=2;
    }
    if(verbose>1) printf("selected ref region := %s\n", temp.c[refi].name);
    /* If ref contains isotope and tis does not, then copy it */
    if(tacGetIsotope(tis)!=ISOTOPE_UNKNOWN) {
      int isotope=tacGetIsotope(&temp);
      if(isotope!=ISOTOPE_UNKNOWN) tacSetIsotope(tis, isotope);
    }
    /* Copy to output struct, if pointer is given; otherwise copy to tis */
    if(ref!=NULL) {
      if(verbose>2) printf("adding %d ref tacs to ref struct\n", temp.tacNr);
      ret=tacDuplicate(&temp, ref);
      if(ret!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, ret);
        tacFree(&temp); return 0;
      }
      if(refIndex!=NULL) *refIndex=refi;
    } else {
      if(verbose>2) printf("adding %d ref tacs to tis struct\n", temp.tacNr);
      for(int i=0; i<tis->tacNr; i++) tis->c[i].sw=0;
      ret=tacAllocateMore(tis, temp.tacNr);
      for(int i=0; i<temp.tacNr && ret==TPCERROR_OK; i++)
        ret=tacCopyTacc(&temp.c[i], &tis->c[tis->tacNr+i], tis->sampleNr);
      if(ret!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, ret);
        tacFree(&temp); return 0;
      }
      if(refIndex!=NULL) *refIndex=tis->tacNr+refi;
      tis->tacNr+=temp.tacNr;
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    n=temp.tacNr; tacFree(&temp);
    if(verbose>1) printf("%d ref region(s) found in file.\n", n);
    return(n);
    // reference file was read and processed
  } else if(verbose>1) {
    printf("  '%s' was not a file\n", reference); fflush(stdout);
  }
 
  /* So 'reference' did not contain valid, accessible filename,
     but is it a region name or number? */
  if(verbose>1) printf("trying to find reference in TAC struct\n");
  int n, refi;
  n=tacSelectTACs(tis, reference, 1, status);
  if(verbose>1) {printf("%d region(s) matched.\n", n); fflush(stdout);}
  if(n<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_REFERENCE);
    return(0);
  }
  if(n==tis->tacNr) {
    if(verbose>2) {printf("... which means all region(s) matched.\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(0);
  }
  refi=tacSelectBestReference(tis);
  if(refi<0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(0);
  }
  tis->c[refi].sw=2; if(refIndex!=NULL) *refIndex=refi;
 
  /* If output struct was not given, then we are ready */
  if(ref==NULL) {
    if(verbose>1) {printf("reference TAC tagged.\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(n);
  }
 
  /* Move reference TACs to output struct */
  if(verbose>1) printf("moving reference TAC(s) to output TAC struct\n");
  int ret;
  ret=tacDuplicate(tis, ref);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return 0;
  }
  int i=ref->tacNr-1; ret=0;
  while(i>=0 && !ret) {
    if(ref->c[i].sw==0) ret=tacDeleteTACC(ref, i); 
    i--;
  }
  i=tis->tacNr-1;
  while(i>=0 && !ret) {
    if(tis->c[i].sw) ret=tacDeleteTACC(tis, i); 
    i--;
  }
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return 0;
  }
  if(refIndex!=NULL) { // set reference index to index inside ref struct
    *refIndex=-1;
    for(i=0; i<ref->tacNr; i++) if(ref->c[i].sw==2) {*refIndex=i; break;}
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(n);
}

tacToPAR()

int tacToPAR ( TAC * tac, PAR * par, TPCSTATUS * status )

extern

Copy the contents of TAC struct into PAR struct.

In addition to data, the header contents such as region names and units are copied if possible.

See also

parInit, parFree, parAllocateWithTAC, parWrite

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

tac	Pointer to the TAC struct to be copied.
par	Pointer to PAR struct; any old contents are deleted.

Precondition

Struct must be initiated with parInit.

Parameters

status

Pointer to status data; enter NULL if not needed

Definition at line 169 of file partac.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(par, tac)\n", __func__);
 
  int ret;
 
  /* Check the input */
  if(tac==NULL || par==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Allocate memory and copy most of headers */
  ret=parAllocateWithTAC(par, tac, tac->sampleNr, status);
  if(ret!=TPCERROR_OK) return(ret);
 
  /* Copy the values */
  for(int i=0; i<tac->tacNr; i++)
    for(int j=0; j<par->parNr; j++)
      par->r[i].p[j]=tac->c[i].y[j];
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}

tacVb()

int tacVb ( TAC * ttac, const int i, TAC * btac, double Vb, const int simVb, const int petVolume, TPCSTATUS * status )

extern

Correct TTACs for vascular blood, or simulate its effect.

See also

tacDelay, tacReadModelingInput, tacInterpolate, imgVb

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

ttac	Pointer to TAC data to process.
i	Index of TAC to be processed; enter <0 to process all.
btac	Pointer to BTAC data to subtract or add; must contain exactly same sample times as TTAC and y values must be in the same units.
Vb	Vb fraction [0,1].
simVb	Switch to either subtract vascular volume (0) or to simulate it (1).
petVolume	Switch to model vascular volume as either 0 : Cpet = (1-Vb)*Ct + Vb*Cb, or 1 : Cpet = Ct + Vb*Cb
status	Pointer to status data; obligatory.

Definition at line 138 of file lisim.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(ttac, %d, btac, %g, %d, %d)\n", __func__, i, Vb, simVb, petVolume);
  if(ttac==NULL || btac==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(ttac->tacNr<1 || btac->tacNr<1 || ttac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(ttac->sampleNr>btac->sampleNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
    return TPCERROR_INCOMPATIBLE_DATA;
  }
  if(!(Vb>=0.0) || !(Vb<=1.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(Vb==1.0 && petVolume==0 && simVb==0) { // combination would lead to divide by zero
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(i>=ttac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
 
  if(simVb==0) {
    if(verbose>1) printf("subtracting blood\n");
    for(int ri=0; ri<ttac->tacNr; ri++) if(i<0 || i==ri) {
      for(int fi=0; fi<ttac->sampleNr; fi++) {
        ttac->c[ri].y[fi]-=Vb*btac->c[0].y[fi];
        if(petVolume==0) ttac->c[ri].y[fi]/=(1.0-Vb);
      }
    }
  } else {
    if(verbose>1) printf("adding blood\n");
    for(int ri=0; ri<ttac->tacNr; ri++) if(i<0 || i==ri) {
      for(int fi=0; fi<ttac->sampleNr; fi++) {
        if(petVolume==0) ttac->c[ri].y[fi]*=(1.0-Vb);
        ttac->c[ri].y[fi]+=Vb*btac->c[0].y[fi];
      }
    }
  }
 
  return(TPCERROR_OK);
}