imgidpvc.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <math.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Iterative deconvolution method for partial volume correction of PET images",
  "in ECAT 6.3, 7.x, and Analyze 7.5 and NIfTI-1 formats.",
  "Not for production use!",
  " ",
  "Usage: @P [Options] image FWHMxy FWHMz pvcimage",
  " ",
  "Options:",
  " -method=<VC|RL>",
  "     Use van Cittert or Richardson-Lucy method; VC by default.",
  " -reblur=<yes|no>",
  "     Include reblurring step; by default yes.",
  " -alpha=<value>",
  "     Converging rate parameter (0.0<alpha<10.0).",
  " -beta=<value>",
  "     Upper limit for deconvolved image as a factor of image maximum value.",
  " -iter=<number>",
  "     Maximum number of iterations.",
  " -term=<value>",
  "     Limit for termination (0.0<term<1.0).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "FWHM can be given in mm (default) or in pixels. If image file does not",
  "contain correct pixel sizes, then FWHM must be entered in pixels.",
  " ",
  "Example: FWHM in mm",
  "     @P i5998dy1.v 2.5 2.65 i5998dy1_pvc.v",
  "Example: FWHM in pixels",
  "     @P i5998dy1.v 2.1pxl 1.7pxl i5998dy1_pvc.v",
  " ",
  "See also: imgfiltg, imgbox, imgmask",
  " ",
  "Keywords: image, PVC",
  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 imgPVCRRL(
  IMG *img1,
  double FWHMxy,
  double FWHMz,
  int reblur,
  double betaFactor,
  double alpha,
  int maxIterNr,
  double termLimit,
  IMG *img2,
  int verbose
) {
  if(verbose>0) {
    printf("%s(i, %g, %g, %d, %g, %g, %d, %g, ...)\n",
           __func__, FWHMxy, FWHMz, reblur, betaFactor, alpha, maxIterNr, termLimit);
    fflush(stdout);
  }
  if(img1==NULL || img2==NULL) return(1);
  if(img1->dimz<1 || img1->dimy<8 || img1->dimx<8 || img1->dimt<1) return(2);
  if(!(FWHMxy>0.0)) return(3);
  if(img1->dimz<2) FWHMz=0.0;
  if(!(betaFactor>0.0)) betaFactor=1.4;
  if(maxIterNr<1) maxIterNr=25;
  if(!(alpha>0.0)) alpha=1.6;
  if(!(termLimit>0.0) || termLimit>=1) termLimit=0.005;
 
  /* Convert FWHM into Gaussian SD */
  float SDxy=FWHMxy/2.355;
  float SDz=FWHMz/2.355;
 
  int x, dimx=img1->dimx;
  int y, dimy=img1->dimy;
  int z, dimz=img1->dimz;
  int t, dimt=img1->dimt;
 
  /* Allocate output image and set headers */
  {
    int ret=imgAllocateWithHeader(img2, dimz, dimy, dimx, dimt, img1); 
    if(ret) return(10+ret);
  }
 
  /* Allocate image frame, for 'true' and correction matrices */
  IMG tkm, cm; imgInit(&tkm); imgInit(&cm);
  {
    int ret=imgAllocateWithHeader(&tkm, dimz, dimy, dimx, 1, img1); 
    if(!ret) ret=imgAllocateWithHeader(&cm, dimz, dimy, dimx, 1, img1); 
    if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(5);}
  }
 
  /*
   *  One image frame at a time
   */
  for(t=0; t<dimt; t++) {
    if(verbose>2 && dimt>1) {printf("  frame %d\n", 1+t); fflush(stdout);}
 
    /* First, copy the original measured contents to "true" image matrix */
    for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
      tkm.m[z][y][x][0]=img1->m[z][y][x][t];
 
    /* Get image frame minimum and maximum for setting beta */
    float omin, omax;
    if(imgMinMax(&tkm, &omin, &omax)) {imgEmpty(&tkm); imgEmpty(&cm); return(6);}
    if(!(omax>omin) || !(omax>0.0)) continue;
    double beta=betaFactor*omax;
    if(verbose>4) {
      printf("  image_min := %g\n  image_max := %g\n", omin, omax);
      printf("  beta := %g\n", beta);
      fflush(stdout);
    }
 
    /* Iterations */
    int ret=0, iterNr=0;
    for(iterNr=0; iterNr<maxIterNr; iterNr++) {
      if(verbose>3) {printf("  iteration %d\n", 1+iterNr); fflush(stdout);}
      /* Convolution of current estimate of true image with PSF (Gaussian filter) */
      for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
        cm.m[z][y][x][0]=tkm.m[z][y][x][0];
      ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
      if(ret) {ret+=100; break;}
      /* Divide the measured image matrix with it */
      for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++) {
        float f=img1->m[z][y][x][t]/cm.m[z][y][x][0];
        if(!isfinite(f)) {
          if(fabs(img1->m[z][y][x][t])<1.0E-06) f=0.0; else f=1.0;
        }
        cm.m[z][y][x][0]=f;
      }
      /* Re-blurring the correction matrix with PSF (Gaussian filter) */
      if(reblur) {
        ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
        if(ret) {ret+=200; break;}
      }
      /* Multiply the previous estimate of true image matrix; */
      /* and calculate factor to check for termination */
      double dsum=0.0, psum=0.0;
      for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++) {
        double prev=tkm.m[z][y][x][0];
        tkm.m[z][y][x][0]*=cm.m[z][y][x][0];
        if(tkm.m[z][y][x][0]<=0.0 && img1->m[z][y][x][t]>0.0)
          tkm.m[z][y][x][0]=0.01*img1->m[z][y][x][t];
        if(tkm.m[z][y][x][0]>beta) tkm.m[z][y][x][0]=beta;
        double d=tkm.m[z][y][x][0]-prev;
        dsum+=d*d; psum+=prev*prev;
      }
      double dssqr=sqrt(dsum)/sqrt(psum); if(verbose>4) printf("  dssqr := %g\n", dssqr);
      if(dssqr<termLimit) break;
    } // next iteration
    if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(ret);}
    if(verbose>3) {printf("  iterNr := %d\n", (iterNr>=maxIterNr?maxIterNr:1+iterNr)); fflush(stdout);}
 
    /* Copy the estimate of true matrix to output image */
    for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
      img2->m[z][y][x][t]=tkm.m[z][y][x][0];
 
  } // next frame
 
  imgEmpty(&tkm); imgEmpty(&cm);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgPVCRVC(
  IMG *img1,
  double FWHMxy,
  double FWHMz,
  int reblur,
  double betaFactor,
  double alpha,
  int maxIterNr,
  double termLimit,
  IMG *img2,
  int verbose
) {
  if(verbose>0) {
    printf("%s(i, %g, %g, %d, %g, %g, %d, %g, ...)\n",
           __func__, FWHMxy, FWHMz, reblur, betaFactor, alpha, maxIterNr, termLimit);
    fflush(stdout);
  }
  if(img1==NULL || img2==NULL) return(1);
  if(img1->dimz<1 || img1->dimy<8 || img1->dimx<8 || img1->dimt<1) return(2);
  if(!(FWHMxy>0.0)) return(3);
  if(img1->dimz<2) FWHMz=0.0;
  if(!(betaFactor>0.0)) betaFactor=1.4;
  if(maxIterNr<1) maxIterNr=25;
  if(!(alpha>0.0)) alpha=1.6;
  if(!(termLimit>0.0) || termLimit>=1) termLimit=0.005;
 
  /* Convert FWHM into Gaussian SD */
  float SDxy=FWHMxy/2.355;
  float SDz=FWHMz/2.355;
 
  int x, dimx=img1->dimx;
  int y, dimy=img1->dimy;
  int z, dimz=img1->dimz;
  int t, dimt=img1->dimt;
 
  /* Allocate output image and set headers */
  {
    int ret=imgAllocateWithHeader(img2, dimz, dimy, dimx, dimt, img1); 
    if(ret) return(10+ret);
  }
 
  /* Allocate image frame, for 'true' and correction matrices */
  IMG tkm, cm; imgInit(&tkm); imgInit(&cm);
  {
    int ret=imgAllocateWithHeader(&tkm, dimz, dimy, dimx, 1, img1); 
    if(!ret) ret=imgAllocateWithHeader(&cm, dimz, dimy, dimx, 1, img1); 
    if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(5);}
  }
 
  /*
   *  One image frame at a time
   */
  for(t=0; t<dimt; t++) {
    if(verbose>2 && dimt>1) {printf("  frame %d\n", 1+t); fflush(stdout);}
 
    /* First, copy the original measured contents to "true" image matrix */
    for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
      tkm.m[z][y][x][0]=img1->m[z][y][x][t];
 
    /* Get image frame minimum and maximum for setting beta */
    float omin, omax;
    if(imgMinMax(&tkm, &omin, &omax)) {imgEmpty(&tkm); imgEmpty(&cm); return(6);}
    if(!(omax>omin) || !(omax>0.0)) continue;
    double beta=betaFactor*omax;
    if(verbose>4) {
      printf("  image_min := %g\n  image_max := %g\n", omin, omax);
      printf("  beta := %g\n", beta);
      fflush(stdout);
    }
 
    /* Iterations */
    int ret=0, iterNr=0;
    for(iterNr=0; iterNr<maxIterNr; iterNr++) {
      if(verbose>3) {printf("  iteration %d\n", 1+iterNr); fflush(stdout);}
      /* Convolution of current estimate of true image with PSF (Gaussian filter) */
      for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
        cm.m[z][y][x][0]=tkm.m[z][y][x][0];
      ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
      if(ret) {ret+=100; break;}
      /* Subtract from measured image matrix */
      for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
        cm.m[z][y][x][0]=img1->m[z][y][x][t]-cm.m[z][y][x][0];
      /* Re-blurring the difference with PSF (Gaussian filter) */
      if(reblur) {
        ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
        if(ret) {ret+=200; break;}
      }
      /* Add to the previous estimate of true image matrix, applying relaxation function; */
      /* and calculate factor to check for termination */
      double dsum=0.0, psum=0.0;
      for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++) {
        double r=alpha*(1.0 - (2.0/beta)*fabs(tkm.m[z][y][x][0] - beta/2.0));
        double d=r*cm.m[z][y][x][0]; dsum+=d*d; psum+=tkm.m[z][y][x][0]*tkm.m[z][y][x][0];
        tkm.m[z][y][x][0]+=r*cm.m[z][y][x][0];
      }
      double dssqr=sqrt(dsum)/sqrt(psum); if(verbose>4) printf("  dssqr := %g\n", dssqr);
      if(dssqr<termLimit) break;
    } // next iteration
    if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(ret);}
    if(verbose>3) {printf("  iterNr := %d\n", (iterNr>=maxIterNr?maxIterNr:1+iterNr)); fflush(stdout);}
 
    /* Copy the estimate of true matrix to output image */
    for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
      img2->m[z][y][x][t]=tkm.m[z][y][x][0];
 
  } // next frame
 
  imgEmpty(&tkm); imgEmpty(&cm);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     petfile[FILENAME_MAX], outfile[FILENAME_MAX];
  double   FWHMxy=nan(""), FWHMz=nan("");
  int      FWHMxyUnit=1, FWHMzUnit=1; // 0=pixels, 1=mm, 2=cm
  int      method=0; // 0=VC; 1=RL
  int      reblur=1;
  double   betaFactor=1.4;
  double   alpha=1.5;
  int      maxIterNr=40;
  double   termLimit=0.01;
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=outfile[0]=(char)0;
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "METHOD=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "VC", 1)==0) {method=0; continue;}
      if(strncasecmp(cptr, "RL", 1)==0) {method=1; continue;}
    } else if(strncasecmp(cptr, "REBLUR=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "YES", 1)==0) {reblur=1; continue;}
      if(strncasecmp(cptr, "NO", 1)==0) {reblur=0; continue;}
    } else if(strncasecmp(cptr, "BETA=", 5)==0 && strlen(cptr)>5) {
      betaFactor=atof_dpi(cptr+5); if(betaFactor>0.0) continue;
    } else if(strncasecmp(cptr, "ALPHA=", 6)==0) {
      alpha=atof_dpi(cptr+6); if(alpha>0.0) continue;
    } else if(strncasecmp(cptr, "TERM=", 5)==0) {
      termLimit=atof_dpi(cptr+5); if(termLimit>0.0 && termLimit<1) continue;
    } else if(strncasecmp(cptr, "ITER=", 5)==0) {
      maxIterNr=atoi(cptr+5); if(maxIterNr>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 */
  if(ai<argc) {strlcpy(petfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    int ret=atof_with_check(argv[ai], &FWHMxy);
    if(ret==0 && FWHMxy>=0.0) {
      /* Check if unit was specified too */
      if(strcasestr(argv[ai], "pxl")!=NULL) FWHMxyUnit=0;
      if(strcasestr(argv[ai], "mm")!=NULL) FWHMxyUnit=1;
      if(strcasestr(argv[ai], "cm")!=NULL) FWHMxyUnit=2;
      ai++;
    } else {
      fprintf(stderr, "Error: invalid x,y FWHM.\n");
      return(1);
    }
  }
  if(ai<argc) {
    int ret=atof_with_check(argv[ai], &FWHMz);
    if(ret==0 && FWHMz>=0.0) {
      /* Check if unit was specified too */
      if(strcasestr(argv[ai], "pxl")!=NULL) FWHMzUnit=0;
      if(strcasestr(argv[ai], "mm")!=NULL) FWHMzUnit=1;
      if(strcasestr(argv[ai], "cm")!=NULL) FWHMzUnit=2;
      ai++;
    } else {
      fprintf(stderr, "Error: invalid z FWHM.\n");
      return(1);
    }
  }
  if(ai<argc) {strlcpy(outfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Did we get all the information that we need? */
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!(FWHMxy>0.0) && !(FWHMz>0.0)) {
    fprintf(stderr, "Error: invalid FWHM.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("FWHMxy := %g\n", FWHMxy);
    printf("FWHMz := %g\n", FWHMz);
    printf("outfile := %s\n", outfile);
    if(method==1) printf("method := RL\n"); else printf("method := VC\n");
    printf("reblur := %d\n", reblur);
    printf("alpha := %g\n", alpha);
    printf("betaFactor := %g\n", betaFactor);
    printf("termLimit := %g\n", termLimit);
    printf("maxIterNr := %d\n", maxIterNr);
  }
  if(verbose>9) IMG_TEST=verbose-10; else IMG_TEST=0;
 
 
  /*
   *  Read image
   */
  if(verbose>0) fprintf(stdout, "reading image %s\n", petfile);
  IMG img; imgInit(&img);
  if(imgRead(petfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg); if(verbose>1) imgInfo(&img);
    return(2);
  }
  if(verbose>1) {
    printf("image dimensions[z y x]: %d %d %d\n", img.dimz, img.dimy, img.dimx);
    printf("image voxel size[z y x]: %g %g %g\n", img.sizez, img.sizey, img.sizex);
    printf("image frame nr: %d\n", img.dimt);
  }
  /* Check if PET data is raw or image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", petfile);
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
  /*
   *  Convert FWHM into pixel units, if necessary
   */
  if((FWHMxyUnit!=0 && (!(img.sizex>0.0) || !(img.sizey>0.0))) || (FWHMzUnit!=0 && !(img.sizez>0.0)))
  {
    fprintf(stderr, "Error: unknown image pixel size.\n");
    imgEmpty(&img); return(3);
  }
  /* Convert cm -> mm */
  if(FWHMxyUnit==2) {FWHMxy*=10.0; FWHMxyUnit=1;}
  if(FWHMzUnit==2) {FWHMz*=10.0; FWHMzUnit=1;}
  /* Convert mm -> pixels */
  if(FWHMxyUnit==1) {FWHMxy/=img.sizex; FWHMxyUnit=0; printf("FWHMxy[pxl] := %g\n", FWHMxy);}
  if(FWHMzUnit==1) {FWHMz/=img.sizez; FWHMzUnit=0; printf("FWHMz[pxl] := %g\n", FWHMz);}
 
 
  /*
   *  PVC
   */
  if(verbose>0) fprintf(stdout, "correcting...\n");
  IMG out; imgInit(&out);
  {
    int ret=0;
    if(method==1)
      ret=imgPVCRRL(&img, FWHMxy, FWHMz, reblur, betaFactor, alpha, maxIterNr, termLimit, 
                    &out, verbose-2);
    else
      ret=imgPVCRVC(&img, FWHMxy, FWHMz, reblur, betaFactor, alpha, maxIterNr, termLimit, 
                    &out, verbose-2);
    if(ret) {
      fprintf(stderr, "Error: cannot correct image.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&out); return(3);
    }
  }
 
  /* Free memory allocated for original image */
  imgEmpty(&img);
 
 
  /*
   *  Save the image
   */
  if(verbose>1) fprintf(stdout, "writing image in %s\n", outfile);
  {
    int ret=imgWrite(outfile, &out);
    if(ret) {
      fprintf(stderr, "Error: %s\n", out.statmsg);
      imgEmpty(&out); return(11);
    }
    if(verbose>0) fprintf(stdout, "Image was written in %s\n", outfile);
  }
 
  /* Free memory allocated for output image */
  imgEmpty(&out);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/