imgmatch.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <float.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
//#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Verify that data in two PET images are similar. ECAT 6.3 and 7, NIfTI-1,",
  "Analyze 7.5 and microPET formats are (partially) supported.",
  "Returns 0 if files match and a nonzero value otherwise.",
  " ",
  "Usage: @P [Options] image1 image2",
  " ",
  "Options:",
//" -a[ll]",
//"     All values in the header are required to match.",
//"     (As default only image sizes and data types are checked.)",
  " -header=<y|N>",
  "     All values in the headers (excluding transformation parameters)",
  "     are required to match (y) or not tested (n, default).",
  " -transform=<y|N>",
  "     Transformation parameters in the headers are required to match (y),",
  "     or not tested (n, default).",
  " -planes=<Y|n>",
  "     Planes numbers are required to match (y, default) or not tested (n).",
  " -frames=<Y|n>",
  "     Frame times are required to match (y, default) or not tested (n).",
  " -abs=<value>",
  "     Absolute result differences must not exceed the specified limit value;",
  "     by default no difference is allowed.",
  " -rel=<value>",
  "     Relative differences must not exceed the specified percent limit.",
//" -roughly",
//"     Pixel values are required to match roughly (99.9%), not exactly.",
//" -around",
//"     Pixel values are required to be around the same (90%), not exactly.",
  " -ss",
  "     Calculates sum-of-squares between image pixel values.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Image1 and image2 must be in the same format, and they must have same",
  "matrix size and plane numbers.",
  " ",
  "See also: tacmatch, resmatch, imgcalc, imgunit, ecatmax, img2tif",
  " ",
  "Keywords: Image, tool, 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;
  double   ROUGH_TEST_LIMIT=0.001;
  double   AROUND_TEST_LIMIT=0.1;
  int      ret;
  int      isHeader=0;
  int      isTransform=0;
  int      isPlanes=1;
  int      isFrames=1;
  int      calcSS=0;
  char     file1[FILENAME_MAX], file2[FILENAME_MAX], *cptr;
  IMG      img1, img2;
  double   accur=0.0;
  double   test_abs=0.0;
  float    maxrel, maxabs;
  VOXEL_4D maxrelvoxel, maxabsvoxel;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  file1[0]=file2[0]=(char)0;
  imgInit(&img1); imgInit(&img2);
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1;
    if(strncasecmp(cptr, "ROUGHLY", 4)==0) { // deprecated
      accur=ROUGH_TEST_LIMIT; continue;
    } else if(strncasecmp(cptr, "AROUND", 3)==0) { // deprecated
      accur=AROUND_TEST_LIMIT; continue;
    } else if(strncasecmp(cptr, "ABS=", 4)==0) {
      cptr+=4; test_abs=atof_dpi(cptr); if(test_abs>=0.0) continue;
    } else if(strncasecmp(cptr, "REL=", 4)==0) {
      cptr+=4; accur=0.01*atof_dpi(cptr); if(accur>=0.0) continue;
    } else if(strncasecmp(cptr, "HEADER=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "YES", 1)==0) {isHeader=1; continue;}
      if(strncasecmp(cptr, "NO", 1)==0) {isHeader=0; continue;}
    } else if(strncasecmp(cptr, "TRANSFORM=", 10)==0) {
      cptr+=10;
      if(strncasecmp(cptr, "YES", 1)==0) {isTransform=1; continue;}
      if(strncasecmp(cptr, "NO", 1)==0) {isTransform=0; continue;}
    } else if(strncasecmp(cptr, "ALL", 1)==0) { // deprecated
      isHeader=1; continue;
    } else if(strncasecmp(cptr, "PLANES=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "YES", 1)==0) {isPlanes=1; continue;}
      if(strncasecmp(cptr, "NO", 1)==0) {isPlanes=0; continue;}
    } else if(strncasecmp(cptr, "FRAMES=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "YES", 1)==0) {isFrames=1; continue;}
      if(strncasecmp(cptr, "NO", 1)==0) {isFrames=0; continue;}
    } else if(strcasecmp(cptr, "SS")==0) {
      calcSS=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);}
 
  /* Arguments */
  if(ai<argc) strlcpy(file1, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(file2, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!file2[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]);
    printf("\n");
    printf("file1 := %s\n", file1);
    printf("file2 := %s\n", file2);
    printf("isHeader := %d\n", isHeader);
    printf("isTransform := %d\n", isTransform);
    printf("isPlanes := %d\n", isPlanes);
    printf("isFrames := %d\n", isFrames);
    printf("accur := %g\n", accur);
    printf("test_abs := %g\n", test_abs);
    printf("calcSS := %d\n", calcSS);
  }
 
#if(0) // Analyze filenames can be funny, do not check
  /* Check that input files do exist */
  if(access(file1, 0) == -1) {
    fprintf(stderr, "Error: file '%s' not found.\n", file1); return(2);}
  if(access(file2, 0) == -1) {
    fprintf(stderr, "Error: file '%s' not found.\n", file2); return(2);}
#endif
 
  /*
   *  Read image files
   */
  if(verbose>1) fprintf(stdout, "reading image %s\n", file1);
  ret=imgRead(file1, &img1); if(verbose>10) imgInfo(&img1);
  if(ret) {fprintf(stderr, "Error: %s\n", img1.statmsg); return(2);}
  if(verbose>1) fprintf(stdout, "reading image %s\n", file2);
  ret=imgRead(file2, &img2); if(verbose>10) imgInfo(&img2);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img2.statmsg);
    imgEmpty(&img1); return(2);
  }
 
  /*
   *  Check that file types, sizes etc are matching
   */
  if(verbose>1) fprintf(stdout, "testing image size etc\n");
  if(img1.type!=img2.type) {
    if(verbose>0) fprintf(stderr, "Mismatching image type.\n");
    imgEmpty(&img1); imgEmpty(&img2); return(3);
  }
  if(img1.dimt!=img2.dimt){
    if(verbose>0) fprintf(stderr, "Mismatching number of frames.\n");
    imgEmpty(&img1); imgEmpty(&img2); return(4);
  }
  if(img1.dimx!=img2.dimx){
    if(verbose>0) fprintf(stderr, "Mismatching number of pixels (dimension x).\n");
    imgEmpty(&img1); imgEmpty(&img2); return(5);
  }
  if(img1.dimy!=img2.dimy){
    if(verbose>0) fprintf(stderr, "Mismatching number of pixels (dimension y).\n");
    imgEmpty(&img1); imgEmpty(&img2); return(6);
  }
  if(img1.dimz!=img2.dimz){
    if(verbose>0) fprintf(stderr, "Mismatching number of planes.\n");
    imgEmpty(&img1); imgEmpty(&img2); return(7);
  }
 
  /*
   *  Match headers if required
   */
  if(isHeader){
    if(verbose>1) fprintf(stdout, "testing header contents\n");
    ret=imgMatchHeader(&img1, &img2);
    if(ret) {
      if(verbose>0) fprintf(stdout, "Image headers do not match.\n");
      imgEmpty(&img1); imgEmpty(&img2); return(11);
    }
  }
 
  /*
   *  Match transformation parameters in headers if required
   */
  if(isTransform){
    if(verbose>1) fprintf(stdout, "testing header transformation parameters\n");
    ret=imgMatchTransform(&img1, &img2);
    if(ret) {
      if(verbose>0) fprintf(stdout, "Image transformation parameters do not match.\n");
      imgEmpty(&img1); imgEmpty(&img2); return(11);
    }
  }
 
  /*
   *  Match plane numbers if required
   */
  if(isPlanes) {
    if(verbose>1) fprintf(stdout, "testing plane numbers\n");
    ret=imgMatchPlanes(&img1, &img2);
    if(ret) {
      if(verbose>0) fprintf(stdout, "Plane numbers do not match.\n");
      imgEmpty(&img1); imgEmpty(&img2); return(21);
    }
  }
 
  /*
   *  Match frame times if required
   */
  if(isFrames) {
    if(verbose>1) fprintf(stdout, "testing frame times\n");
    ret=imgMatchFrames(&img1, &img2);
    if(ret) {
      if(verbose>0) fprintf(stdout, "Frame times do not match.\n");
      imgEmpty(&img1); imgEmpty(&img2); return(31);
    }
  }
 
  /*
   *  Calculate sum-of-squares between image pixel values, if requested
   */
  if(calcSS) {
    double ss=0.0;
    if(imgSS(&img1, &img2, &ss)!=0 || !isfinite(ss)) {
      fprintf(stderr, "Error: cannot calculate sum-of-squares.\n");
      fflush(stderr);
    } else {
      printf("SS := %e\n", ss);
      fflush(stdout);
    }
  }
 
  /*
   *  Check pixel values
   */
  if(verbose>1) fprintf(stdout, "testing pixel values\n");
/*
  ret=imgMatchMatrix(&img1, &img2, accur);
  if(ret) {
    fprintf(stdout, "Images do not match.\n");
    imgEmpty(&img1); imgEmpty(&img2); return(101);
  }
*/
  /* Calculate maximal absolute and relational differences */
  ret=imgMaxDifference(&img1, &img2,
                       &maxabsvoxel, &maxabs, &maxrelvoxel, &maxrel);
  if(ret>0) {
    fprintf(stdout, "Error in testing pixel values (%d).\n", ret);
    imgEmpty(&img1); imgEmpty(&img2); return(101);
  } else if(ret<0) { // no differences
    imgEmpty(&img1); imgEmpty(&img2); 
    if(verbose>0) fprintf(stdout, "Images match.\n"); 
    return(0);
  }
  /* Differences were found; check if below allowed limits */
  if(verbose>1) {
    printf("max_absolute_difference := %g\n", maxabs);
    if(maxabs>0.0)
      printf("max_absolute_difference_at := %d,%d,%d,%d\n",
        maxabsvoxel.z, maxabsvoxel.y, maxabsvoxel.x, maxabsvoxel.t);
    printf("max_relational_difference := %g\n", maxrel);
    if(maxrel>0.0)
      printf("max_relational_difference_at := %d,%d,%d,%d\n",
        maxrelvoxel.z, maxrelvoxel.y, maxrelvoxel.x, maxrelvoxel.t);
  }
  if(maxabs>test_abs && maxrel>accur) {
    if(verbose>0) fprintf(stdout, "Images do not match.\n");
    if(verbose>1) {
      printf("  %g > %g\n", maxabs, test_abs);
      printf("  %g > %g\n", maxrel, accur);
    }
    imgEmpty(&img1); imgEmpty(&img2); return(102);
  }
  /* but difference is below limit */
  if(maxabs>0.0 && maxabs<=test_abs) {
    if(verbose>0) printf("Maximal absolute difference is below the limit.\n");
    if(verbose>1) printf("  %g <= %g\n", maxabs, test_abs);
  }
  if(maxrel>0.0 && maxrel<=accur) {
    if(verbose>0) printf("Maximal relational difference is below the limit.\n");
    if(verbose>1) printf("  %g <= %g\n", maxrel, accur);
  }
 
  imgEmpty(&img1); imgEmpty(&img2); 
  if(verbose>0) fprintf(stdout, "Images are similar within specified limits.\n"); 
  
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/