imgwcbf.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 "tpccm.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).",
  "Partition coefficients for grey and white matter, and blood flow in white",
  "matter, are fixed.",
  " ",
  "Currently only NIfTI format is supported.",
  "Radioactivity concentrations in image and blood data must be in same units.", 
  "Map of CBF (aG*fG + aW*fW) is saved.",
  " ",
  "Usage: @P [options] imgfile btacfile cbfmap",
  " ",
  "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%.",
  " -fG=<filename>",
  "     Map of grey matter blood flow.",
  " -aG=<filename>",
  "     Map of grey matter alpha (PTF).",
  " -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: fit_wcbf, 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], awfile[FILENAME_MAX];
  char vafile[FILENAME_MAX], dtfile[FILENAME_MAX], ssfile[FILENAME_MAX];
  char cbffile[FILENAME_MAX];  // alphaGM*fGM + alphaWM*fWM
  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;
  double fw=0.21; // mL/(min*mL)
  fw/=60.0; // -> mL/(sec*mL)
 
 
  /*
   *  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]=dtfile[0]=ssfile[0]=(char)0;
  llsqfile[0]=kmapfile[0]=cbffile[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, "FG=", 3)==0) {
      strlcpy(fgfile, cptr+3, FILENAME_MAX); if(strlen(fgfile)) continue;
    } else if(strncasecmp(cptr, "AG=", 3)==0) {
      strlcpy(agfile, cptr+3, FILENAME_MAX); if(strlen(agfile)) continue;
    } else if(strncasecmp(cptr, "AW=", 3)==0) {
      strlcpy(awfile, cptr+3, FILENAME_MAX); if(strlen(awfile)) 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(cbffile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
  /* Is something missing? */
  if(!cbffile[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("cbffile := %s\n", cbffile);
    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(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);
    if(fgfile[0]) printf("fgfile := %s\n", fgfile);
    if(agfile[0]) printf("agfile := %s\n", agfile);
    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);
  /* First TAC is used as input */
  if(btac.tacNr>1) btac.tacNr=1;
 
 
  /*
   *  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) {
      ret++; fprintf(stderr, "Error: BTAC does not extend to required end time.\n");}
  }
  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 and its integral, used also to pass frame information to BTAC interpolation */
  TAC ttac; tacInit(&ttac);
  ret=tacAllocateWithIMG(&ttac, &img, 2);
  if(ret) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    tacFree(&btac); imgFree(&img); return(3);
  }
  ttac.sampleNr=frameNr;
  ttac.tacNr=2;
 
 
  /* Make delayed BTACs and integral, matching image frames, for LLSQ fitting */
  TAC inp, inp2; tacInit(&inp); tacInit(&inp2);
  if(tacInputForLLSQ(&btac, &ttac, dtrange[0], dtrange[1], dtstep, &inp, &inp2, NULL, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2);
    return(4);
  } 
  if(verbose>4) {
    printf("writing delayed, interpolated, and integrated BTACs\n");
    FILE *fp=fopen("imgwcbf_input1.tac", "w");
    if(fp!=NULL) tacWrite(&inp, fp, TAC_FORMAT_PMOD, 0, NULL);
    fclose(fp);
    fp=fopen("imgwcbf_input2.tac", "w");
    if(fp!=NULL) tacWrite(&inp2, fp, TAC_FORMAT_PMOD, 0, NULL);
    fclose(fp);
  }
 
 
  /*
   *  Simulate Cw with the known fw, pw, and Ca
   */
  TAC wtac; tacInit(&wtac);
  if(tacDuplicate(&btac, &wtac) || 
     simC1_d(btac.x, btac.c[0].y, btac.sampleNr, fw, fw/pw, wtac.c[0].y)) 
  {
    fprintf(stderr, "Error: cannot simulate Cwm(t).\n");
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2);
    return(4);
  }
  wtac.tacNr=1;
  /* Make delayed WTAC integral, matching image frames, for LLSQ fitting */
  TAC wm2; tacInit(&wm2);
  if(tacInputForLLSQ(&wtac, &ttac, dtrange[0], dtrange[1], dtstep, NULL, &wm2, NULL, &status)) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&wm2); tacFree(&wtac);
    return(4);
  }
  tacFree(&wtac);
  if(verbose>4) {
    printf("writing delayed and integrated WM TACs\n");
    FILE *fp=fopen("imgwcbf_wm.tac", "w");
    if(fp!=NULL) tacWrite(&wm2, 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=4;
  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(&wm2);
    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(&wm2);
    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;
  long long thrsNr=0;
  long long pNr[llsq_n]; for(int i=0; i<llsq_n; i++) pNr[i]=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, NULL, 0)!=0) continue;
      if(liIntegratePET(ttac.x1, ttac.x2, ttac.c[0].y, llsq_m, ttac.c[1].y, NULL, 0)!=0) continue;
 
      if(ttac.c[1].y[llsq_m-1]<thrs) continue;
      thrsNr++;
      /* 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]=ttac.c[0].y[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]=wm2.c[di].y[mi];  // WM TTAC integral
          llsq_mat[mi+3*llsq_m]=-ttac.c[1].y[mi]; // TTAC 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(diBest>=0) {
        for(int pi=0; pi<llsq_n; pi++) if(map.m[zi][yi][xi][pi]>0.0) pNr[pi]+=1;
      }
    }
  }
  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(&wm2);
    imgFree(&map);
    return(7);
  }
  if(verbose>0) {
    printf("  LLSQ successful in %lld/%lld over-threshold pixels\n", okNr, thrsNr);
    for(int pi=0; pi<llsq_n; pi++)
      printf("  p%d non-zero in %lld pixels\n", pi, pNr[pi]);
  }
  /* Free original data */
  tacFree(&btac); imgFree(&img); tacFree(&ttac); tacFree(&inp); tacFree(&inp2); tacFree(&wm2);
 
  /*
   *  Write requested output images
   */
  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);
    imgFree(&map); return(11);
  }
  imgCopyHeader(&map, &pimg);
 
  if(llsqfile[0]) {
    char *pext=filenameGetExtensions(llsqfile);
    char basename[FILENAME_MAX], fname[FILENAME_MAX];
    strcpy(basename, llsqfile);
    filenameRmExtensions(basename);
    for(int i=0; i<llsq_n; i++) {
      if(i==0 && isfinite(fixedVa)) continue;
      sprintf(fname, "%s_p%d%s", basename, 1+i, pext);
      if(verbose>1) printf("  writing %s\n", fname);
      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][i];
      }
      if(imgWrite(&pimg, fname, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
        imgFree(&map); imgFree(&pimg);
        return(12);
      }
      if(verbose>0) {printf("  %s written.\n", fname); fflush(stdout);}
    } 
  }
  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);}
  }
 
  if(cbffile[0]) {
    if(verbose>1) printf("  writing %s\n", cbffile);
    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]=60.*(map.m[zi][yi][xi][1] - map.m[zi][yi][xi][0]*map.m[zi][yi][xi][3]);
      if(!isfinite(pimg.m[zi][yi][xi][0]) || pimg.m[zi][yi][xi][0]<0.0) pimg.m[zi][yi][xi][0]=0.0;
      else if(pimg.m[zi][yi][xi][0]>1.6) pimg.m[zi][yi][xi][0]=1.6;
    }
    if(imgWrite(&pimg, cbffile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr);
      imgFree(&map); imgFree(&pimg);
      return(16);
    }
    if(verbose>0) {printf("  %s written.\n", vafile); fflush(stdout);}
  }
 
  /* Convert LLSQ parameters into parameters fGM, aGM, and aWM */ 
  {
    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 va=p[0];
      double fg=p[3]*pg;
      double aw=0.0; if(p[2]>0.0 && p[3]>(fw/pw)) aw=p[2]/(p[3]-fw/pw);
      if(aw>1.0) aw=1.0;
      double ag=p[1]/fg - aw*fw/fg - p[0]/pg;
      if(!(ag>=0.0)) ag=0.0; else if(ag>1.0) ag=1.0;
      map.m[zi][yi][xi][1]=60.*fg;
      map.m[zi][yi][xi][2]=ag;
      map.m[zi][yi][xi][3]=aw;
    }
  }
  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]=map.m[zi][yi][xi][1];
//      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(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][2];
    }
    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];
    }
    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);
}
/*****************************************************************************/
 
/*****************************************************************************/