imgadd.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[] = {
  "For simulation of simple PET images for software testing;",
  "adds the contents of existing image files as new plane(s) and writes",
  "the combined image as a new image file.",
  " ",
  "Usage: @P [Options] newfile file1 [file2 [file3 ...]]",
  " ",
  "Options:",
  " -x=<N>",
  "     Add the file contents N times; default is 1.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgdim, imgdelpl, ecatcat, dft2img, imgmove, flat2img, simcirc",
  " ",
  "Keywords: image, simulation, software testing",
  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 main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  int    ci, fi, ri, ret, pi, mi, fileNr=0, firstfile=0;
  char   outfile[FILENAME_MAX], file1[FILENAME_MAX], file2[FILENAME_MAX];
  int    multip=1; // How many times file contents are added.
  IMG    img1, img2, img3;
  int    dimz;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  file1[0]=file2[0]=outfile[0]=(char)0;
  imgInit(&img1); imgInit(&img2); imgInit(&img3);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "X=", 2)==0) {
      ret=atoi_with_check(cptr+2, &multip);
      if(ret==0 && multip>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 */
  /* Output filename */
  if(ai<argc) {strlcpy(outfile, argv[ai], FILENAME_MAX); ai++;}
  /* Input filenames */
  for(; ai<argc; ai++) {
    if(firstfile<1) firstfile=ai;
    //if(access(argv[ai], 0) == -1) {
    //  fprintf(stderr, "Error: file '%s' does not exist.\n", argv[ai]);
    //  return(1);
    //}
    /* Check that filename is not the same as the output filename */
    if(strcasecmp(argv[ai], outfile)==0) {
      fprintf(stderr, "Error: the same name for input and output image.\n");
      return(1);
    }
    fileNr++;
  }
 
  /* Is something missing? */
  if(fileNr<1) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("fileNr := %d\n", fileNr);
    printf("outfile := %s\n", outfile);
    printf("multip := %d\n", multip);
  }
 
 
  /*
   *  Read the image dimensions from each image
   *  and verify that images can be combined
   */
  if(verbose>1) printf("reading image headers\n");
  strlcpy(file1, argv[firstfile], FILENAME_MAX);
  if(verbose>2) printf("reading header in %s\n", file1);
  ret=imgReadHeader(file1, &img1, img1._fileFormat);
  if(ret!=STATUS_OK) {
    fprintf(stderr, "Error: %s\n", img1.statmsg); 
    if(verbose>4) imgInfo(&img1);
    imgEmpty(&img1);
    return(2);
  }
  dimz=multip*img1.dimz;
  for(ai=firstfile+1; ai<argc; ai++) {
    strlcpy(file2, argv[ai], FILENAME_MAX);
    imgEmpty(&img2);
    if(verbose>2) printf("reading header in %s\n", file2);
    ret=imgReadHeader(file2, &img2, img2._fileFormat);
    if(ret!=STATUS_OK) {
      fprintf(stderr, "Error: %s\n", img2.statmsg); 
      if(verbose>4) imgInfo(&img2);
      imgEmpty(&img1); imgEmpty(&img2);
      return(2);
    }
    /* Add the nr of planes */
    dimz+=multip*img2.dimz;
    /* Check that images can be combined */
    /* Frame nr */
    if(img1.dimt!=img2.dimt) {
      fprintf(stderr, "Error: different frame number.\n"); 
      imgEmpty(&img1); imgEmpty(&img2);
      return(3);
    }
    /* Data type */
    if(img1.type!=img2.type) {
      fprintf(stderr, "Error: different PET data type.\n"); 
      imgEmpty(&img1); imgEmpty(&img2);
      return(3);
    }
    /* Matrix size */
    if(img1.dimx!=img2.dimx || img1.dimy!=img2.dimy) {
      fprintf(stderr, "Error: different matrix size.\n"); 
      imgEmpty(&img1); imgEmpty(&img2);
      return(3);
    }
    /* Other not that important things */
    if(img1.unit!=img2.unit ||
       img1.isotopeHalflife!=img2.isotopeHalflife ||
       img1.zoom!=img2.zoom ||
       img1.axialFOV!=img2.axialFOV ||
       img1.transaxialFOV!=img2.transaxialFOV) {
      fprintf(stderr, "Warning: different header information.\n");
      if(verbose>1) {
        printf("unit: %d vs %d\n", img1.unit, img2.unit);
        printf("isotopeHalflife: %g vs %g\n",
               img1.isotopeHalflife, img2.isotopeHalflife);
        printf("zoom: %g vs %g\n", img1.zoom, img2.zoom);
        printf("axialFOV: %g vs %g\n", img1.axialFOV, img2.axialFOV);
        printf("transaxialFOV: %g vs %g\n",
               img1.transaxialFOV, img2.transaxialFOV);
      }
    }
  }
  imgEmpty(&img2);
  if(verbose>1) printf("combined_dimz := %d\n", dimz);
 
 
  /*
   *  Read the whole contents of file1
   */
  if(verbose>=0) printf("reading %s\n", file1);
  if(imgRead(file1, &img1)) {
    fprintf(stderr, "Error: %s\n", img1.statmsg); 
    if(verbose>5) imgInfo(&img1);
    imgEmpty(&img1);
    return(4);
  }
 
  /*
   *  Create place for new combined PET data
   */
  if(verbose>1) printf("creating combined data\n");
  ret=imgAllocateWithHeader(&img3, dimz, img1.dimy, img1.dimx, img1.dimt, &img1);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for combined data.\n"); 
    imgEmpty(&img1);
    return(5);
  }
  /* Set plane numbers */
  for(pi=0; pi<img1.dimz; pi++)
    if(img3.planeNumber[pi]<1) img3.planeNumber[pi]=pi+1;
  for(pi=img1.dimz; pi<img3.dimz; pi++)
    img3.planeNumber[pi]=img3.planeNumber[pi-1]+1;
 
  /* Copy pixel values from the first image */
  for(mi=0; mi<multip; mi++)
    for(pi=0; pi<img1.dimz; pi++)
      for(ri=0; ri<img1.dimy; ri++)
        for(ci=0; ci<img1.dimx; ci++)
          for(fi=0; fi<img1.dimt; fi++)
            img3.m[mi*img1.dimz+pi][ri][ci][fi]=img1.m[pi][ri][ci][fi];
  /* Save the current plane number */
  dimz=multip*img1.dimz;
  /* Original image 1 data is no more needed */
  imgEmpty(&img1);
  if(verbose>4) imgInfo(&img3);
 
  /*
   *  Combine the data from the other images
   */
  for(ai=firstfile+1; ai<argc; ai++) {
    strlcpy(file2, argv[ai], FILENAME_MAX);
    if(verbose>=0) printf("reading %s\n", file2);
    ret=imgRead(file2, &img2);
    if(ret!=STATUS_OK) {
      fprintf(stderr, "Error: %s\n", img2.statmsg); 
      if(verbose>5) imgInfo(&img2);
      imgEmpty(&img2); imgEmpty(&img3);
      return(4);
    }
    if(verbose>2) printf("  read ok\n");
    /* Copy pixel values */
    for(mi=0; mi<multip; mi++)
      for(pi=0; pi<img2.dimz; pi++)
        for(ri=0; ri<img2.dimy; ri++)
          for(ci=0; ci<img2.dimx; ci++)
            for(fi=0; fi<img2.dimt; fi++)
              img3.m[dimz+mi*img2.dimz+pi][ri][ci][fi]=img2.m[pi][ri][ci][fi];
    /* Save the current plane number */
    dimz+=multip*img2.dimz;
    imgEmpty(&img2);
  }
 
 
  /*
   *  Write combined PET data
   */
  if(verbose>=0) printf("writing combined data in %s\n", outfile);
  if(imgWrite(outfile, &img3)) {
    fprintf(stderr, "Error: %s\n", img3.statmsg); 
    imgEmpty(&img3);
    return(11);
  }
  imgEmpty(&img3);
  if(verbose>0) printf("done.\n");
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/