RTK  2.1.0
Reconstruction Toolkit
rtkTest.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright RTK Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef rtkTest_h
20 #define rtkTest_h
21 
23 #include <itkImageFileWriter.h>
25 #include "rtkTestConfiguration.h"
26 
27 
28 template <class TImage>
29 #if FAST_TESTS_NO_CHECKS
30 void
31 CheckImageQuality(typename TImage::Pointer itkNotUsed(recon),
32  typename TImage::Pointer itkNotUsed(ref),
33  double itkNotUsed(ErrorPerPixelTolerance),
34  double itkNotUsed(PSNRTolerance),
35  double itkNotUsed(RefValueForPSNR))
36 {}
37 #else
38 void
39 CheckImageQuality(typename TImage::Pointer recon,
40  typename TImage::Pointer ref,
41  double ErrorPerPixelTolerance,
42  double PSNRTolerance,
43  double RefValueForPSNR)
44 {
45  using ImageIteratorType = itk::ImageRegionConstIterator<TImage>;
46  ImageIteratorType itTest(recon, recon->GetBufferedRegion());
47  ImageIteratorType itRef(ref, ref->GetBufferedRegion());
48 
49  using ErrorType = double;
50  ErrorType TestError = 0.;
51  ErrorType EnerError = 0.;
52 
53  itTest.GoToBegin();
54  itRef.GoToBegin();
55 
56  while (!itRef.IsAtEnd())
57  {
58  typename TImage::PixelType TestVal = itTest.Get();
59  typename TImage::PixelType RefVal = itRef.Get();
60  TestError += itk::Math::abs(RefVal - TestVal);
61  EnerError += std::pow(ErrorType(RefVal - TestVal), 2.);
62  ++itTest;
63  ++itRef;
64  }
65  // Error per Pixel
66  ErrorType ErrorPerPixel = TestError / ref->GetBufferedRegion().GetNumberOfPixels();
67  std::cout << "\nError per Pixel = " << ErrorPerPixel << std::endl;
68  // MSE
69  ErrorType MSE = EnerError / ref->GetBufferedRegion().GetNumberOfPixels();
70  std::cout << "MSE = " << MSE << std::endl;
71  // PSNR
72  ErrorType PSNR = 20 * log10(RefValueForPSNR) - 10 * log10(MSE);
73  std::cout << "PSNR = " << PSNR << "dB" << std::endl;
74  // QI
75  ErrorType QI = (RefValueForPSNR - ErrorPerPixel) / RefValueForPSNR;
76  std::cout << "QI = " << QI << std::endl;
77 
78  // // It is often necessary to write the images and look at them
79  // // to understand why a given test fails. This portion of code
80  // // does that. It should be left here but commented out, since
81  // // it is only useful in specific debugging tasks
82  // using FileWriterType = itk::ImageFileWriter<TImage>;
83  // typename FileWriterType::Pointer writer = FileWriterType::New();
84  // writer->SetInput(recon);
85  // writer->SetFileName("Reconstruction.mhd");
86  // writer->Update();
87  // writer->SetInput(ref);
88  // writer->SetFileName("Reference.mhd");
89  // writer->Update();
90  // // End of results writing
91 
92  // Checking results. As a comparison with NaN always returns false,
93  // this design allows to detect NaN results and cause test failure
94  if (!(ErrorPerPixel < ErrorPerPixelTolerance))
95  {
96  std::cerr << "Test Failed, Error per pixel not valid! " << ErrorPerPixel << " instead of " << ErrorPerPixelTolerance
97  << std::endl;
98  exit(EXIT_FAILURE);
99  }
100  if (!(PSNR > PSNRTolerance))
101  {
102  std::cerr << "Test Failed, PSNR not valid! " << PSNR << " instead of " << PSNRTolerance << std::endl;
103  exit(EXIT_FAILURE);
104  }
105 }
106 #endif // FAST_TESTS_NO_CHECKS
107 
108 template <class TImage>
109 #if FAST_TESTS_NO_CHECKS
110 void
111 CheckVectorImageQuality(typename TImage::Pointer itkNotUsed(recon),
112  typename TImage::Pointer itkNotUsed(ref),
113  double itkNotUsed(ErrorPerPixelTolerance),
114  double itkNotUsed(PSNRTolerance),
115  double itkNotUsed(RefValueForPSNR))
116 {}
117 #else
118 void
119 CheckVectorImageQuality(typename TImage::Pointer recon,
120  typename TImage::Pointer ref,
121  double ErrorPerPixelTolerance,
122  double PSNRTolerance,
123  double RefValueForPSNR)
124 {
125  using ImageIteratorType = itk::ImageRegionConstIterator<TImage>;
126  ImageIteratorType itTest(recon, recon->GetBufferedRegion());
127  ImageIteratorType itRef(ref, ref->GetBufferedRegion());
128 
129  using ErrorType = double;
130  ErrorType TestError = 0.;
131  ErrorType EnerError = 0.;
132 
133  itTest.GoToBegin();
134  itRef.GoToBegin();
135 
136  while (!itRef.IsAtEnd())
137  {
138  typename TImage::PixelType TestVal = itTest.Get();
139  typename TImage::PixelType RefVal = itRef.Get();
140  TestError += (RefVal - TestVal).GetNorm();
141  EnerError += (RefVal - TestVal).GetSquaredNorm();
142  ++itTest;
143  ++itRef;
144  }
145  // Error per Pixel
146  ErrorType ErrorPerPixel = TestError / ref->GetBufferedRegion().GetNumberOfPixels();
147  std::cout << "\nError per Pixel = " << ErrorPerPixel << std::endl;
148  // MSE
149  ErrorType MSE = EnerError / ref->GetBufferedRegion().GetNumberOfPixels();
150  std::cout << "MSE = " << MSE << std::endl;
151  // PSNR
152  ErrorType PSNR = 20 * log10(RefValueForPSNR) - 10 * log10(MSE);
153  std::cout << "PSNR = " << PSNR << "dB" << std::endl;
154  // QI
155  ErrorType QI = (RefValueForPSNR - ErrorPerPixel) / RefValueForPSNR;
156  std::cout << "QI = " << QI << std::endl;
157 
158  // // It is often necessary to write the images and look at them
159  // // to understand why a given test fails. This portion of code
160  // // does that. It should be left here but commented out, since
161  // // it is only useful in specific debugging tasks
162  // using FileWriterType = itk::ImageFileWriter<TImage>;
163  // typename FileWriterType::Pointer writer = FileWriterType::New();
164  // writer->SetInput(recon);
165  // writer->SetFileName("Reconstruction.mhd");
166  // writer->Update();
167  // writer->SetInput(ref);
168  // writer->SetFileName("Reference.mhd");
169  // writer->Update();
170  // // End of results writing
171 
172  // Checking results. As a comparison with NaN always returns false,
173  // this design allows to detect NaN results and cause test failure
174  if (!(ErrorPerPixel < ErrorPerPixelTolerance))
175  {
176  std::cerr << "Test Failed, Error per pixel not valid! " << ErrorPerPixel << " instead of " << ErrorPerPixelTolerance
177  << std::endl;
178  exit(EXIT_FAILURE);
179  }
180  if (!(PSNR > PSNRTolerance))
181  {
182  std::cerr << "Test Failed, PSNR not valid! " << PSNR << " instead of " << PSNRTolerance << std::endl;
183  exit(EXIT_FAILURE);
184  }
185 }
186 #endif // FAST_TESTS_NO_CHECKS
187 
188 
189 template <class TImage>
190 #if FAST_TESTS_NO_CHECKS
191 void
192 CheckVariableLengthVectorImageQuality(typename TImage::Pointer itkNotUsed(recon),
193  typename TImage::Pointer itkNotUsed(ref),
194  double itkNotUsed(ErrorPerPixelTolerance),
195  double itkNotUsed(PSNRTolerance),
196  double itkNotUsed(RefValueForPSNR))
197 {}
198 #else
199 void
200 CheckVariableLengthVectorImageQuality(typename TImage::Pointer recon,
201  typename TImage::Pointer ref,
202  double ErrorPerPixelTolerance,
203  double PSNRTolerance,
204  double RefValueForPSNR)
205 {
206  if (!(recon->GetVectorLength() == ref->GetVectorLength()))
207  {
208  std::cerr << "Test Failed, image's vector length is " << recon->GetVectorLength() << " instead of "
209  << ref->GetVectorLength() << std::endl;
210  exit(EXIT_FAILURE);
211  }
212 
213  using ImageIteratorType = itk::ImageRegionConstIterator<TImage>;
214  ImageIteratorType itTest(recon, recon->GetBufferedRegion());
215  ImageIteratorType itRef(ref, ref->GetBufferedRegion());
216 
217  using ErrorType = double;
218  ErrorType TestError = 0.;
219  ErrorType EnerError = 0.;
220 
221  itTest.GoToBegin();
222  itRef.GoToBegin();
223 
224  while (!itRef.IsAtEnd())
225  {
226  typename TImage::PixelType TestVec = itTest.Get();
227  typename TImage::PixelType RefVec = itRef.Get();
228  double accumulatedError = 0;
229  for (unsigned int channel = 0; channel < ref->GetVectorLength(); channel++)
230  {
231  accumulatedError += (RefVec[channel] - TestVec[channel]) * (RefVec[channel] - TestVec[channel]);
232  }
233  TestError += sqrt(accumulatedError);
234  EnerError += accumulatedError;
235  ++itTest;
236  ++itRef;
237  }
238  // Error per Pixel
239  ErrorType ErrorPerPixel = TestError / ref->GetBufferedRegion().GetNumberOfPixels();
240  std::cout << "\nError per Pixel = " << ErrorPerPixel << std::endl;
241  // MSE
242  ErrorType MSE = EnerError / ref->GetBufferedRegion().GetNumberOfPixels();
243  std::cout << "MSE = " << MSE << std::endl;
244  // PSNR
245  ErrorType PSNR = 20 * log10(RefValueForPSNR) - 10 * log10(MSE);
246  std::cout << "PSNR = " << PSNR << "dB" << std::endl;
247  // QI
248  ErrorType QI = (RefValueForPSNR - ErrorPerPixel) / RefValueForPSNR;
249  std::cout << "QI = " << QI << std::endl;
250 
251  // // It is often necessary to write the images and look at them
252  // // to understand why a given test fails. This portion of code
253  // // does that. It should be left here but commented out, since
254  // // it is only useful in specific debugging tasks
255  // using FileWriterType = itk::ImageFileWriter<TImage>;
256  // typename FileWriterType::Pointer writer = FileWriterType::New();
257  // writer->SetInput(recon);
258  // writer->SetFileName("Reconstruction.mhd");
259  // writer->Update();
260  // writer->SetInput(ref);
261  // writer->SetFileName("Reference.mhd");
262  // writer->Update();
263  // // End of results writing
264 
265  // Checking results. As a comparison with NaN always returns false,
266  // this design allows to detect NaN results and cause test failure
267  if (!(ErrorPerPixel < ErrorPerPixelTolerance))
268  {
269  std::cerr << "Test Failed, Error per pixel not valid! " << ErrorPerPixel << " instead of " << ErrorPerPixelTolerance
270  << std::endl;
271  exit(EXIT_FAILURE);
272  }
273  if (!(PSNR > PSNRTolerance))
274  {
275  std::cerr << "Test Failed, PSNR not valid! " << PSNR << " instead of " << PSNRTolerance << std::endl;
276  exit(EXIT_FAILURE);
277  }
278 }
279 #endif // FAST_TESTS_NO_CHECKS
280 
281 void
283 {
284  // // It is often necessary to write the geometries and look at them
285  // // to understand why a given test fails. This portion of code
286  // // does that. It should be left here but commented out, since
287  // // it is only useful in specific debugging tasks
288  // rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter =
289  // rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
290  // xmlWriter->SetFilename("g1.xml");
291  // xmlWriter->SetObject(g1);
292  // xmlWriter->WriteFile();
293  // xmlWriter->SetFilename("g2.xml");
294  // xmlWriter->SetObject(g2);
295  // xmlWriter->WriteFile();
296  // // End of results writing
297 
298  const double e = 1e-10;
299  const unsigned int nproj = g1->GetGantryAngles().size();
300  if (g2->GetGantryAngles().size() != nproj)
301  {
302  std::cerr << "Unequal number of projections in the two geometries" << std::endl;
303  exit(EXIT_FAILURE);
304  }
305  if (e < std::fabs(g1->GetRadiusCylindricalDetector() - g2->GetRadiusCylindricalDetector()))
306  {
307  std::cerr << "Geometries don't have the same cylindrical detector radius" << std::endl;
308  exit(EXIT_FAILURE);
309  }
310 
311  for (unsigned int i = 0; i < nproj; i++)
312  {
313  // std::cout << g1->GetGantryAngles()[i] << " " << g2->GetGantryAngles()[i] << std::endl;
315  std::fabs(g1->GetGantryAngles()[i] - g2->GetGantryAngles()[i])) ||
317  std::fabs(g1->GetOutOfPlaneAngles()[i] - g2->GetOutOfPlaneAngles()[i])) ||
319  std::fabs(g1->GetInPlaneAngles()[i] - g2->GetInPlaneAngles()[i])) ||
320  e < std::fabs(g1->GetSourceToIsocenterDistances()[i] - g2->GetSourceToIsocenterDistances()[i]) ||
321  e < std::fabs(g1->GetSourceOffsetsX()[i] - g2->GetSourceOffsetsX()[i]) ||
322  e < std::fabs(g1->GetSourceOffsetsY()[i] - g2->GetSourceOffsetsY()[i]) ||
323  e < std::fabs(g1->GetSourceToDetectorDistances()[i] - g2->GetSourceToDetectorDistances()[i]) ||
324  e < std::fabs(g1->GetProjectionOffsetsX()[i] - g2->GetProjectionOffsetsX()[i]) ||
325  e < std::fabs(g1->GetProjectionOffsetsY()[i] - g2->GetProjectionOffsetsY()[i]) ||
326  e < std::fabs(g1->GetCollimationUInf()[i] - g2->GetCollimationUInf()[i]) ||
327  e < std::fabs(g1->GetCollimationVInf()[i] - g2->GetCollimationVInf()[i]) ||
328  e < std::fabs(g1->GetCollimationUSup()[i] - g2->GetCollimationUSup()[i]) ||
329  e < std::fabs(g1->GetCollimationVSup()[i] - g2->GetCollimationVSup()[i]))
330  {
331  std::cerr << "Geometry of projection #" << i << " is unvalid." << std::endl;
332  exit(EXIT_FAILURE);
333  }
334  }
335 }
336 
337 
338 template <class TImage1, class TImage2>
339 #if FAST_TESTS_NO_CHECKS
340 void
341 CheckScalarProducts(typename TImage1::Pointer itkNotUsed(im1A),
342  typename TImage1::Pointer itkNotUsed(im1B),
343  typename TImage2::Pointer itkNotUsed(im2A),
344  typename TImage2::Pointer itkNotUsed(im2B))
345 {}
346 #else
347 void
348 CheckScalarProducts(typename TImage1::Pointer im1A,
349  typename TImage1::Pointer im1B,
350  typename TImage2::Pointer im2A,
351  typename TImage2::Pointer im2B)
352 {
353  using Image1IteratorType = itk::ImageRegionConstIterator<TImage1>;
354  using Image2IteratorType = itk::ImageRegionConstIterator<TImage2>;
355  Image1IteratorType itIm1A(im1A, im1A->GetLargestPossibleRegion());
356  Image1IteratorType itIm1B(im1B, im1B->GetLargestPossibleRegion());
357  Image2IteratorType itIm2A(im2A, im2A->GetLargestPossibleRegion());
358  Image2IteratorType itIm2B(im2B, im2B->GetLargestPossibleRegion());
359 
360  typename TImage2::PixelType scalarProductT1, scalarProductT2;
361  scalarProductT1 = 0;
362  scalarProductT2 = 0;
363 
364  while (!itIm1A.IsAtEnd())
365  {
366  scalarProductT1 += itIm1A.Get() * itIm1B.Get();
367  ++itIm1A;
368  ++itIm1B;
369  }
370 
371  while (!itIm2A.IsAtEnd())
372  {
373  scalarProductT2 += itIm2A.Get() * itIm2B.Get();
374  ++itIm2A;
375  ++itIm2B;
376  }
377 
378  // QI
379  double ratio = scalarProductT1 / scalarProductT2;
380  std::cout << "1 - ratio = " << 1 - ratio << std::endl;
381 
382  // // It is often necessary to write the images and look at them
383  // // to understand why a given test fails. This portion of code
384  // // does that. It should be left here but commented out, since
385  // // it is only useful in specific debugging tasks
386  // using FileWriterType1 = itk::ImageFileWriter<TImage1>;
387  // typename FileWriterType1::Pointer writer = FileWriterType1::New();
388  // writer->SetInput(im1A);
389  // writer->SetFileName("im1A.mhd");
390  // writer->Update();
391  // writer->SetInput(im1B);
392  // writer->SetFileName("im1B.mhd");
393  // writer->Update();
394 
395  // using FileWriterType2 = itk::ImageFileWriter<TImage2>;
396  // typename FileWriterType2::Pointer writer2 = FileWriterType2::New();
397  // writer2->SetInput(im2A);
398  // writer2->SetFileName("im2A.mhd");
399  // writer2->Update();
400  // writer2->SetInput(im2B);
401  // writer2->SetFileName("im2B.mhd");
402  // writer2->Update();
403  // // End of results writing
404 
405 
406  // Checking results
407  if (!(itk::Math::abs(ratio - 1) < 0.001))
408  {
409  std::cerr << "Test Failed, ratio not valid! " << ratio << " instead of 1 +/- 0.001" << std::endl;
410  exit(EXIT_FAILURE);
411  }
412 }
413 #endif
414 
415 template <class TImage1, class TImage2>
416 #if FAST_TESTS_NO_CHECKS
417 void
418 CheckVectorScalarProducts(typename TImage1::Pointer im1A,
419  typename TImage1::Pointer im1B,
420  typename TImage2::Pointer im2A,
421  typename TImage2::Pointer im2B)
422 {}
423 #else
424 void
425 CheckVectorScalarProducts(typename TImage1::Pointer im1A,
426  typename TImage1::Pointer im1B,
427  typename TImage2::Pointer im2A,
428  typename TImage2::Pointer im2B)
429 {
430  using Image1IteratorType = itk::ImageRegionConstIterator<TImage1>;
431  using Image2IteratorType = itk::ImageRegionConstIterator<TImage2>;
432  Image1IteratorType itIm1A(im1A, im1A->GetLargestPossibleRegion());
433  Image1IteratorType itIm1B(im1B, im1B->GetLargestPossibleRegion());
434  Image2IteratorType itIm2A(im2A, im2A->GetLargestPossibleRegion());
435  Image2IteratorType itIm2B(im2B, im2B->GetLargestPossibleRegion());
436 
437  typename itk::PixelTraits<typename TImage2::PixelType>::ValueType scalarProductT1, scalarProductT2;
438  scalarProductT1 = 0;
439  scalarProductT2 = 0;
441 
442  while (!itIm1A.IsAtEnd())
443  {
444  for (unsigned int component = 0; component < vectorLength; component++)
445  scalarProductT1 += itIm1A.Get()[component] * itIm1B.Get()[component];
446  ++itIm1A;
447  ++itIm1B;
448  }
449 
450  while (!itIm2A.IsAtEnd())
451  {
452  for (unsigned int component = 0; component < vectorLength; component++)
453  scalarProductT2 += itIm2A.Get()[component] * itIm2B.Get()[component];
454  ++itIm2A;
455  ++itIm2B;
456  }
457 
458  // QI
459  double ratio = scalarProductT1 / scalarProductT2;
460  std::cout << "1 - ratio = " << 1 - ratio << std::endl;
461 
462  // // It is often necessary to write the images and look at them
463  // // to understand why a given test fails. This portion of code
464  // // does that. It should be left here but commented out, since
465  // // it is only useful in specific debugging tasks
466  // using FileWriterType1 = itk::ImageFileWriter<TImage1>;
467  // typename FileWriterType1::Pointer writer = FileWriterType1::New();
468  // writer->SetInput(im1A);
469  // writer->SetFileName("im1A.mhd");
470  // writer->Update();
471  // writer->SetInput(im1B);
472  // writer->SetFileName("im1B.mhd");
473  // writer->Update();
474 
475  // using FileWriterType2 = itk::ImageFileWriter<TImage2>;
476  // typename FileWriterType2::Pointer writer2 = FileWriterType2::New();
477  // writer2->SetInput(im2A);
478  // writer2->SetFileName("im2A.mhd");
479  // writer2->Update();
480  // writer2->SetInput(im2B);
481  // writer2->SetFileName("im2B.mhd");
482  // writer2->Update();
483  // // End of results writing
484 
485 
486  // Checking results
487  if (!(itk::Math::abs(ratio - 1) < 0.001))
488  {
489  std::cerr << "Test Failed, ratio not valid! " << ratio << " instead of 1 +/- 0.001" << std::endl;
490  exit(EXIT_FAILURE);
491  }
492 }
493 #endif
494 
495 #endif // rtkTest_h
void CheckVectorScalarProducts(typename TImage1::Pointer im1A, typename TImage1::Pointer im1B, typename TImage2::Pointer im2A, typename TImage2::Pointer im2B)
Definition: rtkTest.h:425
const std::vector< double > & GetSourceOffsetsY() const
void CheckVectorImageQuality(typename TImage::Pointer recon, typename TImage::Pointer ref, double ErrorPerPixelTolerance, double PSNRTolerance, double RefValueForPSNR)
Definition: rtkTest.h:119
static double ConvertAngleBetween0And2PIRadians(const double a)
void CheckGeometries(const rtk::ThreeDCircularProjectionGeometry *g1, const rtk::ThreeDCircularProjectionGeometry *g2)
Definition: rtkTest.h:282
const std::vector< double > & GetInPlaneAngles() const
typename TPixelType::ValueType ValueType
Projection geometry for a source and a 2-D flat panel.
const std::vector< double > & GetGantryAngles() const
const std::vector< double > & GetSourceToIsocenterDistances() const
const std::vector< double > & GetProjectionOffsetsY() const
const std::vector< double > & GetProjectionOffsetsX() const
void CheckScalarProducts(typename TImage1::Pointer im1A, typename TImage1::Pointer im1B, typename TImage2::Pointer im2A, typename TImage2::Pointer im2B)
Definition: rtkTest.h:348
const std::vector< double > & GetSourceOffsetsX() const
const std::vector< double > & GetCollimationUInf() const
const std::vector< double > & GetSourceToDetectorDistances() const
const std::vector< double > & GetOutOfPlaneAngles() const
const std::vector< double > & GetCollimationVInf() const
const std::vector< double > & GetCollimationVSup() const
const std::vector< double > & GetCollimationUSup() const
void CheckImageQuality(typename TImage::Pointer recon, typename TImage::Pointer ref, double ErrorPerPixelTolerance, double PSNRTolerance, double RefValueForPSNR)
Definition: rtkTest.h:39
void CheckVariableLengthVectorImageQuality(typename TImage::Pointer recon, typename TImage::Pointer ref, double ErrorPerPixelTolerance, double PSNRTolerance, double RefValueForPSNR)
Definition: rtkTest.h:200
virtual double GetRadiusCylindricalDetector() const