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