check_benchmark_pet_dynamic.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cmath>
#ifdef _WIN32
#include <windows.h>
#endif

using namespace std;

void showHelp(int code)
{
  cout << endl;
  cout << "Usage:  check_benchmark  tested_image.hdr" << endl;
  cout << endl;
  exit(code);
}

int main(int argc, char** argv)
{
  // ==========================================================================
  // Parameters
  // ==========================================================================

  if (argc!=3) showHelp(0);
  string hdr_new = (string)argv[1];
  string nb_dyn_imgs_str = (string)argv[2];

  // ==========================================================================
  // Open header and read useful information
  // ==========================================================================

  // Open file
  ifstream fhdr(hdr_new.c_str());
  if (!fhdr)
  {
    cerr << "***** Input header file '" << hdr_new << "' is missing or corrupted !" << endl;
    exit(1);
  }

  // Info to find
  string img_new = "";
  string castor_version = "";

  // Read line after line
  char line[1024];
  fhdr.getline(line,1024);
  while (!fhdr.eof())
  {
    size_t found;
    // Transform it to string to benefit from useful c++ functions
    string test = (string)line;
    // Data file name
    found = test.find("name of data file");
    if (found!=string::npos)
    {
      // Get the file name
      found = test.find("=");
      img_new = test.substr(found+1);
      // Remove blancks and cariage return
      while ( (found=img_new.find(" ")) != string::npos) img_new.replace(found,1,"");
      while ( (found=img_new.find("\r")) != string::npos) img_new.replace(found,1,"");
    }
    // CASToR version
    found = test.find("CASToR version");
    if (found!=string::npos)
    {
      // Get the CASToR version
      found = test.find("=");
      castor_version = test.substr(found+1);
      // Remove blanks and cariage return
      while ( (found=castor_version.find(" ")) != string::npos)  castor_version.replace(found,1,"");
      while ( (found=castor_version.find("\r")) != string::npos) castor_version.replace(found,1,"");
    }
    // Read a new line
    fhdr.getline(line,1024);
  }

  // Close file
  fhdr.close();

  // Check data file name
  if (img_new=="")
  {
    cerr << "***** Did not find the data file name in the provided header file !" << endl;
    exit(1);
  }
  // Check CASToR version
  if (castor_version=="")
  {
    cout << "  --> Did not find the CASToR version in the provided header file. Assume it to be 1.0." << endl;
    castor_version = "1.0";
  }
  else
  {
    cout << "  --> Found CASToR version: " << castor_version << endl;
  }
  // ---------------
  // Truncate any third decimal on the CASToR version, e.g. 2.0.1 becomes 2.0, because we do not expect to be
  // changes in the results between such versions.

  // Truncated CASToR version
  string castor_version_trunc = castor_version;
  // First the position of the first '.'
  size_t first_point = castor_version_trunc.find_first_of(".");
  // Then, while we find a last '.' which is not the same as the first, we truncate the version name
  size_t last_point = string::npos;
  while ( (last_point=castor_version_trunc.find_last_of(".")) != first_point )
  {
    castor_version_trunc = castor_version_trunc.substr(0,last_point);
  }
  // Test if the version is different
  if (castor_version_trunc!=castor_version)
  {
    cout << "  --> Assume same results as CASToR version " << castor_version_trunc << endl;
    castor_version = castor_version_trunc;
  }

  // ---------------
  // Build reference image file name (we replace the substring "challenger" by "reference")
  
  // First find the position of the key word "challenger"
  size_t pos_challenger = img_new.find("challenger");
  if (pos_challenger==string::npos)
  {
    cerr << "***** The name of the image file being tested is supposed to contain the key work \"challenger\" in its name !" << endl;
    exit(1);
  }
  // Copy the file name
  string img_ref = img_new;
  // Set the sub-string that will replace the word "challenger"
  string reference = "reference_v" + castor_version;
  // Then replace "challenger" by the correct reference image and we're done
  img_ref.replace(pos_challenger,10,reference);

  // ==========================================================================
  // Allocate, open and read images
  // ==========================================================================

  // Dimensions
  int nb_dyn_imgs = atoi(nb_dyn_imgs_str.c_str());
  
  int dim_x = 256;
  int dim_y = 256;
  int dim_z = nb_dyn_imgs;
  int dim_xy = dim_x * dim_y;
  int dim_xyz = dim_x * dim_y * dim_z;

  // Allocations
  float* image_ref = (float*)malloc(dim_xyz*sizeof(float));
  float* image_new = (float*)malloc(dim_xyz*sizeof(float));

  // Open and read reference image
  FILE* f_ref = fopen(img_ref.c_str(),"rb");
  if (f_ref==NULL)
  {
    cerr << "***** Input reference image file '" << img_ref << "' is missing or corrupted !" << endl;
    exit(1);
  }
  int nb_data_ref = fread(image_ref,sizeof(float),dim_xyz,f_ref);
  fclose(f_ref);
  if (nb_data_ref!=dim_xyz)
  {
    cerr << "***** Failed to read all data from reference image file '" << img_ref << "' !" << endl;
    exit(1);
  }

  // Open and read test image
  FILE* f_new = fopen(img_new.c_str(),"rb");
  if (f_new==NULL)
  {
    cerr << "***** Input tested image file '" << img_new << "' is missing or corrupted !" << endl;
    exit(1);
  }
  int nb_data_new = fread(image_new,sizeof(float),dim_xyz,f_new);
  fclose(f_new);
  if (nb_data_new!=dim_xyz)
  {
    cerr << "***** Failed to read all data from tested image file '" << img_new << "' !" << endl;
    exit(1);
  }

  // ==========================================================================
  // Compute the difference image
  // ==========================================================================

  // Allocation
  double* image_diff = (double*)malloc(dim_xyz*sizeof(double));

  // Compute the difference
  for (int v=0; v<dim_xyz; v++)
  {
    image_diff[v] = ((double)image_new[v])-((double)image_ref[v]);
  }

  // Number of problems found
  int problems = 0;

  // ==========================================================================
  // Test 1: check that the mean relative difference is below a threshold, for
  //         the entire image, but also slice by slice.
  // ==========================================================================

  // Compute the mean difference per slice, as well as the mean of the reference
  double* mean_difference_per_slice = (double*)calloc(dim_z,sizeof(double));
  double* mean_reference_per_slice = (double*)calloc(dim_z,sizeof(double));
  for (int z=0; z<dim_z; z++)
  {
    for (int xy=0; xy<dim_xy; xy++)
    {
      int v = z * dim_xy + xy;
      mean_difference_per_slice[z] += image_diff[v];
      mean_reference_per_slice[z] += ((double)image_ref[v]);
    }
    mean_difference_per_slice[z] /= ((double)dim_xy);
    mean_reference_per_slice[z] /= ((double)dim_xy);
  }

  // Compute the global mean difference and reference
  double global_mean_difference = 0.;
  double global_mean_reference = 0.;
  for (int z=0; z<dim_z; z++)
  {
    global_mean_difference += mean_difference_per_slice[z];
    global_mean_reference += mean_reference_per_slice[z];
  }
  global_mean_difference /= ((double)dim_z);
  global_mean_reference /= ((double)dim_z);

  // Compute the global mean relative difference
  float global_mean_relative_difference = global_mean_difference / global_mean_reference;

  // Test the value with a 5 per 100 000 threshold
  float global_mean_relative_difference_threshold = 5.e-5;
  if (fabs(global_mean_relative_difference) < global_mean_relative_difference_threshold)
  {
    cout << "  --> Mean relative difference test over the whole image passed (value: " << global_mean_relative_difference << ")." << endl;
  }
  // Test not passed !
  else
  {
    cout << "  --> !!!!! The mean relative difference over the whole image is " << global_mean_relative_difference << " !" << endl;
    cout << "            It is supposed to be between -" << global_mean_relative_difference_threshold << " and "
                                                         << global_mean_relative_difference_threshold << " !" << endl;
    problems++;
  }

  // Test the value of each slice with a 5 per 10 000 threshold
  float slice_relative_difference_threshold = 5.e-4;
  float max_mean_relative_difference_over_slices = 0.;
  int slice_not_passed = 0;
  for (int z=0; z<dim_z; z++)
  {
    float value = fabs(mean_difference_per_slice[z] / mean_reference_per_slice[z]);
    if (value > max_mean_relative_difference_over_slices) max_mean_relative_difference_over_slices = value;
    if (value >= slice_relative_difference_threshold) slice_not_passed++;
  }
  if (slice_not_passed==0)
  {
    cout << "  --> Mean relative difference test for all slices passed (maximum absolute value: " << max_mean_relative_difference_over_slices << ")." << endl;
  }
  // Test not passed
  else
  {
    cout << "  --> !!!!! The mean relative difference for each slice has a maximum absolute value of " << max_mean_relative_difference_over_slices << " !" << endl;
    cout << "            It is supposed to be between -" << slice_relative_difference_threshold << " and " << slice_relative_difference_threshold << " !" << endl;
    cout << "            The test failed for " << slice_not_passed << " slices !" << endl;
    problems++;
  }

  // ==========================================================================
  // Test 2: check that the absolute relative difference is below a threshold
  //         for any voxels away from the border.
  // ==========================================================================

  // Allow some margin around the borders
  int voxel_margin = 5;

  // Compute the maximum absolute relative difference
  float max_abs_relative_difference = 0.;
  for (int z=0; z<dim_z; z++)
  {
    // Skip to next voxel if too close from the border
    if (z<voxel_margin || z>dim_z-voxel_margin-1) continue;
    for (int y=0; y<dim_y; y++)
    {
      // Skip to next voxel if too close from the border
      if (y<voxel_margin || y>dim_y-voxel_margin-1) continue;
      for (int x=0; x<dim_x; x++)
      {
        // Compute 1D cumulative voxel index
        int v = z * dim_y * dim_x + y * dim_x + x;
        // Skip to next voxel if too close from the border
        if (x<voxel_margin || x>dim_x-voxel_margin-1) continue;
        float abs_relative_difference = 0.;
        if (image_ref[v]!=0.) abs_relative_difference = fabs(image_diff[v]/image_ref[v]);
        if (abs_relative_difference > max_abs_relative_difference) max_abs_relative_difference = abs_relative_difference;
      }
    }
  }

  // The the value with a 1 per 1 000 threshold
  float max_abs_relative_difference_threshold = 1.e-3;
  if (max_abs_relative_difference < max_abs_relative_difference_threshold)
  {
    cout << "  --> Maximum absolute relative difference test passed (value: " << max_abs_relative_difference << ")." << endl;
  }
  // Test not passed !
  else
  {
    cout << "  --> !!!!! The maximum absolute relative difference is " << max_abs_relative_difference << " !" << endl;
    cout << "            It is supposed to be below " << max_abs_relative_difference_threshold << " !" << endl;
    problems++;
  }

  // ==========================================================================
  // Summary
  // ==========================================================================

  // Benchmark passed
  if (problems==0)
  {
    cout << "  --> Benchmark successfully passed." << endl;
  }
  // Benchmark not passed !
  else
  {
    cout << endl;
    cout << "  --> !!!!! Benchmark tests failed (" << problems << " problems found) !" << endl;
    cout << "            Double check all what you have done first. Then email the mailing list." << endl;
  }

  // End
  return 0;
}