mkcalhdr.c File Reference

Calculation of quantification headers for the HRRT. More...

Go to the source code of this file.

Functions
int	hrrtMakeCalHdr (char *imgName, int firstplane, int lastplane, int startplane, int endplane, float activityConcentration, int calcSTDV, char *status, int verbose)

Detailed Description

Calculation of quantification headers for the HRRT.

Calculates calibration factor counts/s -> Bq/ml and slice sensitivity coefficients.

The procedure is as follows:

1. get mean value of a circular ROI centered on a HRRT interfile image of a cylindric phantom (for all 207 planes)
2. calculate the average of plane 19-189
   
3. normalise for duration of acquisition (get duration from header)
4. devide by branching factor (get factor from header)
5. apply in-frame decay correction
6. apply deadtime correction (factor from header)
7. get activity concentration in the phantom from command line
8. calculate calibration factor ((Bq/ml)/(counts/s))
9. calculate slice sensitivity coeficients
10. generate output file name according to reconstruction parameters
11. create slice sensitivity file and write output to it.

Copyright

(c) Turku PET Centre

Author

Roman Krais, Jarkko Johansson, Timo Laitinen, Vesa Oikonen

Definition in file mkcalhdr.c.

Function Documentation

hrrtMakeCalHdr()

int hrrtMakeCalHdr	(	char *	imgName,
		int	firstplane,
		int	lastplane,
		int	startplane,
		int	endplane,
		float	activityConcentration,
		int	calcSTDV,
		char *	status,
		int	verbose )

Read transmission scan file and blank scan file and calculate the attenuation correction data and the logarithmic correction data.

Returns

Returns 0 if ok.

Parameters

imgName	Pointer to Interfile name, with or without .i extension.
firstplane	First plane for which the efficient factor is set.
lastplane	Last plane for which the efficient factor is set.
startplane	First plane used for mean calculation.
endplane	Last plane used for mean calculation.
activityConcentration	Activity concentration at scan start time in Bq/mL.
calcSTDV	Calculate (<>0) or do not calculate (0) SD.
status	Pointer to a string (allocated for at least 128 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 105 of file mkcalhdr.c.

  {
  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);
}