check_benchmark_pet_list-mode.cc

Go to the documentation of this file.

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <string>
#include <cmath>
using namespace std;

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

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

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

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

  // Dimensions
  int dim_x = 112;
  int dim_y = 112;
  int dim_z = 109;
  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(file_ref.c_str(),"rb");
  if (f_ref==NULL)
  {
    cerr << "***** Input reference image file 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 !" << endl;
    exit(1);
  }

  // Open and read test image
  FILE* f_new = fopen(file_new.c_str(),"rb");
  if (f_new==NULL)
  {
    cerr << "***** Input tested image file 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 !" << endl;
    exit(1);
  }

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

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

  // Compute the relative difference
  for (int v=0; v<dim_xyz; v++)
  {
    if (image_ref[v]==0.) image_diff[v] = 0.;
    else image_diff[v] = (image_new[v]-image_ref[v]) / 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 relative difference per slice
  float* mean_relative_difference = (float*)calloc(dim_z,sizeof(float));
  for (int z=0; z<dim_z; z++)
  {
    for (int xy=0; xy<dim_xy; xy++)
    {
      int v = z * dim_xy + xy;
      mean_relative_difference[z] += image_diff[v];
    }
    mean_relative_difference[z] /= ((float)dim_xy);
  }

  // Compute the global mean relative difference
  float global_mean_relative_difference = 0.;
  for (int z=0; z<dim_z; z++) global_mean_relative_difference += mean_relative_difference[z];
  global_mean_relative_difference /= ((float)dim_z);

  // Test the value with a 5 per 10 000 000 threshold
  float global_mean_relative_difference_threshold = 5.e-7;
  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 (100 times less critical)
  float mean_relative_difference_threshold = 5.e-5;
  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_relative_difference[z]);
    if (value > max_mean_relative_difference_over_slices) max_mean_relative_difference_over_slices = value;
    if (value >= mean_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 -" << mean_relative_difference_threshold << " and " << mean_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 = fabs(image_diff[v]);
        if (abs_relative_difference > max_abs_relative_difference) max_abs_relative_difference = abs_relative_difference;
      }
    }
  }

  // The the value with a 1 per 10 000 threshold
  float max_abs_relative_difference_threshold = 1.e-4;
  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;
}