![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 #include <iostream> 00002 #include <cstdlib> 00003 #include <cstdio> 00004 #include <string> 00005 #include <cmath> 00006 using namespace std; 00007 00008 void showHelp(int code) 00009 { 00010 cout << endl; 00011 cout << "Usage: check_benchmark reference_image tested_image" << endl; 00012 cout << endl; 00013 exit(code); 00014 } 00015 00016 int main(int argc, char** argv) 00017 { 00018 // ========================================================================== 00019 // Parameters 00020 // ========================================================================== 00021 00022 if (argc!=3) showHelp(0); 00023 string file_ref = (string)argv[1]; 00024 string file_new = (string)argv[2]; 00025 00026 // ========================================================================== 00027 // Allocate, open and read images 00028 // ========================================================================== 00029 00030 // Dimensions 00031 int dim_x = 112; 00032 int dim_y = 112; 00033 int dim_z = 109; 00034 int dim_xy = dim_x * dim_y; 00035 int dim_xyz = dim_x * dim_y * dim_z; 00036 00037 // Allocations 00038 float* image_ref = (float*)malloc(dim_xyz*sizeof(float)); 00039 float* image_new = (float*)malloc(dim_xyz*sizeof(float)); 00040 00041 // Open and read reference image 00042 FILE* f_ref = fopen(file_ref.c_str(),"rb"); 00043 if (f_ref==NULL) 00044 { 00045 cerr << "***** Input reference image file is missing or corrupted !" << endl; 00046 exit(1); 00047 } 00048 int nb_data_ref = fread(image_ref,sizeof(float),dim_xyz,f_ref); 00049 fclose(f_ref); 00050 if (nb_data_ref!=dim_xyz) 00051 { 00052 cerr << "***** Failed to read all data from reference image file !" << endl; 00053 exit(1); 00054 } 00055 00056 // Open and read test image 00057 FILE* f_new = fopen(file_new.c_str(),"rb"); 00058 if (f_new==NULL) 00059 { 00060 cerr << "***** Input tested image file is missing or corrupted !" << endl; 00061 exit(1); 00062 } 00063 int nb_data_new = fread(image_new,sizeof(float),dim_xyz,f_new); 00064 fclose(f_new); 00065 if (nb_data_new!=dim_xyz) 00066 { 00067 cerr << "***** Failed to read all data from tested image file !" << endl; 00068 exit(1); 00069 } 00070 00071 // ========================================================================== 00072 // Compute the relative difference image 00073 // ========================================================================== 00074 00075 // Allocation 00076 float* image_diff = (float*)malloc(dim_xyz*sizeof(float)); 00077 00078 // Compute the relative difference 00079 for (int v=0; v<dim_xyz; v++) 00080 { 00081 if (image_ref[v]==0.) image_diff[v] = 0.; 00082 else image_diff[v] = (image_new[v]-image_ref[v]) / image_ref[v]; 00083 } 00084 00085 // Number of problems found 00086 int problems = 0; 00087 00088 // ========================================================================== 00089 // Test 1: check that the mean relative difference is below a threshold, for 00090 // the entire image, but also slice by slice. 00091 // ========================================================================== 00092 00093 // Compute the mean relative difference per slice 00094 float* mean_relative_difference = (float*)calloc(dim_z,sizeof(float)); 00095 for (int z=0; z<dim_z; z++) 00096 { 00097 for (int xy=0; xy<dim_xy; xy++) 00098 { 00099 int v = z * dim_xy + xy; 00100 mean_relative_difference[z] += image_diff[v]; 00101 } 00102 mean_relative_difference[z] /= ((float)dim_xy); 00103 } 00104 00105 // Compute the global mean relative difference 00106 float global_mean_relative_difference = 0.; 00107 for (int z=0; z<dim_z; z++) global_mean_relative_difference += mean_relative_difference[z]; 00108 global_mean_relative_difference /= ((float)dim_z); 00109 00110 // Test the value with a 5 per 10 000 000 threshold 00111 float global_mean_relative_difference_threshold = 5.e-7; 00112 if (fabs(global_mean_relative_difference) < global_mean_relative_difference_threshold) 00113 { 00114 cout << " --> Mean relative difference test over the whole image passed (value: " << global_mean_relative_difference << ")." << endl; 00115 } 00116 // Test not passed ! 00117 else 00118 { 00119 cout << " --> !!!!! The mean relative difference over the whole image is " << global_mean_relative_difference << " !" << endl; 00120 cout << " It is supposed to be between -" << global_mean_relative_difference_threshold << " and " 00121 << global_mean_relative_difference_threshold << " !" << endl; 00122 problems++; 00123 } 00124 00125 // Test the value of each slice (100 times less critical) 00126 float mean_relative_difference_threshold = 5.e-5; 00127 float max_mean_relative_difference_over_slices = 0.; 00128 int slice_not_passed = 0; 00129 for (int z=0; z<dim_z; z++) 00130 { 00131 float value = fabs(mean_relative_difference[z]); 00132 if (value > max_mean_relative_difference_over_slices) max_mean_relative_difference_over_slices = value; 00133 if (value >= mean_relative_difference_threshold) slice_not_passed++; 00134 } 00135 if (slice_not_passed==0) 00136 { 00137 cout << " --> Mean relative difference test for all slices passed (maximum absolute value: " << max_mean_relative_difference_over_slices << ")." << endl; 00138 } 00139 // Test not passed 00140 else 00141 { 00142 cout << " --> !!!!! The mean relative difference for each slice has a maximum absolute value of " << max_mean_relative_difference_over_slices << " !" << endl; 00143 cout << " It is supposed to be between -" << mean_relative_difference_threshold << " and " << mean_relative_difference_threshold << " !" << endl; 00144 cout << " The test failed for " << slice_not_passed << " slices !" << endl; 00145 problems++; 00146 } 00147 00148 // ========================================================================== 00149 // Test 2: check that the absolute relative difference is below a threshold 00150 // for any voxels in the 92 pixels diameter central cylinder. 00151 // (due to NEGML, large differences can occur around the border of 00152 // the object; also the radial extension of the data has been fairly 00153 // troncated not so far away from the phantom) 00154 // ========================================================================== 00155 00156 // Compute the maximum absolute relative difference 00157 float max_abs_relative_difference = 0.; 00158 for (int z=0; z<dim_z; z++) 00159 { 00160 for (int y=0; y<dim_y; y++) 00161 { 00162 for (int x=0; x<dim_x; x++) 00163 { 00164 // Look if we are inside a cylinder of 45-voxels radius 00165 float y_pos = (((float)(y-dim_y/2))+0.5); 00166 float x_pos = (((float)(x-dim_x/2))+0.5); 00167 if (y_pos*y_pos+x_pos*x_pos>45*45) continue; 00168 // Compute 1D cumulative voxel index 00169 int v = z * dim_y * dim_x + y * dim_x + x; 00170 float abs_relative_difference = fabs(image_diff[v]); 00171 if (abs_relative_difference > max_abs_relative_difference) max_abs_relative_difference = abs_relative_difference; 00172 } 00173 } 00174 } 00175 00176 // The the value with a 1 per 10 000 threshold 00177 float max_abs_relative_difference_threshold = 1.e-4; 00178 if (max_abs_relative_difference < max_abs_relative_difference_threshold) 00179 { 00180 cout << " --> Maximum absolute relative difference test passed (value: " << max_abs_relative_difference << ")." << endl; 00181 } 00182 // Test not passed ! 00183 else 00184 { 00185 cout << " --> !!!!! The maximum absolute relative difference is " << max_abs_relative_difference << " !" << endl; 00186 cout << " It is supposed to be below " << max_abs_relative_difference_threshold << " !" << endl; 00187 problems++; 00188 } 00189 00190 // ========================================================================== 00191 // Summary 00192 // ========================================================================== 00193 00194 // Benchmark passed 00195 if (problems==0) 00196 { 00197 cout << " --> Benchmark successfully passed." << endl; 00198 } 00199 // Benchmark not passed ! 00200 else 00201 { 00202 cout << endl; 00203 cout << " --> !!!!! Benchmark tests failed (" << problems << " problems found) !" << endl; 00204 cout << " Double check all what you have done first. Then email the mailing list." << endl; 00205 } 00206 00207 // End 00208 return 0; 00209 } 00210 00211