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