mkcalhdr.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <ctype.h>
#include <dirent.h>
#include <sys/file.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates HRRT calibration factor, assuming that there is a cylindrical",
  "phantom with 20 cm diameter in the FOV.",
  "Resulting calibration header file is written in file",
  "STUDY_RECMETHOD_iMNsXY_sS.cal.hdr, where STUDY is the part of the calibration",
  "study name preceding delimiter \"_\", \"-\", \".\"  or space, RECMETHOD is",
  "the utilised reconstruction method, iMN is the number of iterations preceded",
  "by char i, sXY is the number of subsets preceded by char s and sS is the span",
  "preceded by char s. If calibration is calculated for 128x128 matrix size,",
  "suffix m128 is appended in header file name.",
  " ",
  "Efficient (sensitivity) factor for each plane is calculated as mean over",
  "planes 19-188 divided by plane average. Efficient factor is set to zero for",
  "planes 1-6 and 201-207 (HRRT has 207 image planes).",
  " ",
  "Usage: @P [Options] -a=<Activity> Interfile(s)",
  " ",
  "Options:",
  " -a=<Activity> or -a <Activity>",
  "     Activity concentration at scan start time in Bq/mL; required option.",
  " -l=<First eff plane> or -l <First eff plane>",
  "     First plane for which the efficient factor is set; default is 7.",
  " -L=<Last eff plane> or -L <Last eff plane>",
  "     Last plane for which the efficient factor is set; default is 200.",
  " -m=<First mean plane> or -m <First mean plane>",
  "     First plane used for mean calculation; default is 19.",
  " -M=<Last mean plane> or M <Last mean plane>",
  "     Last plane used for mean calculation; default is 188.",
  " -sd=<Y|n>",
  "     Calculate standard deviation inside VOI (y, default), or not (n).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P -a=4120 calrx-0004-2006.3.1.15.16.23_em_3D_256i01.i",
  " ",
//  "See also: if2e7, run_if2e7",
//  " ",
  "Keywords: HRRT, Interfile, calibration",
  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 hrrtMakeCalHdr(
  char *imgName,
  int firstplane,
  int lastplane,
  int startplane,
  int endplane,
  float activityConcentration,
  int calcSTDV,
  char *status,
  int verbose
) {
  int i, j, k, ret;
  char buf[512];
 
  /* Check the input */
  if(status!=NULL) sprintf(status, "invalid filename");
  if(imgName==NULL || !imgName[0]) return(1);
  if(status!=NULL) sprintf(status, "invalid efficient plane settings");
  if(firstplane<1 || firstplane>lastplane) return(2);
  if(status!=NULL) sprintf(status, "invalid mean plane settings");
  if(startplane<1 || startplane>endplane) return(3);
  if(status!=NULL) sprintf(status, "invalid activity concentration");
 
  /*
   *  Read Interfile header
   */
  char hdrFile[FILENAME_MAX];
  char imgFile[FILENAME_MAX];
  strlcpy(imgFile, imgName, FILENAME_MAX);
  strlcpy(hdrFile, imgName, FILENAME_MAX-4);
  strcat(hdrFile, ".hdr"); 
  if(verbose>2) printf("reading header in %s\n", hdrFile);
  IFT ift; iftInit(&ift);
  ret=iftRead(&ift, hdrFile, 0, 0);
  /* If error, try to add .i to filenames */
  if(ret!=0) {
    strlcpy(hdrFile, imgName, FILENAME_MAX-4);
    strcat(hdrFile, ".i.hdr"); 
    if(verbose>2) printf("reading header in %s\n", hdrFile);
    ret=iftRead(&ift, hdrFile, 0, 0);
    if(ret==0) {
      /* That was good, therefore fix image filename too */
      strlcpy(imgFile, imgName, FILENAME_MAX);
      strcat(imgFile, ".i"); 
    }
  }
  /* In case of error, stop trying */
  if(ret!=0) {
    if(status!=NULL) sprintf(status, "cannot read file");
    return(11);
  }
 
  /* Matrix size */
  int dimx, dimy, dimz;
  if(iftGetIntValue(&ift, 0, "matrix size [1]", &dimx, 0)<0 || dimx<1) {
    if(status!=NULL) sprintf(status, "cannot read dimx");
    iftEmpty(&ift); return(21);
  }
  if(iftGetIntValue(&ift, 0, "matrix size [2]", &dimy, 0)<0 || dimy<1) {
    if(status!=NULL) sprintf(status, "cannot read dimy");
    iftEmpty(&ift); return(22);
  }
  if(iftGetIntValue(&ift, 0, "matrix size [3]", &dimz, 0)<0 || dimz<1) {
    if(status!=NULL) sprintf(status, "cannot read dimz");
    iftEmpty(&ift); return(23);
  }
  if(verbose>1) printf("  matrix_size := %d x %d x %d\n", dimx, dimy, dimz);
  /* Check plane limits against dimz */
  if(lastplane>dimz) {
    if(status!=NULL) sprintf(status, "invalid efficient plane settings");
    iftEmpty(&ift); return(24);
  }
  if(endplane>dimz) {
    if(status!=NULL) sprintf(status, "invalid mean plane settings");
    iftEmpty(&ift); return(25);
  }
  if(verbose>0) {
    printf("  Efficient factor set to 0 for %d first planes\n", firstplane-1);
    printf("  Efficient factor set to 0 for %d last planes\n", dimz-lastplane);
    printf("  Excluding %d first planes from the mean calculation\n", startplane-1);
    printf("  Excluding %d last planes from the mean calculation\n", dimz-endplane-1);
  }
 
  /* Pixel size */
  double xsize, ysize, zsize;
  if(iftGetDoubleValue(&ift, 0, "scaling factor (mm/pixel) [1]", &xsize, 0)<0 || xsize<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read x size");
    iftEmpty(&ift); return(31);
  }
  if(iftGetDoubleValue(&ift, 0, "scaling factor (mm/pixel) [2]", &ysize, 0)<0 || ysize<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read y size");
    iftEmpty(&ift); return(32);
  }
  if(iftGetDoubleValue(&ift, 0, "scaling factor (mm/pixel) [3]", &zsize, 0)<0 || zsize<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read z size");
    iftEmpty(&ift); return(33);
  }
  if(verbose>1) printf("  pixel_size := %g x %g x %g\n", xsize, ysize, zsize);
  if(fabs(xsize-ysize)>0.01) {
    if(status!=NULL) sprintf(status, "invalid x and y pixel sizes");
    iftEmpty(&ift); return(34);
  }
  /* Calculate ROI radius in pixels */
  float ROIradius;
  ROIradius=75.0/xsize; /* 15 cm diameter ROI = 75 mm Radius */
  if(verbose>0) printf("  ROIradius := %g\n", ROIradius);
 
  /* Scan duration */
  double scanDuration;
  if(iftGetDoubleValue(&ift, 0, "image duration", &scanDuration, 0)<0 || scanDuration<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read scan duration");
    iftEmpty(&ift); return(41);
  }
  if(verbose>1) printf("  scanDuration := %g\n", scanDuration);
 
  /* Branching fraction */
  double branchingFraction;
  if(iftGetDoubleValue(&ift, 0, "branching factor", &branchingFraction, 0)<0 || branchingFraction<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read branching fraction");
    iftEmpty(&ift); return(42);
  }
  if(verbose>1) printf("  branchingFraction := %g\n", branchingFraction);
 
  /* Halflife */
  double halflife;
  if(iftGetDoubleValue(&ift, 0, "isotope halflife", &halflife, 0)<0 || halflife<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read halflife");
    iftEmpty(&ift); return(43);
  }
  if(verbose>1) printf("  halflife := %g\n", halflife);
 
  /* Dead time */
  double deadtimeCorFactor;
  if(iftGetDoubleValue(&ift, 0, "Dead time correction factor", &deadtimeCorFactor, 0)<0
     || deadtimeCorFactor<=0.0) {
    if(status!=NULL) sprintf(status, "cannot read dead time");
    iftEmpty(&ift); return(44);
  }
  if(verbose>1) printf("  deadtimeCorFactor := %g\n", deadtimeCorFactor);
 
  /* Span */
  int span;
  if(iftGetIntValue(&ift, 0, "axial compression", &span, 0)<0) {
    if(status!=NULL) sprintf(status, "cannot read span");
    iftEmpty(&ift); return(45);
  }
  if(verbose>1) printf("  span := %d\n", span);
 
  /* Reconstruction method */
  char recMethod[256];
  i=iftGetFrom(&ift, 0, "method of reconstruction used", 0);
  if(i<0) i=iftGetFrom(&ift, 0, "reconstruction method", 0);
  if(i<0) {
    if(status!=NULL) sprintf(status, "cannot read reconstruction method");
    iftEmpty(&ift); return(46);
  }
  strlcpy(buf, ift.item[i].value, 256);
  for(i=j=0; i<(int)strlen(buf); i++) {
    if(buf[i]!='/' && buf[i]!=' ') recMethod[j++]=buf[i];
  }
  recMethod[j]=(char)0;
  if(verbose>1) printf("  reconstruction method := %s\n", recMethod);
 
  /* Nr of iterations (not always present) */
  int numIter;
  if(iftGetIntValue(&ift, 0, "number of iterations", &numIter, 0)<0) {
    numIter=0;
  }
  if(verbose>1) printf("  numIter := %d\n", numIter);
 
  /* Nr of subsets (not always present) */
  int numSubsets;
  if(iftGetIntValue(&ift, 0, "number of subsets", &numSubsets, 0)<0) {
    numSubsets=0;
  }
  if(verbose>1) printf("  numSubsets := %d\n", numSubsets);
 
  /* Check that image data is saved as 4-byte floats */
  i=iftGetFrom(&ift, 0, "number format", 0);
  if(i<0 || strcasecmp(ift.item[i].value, "float")!=0) {
    if(status!=NULL) sprintf(status, "data is not stored as floats");
    iftEmpty(&ift); return(51);
  }
  i=iftGetFrom(&ift, 0, "number of bytes per pixel", 0);
  if(i<0 || strcasecmp(ift.item[i].value, "4")!=0) {
    if(status!=NULL) sprintf(status, "data is not stored as 4-byte floats");
    iftEmpty(&ift); return(52);
  }
 
  /* free header data */
  iftEmpty(&ift);
 
 
  /*
   *  Read Interfile image data
   */
  if(verbose>2) printf("reading image data in %s\n", imgFile);
  float *img;
  unsigned int pxlNr, n;
  pxlNr=dimx*dimy*dimz;
  img=(float*)malloc(pxlNr*sizeof(float));
  if(img==NULL) {
    if(status!=NULL) sprintf(status, "out of memory");
    return(100);
  }
  FILE *fp;
  fp=fopen(imgFile, "rb");
  if(fp==NULL) {
    if(status!=NULL) sprintf(status, "cannot open image file");
    free(img); return(101);
  }
  n=fread(img, sizeof(float), pxlNr, fp);
  fclose(fp);
  if(n<pxlNr) {
    if(status!=NULL) sprintf(status, "cannot read image");
    free(img); return(102);
  }
  /* Swap bytes if necessary */
  if(little_endian()==0) {
    if(verbose>1) printf("  swapping bytes\n");
    for(n=0; n<pxlNr; n++) swawbip(&img[n], sizeof(float));
  }
 
  /*
   *  ROI evaluation
   */
  if(verbose>2) printf("ROI evaluation (circular, centered, 15 cm diameter)\n");
  int numpix=0, planenumpix, x, y, z;
  double sqrSum=0.0;
  double VoISTDV=0.0;      
  float ROIsum=0.0, ROIavg;
  float planesum, xsum, ysum, x0, y0;
  float planeROIavg[dimz], planeROIsum;
  for(i=0; i<dimz; i++) planeROIavg[i]=0.0;
  j=k=0; //dimx*dimy*(startplane-1);
  //for(z=startplane-1; z<endplane; z++) {
  for(z=0; z<dimz; z++) {
    planesum=xsum=ysum=0.0;
    /* centre of gravity */
    for(y=0; y<dimy; y++) {
      for(x=0; x<dimx; x++) {
        planesum+=img[j];
        xsum+=(float)x*img[j];
        ysum+=(float)y*img[j];
        j++;    
      }
    }
    x0=xsum/planesum; 
    y0=ysum/planesum; 
    if(!isnormal(x0) || (int)x0<=0 || (int)x0>=dimx) x0=0.5*dimx;
    if(!isnormal(y0) || (int)y0<=0 || (int)y0>=dimy) y0=0.5*dimy;
    if(verbose>3) printf("  ROIcenter[%d] = (%f/%f)\n", z, x0, y0);
    planeROIsum=0.0;
    planenumpix=0;
    for(y=0; y<dimy; y++) {
      for(x=0; x<dimx; x++) {
        if(hypotf((float)x-x0, (float)y-y0) <= ROIradius) {
          planenumpix++;
          planeROIsum+=img[k];
          if((z+1)>=startplane && (z+1)<=endplane) {
            /* Only these are included in the VOI mean */
            numpix++;
            ROIsum+=img[k];
          }
        }
        k++;
      }
    }
    planeROIavg[z]=planeROIsum/(float)planenumpix;
    if(verbose>0) printf("  VOIaverage[%d] := %g\n", z, planeROIavg[z]);
  }
  ROIavg=ROIsum/(float)numpix;
  if(verbose>2) printf("  ROIavg := %g\n", ROIavg);
 
  /* Calculate standard deviation inside VOI, if requested */
  if(calcSTDV!=0) {
    j=k=dimx*dimy*(startplane-1);
    for(z=startplane; z<endplane; z++) {
      planesum=xsum=ysum=0.0;
      /* centre of gravity */
      for(y=0; y<dimy; y++) {
        for(x=0; x<dimx; x++) {
          planesum+=img[j];
          xsum+=(float)x*img[j];
          ysum+=(float)y*img[j];
          j++;         
        }
      }
      x0=xsum/planesum;
      y0=ysum/planesum;
      for(y=0; y<dimy; y++) {
        for(x=0; x<dimx; x++) {
          if(hypotf((float)x-x0, (float)y-y0) <= ROIradius) {
            sqrSum+=(img[k]-ROIavg)*(img[k]-ROIavg);
          }
          k++;
        }
      }
    }
    VoISTDV=sqrt(sqrSum/(double)numpix);
    if(verbose>2) {
      fprintf(stdout, "  standard deviation in VoI := %f\n", VoISTDV);
      fflush(stdout);
    }
  }
 
  free(img);
 
 
  /*
   *  Calculate calibration factor
   */
  if(verbose>2) printf("Calculating calibration factor\n");
  /* normalise mean for duration of acquisition */
  double corROIavg=ROIavg/scanDuration;
  /* divide by branching fraction */
  corROIavg=corROIavg/branchingFraction;
  /* in frame decay correction */
  double dfrac=scanDuration/halflife*log(2);
  dfrac/=(1.0-exp(-dfrac));
  if(verbose>0) printf("  decay correction := %f\n", dfrac);
  corROIavg*=dfrac;
  /* deadtime correction */
  corROIavg*=deadtimeCorFactor;
  /* calculate calibration factor */
  float calibFactor=activityConcentration/corROIavg;
  if(verbose>2) {
    fprintf(stdout, "  calibration factor := %e (counts/s)/(Bq/ml)\n",
            calibFactor);
    fflush(stdout);
  }
 
  /*
   *  Write output file
   */
  if(verbose>2) printf("Writing output file\n");
  /* Construct output file name */
  char *cptr, calHdrName[512];
  strlcpy(buf, imgFile, 512);
  if((cptr=strpbrk(buf, "_-. "))!=NULL) *cptr='\0';
  strcpy(calHdrName, buf);
  strlcat(calHdrName, "_", 512); strlcat(calHdrName, recMethod, 512);
  if(numIter>0) {
    sprintf(buf, "_i%02d", numIter); strlcat(calHdrName, buf, 512);
  } else {
    strlcat(calHdrName, "_i0", 512);
  }
  sprintf(buf, "s%d", numSubsets); strlcat(calHdrName, buf, 512);
  sprintf(buf, "_s%d", span); strlcat(calHdrName, buf, 512);
  if(dimx!=256) {sprintf(buf, "_m%d", dimx); strlcat(calHdrName, buf, 512);}
  strlcat(calHdrName, ".cal.hdr", 512); 
  /* Fill IFT struct */
  ret=iftPut(&ift, NULL, imgName, ";", 0);
  if(ret==0)
    ret=iftPutDouble(&ift, "activity concentration used (Bq/ml)", 
                     (double)activityConcentration, ";", 0);
  if(ret==0) {
    sprintf(buf, "VOI mean from planes %d - %d", startplane, endplane);
    ret=iftPut(&ift, NULL, buf, ";", 0);
  }
  if(ret==0) ret=iftPutDouble(&ift, "mean activity concentration in VOI", (double)ROIavg, ";", 0);
  if(calcSTDV!=0) {
    if(ret==0) ret=iftPutDouble(&ift, "standard deviation in VOI", (double)VoISTDV, ";", 0);
    if(ret==0) ret=iftPutDouble(&ift, "CV in VOI", (double)VoISTDV/ROIavg, ";", 0);
  }
  if(ret==0)
    ret=iftPutDouble(&ift, "corrected VOI activity concentration", (double)corROIavg, ";", 0);
  if(ret==0) ret=iftPut(&ift, NULL, " ", NULL, 0);
  if(ret==0) ret=iftPut(&ift, NULL, "calibration factor in (Bq/cc)/(counts/s)", ";", 0);
  if(ret==0) ret=iftPutDouble(&ift, "calibration factor", (double)calibFactor, NULL, 0);
  if(ret==0) ret=iftPut(&ift, NULL, " ", NULL, 0);
  for(j=1; j<=dimz && ret==0; j++) {
    sprintf(buf, "efficient factor for plane %d", j-1);
    if(j<firstplane || j>lastplane)
      ret=iftPutDouble(&ift, buf, 0.0, NULL, 0);
    else
      ret=iftPutDouble(&ift, buf, ROIavg/planeROIavg[j-1], NULL, 0);
  }
  if(ret!=0) {
    if(status!=NULL) sprintf(status, "cannot make calibration file");
    iftEmpty(&ift); return(201);
  }
  if(verbose>0) {
    printf("\n--------------------------------------------------\n");
    printf("%s :\n", calHdrName);
    printf("--------------------------------------------------\n");
    iftWrite(&ift, "stdout", 0);
  }
  /* Write the file */
  ret=iftWrite(&ift, calHdrName, 0);
  iftEmpty(&ift);
  if(ret!=0) {
    if(status!=NULL) sprintf(status, "cannot write calibration file");
    return(202);
  }
  fprintf(stdout, "\ncalibrations written in %s\n\n", calHdrName);
  fflush(stdout);
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  int     ret, fileNr=0, firstfile=0;
  int     calcSTDV=1; // 0=SD not calculated; 1=SD calculated
  float   activityConcentration=0.0;
  int     firstplane, lastplane, startplane, endplane;
  char    imgfile[FILENAME_MAX];
  char   *cptr;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  firstplane=7; lastplane=200; startplane=19; endplane=188;
  imgfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "SD=", 3)==0) {
      cptr+=3;
      if(strncasecmp(cptr, "YES", 1)==0 || strcasecmp(cptr, "ON")==0) {
        calcSTDV=1; continue;
      }
      if(strncasecmp(cptr, "NO", 1)==0 || strcasecmp(cptr, "OFF")==0) {
        calcSTDV=0; continue;
      }
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      double v;
      ret=atof_with_check(cptr+2, &v);
      if(ret==0 && v>0.0) {activityConcentration=v; continue;}
    } else if(strcasecmp(cptr, "A")==0 && ai<argc-1) {
      ai++; cptr=argv[ai];
      double v;
      ret=atof_with_check(cptr, &v);
      if(ret==0 && v>0.0) {activityConcentration=v; continue;}
    } else if(strncmp(cptr, "l=", 2)==0) {
      ret=atoi_with_check(cptr+2, &firstplane);
      if(ret==0 && firstplane>0) {continue;}
    } else if(strcmp(cptr, "l")==0 && ai<argc-1) {
      ai++; cptr=argv[ai];
      ret=atoi_with_check(cptr, &firstplane);
      if(ret==0 && firstplane>0) {continue;}
    } else if(strncmp(cptr, "L=", 2)==0) {
      ret=atoi_with_check(cptr+2, &lastplane);
      if(ret==0 && lastplane>0) {continue;}
    } else if(strcmp(cptr, "L")==0 && ai<argc-1) {
      ai++; cptr=argv[ai];
      ret=atoi_with_check(cptr, &lastplane);
      if(ret==0 && lastplane>0) {continue;}
    } else if(strncmp(cptr, "m=", 2)==0) {
      ret=atoi_with_check(cptr+2, &startplane);
      if(ret==0 && startplane>0) {continue;}
    } else if(strcmp(cptr, "m")==0 && ai<argc-1) {
      ai++; cptr=argv[ai];
      ret=atoi_with_check(cptr, &startplane);
      if(ret==0 && startplane>0) {continue;}
    } else if(strncmp(cptr, "M=", 2)==0) {
      ret=atoi_with_check(cptr+2, &endplane);
      if(ret==0 && endplane>0) {continue;}
    } else if(strcmp(cptr, "M")==0 && ai<argc-1) {
      ai++; cptr=argv[ai];
      ret=atoi_with_check(cptr, &endplane);
      if(ret==0 && endplane>0) {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 */
  for(; ai<argc; ai++) {
    if(firstfile==0) firstfile=ai;
    fileNr++;
  }
  /* Did we get all the information that we need? */
  if(fileNr==0) {tpcPrintUsage(argv[0], info, stderr); return(1);}
 
  /* Are options set to reasonable values? */
  if(firstplane>=lastplane) {
    fprintf(stderr, "Error: invalid efficient plane settings.\n");
    return(1);
  }
  if(startplane>=endplane) {
    fprintf(stderr, "Error: invalid mean plane settings.\n");
    return(1);
  }
  if(activityConcentration<=0.0) {
    fprintf(stderr, "Error: missing activity concentration.\n");
    return(1);
  }
  if(activityConcentration<=50.0) {
    fprintf(stderr, "Warning: activity concentration must be given in Bq/mL.\n");
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>3) {
    printf("activityConcentration := %g\n", activityConcentration);
    printf("firstplane := %d\n", firstplane);
    printf("lastplane := %d\n", lastplane);
    printf("startplane := %d\n", startplane);
    printf("endplane := %d\n", endplane);
    printf("fileNr := %d\n", fileNr);
  }
 
 
  /*
   *  Process each file
   */
  if(verbose>0) printf("processing %d file(s)...\n", fileNr);
  char buf[256];
  for(ai=firstfile; ai<argc; ai++) {
    strlcpy(imgfile, argv[ai], FILENAME_MAX); 
    if(verbose>1) printf("processing %s\n", imgfile);
    ret=hrrtMakeCalHdr(imgfile, firstplane, lastplane, startplane, endplane,
                       activityConcentration, calcSTDV, buf, verbose-1);
    if(ret!=0) {
      fprintf(stderr, "Error: %s\n", buf);
      return(100+ret);
    }
  }
 
  if(verbose>0) printf("... done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/