mapdiff.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "M3.h"
#include "M3_system.h"

int main( int argc, char **argv)
{
  FILE *f1, *f2, *f3;
  M3_MapEl *buffer1;
  M3_MapEl *buffer2;
  M3_Interval interval, bufferInterval;
  int64_t numReads = 32, n, loopCount;
  int32_t header1[2];
  int32_t header2[2];
  int64_t bufferLength, i, j;
  double norm = 0;
  double mean = 0;
  static int64_t diskBufferSize = -1;
  char *tempPtr;

  if( diskBufferSize == -1 )
  {
    tempPtr = getenv("M3_DISK_BUFFER_SIZE");
    if( tempPtr )
      diskBufferSize = strtoll( tempPtr, NULL, 10);
    else
      diskBufferSize = M3_DISK_BUFFER_SIZE;
  }

  if( argc < 4 )
  {
    fprintf(stderr, "Usage:  %s map1 map2 diffMap [relDiffMap]\n\n", argv[0] );
    return 1;
  }

  f1 = fopen( argv[1], "r");
  M3_ErrorCheck( -1, argv[1], (f1?1:0), M3_EFLAG_FOPEN_READ );
  f2 = fopen( argv[2], "r");
  M3_ErrorCheck( -1, argv[1], (f2?1:0), M3_EFLAG_FOPEN_READ );

  fread( header1, sizeof(int), 2, f1);
  fread( header2, sizeof(int), 2, f2);

  if( header1[0] != header2[0] || 
      header1[1] != header2[1] )
  {
    fprintf(stderr, "ERORR: File headers do not match.\n");
    return(1);
  }

  buffer1 = malloc(diskBufferSize);
  buffer2 = malloc(diskBufferSize);
  M3_ErrorCheck(-1, "main()", (buffer1 && buffer2), M3_EFLAG_MALLOC_FSTRING );


  f3 = fopen( argv[3], "w");
  fwrite(header1, sizeof(int), 2, f3);
  M3_ErrorCheck(-1, argv[3], (f3?1:0), M3_EFLAG_FOPEN_READ );

  if( argc > 4 )
    loopCount = 2;
  else 
    loopCount = 1;

  while( loopCount )
  {
    if( loopCount == 2 )
      f3 = fopen( argv[4], "w");
    else 
      f3 = fopen( argv[3], "w");
    M3_ErrorCheck(-1, argv[3], (f3?1:0), M3_EFLAG_FOPEN_READ );

    fwrite(header1, sizeof(int), 2, f3);

    interval.firstSample = 0;
    interval.lastSample = header1[1]-1;

    M3_bufferIntervalStart( interval, diskBufferSize, sizeof(M3_MapEl), &bufferLength, &numReads );
    for( i = 0; i < numReads; i++ )
    {
      M3_bufferIntervalStep( interval, bufferLength, i, &bufferInterval );
      n = bufferInterval.lastSample - bufferInterval.firstSample;

      fread(buffer1, sizeof(M3_MapEl), n, f1);
      fread(buffer2, sizeof(M3_MapEl), n, f2);

      for( j = 0; j < n; j++ )
      {
        if( buffer1[j].pixel != buffer2[j].pixel )
        {
          fprintf( stderr, "ERROR:  Observed pixels do not match.\n");
          return(1);
        }
        if( loopCount == 2 )
          buffer1[j].value = 2.0*(buffer1[j].value - buffer2[j].value)/(buffer1[j].value + buffer2[j].value);
        else
          buffer1[j].value -= buffer2[j].value;

        norm += buffer1[j].value*buffer1[j].value;
        mean += buffer1[j].value;
      }
    
      fwrite(buffer1, sizeof(M3_MapEl), n, f3);
    }
    norm = sqrt((norm - mean*mean/header1[1])/(header1[1]-1));
    mean = mean/header1[1];
    
    if( loopCount == 2 )
      fprintf(stdout, "Relative difference:\n");
    else 
      fprintf(stdout, "Absolute difference:\n");
    fprintf( stdout, "  rms = %e, mean = %e\n", norm, mean);

    fclose(f3);
    fseek(f1, SEEK_SET, 2*sizeof(int));
    fseek(f2, SEEK_SET, 2*sizeof(int));
    loopCount--;
  }

  fclose(f2);
  fclose(f1);

  free(buffer1);
  free(buffer2);

}

---