bf_1tcm.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpctac.h"
#include "tpctacmod.h"
#include "tpccm.h"
/*****************************************************************************/
#include "tpcbfm.h"
/*****************************************************************************/
 
/*****************************************************************************/
int bfm1TCM(
  TAC *input,
  TAC *tissue,
  int bfNr,
  const double k2min,
  const double k2max,
  const int distr, 
  TAC *bf,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(itac, ttac, %d, %g, %g, %d, bf)\n", __func__, bfNr, k2min, k2max, distr);
  if(input==NULL || tissue==NULL || bf==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(input->sampleNr<1 || input->tacNr<1 || tissue->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(input->sampleNr<3 || bfNr<2) {
    if(verbose>1) printf("invalid sample or BF number\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return TPCERROR_TOO_FEW;
  }
  if(fabs(k2min)<1.0E-10 || fabs(k2min)>=k2max) { // range cannot be computed otherwise
    if(verbose>1) printf("invalid k2 range\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  int ret;
 
  /* Check units */
  if(status==NULL || !status->forgiving) {
    if(!unitIsTime(input->tunit) || input->tunit!=tissue->tunit) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_UNIT);
      return TPCERROR_UNKNOWN_UNIT;
    }
  }
 
  /* Allocate memory for basis functions */
  if(verbose>1) printf("allocating memory for basis functions\n");
  tacFree(bf);
  ret=tacAllocate(bf, tissue->sampleNr, bfNr);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return ret;
  }
  /* Copy and set information fields */
  bf->tacNr=bfNr; bf->sampleNr=tissue->sampleNr;
  bf->format=TAC_FORMAT_TSV_UK;
  bf->isframe=tissue->isframe; tacXCopy(tissue, bf, 0, tissue->sampleNr-1);
  bf->tunit=tissue->tunit; bf->cunit=tissue->cunit;
  for(int bi=0; bi<bf->tacNr; bi++) sprintf(bf->c[bi].name, "B%5.5d", bi+1);
  /* Compute the range of k2 values to size fields */
  if(verbose>1) printf("computing k2 values\n");
  int zero=0; if(k2min<0.0) {zero=1; bf->c[0].size=0.0;}
  if(distr==0) { // printf("  log-based distribution\n");
    double a, b, c;
    a=log10(fabs(k2min)); b=log10(k2max); c=(b-a)/(double)(bfNr-1-zero);
    for(int bi=zero; bi<bf->tacNr; bi++) {
      bf->c[bi].size=pow(10.0, (double)(bi-zero)*c+a);
    }
  } else { // printf("  even distribution\n");
    double c=(k2max-fabs(k2min))/(double)(bfNr-1-zero);
    for(int bi=zero; bi<bf->tacNr; bi++) {
      bf->c[bi].size=fabs(k2min)+(double)(bi-zero)*c;
    }
  }
 
 
#if(0) // The replasment works better with image-derived input with long frames
  
  /* Allocate memory for simulated TAC */
  double *sim;
  sim=(double*)malloc(input->sampleNr*sizeof(double));
  if(sim==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    tacFree(bf); return TPCERROR_OUT_OF_MEMORY;
  }
 
  /* Simulate the basis functions at input time points, and interpolate to tissue sample times */
  double *cbi=NULL;
  if(input->isframe!=0) {
    // Time frames are available:
    // Integrate BTAC to frame middle times and use it as input */ 
    if(verbose>1) printf("using BTAC integral as simulation input\n");
    cbi=(double*)malloc(input->sampleNr*sizeof(double));
    if(cbi==NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      tacFree(bf); free(sim); return TPCERROR_OUT_OF_MEMORY;
    }
    ret=liIntegratePET(input->x1, input->x2, input->c[0].y, input->sampleNr, cbi, NULL, verbose-10);
    if(ret) {
      if(verbose>1) printf("liIntegratePET() = %d\n", ret);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      tacFree(bf); free(sim); free(cbi); return TPCERROR_INVALID_VALUE;
    }
  }
  if(verbose>1) printf("computing basis functions\n");
  ret=0;
  for(int bi=0; bi<bf->tacNr && !ret; bi++) {
    if(input->isframe==0) { // time frames not available
      ret=simC1(input->x, input->c[0].y, input->sampleNr, 1.0, bf->c[bi].size, sim);
    } else { // time frames are available
      ret=simC1_i(input->x, cbi, input->sampleNr, 1.0, bf->c[bi].size, sim);
    }
    if(ret) break;
    if(verbose>100) {
      printf("\nk2 := %g\n", bf->c[bi].size);
      printf("simulated TAC:\n");
      for(int i=0; i<input->sampleNr; i++) printf("  %12.6f  %12.3f\n", input->x[i], sim[i]);
    }
    /* interpolate to PET time frames */
    if(tissue->isframe)
      ret=liInterpolateForPET(input->x, sim, input->sampleNr, 
               tissue->x1, tissue->x2, bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
    else
      ret=liInterpolate(input->x, sim, input->sampleNr, tissue->x,
               bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
    if(ret) break;
  } // next basis function
  free(sim); free(cbi);
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
#else
 
  /* Interpolate the input function to times based on both input and tissue data */
  /* When input TAC has long frames, the last simulated tissue frame will be biased without this. */
  TAC inp; tacInit(&inp);
  tacInput2sim(input, tissue, &inp, status);
  
  /* Allocate memory for simulated TAC */
  double *sim;
  sim=(double*)malloc(inp.sampleNr*sizeof(double));
  if(sim==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    tacFree(bf); return TPCERROR_OUT_OF_MEMORY;
  }
 
  if(verbose>1) printf("computing basis functions\n");
  ret=0;
  for(int bi=0; bi<bf->tacNr && !ret; bi++) {
    ret=simC1(inp.x, inp.c[0].y, inp.sampleNr, 1.0, bf->c[bi].size, sim);
    if(ret) break;
    if(verbose>100) {
      printf("\nk2 := %g\n", bf->c[bi].size);
      printf("simulated TAC:\n");
      for(int i=0; i<inp.sampleNr; i++) printf("  %12.6f  %12.3f\n", inp.x[i], sim[i]);
    }
    /* interpolate to PET time frames */
    if(tissue->isframe)
      ret=liInterpolateForPET(inp.x, sim, inp.sampleNr, 
               tissue->x1, tissue->x2, bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
    else
      ret=liInterpolate(inp.x, sim, inp.sampleNr, tissue->x,
               bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
    if(ret) break;
  } // next basis function
  free(sim);
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  tacFree(&inp);
#endif
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/