imgbkgrm.c File Reference

Set pixels in dynamic PET image to zero if their TACs seem to contain only noise. More...

Go to the source code of this file.

Functions
int	imgBorderAverageTAC (IMG *img, float *tac, float *sd, float *med)
int	imgNoiseTemplate (IMG *img, int keepNegat, int border, IMG *timg, int verbose, char *status)

Detailed Description

Set pixels in dynamic PET image to zero if their TACs seem to contain only noise.

Remarks

Application name was previously ebkgrm.

Copyright

(c) Turku PET Centre

Author

Vesa Oikonen

Definition in file imgbkgrm.c.

Function Documentation

imgBorderAverageTAC()

int imgBorderAverageTAC	(	IMG *	img,
		float *	tac,
		float *	sd,
		float *	med )

Calculates an average time-activity curve of all pixels at the x,y-border in the specified IMG data.

Border values can often be assumed to consist of background noise only. Z-borders are not included because imaging target usually extends over z-borders.

See also

imgAverageTAC, imgAverageMaskTAC

Returns

Returns 0, if ok.

Parameters

img	Pointer to the image data.
tac	Pointer to an array, allocated for the size of img->dimt, where border pixel mean values will be written; enter NULL if not needed.
sd	Pointer to an array, allocated for the size of img->dimt, where border pixel s.d. values will be written; enter NULL if not needed.
med	Pointer to an array, allocated for the size of img->dimt, where border pixel median values will be written; enter NULL if not needed.

Definition at line 85 of file imgbkgrm.c.

  {
  if(img==NULL || img->status<IMG_STATUS_OCCUPIED) return(1);
  if(tac==NULL && sd==NULL && med==NULL) return(0);
 
  /* compute the TAC */
  int pxlNr=img->dimz*(2*img->dimx+2*(img->dimy-2));
  float a[pxlNr], fm, fs;
  for(int fi=0; fi<img->dimt; fi++) {
    int i=0;
    for(int zi=0; zi<img->dimz; zi++) {
      int yi, xi;
      yi=0;
      for(xi=0; xi<img->dimx; xi++) a[i++]=img->m[zi][yi][xi][fi];
      yi=img->dimy-1;
      for(xi=0; xi<img->dimx; xi++) a[i++]=img->m[zi][yi][xi][fi];
      xi=0;
      for(yi=1; yi<img->dimy-1; yi++) a[i++]=img->m[zi][yi][xi][fi];
      xi=img->dimx-1;
      for(yi=1; yi<img->dimy-1; yi++) a[i++]=img->m[zi][yi][xi][fi];
    }
    if(tac!=NULL || sd!=NULL) {
      fm=fmean(a, pxlNr, &fs);
      if(tac!=NULL) tac[fi]=fm;
      if(sd!=NULL) sd[fi]=fs;
    }
    if(med!=NULL) med[fi]=fmedian(a, pxlNr);
  }
  return(0);
}

Referenced by imgNoiseTemplate().

imgNoiseTemplate()

int imgNoiseTemplate	(	IMG *	img,
		int	keepNegat,
		int	border,
		IMG *	timg,
		int	verbose,
		char *	status )

Creates a mask (template) image based on whether pixel TAC is closer to scaled mean TAC of the image data than to the zero, based on AICs.

Pixel considered as noise will lead to mask value 0, otherwise 1.

See also

imgThresholdByMask, imgRawCountsPerTime

Returns

Returns 0 if ok.

Parameters

img	Original dynamic image; if sinogram data, then divide by frame lengths first. If data does not contain frame times, then equal frame times are assumed.
keepNegat	Process pixels with AUC<0 normally (1), or consider as noise (0).
border	Set to 0 to compare pixels TACs to zero, or set to <>0 to compare to background TAC calculated as median of image x,y-border pixels.
timg	Mask image; if empty, then it will be allocated here; if pre-allocated, then mask pixel value changed to 0 when necessary, but 0 is never changed to 1.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.

Definition at line 137 of file imgbkgrm.c.

  {
  if(verbose>0) printf("imgNoiseTemplate(*img, %d, %d, *timg, %d, status)\n",
                       keepNegat, border, verbose);
  if(status!=NULL) strcpy(status, "fault in calling routine");
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(timg==NULL) return(1);
  if(timg->status==IMG_STATUS_OCCUPIED) {
    if(img->dimx!=timg->dimx) return(2);
    if(img->dimy!=timg->dimy) return(2);
    if(img->dimz!=timg->dimz) return(2);
  }
  if(img->dimt<2) {
    if(status!=NULL) strcpy(status, "method works only for dynamic data");
    return(3);
  }
 
  int ret, n;
 
  /*
   *  Compute the mean TAC of all pixels
   */
  n=img->dimt;
  float fdur[n], meany[n], bkgtac[n];
 
  if(verbose>1) fprintf(stdout, "computing mean of pixels\n");
  if(imgExistentTimes(img)) {
    for(int i=0; i<n; i++) fdur[i]=img->end[i]-img->start[i];
  } else {
    for(int i=0; i<n; i++) fdur[i]=1.0;
  }
  ret=imgAverageTAC(img, meany);
  if(ret) {
    if(status!=NULL) strcpy(status, "cannot calculate average TAC");
    if(verbose>1) printf("ret := %d\n", ret);
    return(4);
  }
 
  if(verbose>1) fprintf(stdout, "computing AUC of the mean TAC\n");
  float mean_tac_auc=0.0;
  for(int i=0; i<n; i++) mean_tac_auc+=fdur[i]*meany[i];
  if(verbose>2) printf("mean_tac_auc := %g\n", mean_tac_auc);
  if(fabs(mean_tac_auc)<1.0E-10) {
    if(status!=NULL) strcpy(status, "very small TAC average");
    return(4);
  }
 
  /* If requested, compute the background TAC as the mean of x,y-border pixels */
  double bkg_tac_auc=0.0;
  if(border==0) {
    /* or just zeroes */
    for(int i=0; i<n; i++) bkgtac[i]=0.0;
  } else {
    if(verbose>1) fprintf(stdout, "computing background from border pixels\n");
    if(imgBorderAverageTAC(img, NULL, NULL, bkgtac)!=0) {
      if(status!=NULL) strcpy(status, "cannot calculate background TAC");
      return(5);
    }
    /* Compute AUC for background TAC */
    for(int i=0; i<n; i++) bkg_tac_auc+=fdur[i]*bkgtac[i];
  }
  if(verbose>2) printf("bkg_tac_auc := %g\n", bkg_tac_auc);
 
  /*
   *  Allocate the mask image, if not allocated already
   */
  int maskExisted=0;
  if(timg->status==IMG_STATUS_OCCUPIED) {
    /* Pre-existing mask contents are used later */ 
    maskExisted=1;
  } else {
    ret=imgAllocateWithHeader(timg, img->dimz, img->dimy, img->dimx, 1, img);
    if(ret) {
      if(status!=NULL) strcpy(status, "cannot allocate memory for mask data");
      return(6);
    }
  }
 
  /*
   *  Go through all the pixels
   */
  if(verbose>1) fprintf(stdout, "processing pixels...\n");
  float wsum=0.0; for(int i=0; i<n; i++) wsum+=fdur[i]; if(verbose>4) printf("wsum=%g\n", wsum);
  int zi, yi, xi;
  for(zi=0; zi<img->dimz; zi++) {
    for(yi=0; yi<img->dimy; yi++) {
      for(xi=0; xi<img->dimx; xi++) {
        /* If pixel was already masked out, then keep it that way */
        if(maskExisted && timg->m[zi][yi][xi][0]<0.1) {
          timg->m[zi][yi][xi][0]=0.0; continue;
        }
        /* Compute AUC for this pixel */
        double pxl_tac_auc=0.0;
        for(int i=0; i<n; i++) pxl_tac_auc+=fdur[i]*img->m[zi][yi][xi][i];
        /* If tacAUC<AUCbkg and user did not want to keep negatives, then remove it */
        if(keepNegat==0 && pxl_tac_auc<=bkg_tac_auc) {
          timg->m[zi][yi][xi][0]=0.0; continue;
        }
 
        /* Calculate WSS between this TAC and zero */
        double wss0=0.0;
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i];
          wss0+=fdur[i]*g*g;
        }
 
        /* Calculate WSS between this TAC and background TAC */
        double wss1=0.0;
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i]-bkgtac[i];
          wss1+=fdur[i]*g*g;
        }
 
        /* Calculate WSS between this TAC and scaled background TAC */
        double wss2=0.0; double sf2=1.0;
        if(bkg_tac_auc>1.0E-10) sf2=pxl_tac_auc/bkg_tac_auc;
        if(sf2>2.0) sf2=2.0; // Limit for the scale
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i]-sf2*bkgtac[i];
          wss2+=fdur[i]*g*g;
        }
 
        /* Calculate WSS between this TAC and scaled mean TAC */
        double wss3=0.0; double sf3=pxl_tac_auc/mean_tac_auc;
        if(sf3<0.5) sf3=0.5; // Limit for the scale
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i]-sf3*meany[i];
          wss3+=fdur[i]*g*g;
        }
 
        /* Compute AICs */
        double aic0=aicSS(wss0, n, 0);
        double aic1=aicSS(wss1, n, 0);
        double aic2=aicSS(wss2, n, 1);
        double aic3=aicSS(wss3, n, 1);
        /* Conserve pixel value if aic3 is better than others */
        if(aic3<aic0 && aic3<aic1 && aic3<aic2) {
          timg->m[zi][yi][xi][0]=1.0;
        } else  {
          timg->m[zi][yi][xi][0]=0.0;
        }
        if(verbose>10) {
          printf(" p[%d][%d][%d]: wss0=%g wss1=%g wss2=%g wss3=%g;  aic0=%g aic1=%g aic2=%g aic3=%g\n",
                 zi, yi, xi, wss0, wss1, wss2, wss3, aic0, aic1, aic2, aic3);
          printf("   -> %g\n", timg->m[zi][yi][xi][0]);
        } else if(verbose>1 && yi==img->dimy-1 && xi==img->dimx-1 && timg->m[zi][yi][xi][0]==1.0) {
          printf(" p[%d][%d][%d]: wss0=%g wss1=%g wss2=%g wss3=%g;  aic0=%g aic1=%g aic2=%g aic3=%g\n",
                 zi, yi, xi, wss0, wss1, wss2, wss3, aic0, aic1, aic2, aic3);
          printf("   sf2=%g sf3=%g\n", sf2, sf3);
          printf("   -> %g\n", timg->m[zi][yi][xi][0]);
        }
      }
    }
  } // next zi
  if(verbose>1) fprintf(stdout, "... done.\n");
 
 
  if(status!=NULL) strcpy(status, "ok");
  return(0);
}