seems to work, but not yet tested

v2/pk/fit2auc.c

@@ -2,9 +2,6 @@
     @brief Calculation of PTAC AUC from specified time to infinity.
     @copyright (c) Turku PET Centre
     @author Vesa Oikonen
-
-    @note Work-in-progress.
-
  */
 /// @cond
 /*****************************************************************************/
@@ -110,8 +107,86 @@ int main(int argc, char **argv)
   }
 
 
+  /* 
+   *  Read the file
+   */
+  if(verbose>1) printf("reading %s\n", parfile);
+  PAR par; parInit(&par);
+  if(parRead(&par, parfile, &status)!=TPCERROR_OK) {
+    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
+    parFree(&par); return(2);
+  }
+  if(verbose>2) {
+    printf("fileformat := %s\n", parFormattxt(par.format));
+    printf("parNr := %d\n", par.parNr);
+    printf("tacNr := %d\n", par.tacNr);
+  }
+
+
+  /*
+   *  Allocate place for AUCs
+   */
+  PAR auc; parInit(&auc);
+  if(parAllocate(&auc, 1, par.tacNr)!=TPCERROR_OK) {
+    fprintf(stderr, "Error: cannot allocate memory.\n");
+    parFree(&par); return(5);
+  }
+  auc.parNr=1; auc.tacNr=0;
+  /* Set the file format */
+  auc.format=PAR_FORMAT_UNKNOWN;
+  if(aucfile[0]) auc.format=parFormatFromExtension(aucfile);
+  if(auc.format==PAR_FORMAT_UNKNOWN) auc.format=PAR_FORMAT_TSV_UK;
+  /* Set parameter name */
+  strcpy(auc.n[0].name, "AUC");
+
+
+  /*
+   *  Compute the AUCs (if possible)
+   */
+  for(int ti=0; ti<par.tacNr; ti++) {
+    if(verbose>2) printf("%s\n", par.r[ti].name);
+    double a;
+    if(mfEvalIntToInf(modelCode(par.r[ti].model), par.parNr, par.r[ti].p, starttime, &a, verbose-2)) {
+      if(verbose>0) fprintf(stderr, "Error: model %s not supported.\n", modelCode(par.r[ti].model)); 
+      continue;
+    }
+    if(verbose>1) printf("AUC for %s := %g\n", par.r[ti].name, a);
+    strcpy(auc.r[auc.tacNr].name, par.r[ti].name);
+    auc.r[auc.tacNr].p[0]=a;
+    auc.tacNr++;
+  } // next curve
+  if(auc.tacNr==0) {
+    fprintf(stderr, "Error: model not supported.\n"); 
+    parFree(&par); parFree(&auc);
+    return(3);
+  }
+
+
+  /*
+   *  Print and save results
+   */
+  if(verbose>0 || !aucfile[0]) parWrite(&auc, stdout, PAR_FORMAT_TSV_UK, 0, NULL);
+  if(aucfile[0]) {
+    /* Save file */
+    if(verbose>1) printf("  saving %s\n", aucfile);
+    FILE *fp;
+    fp=fopen(aucfile, "w");
+    if(fp==NULL) {
+      fprintf(stderr, "Error: cannot open file for writing.\n");
+      parFree(&par); parFree(&auc); return(11);
+    }
+    int ret=parWrite(&auc, fp, PAR_FORMAT_UNKNOWN, 1, &status);
+    fclose(fp);
+    if(ret!=TPCERROR_OK) {
+      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
+      parFree(&par); parFree(&auc); return(12);
+    }
+    if(verbose>0) printf("AUC(s) saved in %s\n", aucfile);
+  }
 
 
+  /* Free memory */
+  parFree(&par); parFree(&auc);
 
   return(0);
 }