imgflowb.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpccsv.h"
#include "tpctac.h"
#include "tpcimage.h"
#include "tpcli.h"
#include "tpclinopt.h"
#include "tpctacmod.h"
#include "tpctacimg.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Make a map of cerebral blood flow from dynamic radiowater PET image,.",
  "correcting for the PVE (Iida et al., 2000).",
  " ",
  "Currently only NIfTI format is supported.",
  "Radioactivity concentrations in image and blood data must be in same units.", 
  "Maps of grey matter blood flow and alpha (PTF) are saved.",
  " ",
  "Usage: @P [options] imgfile btacfile fgmap agmap",
  " ",
  "Options:",
  " -end=<Fit end time (sec)>",
  "     Use data from 0 to end time; by default, 300 s.",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x BTAC AUC) are set to zero.",
  "     Default is 1%.",
  " -aW=<filename>",
  "     Map of white matter alpha (PTF).",
  " -Va=<filename | 0>",
  "     Va map is saved in units mL blood/mL PET volume, or, Va is fixed to 0.",
  " -dT=<filename | delay>",
  "     Time delay map is saved, or, delay is fixed to given time in sec.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgflow, imgdelay",
  " ",
  "Keywords: image, blood flow, radiowater, brain",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/ 
int tacInputForLLSQ( 
  TAC *btac,
  TAC *ttac,
  double dtmin,
  double dtmax,
  double dtstep,
  TAC *b1,
  TAC *b2,
  TAC *b3,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>2) {printf("%s(inp, tis, %g, %g, %g, ...)\n", __func__, dtmin, dtmax, dtstep); fflush(stdout);}
 
  /* Check function input */
  if(btac==NULL || ttac==NULL || btac->sampleNr<1 || ttac->sampleNr<1 || btac->tacNr<1 || ttac->tacNr<1 ) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  if(verbose>3) {
    printf(" %d BTAC ", btac->sampleNr); if(btac->isframe) printf("frames\n"); else printf("samples\n");
    printf(" %d TTAC ", ttac->sampleNr); if(ttac->isframe) printf("frames\n"); else printf("samples\n");
  }
  if(btac->sampleNr<5 || ttac->sampleNr<5) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return(TPCERROR_TOO_FEW);
  }
  double brange[2], trange[2];
  if(tacXRange(ttac, &trange[0], &trange[1])!=0 || tacXRange(btac, &brange[0], &brange[1])!=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",
           trange[0], trange[1], brange[0], brange[1]);
  }
  int dtNr=1;
  if(isfinite(dtstep) && dtstep>0.0) {
    double r=dtmax-dtmin;
    if(!(r>5.0*dtstep)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
      return(TPCERROR_FAIL);
    }
    dtNr=1+r/dtstep;
    if(dtNr>5000) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
      return(TPCERROR_FAIL);
    }
    /* Check that input curve is long enough */
    if((brange[1]+dtmin)<0.95*trange[1]) {
      if(verbose>2) printf("BTAC tmax=%g TTAC tmax=%g dtmin=%g\n", brange[1], trange[1], dtmin);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
  }
 
  /* Check if there is anything to do */
  if(b1==NULL && b2==NULL && b3==NULL) {
    if(verbose>2) {printf("  no output required\n"); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
 
  /* Integrate BTAC */
  double bi[btac->sampleNr], bii[btac->sampleNr], *ix;
  int ret=0;
  if(btac->isframe==0) {
    if(verbose>3) {printf("  integrating BTAC samples\n"); fflush(stdout);}
    ret=liIntegrate(btac->x, btac->c[0].y, btac->sampleNr, bi, 3, 0);
    if(!ret) liIntegrate(btac->x, bi, btac->sampleNr, bii, 3, 0);
    ix=btac->x;
  } else {
    if(verbose>3) {printf("  integrating BTAC frames\n"); fflush(stdout);}
    ret=liIntegrateFE(btac->x1, btac->x2, btac->c[0].y, btac->sampleNr, bi, bii, 0);
    ix=btac->x2;
  }
  if(ret) {
    if(verbose>2) {printf("  -> returns %d\n", ret); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(TPCERROR_INVALID_VALUE);
  }
 
  /* Set space for result BTACs */
  if(verbose>3) {printf("  room for output\n"); fflush(stdout);}
  if(b1!=NULL) tacFree(b1);
  if(b2!=NULL) tacFree(b2);
  if(b3!=NULL) tacFree(b3);
  if(b1!=NULL) ret=tacDuplicate(ttac, b1);
  if(b2!=NULL && !ret) ret=tacDuplicate(ttac, b2);
  if(b3!=NULL && !ret) ret=tacDuplicate(ttac, b3);
  if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  if(dtNr>1) {
    if(b1!=NULL) {ret=tacAllocateMore(b1, dtNr); b1->tacNr=dtNr;}
    if(b2!=NULL && !ret) {ret=tacAllocateMore(b2, dtNr); b2->tacNr=dtNr;}
    if(b3!=NULL && !ret) {ret=tacAllocateMore(b3, dtNr); b3->tacNr=dtNr;}
    if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  }
 
  /* Interpolate to TTAC times */
  if(b1!=NULL) {
    if(verbose>3) {printf("  interpolating BTAC\n"); fflush(stdout);}
    double dx[btac->sampleNr];
    for(int ci=0; ci<b1->tacNr; ci++) {
      double dt=0.0; if(dtNr>1) dt=dtmin+dtstep*(double)ci;
      b1->c[ci].size=dt;
      for(int ti=0; ti<btac->sampleNr; ti++) dx[ti]=btac->x[ti]+dt;
      if(ttac->isframe==0)
        ret=liInterpolate(dx, btac->c[0].y, btac->sampleNr, b1->x, b1->c[ci].y, 
                          NULL, NULL, b1->sampleNr, 3, 1, 0);
      else
        ret=liInterpolateForPET(dx, btac->c[0].y, btac->sampleNr, b1->x1, b1->x2, b1->c[ci].y, 
                                NULL, NULL, b1->sampleNr, 3, 1, 0);
      if(ret) break;
    }
  }
  if(b2!=NULL && !ret) {
    if(verbose>3) {printf("  interpolating BTAC integral\n"); fflush(stdout);}
    double dx[btac->sampleNr];
    for(int ci=0; ci<b2->tacNr; ci++) {
      double dt=0.0; if(dtNr>1) dt=dtmin+dtstep*(double)ci;
      b2->c[ci].size=dt;
      for(int ti=0; ti<btac->sampleNr; ti++) dx[ti]=ix[ti]+dt;
      if(ttac->isframe==0)
        ret=liInterpolate(dx, bi, btac->sampleNr, b2->x, b2->c[ci].y, 
                          NULL, NULL, b2->sampleNr, 3, 1, 0);
      else
        ret=liInterpolateForPET(dx, bi, btac->sampleNr, b2->x1, b2->x2, b2->c[ci].y, 
                                NULL, NULL, b2->sampleNr, 3, 1, 0);
      if(ret) break;
    }
  }
  if(b3!=NULL && !ret) {
    if(verbose>3) {printf("  interpolating BTAC 2nd integral\n"); fflush(stdout);}
    double dx[btac->sampleNr];
    for(int ci=0; ci<b3->tacNr; ci++) {
      double dt=0.0; if(dtNr>1) dt=dtmin+dtstep*(double)ci;
      b3->c[ci].size=dt;
      for(int ti=0; ti<btac->sampleNr; ti++) dx[ti]=ix[ti]+dt;
      if(ttac->isframe==0)
        ret=liInterpolate(dx, bii, btac->sampleNr, b3->x, b3->c[ci].y, 
                          NULL, NULL, b3->sampleNr, 3, 1, 0);
      else
        ret=liInterpolateForPET(dx, bii, btac->sampleNr, b3->x1, b3->x2, b3->c[ci].y, 
                                NULL, NULL, b3->sampleNr, 3, 1, 0);
      if(ret) break;
    }
  }
  if(ret) {
    if(verbose>2) {printf("  -> returns %d\n", ret); fflush(stdout);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(TPCERROR_INVALID_VALUE);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int ai, help=0, version=0, verbose=1;
  int ret;
  char imgfile[FILENAME_MAX], btacfile[FILENAME_MAX];
  char fgfile[FILENAME_MAX], agfile[FILENAME_MAX], fwfile[FILENAME_MAX], awfile[FILENAME_MAX];
  char vafile[FILENAME_MAX], dtfile[FILENAME_MAX], ssfile[FILENAME_MAX];
  char llsqfile[FILENAME_MAX]; // Results from LLSQ: p0-4, dT, sum-of-squares
  char kmapfile[FILENAME_MAX]; // Results converted to CM: Va, K1a, k2a, K1b, k2b, dT, SS; 1/min 
  double fixedVa=nan("");
  double fixeddT=nan("");
  double endtime=300.0; // fit end time in seconds
  double endtimemin=120.0; // shortest fit end time allowed
  double dtrange[2]={-10.,+50.};
  double dtstep=1.0;
  float calcThreshold=0.01;
 
  double pg=1.03;
  double pw=0.86;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=btacfile[0]=(char)0;
  fgfile[0]=agfile[0]=vafile[0]=awfile[0]=fwfile[0]=dtfile[0]=ssfile[0]=(char)0;
  llsqfile[0]=kmapfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "END=", 4)==0) {
      ret=atofCheck(cptr+4, &endtime); if(!ret && endtime>0.0) continue;
    } else if(strncasecmp(cptr, "MIN=", 4)==0) {
      ret=atofCheck(cptr+4, &dtrange[0]); if(!ret && isfinite(dtrange[0])) continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      ret=atofCheck(cptr+4, &dtrange[1]); if(!ret && isfinite(dtrange[1])) continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      double v; ret=atofCheck(cptr+4, &v);
      if(!ret && v<100.0) {calcThreshold=(float)(0.01*v); continue;}
    } else if(strncasecmp(cptr, "AW=", 3)==0) {
      strlcpy(awfile, cptr+3, FILENAME_MAX); if(strlen(awfile)) continue;
    } else if(strncasecmp(cptr, "FW=", 3)==0) {
      strlcpy(fwfile, cptr+3, FILENAME_MAX); if(strlen(fwfile)) continue;
    } else if(strcasecmp(cptr, "VA=0")==0) {
      fixedVa=0.0; continue;
    } else if(strncasecmp(cptr, "VA=", 3)==0) {
      strlcpy(vafile, cptr+3, FILENAME_MAX); if(strlen(vafile)) continue;
    } else if(strncasecmp(cptr, "DT=", 3)==0) {
      double v; ret=atofCheck(cptr+3, &v);
      if(!ret && v>=-50.0 && v<=50.0) {fixeddT=v; continue;}
      strlcpy(dtfile, cptr+3, FILENAME_MAX); if(strlen(dtfile)) continue;
    } else if(strncasecmp(cptr, "SS=", 3)==0) {
      strlcpy(ssfile, cptr+3, FILENAME_MAX); if(strlen(ssfile)) continue;
    } else if(strncasecmp(cptr, "LLSQ=", 5)==0) {
      strlcpy(llsqfile, cptr+5, FILENAME_MAX); if(strlen(llsqfile)) continue;
    } else if(strncasecmp(cptr, "KMAP=", 5)==0) {
      strlcpy(kmapfile, cptr+5, FILENAME_MAX); if(strlen(kmapfile)) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(imgfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(fgfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(agfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
  /* Is something missing? */
  if(!agfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  /* Check options */
  if(endtime<endtimemin) {
    fprintf(stderr, "Error: too short fit time set with option -end.\n");
    return(1);
  }
  if(isfinite(fixeddT)) {
    dtrange[0]=dtrange[1]=fixeddT;
    dtstep=nan("");
  } else {
    if(dtrange[0]>dtrange[1]) {int i=dtrange[0]; dtrange[0]=dtrange[1]; dtrange[1]=i;}
    if((dtrange[1]-dtrange[0])>180.) {
      fprintf(stderr, "Error: too wide delay range.\n");
      return(1);
    } else if((dtrange[1]-dtrange[0])>120.) {
      fprintf(stderr, "Warning: large delay range.\n");
    } else if((dtrange[1]-dtrange[0])<4.) {
      fprintf(stderr, "Error: too short delay range.\n");
      return(1);
    }
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    printf("btacfile := %s\n", btacfile);
    printf("fgfile := %s\n", fgfile);
    printf("agfile := %s\n", agfile);
    printf("endtime := %g s\n", endtime);
    printf("delay_range := %g - %g s\n", dtrange[0], dtrange[1]);
    printf("delay_step_size := %g s\n", dtstep);
    printf("threshold := %g\n", calcThreshold);
    if(awfile[0]) printf("awfile := %s\n", awfile);
    if(fwfile[0]) printf("fwfile := %s\n", fwfile);
    if(vafile[0]) printf("vafile := %s\n", vafile);
    if(isfinite(fixedVa)) printf("fixed_Va := %g\n", fixedVa);
    if(dtfile[0]) printf("dtfile := %s\n", dtfile);
    if(isfinite(fixeddT)) printf("fixed_dT := %g\n", fixeddT);
    if(ssfile[0]) printf("ssfile := %s\n", ssfile);
    if(llsqfile[0]) printf("llsqfile := %s\n", llsqfile);
    if(kmapfile[0]) printf("kmapfile := %s\n", kmapfile);
    fflush(stdout);
  }
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
 
 
  /*
   *  Read BTAC
   */
  if(verbose>1) printf("reading %s\n", btacfile);
  TAC btac; tacInit(&btac);
  ret=tacRead(&btac, btacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), btacfile);
    tacFree(&btac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(btac.format));
    printf("tacNr := %d\n", btac.tacNr);
    printf("sampleNr := %d\n", btac.sampleNr);
    printf("xunit := %s\n", unitName(btac.tunit));
    printf("yunit := %s\n", unitName(btac.cunit));
    printf("isframe := %d\n", btac.isframe);
    fflush(stdout);
  }
  /* Make sure that BTAC frame mid times are set */
  tacSetX(&btac, NULL);
 
 
  /*
   *  Read image data
   */
  if(verbose>1) {printf("reading %s\n", imgfile); fflush(stdout);}
  IMG img; imgInit(&img);
  ret=imgRead(&img, imgfile, &status);
  if(ret!=TPCERROR_OK) { // error
    tacFree(&btac); imgFree(&img);
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), imgfile); fflush(stderr);
    return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  if(verbose>2) imgContents(&img, stdout);
  double ixmin, ixmax;
  if(imgXRange(&img, &ixmin, &ixmax)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times in %s\n", imgfile); fflush(stderr);
    tacFree(&btac); imgFree(&img); return(2);
  }
  if(verbose>1) printf("Image range: %g - %g\n", ixmin, ixmax);
  /* Refine fit end time */
  if(ixmax<endtimemin) {
    fprintf(stderr, "Error: image time range is too short.\n");
    tacFree(&btac); imgFree(&img); return(2);
  }
  if(ixmax<endtime) endtime=ixmax;
 
  /* If BTAC has no frame start and end times, check if those can be obtained from image */
  if(tacimgXMatch(&btac, &img)) tacimgXCopy(&btac, &img);
 
  /* Convert BTAC time units to same as in the image */
  if(btac.tunit==UNIT_UNKNOWN) {
    if(verbose>0) fprintf(stderr, "Warning: missing BTAC time units.\n");
    double bxmin, bxmax;
    if(tacXRange(&btac, &bxmin, &bxmax)!=TPCERROR_OK) {
      fprintf(stderr, "Error: invalid sample times in %s\n", btacfile); fflush(stderr);
      tacFree(&btac); imgFree(&img); return(2);
    }
    unit guessedUnit=img.tunit;
    if(img.tunit==UNIT_SEC && bxmax<0.2*ixmax) guessedUnit=UNIT_MIN;
    if(img.tunit==UNIT_MIN && bxmax>20.*ixmax) guessedUnit=UNIT_SEC;
    fprintf(stderr, "Warning: assuming BTAC sample times are in units %s.\n", unitName(guessedUnit));
    btac.tunit=guessedUnit;
  }
  if(tacXUnitConvert(&btac, img.tunit, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), btacfile); fflush(stderr);
    tacFree(&btac); imgFree(&img); return(2);
  }
  /* If times are in minutes, convert to sec */
  if(btac.tunit==UNIT_MIN) {
    tacXUnitConvert(&btac, UNIT_SEC, NULL);
    imgXUnitConvert(&img, UNIT_SEC);
  }
  double bxmin, bxmax;
  if(tacXRange(&btac, &bxmin, &bxmax)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times in %s\n", btacfile); fflush(stderr);
    tacFree(&btac); imgFree(&img); return(2);
  }
  if(verbose>1) printf("BTAC range: %g - %g\n", bxmin, bxmax);
 
  /* Check that concentration units match in image and BTAC */
  if(img.cunit==UNIT_UNKNOWN) {
    if(btac.cunit!=UNIT_UNKNOWN) {
      img.cunit=btac.cunit;
      fprintf(stderr, "Warning: assuming image concentrations are in units %s.\n", unitName(img.cunit));
    }
  } else if(btac.cunit==UNIT_UNKNOWN) {
    btac.cunit=img.cunit;
    fprintf(stderr, "Warning: assuming BTAC concentrations are in units %s.\n", unitName(img.cunit));
  } else {
    if(tacYUnitConvert(&btac, img.cunit, &status)) { 
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&btac); imgFree(&img); return(2);
    }
  }
 
 
  /*
   *  Check time ranges
   */
  if(verbose>1) {printf("checking time range\n"); fflush(stdout);}
  ret=0;
  if(btac.sampleNr<5) {fprintf(stderr, "Error: BTAC has too few samples.\n"); ret++;}
  if(img.dimt<5) {fprintf(stderr, "Error: image has too few frames.\n"); ret++;}
  if(bxmax<endtimemin) {fprintf(stderr, "Error: BTAC time range is too short.\n"); ret++;}
  if(ixmax<endtimemin) {fprintf(stderr, "Error: image time range is too short.\n"); ret++;}
  if(dtrange[0]<0.0 && bxmax<(endtime+dtrange[0])) {
    endtime+=dtrange[0];
    if(endtime<endtimemin) {
      fprintf(stderr, "Error: BTAC does not extend to required end time.\n"); ret++;}
  }
  int frameNr=0;
  for(int i=0; i<img.dimt; i++) if(img.x[i]<=endtime) frameNr++;
  if(frameNr<5) {fprintf(stderr, "Error: image has too few frames in fit time range.\n"); ret++;}
  if(ret>0) {tacFree(&btac); imgFree(&img); return(2);}
 
  /* Make room for TTAC integrals, used also to pass frame information to BTAC interpolation */
  TAC ttac; tacInit(&ttac);
  ret=tacAllocateWithIMG(&ttac, &img, 3);
  if(ret) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    tacFree(&btac); imgFree(&img); return(3);
  }
  ttac.sampleNr=frameNr;
  ttac.tacNr=3;
 
  /* Make delayed BTACs and integrals, matching image frames, for LLSQ fitting */
  TAC inp, inp2, inp3; tacInit(&inp); tacInit(&inp2); tacInit(&inp3);
  if(tacInputForLLSQ(&btac, &ttac, dtrange[0], dtrange[1], dtstep, &inp, &inp2, &inp3, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&inp3);
    return(4);
  } 
  if(verbose>4) {
    printf("writing delayed, interpolated, and integrated BTACs\n");
    FILE *fp=fopen("imgflowb_input1.tac", "w");
    if(fp!=NULL) tacWrite(&inp, fp, TAC_FORMAT_PMOD, 0, NULL);
    fclose(fp);
    fp=fopen("imgflowb_input2.tac", "w");
    if(fp!=NULL) tacWrite(&inp2, fp, TAC_FORMAT_PMOD, 0, NULL);
    fclose(fp);
    fp=fopen("imgflowb_input3.tac", "w");
    if(fp!=NULL) tacWrite(&inp3, fp, TAC_FORMAT_PMOD, 0, NULL);
    fclose(fp);
  }
 
 
  /* Allocate matrices for LLSQ */
  if(verbose>1) {printf("allocating memory for LLSQ\n"); fflush(stdout);}
  int llsq_m=ttac.sampleNr;
  int llsq_n=5;
  double **llsq_a=(double**)malloc(llsq_n*sizeof(double*));
  double *llsq_mat=(double*)malloc((llsq_n*llsq_m)*sizeof(double));
  if(llsq_a==NULL || llsq_mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for LLSQ.\n");
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&inp3);
    return(5);
  }
  for(int ni=0; ni<llsq_n; ni++) llsq_a[ni]=llsq_mat+ni*llsq_m;
  double llsq_b[llsq_m], llsq_x[llsq_n], llsq_wp[llsq_n], llsq_zz[llsq_m];
  int indexp[llsq_n];
 
 
  /*
   *  Allocate memory for LLSQ parameter images (llsq_n & delay & SS)
   */
  if(verbose>1) {printf("preparing LLSQ parameter map\n"); fflush(stdout);}
  IMG map; imgInit(&map);
  if(imgAllocate(&map, img.dimz, img.dimy, img.dimx, llsq_n+2, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory\n"); fflush(stderr);
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&inp3);
    free(llsq_a); free(llsq_mat);
    return(6);
  }
  imgCopyHeader(&img, &map);
  map.cunit=UNIT_UNKNOWN;
 
 
  /* Calculate threshold value */
  double thrs=tacAUC(&btac, 0, 0.0, img.x2[frameNr-1], NULL);
  thrs*=calcThreshold;
 
  /*
   *  Pixel-by-pixel LLSQ
   */
  long long pxlNr=img.dimx*img.dimy*img.dimz;
  long long okNr=0;
  if(verbose>0) {printf("processing %llu image pixels\n", pxlNr); fflush(stdout);}
  for(int zi=0; zi<img.dimz; zi++) {
    if(img.dimz>2 && verbose>0) {fprintf(stdout, "."); fflush(stdout);}
    for(int yi=0; yi<img.dimy; yi++) for(int xi=0; xi<img.dimx; xi++) {
      for(int pi=0; pi<map.dimt; pi++) map.m[zi][yi][xi][pi]=0.0;
      for(int mi=0; mi<llsq_m; mi++) ttac.c[0].y[mi]=img.m[zi][yi][xi][mi];
//      if(liIntegrateFE(ttac.x1, ttac.x2, ttac.c[0].y, llsq_m, ttac.c[1].y, ttac.c[2].y, 0)!=0) continue;
      if(liIntegratePET(ttac.x1, ttac.x2, ttac.c[0].y, llsq_m, ttac.c[1].y, ttac.c[2].y, 0)!=0) continue;
 
      if(ttac.c[1].y[llsq_m-1]<thrs) continue;
      /* Try with each delayed BTAC */
      int diBest=-1;
      double r2Best=nan("");
      for(int di=0; di<inp.tacNr; di++) {
        /* Setup data matrix A and vector B */
        for(int mi=0; mi<llsq_m; mi++)
          llsq_b[mi]=img.m[zi][yi][xi][mi];       // TTAC
        for(int mi=0; mi<llsq_m; mi++) {
          if(isfinite(fixedVa))
            llsq_mat[mi]=0.0;
          else
            llsq_mat[mi]=inp.c[di].y[mi];         // BTAC
          llsq_mat[mi+1*llsq_m]=inp2.c[di].y[mi]; // BTAC integral
          llsq_mat[mi+2*llsq_m]=-ttac.c[1].y[mi]; // TTAC integral
          llsq_mat[mi+3*llsq_m]=inp3.c[di].y[mi]; // BTAC 2nd integral
          llsq_mat[mi+4*llsq_m]=-ttac.c[2].y[mi]; // TTAC 2nd integral
        }
        /* Compute NNLS */
        double r2=nan("");
        int ret=nnls(llsq_a, llsq_m, llsq_n, llsq_b, llsq_x, &r2, llsq_wp, llsq_zz, indexp);
        if(ret>1) continue;
        if((isfinite(r2) && !isfinite(r2Best)) || r2<r2Best) {
          r2Best=r2;
          diBest=di; 
          for(int pi=0; pi<llsq_n; pi++) map.m[zi][yi][xi][pi]=llsq_x[pi];
          map.m[zi][yi][xi][llsq_n]=inp.c[di].size; // delay time was stored there
          map.m[zi][yi][xi][llsq_n+1]=r2;
        }
      } // next delay
      if(diBest>=0) okNr++;
    }
  }
  if(img.dimz>2 && verbose>0) {fprintf(stdout, "\n"); fflush(stdout);}
  free(llsq_a); free(llsq_mat); /* Free NNLS matrix data */
  if(okNr==0) {
    fprintf(stderr, "Error: no pixels fitted successfully\n"); fflush(stderr);
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&inp3);
    imgFree(&map);
    return(7);
  }
  if(verbose>0) printf("  LLSQ successful in %lld/%lld pixels\n", okNr, pxlNr);
  /* Free original data */
  tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&inp3);
 
  /*
   *  Write requested output images
   */
  if(llsqfile[0]) {
    if(verbose>1) printf("  writing %s\n", llsqfile);
    if(imgWrite(&map, llsqfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map);
      return(11);
    }
    if(verbose>0) {printf("  %s written.\n", llsqfile); fflush(stdout);}
  }
  IMG pimg; imgInit(&pimg); // one-frame image for one outcome parameter
  if(imgAllocate(&pimg, map.dimz, map.dimy, map.dimx, 1, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory\n"); fflush(stderr);
    return(12);
  }
  imgCopyHeader(&map, &pimg);
  if(dtfile[0]) {
    if(verbose>1) printf("  writing %s\n", dtfile);
    pimg.cunit=UNIT_SEC;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=map.m[zi][yi][xi][map.dimt-2];
    }
    if(imgWrite(&pimg, dtfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(13);
    }
    if(verbose>0) {printf("  %s written.\n", dtfile); fflush(stdout);}
  }
  if(ssfile[0]) {
    if(verbose>1) printf("  writing %s\n", ssfile);
    pimg.cunit=UNIT_UNKNOWN;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=map.m[zi][yi][xi][map.dimt-1];
    }
    if(imgWrite(&pimg, ssfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(14);
    }
    if(verbose>0) {printf("  %s written.\n", ssfile); fflush(stdout);}
  }
  if(vafile[0]) {
    if(verbose>1) printf("  writing %s\n", vafile);
    pimg.cunit=UNIT_ML_PER_ML;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=map.m[zi][yi][xi][0];
    }
    if(imgWrite(&pimg, vafile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(15);
    }
    if(verbose>0) {printf("  %s written.\n", vafile); fflush(stdout);}
  }
 
  /* Convert LLSQ parameters into parameters K1a, k2a, K1b, and k2b */ 
  {
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      double p[llsq_n]; for(int i=0; i<llsq_n; i++) p[i]=map.m[zi][yi][xi][i];
      double k2a=(p[2]+sqrt(p[2]*p[2] - 4.0*p[4]))/2.0;
      double k2b=(p[2]-sqrt(p[2]*p[2] - 4.0*p[4]))/2.0;
      double K1a=(k2a*(p[1]-p[0]*p[2]) - p[3] + p[0]*p[4])/sqrt(p[2]*p[2]-4.0*p[4]);
      double K1b=(p[3] - p[0]*p[4] - k2b*(p[1]-p[0]*p[2]))/sqrt(p[2]*p[2]-4.0*p[4]);
//printf("p[] = %g %g %g %g %g\n", p[0], p[1], p[2], p[3], p[4]);
//printf("K1a=%g k2a=%g K1b=%g k2b=%g\n", K1a, k2a, K1b, k2b);
      map.m[zi][yi][xi][1]=60.*K1a;
      map.m[zi][yi][xi][2]=60.*k2a;
      map.m[zi][yi][xi][3]=60.*K1b;
      map.m[zi][yi][xi][4]=60.*k2b;
    }
  }
  if(kmapfile[0]) {
    if(verbose>1) printf("  writing %s\n", kmapfile);
    if(imgWrite(&map, kmapfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(21);
    }
    if(verbose>0) {printf("  %s written.\n", kmapfile); fflush(stdout);}
  }
 
 
  if(fgfile[0]) {
    if(verbose>1) printf("  writing %s\n", fgfile);
    pimg.cunit=UNIT_ML_PER_ML_MIN;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=pg*map.m[zi][yi][xi][2];
      if(pimg.m[zi][yi][xi][0]>1.0) pimg.m[zi][yi][xi][0]=1.0;
    }
    if(imgWrite(&pimg, fgfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(31);
    }
    if(verbose>0) {printf("  %s written.\n", fgfile); fflush(stdout);}
  }
  if(fwfile[0]) {
    if(verbose>1) printf("  writing %s\n", fwfile);
    pimg.cunit=UNIT_ML_PER_ML_MIN;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=pw*map.m[zi][yi][xi][4];
      if(pimg.m[zi][yi][xi][0]>1.0) pimg.m[zi][yi][xi][0]=1.0;
    }
    if(imgWrite(&pimg, fwfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(32);
    }
    if(verbose>0) {printf("  %s written.\n", fwfile); fflush(stdout);}
  }
 
  if(agfile[0]) {
    if(verbose>1) printf("  writing %s\n", agfile);
    pimg.cunit=UNIT_ML_PER_ML;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=(map.m[zi][yi][xi][1]/map.m[zi][yi][xi][2])/pg;
      if(!isfinite(pimg.m[zi][yi][xi][0])) pimg.m[zi][yi][xi][0]=0.0;
      if(pimg.m[zi][yi][xi][0]<0.0) pimg.m[zi][yi][xi][0]=0.0;
      if(pimg.m[zi][yi][xi][0]>1.0) pimg.m[zi][yi][xi][0]=1.0;
    }
    if(imgWrite(&pimg, agfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(41);
    }
    if(verbose>0) {printf("  %s written.\n", agfile); fflush(stdout);}
  }
  if(awfile[0]) {
    if(verbose>1) printf("  writing %s\n", awfile);
    pimg.cunit=UNIT_ML_PER_ML;
    for(int zi=0; zi<pimg.dimz; zi++) for(int yi=0; yi<pimg.dimy; yi++) for(int xi=0; xi<pimg.dimx; xi++) {
      pimg.m[zi][yi][xi][0]=(map.m[zi][yi][xi][3]/map.m[zi][yi][xi][4])/pw;
      if(!isfinite(pimg.m[zi][yi][xi][0])) pimg.m[zi][yi][xi][0]=0.0;
      if(pimg.m[zi][yi][xi][0]<0.0) pimg.m[zi][yi][xi][0]=0.0;
      if(pimg.m[zi][yi][xi][0]>1.0) pimg.m[zi][yi][xi][0]=1.0;
    }
    if(imgWrite(&pimg, awfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(42);
    }
    if(verbose>0) {printf("  %s written.\n", awfile); fflush(stdout);}
  }
 
 
  imgFree(&map); imgFree(&pimg);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/