OTB  5.9.0
Orfeo Toolbox
otbStreamingCompareImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef otbStreamingCompareImageFilter_txx
19 #define otbStreamingCompareImageFilter_txx
21 
22 #include "itkImageRegionIterator.h"
23 #include "itkProgressReporter.h"
24 #include "itkMath.h"
25 #include "otbMacro.h"
26 
27 namespace otb
28 {
29 
30 template<class TInputImage>
32 ::PersistentCompareImageFilter() : m_SquareOfDifferences(1), m_AbsoluteValueOfDifferences(1),
33  m_ThreadMinRef(1), m_ThreadMaxRef(1), m_Count(1), m_DiffCount(1), m_PhysicalSpaceCheck(true)
34 {
35  this->SetNumberOfRequiredInputs( 2 );
36  // first output is a copy of the image, DataObject created by
37  // superclass
38 
39  // allocate the data objects for the outputs which are
40  // just decorators around real types
41  for (int i = 1; i < 5; ++i)
42  {
43  typename RealObjectType::Pointer output
44  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
46  }
47 
48  this->GetPSNROutput()->Set(itk::NumericTraits<RealType>::max());
49  this->GetMSEOutput()->Set(itk::NumericTraits<RealType>::max());
50  this->GetMAEOutput()->Set(itk::NumericTraits<RealType>::max());
51  this->GetDiffCountOutput()->Set(itk::NumericTraits<RealType>::Zero);
52 
53  this->Reset();
54 }
55 
59 template<class TInputImage>
60 void
62 ::SetInput1(const TInputImage *image)
63 {
64  // The ProcessObject is not const-correct so the const_cast is required here
65  this->SetNthInput(0, const_cast<TInputImage *>(image));
66 }
67 
71 template<class TInputImage>
72 void
74 ::SetInput2(const TInputImage *image)
75 {
76  // The ProcessObject is not const-correct so the const_cast is required here
77  this->SetNthInput(1, const_cast<TInputImage *>(image));
78 }
79 
80 template<class TInputImage>
81 const TInputImage *
84 {
85  if (this->GetNumberOfInputs()<1)
86  {
87  return 0;
88  }
89  return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(0));
90 }
91 
92 template<class TInputImage>
93 const TInputImage *
96 {
97  if (this->GetNumberOfInputs()<2)
98  {
99  return 0;
100  }
101  return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(1));
102 }
103 
104 template<class TInputImage>
108 {
109  switch (output)
110  {
111  case 0:
112  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
113  break;
114  case 1:
115  case 2:
116  case 3:
117  case 4:
118  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
119  break;
120  default:
121  // might as well make an image
122  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
123  break;
124  }
125 }
126 
127 template<class TInputImage>
131 {
132  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
133 }
134 
135 template<class TInputImage>
139 {
140  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
141 }
142 
143 template<class TInputImage>
147 {
148  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
149 }
150 
151 template<class TInputImage>
155 {
156  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
157 }
158 
159 template<class TInputImage>
163 {
164  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
165 }
166 
167 template<class TInputImage>
171 {
172  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
173 }
174 
175 template<class TInputImage>
179 {
180  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
181 }
182 
183 template<class TInputImage>
187 {
188  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
189 }
190 
191 template<class TInputImage>
192 void
195 {
196  Superclass::GenerateOutputInformation();
197  if (this->GetInput())
198  {
199  this->GetOutput()->CopyInformation(this->GetInput());
200  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
201 
202  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
203  {
204  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
205  }
206  }
207 }
208 template<class TInputImage>
209 void
212 {
213  // This is commented to prevent the streaming of the whole image for the first stream strip
214  // It shall not cause any problem because the output image of this filter is not intended to be used.
215  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
216  //this->GraftOutput( image );
217  // Nothing that needs to be allocated for the remaining outputs
218 }
219 
220 template<class TInputImage>
221 void
224 {
225  int i;
226  unsigned long count;
227  unsigned long diffCount;
228  RealType squareOfDifferences, absoluteValueOfDifferences;
229 
230  int numberOfThreads = this->GetNumberOfThreads();
231 
232  PixelType minimumRef, maximumRef;
233  RealType mse;
234  RealType mae;
235  RealType psnr;
236 
237  squareOfDifferences = absoluteValueOfDifferences = itk::NumericTraits<RealType>::Zero;
238  count = 0;
239  diffCount=0;
240 
241  // Find min/max and the accumulate count and difference of squares over all threads
242  minimumRef = itk::NumericTraits<PixelType>::max();
244 
245  for (i = 0; i < numberOfThreads; ++i)
246  {
247  count += m_Count[i];
248  diffCount += m_DiffCount[i];
249  squareOfDifferences += m_SquareOfDifferences[i];
250  absoluteValueOfDifferences += m_AbsoluteValueOfDifferences[i];
251 
252  if (m_ThreadMinRef[i] < minimumRef)
253  {
254  minimumRef = m_ThreadMinRef[i];
255  }
256  if (m_ThreadMaxRef[i] > maximumRef)
257  {
258  maximumRef = m_ThreadMaxRef[i];
259  }
260  }
261 
262  // compute mse
263  mse = (count != 0) ? squareOfDifferences / static_cast<RealType>(count) : 0.;
264  // compute mse
265  mae = (count != 0) ? absoluteValueOfDifferences / static_cast<RealType>(count) : 0.;
266 
267  //compute psnr
268  psnr = (vcl_abs(mse)>0.0000000001 && (maximumRef-minimumRef)>0.0000000001) ? 10.*vcl_log10(((maximumRef-minimumRef)*(maximumRef-minimumRef))/mse):0.;
269  // Set the output
270  this->GetMSEOutput()->Set(mse);
271  this->GetMAEOutput()->Set(mae);
272  this->GetPSNROutput()->Set(psnr);
273  this->GetDiffCountOutput()->Set(static_cast<RealType>(diffCount));
274 }
275 
276 template<class TInputImage>
277 void
280 {
281  int numberOfThreads = this->GetNumberOfThreads();
282 
283  // Resize the thread temporaries
284  m_Count.SetSize(numberOfThreads);
285  m_DiffCount.SetSize(numberOfThreads);
286  m_SquareOfDifferences.SetSize(numberOfThreads);
287  m_AbsoluteValueOfDifferences.SetSize(numberOfThreads);
288 
289  m_ThreadMinRef.SetSize(numberOfThreads);
290  m_ThreadMaxRef.SetSize(numberOfThreads);
291 
292  // Initialize the temporaries
293  m_Count.Fill(itk::NumericTraits<long>::Zero);
294  m_DiffCount.Fill(itk::NumericTraits<long>::Zero);
295  m_SquareOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
296  m_AbsoluteValueOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
297  m_ThreadMinRef.Fill(itk::NumericTraits<PixelType>::max());
298  m_ThreadMaxRef.Fill(itk::NumericTraits<PixelType>::NonpositiveMin());
299 }
300 
301 template<class TInputImage>
302 void
305 {
306  if (m_PhysicalSpaceCheck)
307  Superclass::VerifyInputInformation();
308 }
309 
310 template<class TInputImage>
311 void
313 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
314  itk::ThreadIdType threadId)
315 {
319  InputImagePointer inputPtr1 = const_cast<TInputImage *>(this->GetInput(0));
320  InputImagePointer inputPtr2 = const_cast<TInputImage *>(this->GetInput(1));
322 
323  // support progress methods/callbacks
324  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
325 
326  RealType realValue1, realValue2;
327  PixelType value1, value2;
328 
329  itk::ImageRegionConstIterator<TInputImage> it1(inputPtr1, outputRegionForThread);
330  itk::ImageRegionConstIterator<TInputImage> it2(inputPtr2, outputRegionForThread);
331 
332  it1.GoToBegin();
333  it2.GoToBegin();
334  // do the work
335  while (!it1.IsAtEnd() && !it2.IsAtEnd())
336  {
337  value1 = it1.Get();
338  realValue1 = static_cast<RealType>(value1);
339 
340  value2 = it2.Get();
341  realValue2 = static_cast<RealType>(value2);
342 
343  if (value1 < m_ThreadMinRef[threadId])
344  {
345  m_ThreadMinRef[threadId] = value1;
346  }
347  if (value1 > m_ThreadMaxRef[threadId])
348  {
349  m_ThreadMaxRef[threadId] = value1;
350  }
351 
352  RealType diffVal = realValue1 - realValue2;
353  m_SquareOfDifferences[threadId] += diffVal * diffVal;
354  m_AbsoluteValueOfDifferences[threadId] += vcl_abs( diffVal );
355  if (! itk::Math::FloatAlmostEqual(realValue1, realValue2))
356  {
357  m_DiffCount[threadId]++;
358  }
359  m_Count[threadId]++;
360  ++it1;
361  ++it2;
362  progress.CompletedPixel();
363  }
364 }
365 template <class TImage>
366 void
368 ::PrintSelf(std::ostream& os, itk::Indent indent) const
369 {
370  Superclass::PrintSelf(os, indent);
371 
372  os << indent << "PSNR: " << this->GetPSNR() << std::endl;
373  os << indent << "MSE: " << this->GetMSE() << std::endl;
374  os << indent << "MAE: " << this->GetMAE() << std::endl;
375  os << indent << "Count: " << this->GetDiffCount() << std::endl;
376 }
377 } // end namespace otb
378 #endif
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE
ObjectType * GetPointer() const
static ITK_CONSTEXPR_FUNC T max(const T &)
unsigned int ThreadIdType
DataObject * GetInput(const DataObjectIdentifierType &key)
static ITK_CONSTEXPR_FUNC T NonpositiveMin()
bool FloatAlmostEqual(T x1, T x2, typename Detail::FloatIEEE< T >::IntType maxUlps=4, typename Detail::FloatIEEE< T >::FloatType maxAbsoluteDifference=0.1 *itk::NumericTraits< T >::epsilon())
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) ITK_OVERRIDE
itk::NumericTraits< PixelType >::RealType RealType
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
bool IsAtEnd(void) const
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
DataObject * GetOutput(const DataObjectIdentifierType &key)