Select Git revision
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
absscal.c 46.07 KiB
/** @file absscal.c
* @brief Extraction and calibration of blood sampler data (previously blo2kbq).
* @copyright (c) Turku PET Centre
* @author Vesa Oikonen
*/
/// @cond
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
/*****************************************************************************/
/*****************************************************************************/
static char *info[] = {
"Extraction and calibration of blood TACs measured using",
"automatic blood sampling systems (ABSS, \"blood pump\") (1).",
"Application name was previously blo2kbq, version 3.8.1.",
" ",
"Usage: @P [Options] ABSS_file [BTAC_file]",
" ",
"Options:",
" -c=<filename>",
" Calibration coefficients for ABSS and well-counter should be specified",
" with calibration file (2), given with option -c.",
" -i=<isotope>",
" If ABSS data file does not contain the isotope code or halflife, then",
" it can be given with this option, using codes",
" C-11, F-18, Ga-68, Ge-68, O-15, ...",
" -decay=<on|off>",
" Blood data is corrected for radioactive decay (on, default), or not",
" corrected (off).",
//" Possible background subtraction is done before decay correction.",
" -start=<t>",
" If ABSS measurement was not started at the time of injection,",
" the delay t (sec) must be entered with option -start=<t>. If data",
" collection was started before injection, enter a negative value for t;",
" activities before that time are removed from the output data.",
" -density[=<blood density>]",
" With this option blood radioactivity concentrations (kBq/mL) are",
" divided by the given blood density (by default 1.06)",
" to get concentrations in units kBq/g.",
" -min",
" Sample times are written in minutes (in seconds by default).",
" -mean, -left, -right, -both",
" The counts or radioactivity concentrations are by default calculated",
" as the mean of the channels (-mean). With options -left and -right ",
" only the left (ch2) or right (ch1) channel is used.",
" With option -both the results from both channels are saved.",
" Only applicable to Scanditronics and GEMS data.",
/*
" -z",
" For reading of the measurement date, local time zone can be used with",
" option -z instead of Greenwich main time (default).",
*/
" -abss=<device>",
" ABSS produces device-specific file formats, and detectors have",
" calibration factors. This option is needed if ABSS can not be correctly",
" determined from the filename (*.lis, *.bld, *.alg, or *.txt).",
" Accepted devices are Scanditronics, GEMS, Allogg1, and Allogg2, or",
" numbers 1-4, accordingly.",
" -stdoptions", // List standard options like --help, -v, etc
" ", "Filename for corrected blood TAC data can be entered as the last command-line",
"argument. If not given, then file is written to the directory where ABSS file",
"is; calibrated file is named *.kbq and non-calibrated *.dat by default.",
" ",
"Example:",
" @P -i=O-15 -c=S:/Lab/plasma/bsampler_calibration/pump_cal.dat sr452_blo.bld",
" ",
//"In case of the GEMS online sampler (GE Advange/PET CT), the injection",
//"time is always taken from the first sample time in the *.bld file,",
//"not from the header of the file, like with Scanditronics data.",
//" ",
"References:",
"1. http://www.turkupetcentre.net/petanalysis/input_abss.html",
"2. http://www.turkupetcentre.net/petanalysis/input_abss_calibration.html",
" ",
"See also: absszero, absstime, abssexam, abssfch, dftcat, disp4dft, fitdelay",
" ",
"Keywords: blood, input, ABSS, calibration",
0};
/*****************************************************************************/
/*****************************************************************************/
/* Types and defaults */
#define MAX_PUMPS 4
#define BLOOD_DENSITY 1.06
typedef enum {SCANDITRONICS,ALLOGG,ALLOGG2} Format;
static int both_cols=0;
static int gm_timez=1; /*time zone = Greenwich main time*/
static int pump=0;
char isotope[10]="UNKNOWN";
double halflife=-1.0;
char studynr[MAX_STUDYNR_LEN+1];
/*****************************************************************************/
/* Results data */
struct BDataResults {
double t, ch1, ch2;
};
/*****************************************************************************/
/* Scanditronics data */
struct BDataS {
struct BDataResults results;
double time_of_day;
double time_since_start;
double meas_interval;
int pair1_coincid;
int pair1_det1_cnts;
int pair1_det2_cnts;
int pair2_coincid;
int pair2_det1_cnts;
int pair2_det2_cnts;
int aux_cnts;
};
/* Allogg data */
struct BDataA {
struct BDataResults results;
double time_since_start;
double meas_interval; /* This is calculated from the data */
int data[2];
};
struct DataList {
union {
struct BDataS scanditronics;
struct BDataA allogg;
struct BDataResults results;
};
struct DataList *prev, *next;
};
struct Blood {
/** Date acquired from the file header (comments) */
struct tm date;
/** Date in time_t format. Possibly altered by time data from actual data */
time_t time;
/** When was the background radiation measured */
char bgdate[11];
/** Background radiation */
float background;
/** List of datapoints */
struct DataList *data;
};
/*****************************************************************************/
/* Local functions */
int read_scanditronics_data(
const char *inputfile, struct Blood *blood, int verbose);
int read_allogg_data(
const char *inputfile, struct Blood *blood, int verbose);
int read_allogg2_data(
char *inputfile, struct Blood *blood, int verbose, char *msg);
void correct_scanditronics(struct Blood *blood, int verbose);
int correct_allogg(struct Blood *blood);
int correct_allogg2(struct Blood *blood, int verbose);
void subtract_background(struct Blood *blood);
int get_calfactors(const char *file,time_t bloodt, double *wcf,
double *pumpf, int pump, char *isotope, char *caldate, int verbose);
int find_date(char *str, struct Blood *blood, int verbose);
int find_halflife(char *str, int verbose);
/*****************************************************************************/
/*****************************************************************************/
/* 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;
/*****************************************************************************/
/*****************************************************************************/
/**
* Main
*/
int main(int argc, char *argv[])
{
int ai, help=0, version=0, verbose=1;
double startTime=0.0, density=1.0, wc_factor=-1.0, pump_factor=-1.0;
int decay_correct=1; // 0=no, 1=yes
char pumpfile[FILENAME_MAX], calfile[FILENAME_MAX], outfile[FILENAME_MAX];
char tmp[FILENAME_MAX];
float sum1, sum2;
Format format=SCANDITRONICS;
struct DataList *ptr;
int timeInMinutes=0;
struct Blood blood;
char caldate[25];
int column=0; // use mean (0), left (1), or right (2) channel
char *cptr;
FILE *fp;
int r;
int zero1=1, zero2=1;
time_t t_injection, t_samp;
struct tm tm_injection, tm_samp;
/*
* Get arguments */
if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
pumpfile[0]=calfile[0]=outfile[0]=(char)0;
studynr[0]=(char)0;
/* Options */
for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
//printf("argv[%d] := '%s'\n", ai, argv[ai]);
cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
if(strncasecmp(cptr, "C=", 2)==0 && strlen(cptr)>2) {
strlcpy(calfile, cptr+2, FILENAME_MAX); continue;
} else if(strncasecmp(cptr, "O=", 2)==0 && strlen(cptr)>2) { // deprecated
strlcpy(outfile, cptr+2, FILENAME_MAX); continue;
} else if(strncasecmp(cptr, "I=", 2)==0 && strlen(cptr)>2) {
if(hlCorrectIsotopeCode(cptr+2)!=NULL) {
strcpy(isotope, hlCorrectIsotopeCode(cptr+2));
halflife=hlFromIsotope(isotope);
continue;
}
printf("Error: Incorrect isotope code '%s'\n", cptr+2);
return(1);
} else if(strncasecmp(cptr, "DECAY=", 6)==0 && strlen(cptr)>6) {
cptr+=6;
if(strcasecmp(cptr, "OFF")==0) {decay_correct=0; continue;}
if(strcasecmp(cptr, "ON")==0) {decay_correct=1; continue;}
} else if(strcasecmp(cptr, "DENSITY")==0) {
density=BLOOD_DENSITY; continue;
} else if(strncasecmp(cptr, "DENSITY=", 8)==0) {
cptr+=8; density=atof_dpi(cptr); if(density>0.0) continue;
} else if(strncasecmp(cptr, "MIN", 1)==0) {
timeInMinutes=1; continue;
} else if(strcasecmp(cptr, "SEC")==0) {
timeInMinutes=0; continue;
} else if(strncasecmp(cptr, "Z", 1)==0) {
gm_timez=0; continue;
} else if(strncasecmp(cptr, "LEFT", 1)==0) {
column=1; continue;
} else if(strncasecmp(cptr, "RIGHT", 1)==0) {
column=2; continue;
} else if(strcasecmp(cptr, "BOTH")==0) {
both_cols=1; continue;
} else if(strncasecmp(cptr, "PUMP=", 5)==0 && strlen(cptr)>5) {
cptr+=5;
if(strncasecmp(cptr, "SCANDITRONICS", 1)==0) {format=SCANDITRONICS; pump=1; continue;}
if(strncasecmp(cptr, "GEMS", 1)==0) {format=SCANDITRONICS; pump=2; continue;}
if(strncasecmp(cptr, "ALLOGG1", 1)==0) {format=ALLOGG; pump=3; continue;}
if(strcasecmp(cptr, "ALLOGG2")==0) {format=ALLOGG2; pump=4; continue;}
pump=atoi(cptr);
if(pump==1 || pump==2) {format=SCANDITRONICS; continue;}
if(pump==3) {format=ALLOGG; continue;}
if(pump==4) {format=ALLOGG2; continue;}
} else if(strncasecmp(cptr, "START=", 6)==0 && strlen(cptr)>6) {
startTime=atof_dpi(cptr+6); 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(pumpfile, argv[ai++], FILENAME_MAX);
if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
/* Is something missing? */ if(!pumpfile[0]) {
fprintf(stderr, "Error: missing ABSS file name.\n");
return(1);
}
/* In verbose mode print arguments and options */
if(verbose>1) {
printf("pumpfile := %s\n", pumpfile);
printf("isotope := %s\n", isotope);
printf("halflife := %g\n", halflife);
printf("decay_correct := %d\n", decay_correct);
printf("startTime := %g\n", startTime);
printf("density := %g\n", density);
printf("timeInMinutes := %d\n", timeInMinutes);
printf("column := %d\n", column);
printf("both_cols := %d\n", both_cols);
fflush(stdout);
}
/* Decide which pump we are dealing with, if not given by user */
if(pump==0) {
cptr=pumpfile+strlen(pumpfile)-4;
if(strcasecmp(cptr, ".lis")==0) {pump=1; format=SCANDITRONICS;}
else if(strcasecmp(cptr, ".bld")==0) {pump=2; format=SCANDITRONICS;}
else if(strcmp(cptr, ".alg")==0) {pump=3; format=ALLOGG;}
else if(strcmp(cptr, ".txt")==0) {pump=4; format=ALLOGG2;}
else {
fprintf(stderr, "Error: cannot identify the pump;\n");
fprintf(stderr, "please rename the file to *.lis, *.bld, *.alg or *.txt\n");
fprintf(stderr, "or use option -abss=<id>.\n");
return(1);
}
}
if(verbose>1) {
printf("pump := %d\n", pump);
printf("format := %d\n", format);
}
/*
* Read in data
*/
if(verbose>0) printf("reading %s\n", pumpfile);
switch(format) {
case SCANDITRONICS:
if(read_scanditronics_data(pumpfile, &blood, verbose-2)) return 2;
break;
case ALLOGG:
if(read_allogg_data(pumpfile, &blood, verbose-2)) return 2;
break;
case ALLOGG2:
if(read_allogg2_data(pumpfile, &blood, verbose-2, tmp)) {
fprintf(stderr, "Error in reading %s: %s.\n", pumpfile, tmp);
return 2;
}
break;
}
/* Check that isotope is specified with decay correction */
if(decay_correct!=0) {
if(halflife<=0.0 || strcasecmp(isotope, "UNKNOWN")==0) {
fprintf(stderr,"Error: isotope must be specified for decay correction.\n");
return 1;
}
}
/* Check that isotope is specified with calfile */
if(calfile[0] && strcasecmp(isotope, "UNKNOWN")==0) {
fprintf(stderr,"Error: isotope must be specified for calibration.\n");
return 1;
}
/* Process the data */
if(verbose>0) printf("processing the data\n");
switch(format) {
case SCANDITRONICS:
correct_scanditronics(&blood, verbose-2);
break;
case ALLOGG:
if(correct_allogg(&blood) == 1) return 3;
break;
case ALLOGG2:
if(correct_allogg2(&blood, verbose-2) == 1) return 3;
break;
}
/* Check that neither of the channels has all values zero*/
if(verbose>1) printf("verifying that data is not all zeroes\n");
zero1=zero2=1;
if(format==SCANDITRONICS) {
ptr=blood.data;
while(ptr) {
if(fabs(ptr->results.ch1) > 1.0e-10) zero1=0;
if(fabs(ptr->results.ch2) > 1.0e-10) zero2=0;
ptr=ptr->next;
}
if(zero1){
fprintf(stderr, "Error: All values from channel 1 are zeros. Contact the physicist!\n");
if(column!=2) return 4;
}
if(zero2){
fprintf(stderr, "Error: All values from channel 2 are zeros. Contact the physicist!\n");
if(column!=1) return 4;
}
}
/* Subtract background radiation; must be done before decay correction */
if(blood.background>0.0) {
if(verbose>=0) printf("subtracting background\n");
subtract_background(&blood);
} else if(blood.background<-1.0) {
if(format==ALLOGG || format==ALLOGG2) {
if(calfile[0]) {
fprintf(stderr, "Error: %s does not contain background.\n", pumpfile);
return(5);
} else
fprintf(stderr, "Warning: %s does not contain background.\n", pumpfile);
}
}
/*
* Calibration
*/
if(calfile[0]) {
/* Check that we have the measurement date */
if(blood.date.tm_year==0) {
fprintf(stderr, "\nError: Blood file must contain the measurement date!\n");
fprintf(stderr, "Open the blood file in a text editor and add the date to\n");
fprintf(stderr, "the beginning as its own comment line, e.g. # 1999-11-31\n\n");
return(6);
}
if(verbose>0) printf("calibrating\n");
ptr=blood.data;
if(get_calfactors(calfile, timegm(&blood.date), &wc_factor, &pump_factor,
pump, isotope, caldate, verbose-1))
// if(get_calfactors(calfile,mktime(&blood.date),&wc_factor,&pump_factor,
// pump,isotope,caldate))
return 7;
if(verbose>=0) printf(" PC := %g\n WC/BR := %g\n calibration_date := %s\n",
pump_factor, wc_factor, caldate);
while(ptr) {
ptr->results.ch1*=wc_factor*pump_factor;
ptr->results.ch2*=wc_factor*pump_factor;
ptr=ptr->next;
}
} else if(verbose>=0) {
printf("Warning: no calibration file was entered; data is not calibrated.\n");
}
/* Correct the start time for decay correction */
/* If negative, then correct them back again later */
/* Set earlier activities to 0 */
if(startTime!=0.0) {
if(verbose>1) printf("correcting the start time\n");
ptr=blood.data;
while(ptr) {
ptr->results.t+=startTime;
if(format==SCANDITRONICS)
ptr->scanditronics.time_since_start+=startTime;
else
ptr->allogg.time_since_start+=startTime;
if(ptr->results.t<0.0) {
ptr->results.ch1=0;
ptr->results.ch2=0;
}
ptr=ptr->next;
}
}
/* Print the time of the first sample */
if(format==SCANDITRONICS) {
ptr=blood.data; t_samp=ptr->scanditronics.time_of_day;
if(gmtime_r(&t_samp, &tm_samp)==NULL) {
fprintf(stderr, "Warning: invalid sample time.\n");
} else {
tm_samp.tm_year=blood.date.tm_year;
tm_samp.tm_mon=blood.date.tm_mon;
tm_samp.tm_mday=blood.date.tm_mday;
t_samp=timegm(&tm_samp);
char buf[32];
if(!ctime_r_int(&t_samp, buf)) {
fprintf(stderr, "Warning: invalid sample time.\n");
} else {
fprintf(stdout, "Time point of the first sample is: %s\n", buf);
}
}
}
/* Decay correction */
if(decay_correct!=0) {
double lambda=hl2lambda(halflife*60.0);
double dc=0;
ptr=blood.data;
if(verbose>=0) printf("decay correction with half-life %.2f min.\n", halflife);
while(ptr) {
dc=1.0;
if(format==SCANDITRONICS) {
if(ptr->scanditronics.time_since_start>=0.0) {
dc = lambda2factor(lambda,
ptr->scanditronics.time_since_start, ptr->scanditronics.meas_interval);
}
} else {
if(ptr->allogg.time_since_start>=0.0) {
dc = lambda2factor(lambda, ptr->allogg.time_since_start, ptr->allogg.meas_interval);
if(verbose>5) printf("dc:%f %f\n",dc,ptr->allogg.time_since_start);
} }
ptr->results.ch1*=dc;
ptr->results.ch2*=dc;
ptr=ptr->next;
}
}
#if(0) // REMOVED IN BLO2KBQ VERSION 3.4.3
/* If start time is negative, then set orig times */
if(startTime<0.0) {
ptr=blood.data;
while(ptr) {
ptr->results.t-=startTime;
ptr=ptr->next;
}
}
#endif
/* Blood density correction */
if(density!=1.0) {
if(verbose>0)
printf("kBq/mL to kBq/g correction with blood density of %.2f g/mL.\n", density);
ptr=blood.data;
while(ptr) {
ptr->results.ch1/=density;
ptr->results.ch2/=density;
ptr=ptr->next;
}
}
/* Check that channels are not too different */
ptr=blood.data; sum1=sum2=0.0;
while(ptr) {
if(format==ALLOGG || format==ALLOGG2) {
if(ptr->results.ch1 < -2.0)
printf("Warning: channel below zero:%lf\n", ptr->results.ch1);
} else {
if(ptr->results.ch1 < 0)
printf("Warning: channel one below zero:%lf\n", ptr->results.ch1);
if(ptr->results.ch2 < 0)
printf("Warning: channel two below zero:%lf\n", ptr->results.ch2);
}
sum1+=ptr->results.ch1;
sum2+=ptr->results.ch2;
ptr=ptr->next;
}
if(column!=1 && column!=2 && (fabs(sum1-sum2)/(1.0+sum1+sum2))>0.3)
fprintf(stderr, "Warning: Large difference between channels of online blood sampler!\n");
/*
* Write output
*/
/* Create output filename, if not given */
if(!outfile[0]) {
if(studynr_validity_check2(studynr, 1)) strcpy(outfile, studynr);
else strcpy(outfile, pumpfile);
cptr=strrchr(outfile,'.'); if(cptr!=NULL) *cptr='\0';
if(calfile[0]) strcat(outfile,".kbq"); else strcat(outfile,".dat");
if(verbose>1) printf("outputfile := %s\n", outfile);
}
/* Check if output file exists; backup, if necessary */
if(backupExistingFile(outfile, NULL, tmp)) {
fprintf(stderr, "Error: %s\n", tmp);
}
/* Open the file */
if(verbose>=0) printf("writing output in file %s\n", outfile);
fp=fopen(outfile, "w");
if(fp==NULL) { fprintf(stderr, "Error: cannot write %s\n", outfile);
perror("fopen");
return 11;
}
/*
* Write header information
*/
if(verbose>2) printf("writing 'header'\n");
/* measurement start time */
if(gmtime_r(&t_injection, &tm_samp)==NULL) {
fprintf(stderr, "Warning: invalid measurement start time.\n");
}
#if(0)
} else {
/* Greenwich main time->Hki time zone*/
if(blood.date.tm_hour==tm_samp.tm_hour-2)
blood.date.tm_hour+=2;
else if(blood.date.tm_hour==tm_samp.tm_hour-3)
blood.date.tm_hour+=3;
}
#endif
fprintf(fp, "# measurement_start_time := %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
blood.date.tm_year+1900,blood.date.tm_mon+1,blood.date.tm_mday,
blood.date.tm_hour,blood.date.tm_min,blood.date.tm_sec);
/* detector id */
fprintf(fp, "# detector := %d\n", pump);
/* sample time unit */
if(timeInMinutes==1) fprintf(fp, "# time_unit := min\n");
else fprintf(fp, "# time_unit := sec\n");
/* calibration */
if(calfile[0]) {
fprintf(fp, "# calibrated := yes\n");
fprintf(fp, "# unit := kBq/ml\n");
fprintf(fp, "# calibration_date := %s\n", caldate);
fprintf(fp, "# pump_calibration_factor := %g\n", pump_factor);
fprintf(fp, "# wc_calibration_factor_per_br := %g\n", wc_factor);
fprintf(fp, "# calibration_file := %s\n", calfile);
} else {
fprintf(fp, "# calibrated := no\n");
fprintf(fp, "# unit := cps\n");
}
/* isotope */
fprintf(fp, "# isotope := %s\n", isotope);
/* decay correction */
if(decay_correct!=0) {
fprintf(fp, "# decay_correction := yes\n");
fprintf(fp, "# halflife := %g [min]\n", halflife);
} else {
fprintf(fp, "# decay_correction := no\n");
}
/* injection time */
r=3600*blood.date.tm_hour + 60*blood.date.tm_min + blood.date.tm_sec;
if(startTime==0.0) { /* injection time = measurement start time */
fprintf(fp, "# injection_time := %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
blood.date.tm_year+1900, blood.date.tm_mon+1, blood.date.tm_mday,
blood.date.tm_hour, blood.date.tm_min, blood.date.tm_sec);
} else if(r==0) {
fprintf(fp, "# injection_time := %g [sec]\n", startTime);
} else {
t_injection=timegm(&blood.date);
if(t_injection==-1) {
fprintf(stderr, "Error: invalid injection time.\n");
fclose(fp); return(12);
} else {
t_injection-=startTime;
if(gmtime_r(&t_injection, &tm_injection)==NULL) {
fprintf(stderr, "Error: invalid injection time.\n");
fclose(fp); return(12); }
}
fprintf(fp, "# injection_time := %.4d-%.2d-%.2d %.2d:%.2d:%.2d\n",
tm_injection.tm_year+1900, tm_injection.tm_mon+1, tm_injection.tm_mday,
tm_injection.tm_hour, tm_injection.tm_min, tm_injection.tm_sec);
}
/* background activity */
if(blood.background>0.0) {
if(blood.bgdate[0])
fprintf(fp,"# background_date := %s\n", blood.bgdate);
fprintf(fp,"# background := %g [cps]\n", blood.background);
}
/* Study number */
if(studynr_validity_check2(studynr, 1)) {
fprintf(fp, "# studynr := %s\n", studynr);
}
/*
* Write the data lines
*/
if(verbose>2) printf("Writing data lines\n");
ptr=blood.data;
while(ptr) {
/* Don't print if sample time is negative */
if(ptr->results.t<0.0) {ptr=ptr->next; continue;}
/* Write the sample time */
if(timeInMinutes) fprintf(fp, "%12.4f", ptr->results.t/60.0);
else fprintf(fp, "%12.3f", ptr->results.t);
/* Write the channel avg, or both separately, or only one */
if(both_cols) {
fprintf(fp, " %12.3f %12.3f\n", ptr->results.ch1, ptr->results.ch2);
} else if(column==1) {
fprintf(fp, " %12.3f\n", ptr->results.ch1);
} else if(column==2) {
fprintf(fp, " %12.3f\n", ptr->results.ch2);
} else {
fprintf(fp, " %12.3f\n", 0.5*(ptr->results.ch1+ptr->results.ch2));
}
ptr=ptr->next;
}
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/** Try to find the date from a comment line.
Works on both Scanditronics and Allogg formats.
@return Returns 1 if not found
*/
int find_date(
/** String to process. */
char *str,
/** Data structure to write date and time into. */
struct Blood *blood,
/** Verbose prints; set to 0 to print nothing in stdout. */
int verbose
) {
int c,yyyy,MM,dd,hh=0,mm=0,ss=0;
if(str[6]=='-' && (str[8]=='-'||str[9]=='-')) {
c=sscanf(str,"# %d-%d-%d %d:%d:%d",&yyyy,&MM,&dd,&hh,&mm,&ss);
/* Found all fields */
if(c==6 || c==3) {
blood->date.tm_year=yyyy-1900;
blood->date.tm_mday=dd;
blood->date.tm_mon=MM-1;
blood->date.tm_hour=hh;
blood->date.tm_min=mm;
blood->date.tm_sec=ss; blood->date.tm_isdst=-1;
if(verbose>0)
printf("find_date(): %04d-%02d-%02d %02d:%02d:%02d\n", yyyy, MM, dd, hh, mm, ss);
return 0;
}
}
return 1;
}
/*****************************************************************************/
/*****************************************************************************/
/** Try to find the isotope half-life from a comment line.
Works at least on the new Scanditronics format.
@return Returns 1 if not found
*/
int find_halflife(
/** String to process. */
char *str,
/** Verbose prints; set to 0 to print nothing in stdout. */
int verbose
) {
char *cptr=strcasestr(str, "half-life: ");
if(cptr!=NULL) {
cptr+=11;
} else {
cptr=strcasestr(str, "halflife: ");
if(cptr!=NULL) {
cptr+=10;
} else {
return(1);
}
}
if(verbose>0) printf("%s(): %s\n", __func__, cptr);
double f=atof_dpi(cptr); if(!isnormal(f)) return(1);
if(!isnan(halflife) && halflife>0.0 && !doubleMatch(halflife, f, 0.1)) {
fprintf(stderr, "Error: different isotope in file and given by user.\n");
exit(1);
}
halflife=f;
int c=hlIsotopeFromHalflife(halflife); if(c<0) return(1);
strcpy(isotope, hlIsotopeCode(c));
return(1);
}
/*****************************************************************************/
/*****************************************************************************/
/** Try to find the background radiation comment line.
@return Returns 1 if not found
*/
int find_background(
/** String to process. */
const char *str,
/** Data structure to write result into. */
struct Blood *blood
) {
char tmp[11];
while(isspace((int)*str) || *str=='#') str++;
if(*str) {
if(strncasecmp(str, "Background:", 11)==0) {
str+=11;
if(sscanf(str, "%f", &blood->background) == 1) {
return 0;
}
} else if(strncasecmp(str, "Background on", 13)==0) {
str+=13;
if(sscanf(str,"%8s: %f", tmp, &blood->background) == 2) {
sprintf(blood->bgdate, "%4.4s-%2.2s-%2.2s", tmp, tmp+4, tmp+6);
return 0; } else if(sscanf(str, "%10s: %f", blood->bgdate, &blood->background)==2) {
return 0;
}
}
}
return 1;
}
/*****************************************************************************/
/*****************************************************************************/
/** Read and parse datafile in the format used by the Scanditronics Automated
Blood Measurement system. The program 'bsampler' outputs this format.
*/
int read_scanditronics_data(
/** Pointer to file name. */
const char *inputfile,
/** Data structure to write results into. */
struct Blood *blood,
/** Verbose prints; set to 0 to print nothing in stdout. */
int verbose
) {
char line[512],*str;
struct DataList *points;
struct BDataS *bdata;
int datenotfound=1;
int hlnotfound=1;
FILE *fp;
int n;
//time_t t;
if(verbose>0) printf("%s(%s)\n", __func__, inputfile);
blood->bgdate[0]='\0';
blood->background = 0.0;
fp=fopen(inputfile,"r");
if(!fp) {
fprintf(stderr,"Error: cannot open \"%s\".\n",inputfile);
perror("fopen");
return 1;
}
/*
t=0; if(gm_timez) blood->date=*gmtime(&t); else blood->date=*localtime(&t);
*/
points=malloc(sizeof(struct DataList));
points->prev=NULL; points->next=NULL;
bdata=&points->scanditronics;
memset(bdata,0,sizeof(struct BDataS));
blood->data=points;
/* Read data */
while(fgets(line,sizeof(line),fp)) {
str=line;
while(*str && isspace((unsigned char)*str)) str++;
if(*str=='#') { /* Comment line; find date */
if(datenotfound) datenotfound=find_date(str, blood, verbose);
if(hlnotfound) hlnotfound=find_halflife(str, verbose);
continue;
}
if(datenotfound) {
fprintf(stderr, "Error: valid measurement date was not found in header.\n");
fclose(fp); return 5;
}
if(*str=='*' || *str=='\0') continue;
n=sscanf(str, "%lf %lf %lf %i %i %i %i %i %i %i",
&bdata->time_of_day, &bdata->time_since_start, &bdata->meas_interval,
&bdata->pair1_coincid, &bdata->pair1_det1_cnts, &bdata->pair1_det2_cnts,
&bdata->pair2_coincid, &bdata->pair2_det1_cnts, &bdata->pair2_det2_cnts,
&bdata->aux_cnts);
if(n==10) {
points->next=malloc(sizeof(struct DataList));
points->next->prev=points; points=points->next;
points->next=NULL;
bdata=&points->scanditronics;
memset(bdata,0,sizeof(struct BDataS));
} else {
if(n>0 && n!=7) fprintf(stderr,"Warning: malformed line!:{%s}\n",str);
}
}
fclose(fp);
if(bdata->time_of_day==0 && bdata->meas_interval==0) {
/* Remove extra entry if any */
struct DataList *prev=points->prev;
free(points);
points=prev;
if(points) {
points->next=NULL;
} else {
blood->data=NULL;
}
}
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/** Read and parse datafile in the format used by the Allogg Automated Blood Sampling system.
The program 'bloodsampler' outputs this format.
@return Returns 0 if successful and >0 in case of an error.
*/
int read_allogg_data(
/** Pointer to file name. */
const char *inputfile,
/** Data structure to write results into. */
struct Blood *blood,
/** Verbose prints; set to 0 to print nothing in stdout. */
int verbose
) {
char line[512],*str;
struct DataList *points;
struct BDataA *bdata;
int datenotfound=1, bgnotfound=1, timenotfound=1;
FILE *fp;
int n;
//time_t t;
if(verbose>0) printf("%s(%s)\n", __func__, inputfile);
blood->bgdate[0] = '\0';
blood->background = -1.0E20;
fp=fopen(inputfile,"r");
if(!fp) {
fprintf(stderr,"Error: cannot open \"%s\".\n",inputfile);
perror("fopen");
return 1;
}
/*
t=0; if(gm_timez) blood->date=*gmtime(&t); else blood->date=*localtime(&t);
*/
memset(&blood->date,0,sizeof(struct tm));
points=malloc(sizeof(struct DataList));
points->prev=NULL; points->next=NULL;
bdata=&points->allogg;
memset(bdata,0,sizeof(struct BDataA));
blood->data=points;
/* Read data */
while(fgets(line,sizeof(line),fp)) {
str=line; while(isspace((unsigned char)*str)) str++;
if(*str=='#') { /* Comment line; find date */
if(datenotfound) datenotfound=find_date(str, blood, verbose-2);
if(bgnotfound) bgnotfound=find_background(str, blood);
continue;
}
if(*str=='\0') continue;
n=sscanf(str,"%lf %d %d", &bdata->time_since_start,&bdata->data[0],&bdata->data[1]);
if(timenotfound && strlen(str)>6&&strlen(str)<9) {
blood->date.tm_hour=(str[0]-'0')*10+(str[1]-'0');
blood->date.tm_min=(str[2]-'0')*10+(str[3]-'0');
blood->date.tm_sec=(str[4]-'0')*10+(str[5]-'0');
blood->date.tm_isdst=-1; timenotfound=0;
} else if(n<2) {
fprintf(stderr,"Error: malformed line %sCorrect syntax is hhmmss.\n",str);
return 2;
} else {
points->next=malloc(sizeof(struct DataList));
points->next->prev=points;
points=points->next;
points->next=NULL;
bdata=&points->allogg;
memset(bdata,0,sizeof(struct BDataA));
}
}
fclose(fp);
if(bdata->time_since_start==0) { /* Remove extra entry if any */
struct DataList *prev=points->prev;
free(points);
points=prev;
if(points) {
points->next=NULL;
} else {
blood->data=NULL;
}
}
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/** Read and parse datafile in the format used by the new Allogg Automated Blood Sampling system.
@return Returns 0 if successful and >0 in case of an error.
*/
int read_allogg2_data(
/** Name of file to be read */
char *inputfile,
/** Pointer to struct where data is written */
struct Blood *blood,
/** Verbose level */
int verbose,
/** Pointer to string where error message is written; NULL if not needed. */
char *msg
) {
IFT ift;
int ii, n, ret;
int singles, coincidents;
char tmp[1024];
double f, tas;
//time_t t;
struct DataList *points;
struct BDataA *bdata;
if(verbose>0) printf("%s(%s)\n", __func__, inputfile);
/* Check the arguments */
if(inputfile==NULL || blood==NULL || strlen(inputfile)<1) {
if(msg!=NULL) strcpy(msg, "program error");
return(1); }
blood->bgdate[0] = '\0';
blood->background = -1.0E20;
/* Read Allogg2 file into IFT struct */
iftInit(&ift); ret=iftRead(&ift, inputfile, 0, 0);
if(ret) {if(msg!=NULL) strcpy(msg, ift.status); iftEmpty(&ift); return(2);}
if(verbose>=12) {iftWrite(&ift, "stdout", 0); fflush(stdout);}
/* Check that this is TPC Allogg2 */
strcpy(tmp, "System ID"); ii=iftGet(&ift, tmp, 0);
if(ii<0) {
if(msg!=NULL) strcpy(msg, "wrong format");
iftEmpty(&ift); return(3);
}
if(strcmp(ift.item[ii].value, "ABSS09282")!=0) {
if(msg!=NULL) strcpy(msg, "wrong sampler ID");
iftEmpty(&ift); return(4);
}
/* Read the background */
strcpy(tmp, "//Average background counts"); ii=iftFindNthKey(&ift, tmp, 1, 0);
if(ii>=0) {
if(verbose>3) printf("key := '%s'\n", ift.item[ii].key);
/* Get background date from the key */
ret=isdate(ift.item[ii].key+28);
if(ret!=0) ret=isdate2(ift.item[ii].key+28, tmp);
if(ret!=0) ret=isdate3(ift.item[ii].key+28, tmp);
if(ret==0) {
if(verbose>1) printf("background_date := %s\n", tmp);
strcpy(blood->bgdate, tmp);
} else if(verbose>0) printf("background_date not valid\n");
/* Get background counts/sec from value */
if(verbose>1) printf("background_value := %s\n", ift.item[ii].value);
blood->background=atof(ift.item[ii].value);
} else if(verbose>0) printf("background not found\n");
/* Read the isotope (half-life) */
strcpy(tmp, "HalfTime"); ii=iftGet(&ift, tmp, 0);
if(ii>=0) {
f=atof(ift.item[ii].value);
ii=hlIsotopeFromHalflife(f);
/* Use it only if user did not specify the isotope */
if(strcmp(isotope, "UNKNOWN")==0) {
strcpy(isotope, hlIsotopeCode(ii));
halflife=hlFromIsotope(isotope);
if(verbose) printf("half-life := %g min\n", halflife);
}
}
/* Read study number */
strcpy(tmp, "Run number"); ii=iftGet(&ift, tmp, 0);
if(ii>=0) {
#if(1)
studynr_from_fname2(ift.item[ii].value, studynr, 0);
#else
int len, i, j;
strncpy(studynr, ift.item[ii].value, MAX_STUDYNR_LEN);
/* Remove possible extra tail */
len=strlen(studynr);
for(i=0; i<len; i++) if(!isalnum((int)studynr[i])) {strcpy(studynr, ""); len=0;}
for(i=0; i<len; i++) if(!isalpha((int)studynr[i])) break;
if(i<1 || i>5) {strcpy(studynr, ""); len=0;}
for(j=0; (i+j)<len; j++) if(!isdigit((int)studynr[i+j])) {studynr[i+j]=(char)0; break;}
if(j<1 || j>5) {strcpy(studynr, ""); len=0;}
#endif
if(verbose>0) printf("studynr := %s\n", studynr);
}
/* Find the first line which starts with date and time, that is, */
/* the first count data line */
for(ii=0; ii<ift.keyNr; ii++) if(!ift.item[ii].key[0]) {
if(isdatetime(ift.item[ii].value, NULL)!=0) continue;
if(verbose>7) printf("%s\n", ift.item[ii].value);
break;
}
/* Get start time from it */
/*
t=0; if(gm_timez) blood->date=*gmtime(&t); else blood->date=*localtime(&t);
*/
//memset(&blood->date,0,sizeof(struct tm));
//find_date(ift.item[ii].value, blood);
ret=get_datetime(ift.item[ii].value, &(blood->date), verbose-4);
if(ret!=0 && verbose>1) {
fprintf(stderr, "Error %d in reading date and time in %s\n", ret, inputfile);
}
/* Read the count data lines */
if(verbose>1) printf("reading count data\n");
/* Initiate the list */
points=malloc(sizeof(struct DataList));
points->prev=NULL; points->next=NULL;
bdata=&points->allogg;
memset(bdata,0,sizeof(struct BDataA));
blood->data=points;
n=0;
for(; ii<ift.keyNr; ii++) {
/* Stop if item key is found */
if(ift.item[ii].key[0]) break;
/* Ignore comment lines */
if(strncmp(ift.item[ii].value, "//", 2)==0) continue;
if(strncmp(ift.item[ii].value, "#", 1)==0) continue;
if(strncmp(ift.item[ii].value, ";", 1)==0) continue;
/* Stop if line does not start with valid date and time */
if(isdatetime(ift.item[ii].value, NULL)!=0) break;
/* Read time after start, singles and coincidents */
ret=sscanf(ift.item[ii].value+19, "%lf %d %d", &tas, &singles, &coincidents);
if(ret!=3) continue;
if(verbose>10 && n<4) printf(" %f %d\n", tas, coincidents);
/* Increase the list */
if(n>0) {
points->next=malloc(sizeof(struct DataList));
points->next->prev=points;
points=points->next;
points->next=NULL;
bdata=&points->allogg;
memset(bdata,0,sizeof(struct BDataA));
}
/* Copy to list */
bdata->time_since_start=tas;
bdata->data[0]=coincidents;
bdata->data[1]=0;
n++;
}
iftEmpty(&ift);
if(n<1) {if(msg!=NULL) strcpy(msg, "no valid count data"); return(9);}
if(msg!=NULL) strcpy(msg, "ok");
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/** Calculate mid times and duration corrected activity values in Scanditronics format.
*/
void correct_scanditronics(
/** Pointer to data structure. */
struct Blood *blood,
/** Verbose level. */ int verbose
) {
struct DataList *data;
struct tm date; //*date;
time_t tim;
if(verbose>0) printf("%s()\n", __func__);
data=blood->data;
tim=data->scanditronics.time_of_day;
if(verbose>2) {
char tmp[32];
if(!ctime_r_int(&tim, tmp)) strcpy(tmp, "1900-01-01 00:00:00");
printf(" %s\n", tmp);
}
if(verbose>3) {
char buf[32];
strftime(buf, 32, "%Y-%m-%d %H:%M:%S", &blood->date);
printf(" blood->date: %s\n", buf);
}
/* Set time if it was not properly set in the comments */
if(pump==2 || (blood->date.tm_hour==0 && blood->date.tm_min==0 && blood->date.tm_sec==0)) {
//if(gm_timez) date=gmtime(&tim); else date=localtime(&tim);
if(gmtime_r(&tim, &date)!=NULL) {
blood->date.tm_hour=date.tm_hour;
blood->date.tm_min=date.tm_min;
blood->date.tm_sec=date.tm_sec;
} else {
fprintf(stderr, "Warning: invalid date in file.\n");
}
}
while(data) {
struct BDataS *d=&data->scanditronics;
d->results.t=d->time_since_start+0.5*d->meas_interval;
d->results.ch1=(double)d->pair1_coincid;
d->results.ch2=(double)d->pair2_coincid;
if(d->meas_interval>0) {
d->results.ch1/=(double)d->meas_interval;
d->results.ch2/=(double)d->meas_interval;
}
data=data->next;
}
}
/*****************************************************************************/
/*****************************************************************************/
/** Calculate measurement interval for Allogg data and set result values.
@returns 0 on success, 1 otherwise.
*/
int correct_allogg(
/** Data structure for correction. */
struct Blood *blood
) {
struct DataList *data;
double meas_interval=0.0;
data=blood->data;
if(data->next)
meas_interval=data->next->allogg.time_since_start-data->allogg.time_since_start;
else {
printf("Error: Only one blood data line found, cannot calculate frame time.\n");
return 1;
}
while(data) {
struct BDataA *d=&data->allogg;
d->meas_interval=meas_interval;
d->results.t=d->time_since_start-0.5*meas_interval; /* '+'->'-' by VO */
#if 0
d->results.ch1=(d->data[0]-d->data[1])/d->meas_interval;
d->results.ch2=(d->data[0]-d->data[1])/d->meas_interval;#else /* Use only second discriminator */
if(meas_interval == 0) {
printf("Error: zero measuring interval found, would lead to divide by zero\n");
return 1;
}
d->results.ch1=d->data[1]/d->meas_interval;
d->results.ch2=d->data[1]/d->meas_interval;
#endif
data=data->next;
}
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/** Calculate measurement interval for allogg2 data and set result values.
@return Returns 0 on success, <>0 otherwise
*/
int correct_allogg2(
/** Blood data structure for correction. */
struct Blood *blood,
/** Verbose level. */
int verbose
) {
if(verbose>0) printf("%s()\n", __func__);
struct DataList *data;
double meas_interval=0.0;
data=blood->data;
if(data->next) {
meas_interval=data->next->allogg.time_since_start-data->allogg.time_since_start;
if(verbose>1) {
printf("time_since_start := %g\n", data->allogg.time_since_start);
printf("time_since_start_next := %g\n", data->next->allogg.time_since_start);
printf("data0 := %d\n", data->allogg.data[0]);
printf("data0 next := %d\n", data->next->allogg.data[0]);
}
} else {
printf("Error: Only one blood data line found, cannot calculate frame time.\n");
return 1;
}
while(data) {
struct BDataA *d=&data->allogg;
d->meas_interval=meas_interval;
if(meas_interval == 0) {
printf("Error: zero measuring interval found, would lead to divide by zero\n");
if(verbose>1) {
printf("time_since_start := %g\n", data->allogg.time_since_start);
printf("time_since_start_next := %g\n", data->next->allogg.time_since_start);
}
return 1;
}
d->results.t=d->time_since_start+0.5*meas_interval;
d->results.ch1=d->data[0]/d->meas_interval;
d->results.ch2=d->data[0]/d->meas_interval;
data=data->next;
}
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/** Subtract background radiation values */
void subtract_background(
/** Blood data structure for correction. */
struct Blood *blood
) {
struct DataList *data=blood->data;
while(data) { struct BDataResults *res=&data->results;
res->ch1 -= blood->background;
res->ch2 -= blood->background;
data=data->next;
}
}
/*****************************************************************************/
/*****************************************************************************/
/** Reads calibration factors from file
@param file absolute or relative path to file
@param bloodt blood data measurement date
@param wcf for writing of well-counter/branching-ratio factor
@param pumpf for writing of pump calibration factor
@param pump pump number, in [1..3]
@param isotope used isotope
@param caldate for writing of calibration date
@param verbose verbose printing
@return Returns 0 on success, <>0 otherwise
*/
int get_calfactors(
const char *file,
time_t bloodt,
double *wcf,
double *pumpf,
int pump,
char *isotope,
char *caldate,
int verbose
) {
if(verbose>0) printf("%s(%s, ..., %d, %s, ...)\n", __func__, file, pump, isotope);
const time_t too_old = 60*60*24*31*6; /* 6 months */
struct tm calibrdate; //*calibrdate;
char line[512],*str;
int i,yyyy,mm,dd;
time_t calibrt, last_date=0;
char sdate[25];
FILE *fp;
strcpy(sdate, "");
if(pump<1) {
fprintf(stderr,"Warning: pump number %d too low, set to 1\n",pump);
pump=1;
} else if(pump>MAX_PUMPS) {
fprintf(stderr,"Warning: pump number %d too high, set to %d\n", pump, MAX_PUMPS);
pump=MAX_PUMPS;
}
fp=fopen(file,"r");
if(!fp) {
fprintf(stderr,"Error: cannot open file %s\n",file);
perror("fopen");
return 1;
}
while(fgets(line,sizeof(line),fp)) {
str=line;
while(*str && isspace((unsigned char)*str)) str++;
if(*str=='#' || *str=='\0') continue;
double p[MAX_PUMPS]; /* Pumps */
double w; /* Well-counter */
i=sscanf(str,"%d-%d-%d %lf %lf %lf %lf %lf",
&yyyy, &mm, &dd, &p[0], &p[1], &p[2], &p[3], &w);
if(i<5) {
printf("Error: Malformed line in calibration file, cannot continue\n");
fclose(fp); return 1; }
if(pump>i-4) {
printf("Error: no calibration for this blood sampler in %s\n", file);
fclose(fp); return 2;
}
if(i==5) {w=p[1];} else if(i==6) {w=p[2];} else if(i==7) {w=p[3];}
calibrdate.tm_mday=dd; calibrdate.tm_mon=mm-1;
calibrdate.tm_year=yyyy-1900;
calibrdate.tm_hour=calibrdate.tm_min=calibrdate.tm_sec=0;
calibrdate.tm_isdst=-1;
calibrt=timegm(&calibrdate);
if(calibrt==-1) calibrt=0;
if(calibrt>=last_date) last_date=calibrt;
/*
calibrt=time(NULL); calibrdate=localtime(&calibrt);
calibrdate->tm_mday=dd; calibrdate->tm_mon=mm-1;
calibrdate->tm_year=yyyy-1900;
calibrdate->tm_isdst=-1;
calibrt=mktime(calibrdate);
if(calibrt>last_date) last_date=calibrt;
*/
/* Check that calibration date is not later than measurement date */
if(calibrt>bloodt) break;
//if(calibrt!=(time_t)-1 && difftime(bloodt, calibrt)<0.0) break;
strftime(sdate, 24, "%Y-%m-%d", &calibrdate);
/* It was not, so set calibration factors and try the next line */
*pumpf=p[pump-1];
*wcf=w;
}
if(verbose>1) {
printf(" pump_factor: %g\n", *pumpf);
printf(" wc_factor: %g\n", *wcf);
}
/* Correct for branching ratio */
if(strcasecmp(isotope, "C-11")==0) *wcf/=BRANCHING_C;
else if(strcasecmp(isotope, "F-18")==0) *wcf/=BRANCHING_F;
else if(strcasecmp(isotope, "O-15")==0) *wcf/=BRANCHING_O;
else if(strcasecmp(isotope, "GE-68")==0) *wcf/=BRANCHING_Ge;
else if(strcasecmp(isotope, "GA-68")==0) *wcf/=BRANCHING_Ga;
else if(strcasecmp(isotope, "Cu-64")==0) *wcf/=BRANCHING_Cu64;
else fprintf(stderr, "Warning: branching factor not included in calibration.\n");
if(verbose>1) {
printf(" wc_factor_including_branching_correction: %g\n", *wcf);
}
fclose(fp);
/* Warn user if latest calibration date is older than 6 months */
if(bloodt>last_date && bloodt-last_date > too_old) {
fprintf(stderr, "Warning: check that %s contains the latest calibration coefficients.\n", file);
}
/* Set return value for calibration date */
if(caldate!=NULL) strcpy(caldate, sdate);
return 0;
}
/*****************************************************************************/
/*****************************************************************************/
/// @endcond