imgcalc.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 <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simple arithmetic calculations on PET images.",
  " ",
  "Usage: @P [Options] file1 operation constant|file2 outputfile",
  " ",
  ", where",
  "    Operation = +, -, x , :",
  "    Constant = Image 1 is operated with this value.",
  " ",
  "Options:",
  " -frames",
  "     All frames in image 1 are processed with the 1st time frame of image 2.",
  " -max=<value>",
  "     Upper limit for output pixels (may be necessary in division).",
//" -scale = image2 data is scaled to image1 maximum before operation",
  " -abs[olute]",
  "     Save absolute voxel values.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Image 1 and 2 must be in the same format, and they must have same",
  "matrix size and plane numbers.",
  "Data is not interpolated or corrected for physical decay, and thus it is",
  "on the responsibility of the user to check that PET time frames are similar",
  "and that physical decay and other corrections are appropriately considered.",
  " ",
  "Constant can be given directly as value on the command line, or as name of",
  "file containing only the value, or value with key 'constant'.",
  " ",
  "Input files must not be overwritten with the output file.",
  " ",
  "Examples:",
  "     @P a999.flow.img - a1000.flow.img subtr.img",
//"     @P -scale s1444dy1.img - s1445dy1.img s1444.subtr.img",
  "     @P o657dy1.v : 2.34 o657rda.v",
//"     @P -frames -max=0.25 l333dy1.scn x l333.nrm l333dy1.corr.scn",
  " ",
  "See also: imgdecay, imgunit, eframe, imginv, esplit, imginteg, imgthrs",
  " ",
  "Keywords: image, modelling, simulation, tool",
  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 ecat7CopyHeadersNoQuant(
  const char *file1,
  const char *file2,
  int verbose
) {
  int ii, mi, ret;
  char tmp[1024];
  ECAT_HEADERS ehdr;
 
  if(verbose>1) printf("ecatCopyHeadersNoQuant(%s, %s)\n", file1, file2);  
  if(file1==NULL || file2==NULL) return(1);
  
  
  /* Read headers */
  ehdrInitiate(&ehdr);
  ret=ecat7ReadHeaders(file1, &ehdr, verbose-1); if(ret!=0) return(10);
  /* Remove main header information that should NOT be copied */
  strcpy(tmp, "user_process_code");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
  strcpy(tmp, "ecat_calibration_factor");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
 
  /* Remove subheader information that should NOT be copied */
  for(mi=0; mi<ehdr.nr; mi++) {
    strcpy(tmp, "scale_factor");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "image_min");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "image_max");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "scan_min");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "scan_max");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
  }
  
  /* Write headers */
  ret=ecat7WriteHeaders(file2, &ehdr, verbose-1);
  ehdrEmpty(&ehdr);
  if(ret!=0) return(10+ret);
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  int     ret, fi;
  int     allframes_with_one=0;
  int     frame_nr=0;
  char    file1[FILENAME_MAX], file2[FILENAME_MAX], file3[FILENAME_MAX];
  char   *cptr;
  char    operation=0; // '+', '-', '*', '/'
  int     isvalue=0;
  IMG     img1, img2;
  float   opconst=nan("");
  float   ulimit=0.0;
  int     absolutes=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  file1[0]=file2[0]=file3[0]=(char)0;
  imgInit(&img1); imgInit(&img2);
  /* 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(strcasecmp(cptr, "FRAMES")==0) {
      allframes_with_one=1; continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      ulimit=atof_dpi(cptr+4); if(ulimit>0.0) continue;
    } else if(strncasecmp(cptr, "ABSOLUTES", 3)==0) {
      absolutes=1; 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(file1, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    if(strncasecmp(argv[ai], "DIVISION", 3)==0) {
      operation='/';
    } else if(strncasecmp(argv[ai], "MULTIPLICATION", 3)==0) {
      operation='*';
    } else if(strncasecmp(argv[ai], "MINUS", 3)==0) {
      operation='-';
    } else if(strncasecmp(argv[ai], "PLUS", 3)==0) {
      operation='+';
    } else if(strlen(argv[ai])==1) switch(*argv[ai]) {
       case '+' : operation='+'; break;
       case '-' : operation='-'; break;
       case '*' : // this might be hard to use, but accept it anyway
       case '.' :
       case 'X' :
       case 'x' : operation='*'; break;
       case '/' : // MinGW does not work with '/', it requires '//'
       case ':' : operation='/'; break;
    }
    if(operation==0) {fprintf(stderr, "Error: invalid operator.\n"); return(1);}
    ai++;
  }
  if(ai<argc) {
    double v;
    char *s=iftReadValue(argv[ai], "constant", 0);
    if(s!=NULL) {ret=atof_with_check(s, &v); free(s);}
    else ret=atof_with_check(argv[ai], &v);
    if(ret==0) {isvalue=1; opconst=v; allframes_with_one=1;}
    else strlcpy(file2, argv[ai], FILENAME_MAX);
    ai++;
  }
  if(ai<argc) {strlcpy(file3, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing? */
  if(!file3[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    fflush(stdout); 
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]);
    printf("\n");
    printf("file1 := %s\n", file1);
    printf("operation := %c\n", operation);
    if(isvalue) printf("value := %g\n", opconst);
    else printf("file2 := %s\n", file2);
    printf("ulimit := %g\n", ulimit);
    printf("absolute_values := %d\n", absolutes);
    printf("allframes_with_one := %d\n", allframes_with_one);
  }
  if(verbose>1) IMG_TEST=verbose-1;
 
  /* Check that output file is not the same as either of input files */
  if(strcasecmp(file3, file1)==0 || strcasecmp(file3, file2)==0) {
    fprintf(stderr, "Error: original file must not be overwritten.\n");
    return(1);
  }
 
 
  /*
   *  Read information about input files
   */
  if(verbose>0) printf("reading image(s)...\n");
  ret=imgReadHeader(file1, &img1, IMG_UNKNOWN); if(verbose>10) imgInfo(&img1);
  if(ret) {
    fprintf(stderr, "Error in reading %s: %s\n", file1, imgStatus(ret));
    return(2);
  }
  if(*file2) {
    ret=imgReadHeader(file2, &img2, IMG_UNKNOWN); if(verbose>10) imgInfo(&img2);
    if(ret) {
      fprintf(stderr, "Error in reading %s: %s\n", file2, imgStatus(ret));
      imgEmpty(&img1); return(2);
    }
    /* Check that file types, sizes etc are matching */
    if(img1.type!=img2.type || img1.dimx!=img2.dimx || img1.dimy!=img2.dimy) {
      fprintf(stderr, "Error: non-matching files.\n");
      imgEmpty(&img1); imgEmpty(&img2); return(2);
    }
    if(img1.zoom!=img2.zoom) fprintf(stderr, "Warning: different zoom.\n");
    /* If file1 has several frames and file2 has only one frame,     */
    /* set allframes_with_one=1                                      */
    if(allframes_with_one==0 && img1.dimt>1 && img2.dimt==1) {
      if(verbose>0) {
        fprintf(stderr, "Warning: automatically applying option -frames\n"); fflush(stderr);}
      allframes_with_one=1;
    }
    /* Error, if frame nr is different but allframes_with_one==0 */
    if(allframes_with_one==0 && img1.dimt!=img2.dimt) {
      fprintf(stderr, "Error: non-matching frame number.\n");
      imgEmpty(&img1); imgEmpty(&img2); return(2);
    }
  }
 
 
  /*
   *  Operation, one frame at a time
   */
  if(verbose>0) printf("computing");
  frame_nr=img1.dimt;
  imgEmpty(&img1); imgEmpty(&img2);
  for(fi=0; fi<frame_nr; fi++) {
    if(verbose>1) printf("\n frame %d\n", fi+1);
    /* Read first/next frame in file1 */
    ret=imgReadFrame(file1, fi+1, &img1, 0);
    if(ret) {
      fprintf(stderr, "\nError in reading %s: %s\n", file1, imgStatus(ret));
      imgEmpty(&img1); imgEmpty(&img2); (void)remove(file3);
      return(3);
    }
    if(*file2 && (fi==0 || allframes_with_one==0)) {
      if(verbose>2) printf("\n reading frame %d in file2\n", fi+1);
      /* Read first/next frame in file2 */
      ret=imgReadFrame(file2, fi+1, &img2, 0);
      if(ret) {
        fprintf(stderr, "\nError in reading %s: %s\n", file2, imgStatus(ret));
        imgEmpty(&img1); imgEmpty(&img2); (void)remove(file3);
        return(4);
      }
    }
    /* Do the operation */
    if(verbose>2) printf("\noperating frame %d\n", fi+1);
    if(*file2) {
      if(allframes_with_one==0)
        ret=imgArithm(&img1, &img2, operation, ulimit, verbose-3);
      else
        ret=imgArithmFrame(&img1, &img2, operation, ulimit, verbose-3);
    } else {
      ret=imgArithmConst(&img1, opconst, operation, ulimit, verbose-3);
    }
    if(ret) {
      fprintf(stderr, "Error %d in calculation.\n", ret);
      imgEmpty(&img1); imgEmpty(&img2); (void)remove(file3); return(6);
    }
    /* Convert result pixel values into absolute values, if requested */
    if(absolutes && imgAbs(&img1)!=0) {
      fprintf(stderr, "Error in calculation.\n");
      imgEmpty(&img1); imgEmpty(&img2); (void)remove(file3); return(7);
    }
    /* Remove previous output image */
    if(fi==0) (void)remove(file3);
    /* Attenuation format cannot currently be saved from IMG */
    if(img1.type==IMG_TYPE_ATTN) img1.type=IMG_TYPE_RAW;
    /* Save the frame */
    ret=imgWriteFrame(file3, fi+1, &img1, 0);
    if(ret) {
      fprintf(stderr, "\nError in writing %s: %s\n", file3, imgStatus(ret));
      imgEmpty(&img1); imgEmpty(&img2); (void)remove(file3);
      return(11);
    }
    if(verbose==1) {fprintf(stdout, "."); fflush(stdout);}
  }
  if(verbose==1) printf("\n");
 
  if(img1._fileFormat==IMG_E7 || img1._fileFormat==IMG_E7_2D) {
    if(verbose>1) printf("copying header information\n");
    ret=ecat7CopyHeadersNoQuant(file1, file3, verbose-1);
    if(verbose>0 && ret!=0) {printf("copying headers not successful (%d).\n", ret); fflush(stdout);}
  }
 
  imgEmpty(&img1); imgEmpty(&img2);
  if(verbose>0) printf("done.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/