OTB  6.1.0
Orfeo Toolbox
otbStreamingCompareImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbStreamingCompareImageFilter_txx
22 #define otbStreamingCompareImageFilter_txx
24 
25 #include "itkImageRegionIterator.h"
26 #include "itkProgressReporter.h"
27 #include "itkMath.h"
28 #include "otbMacro.h"
29 
30 namespace otb
31 {
32 
33 template<class TInputImage>
35 ::PersistentCompareImageFilter() : m_SquareOfDifferences(1), m_AbsoluteValueOfDifferences(1),
36  m_ThreadMinRef(1), m_ThreadMaxRef(1), m_Count(1), m_DiffCount(1), m_PhysicalSpaceCheck(true)
37 {
38  this->SetNumberOfRequiredInputs( 2 );
39  // first output is a copy of the image, DataObject created by
40  // superclass
41 
42  // allocate the data objects for the outputs which are
43  // just decorators around real types
44  for (int i = 1; i < 5; ++i)
45  {
46  typename RealObjectType::Pointer output
47  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
49  }
50 
51  this->GetPSNROutput()->Set(itk::NumericTraits<RealType>::max());
52  this->GetMSEOutput()->Set(itk::NumericTraits<RealType>::max());
53  this->GetMAEOutput()->Set(itk::NumericTraits<RealType>::max());
54  this->GetDiffCountOutput()->Set(itk::NumericTraits<RealType>::Zero);
55 
56  this->Reset();
57 }
58 
62 template<class TInputImage>
63 void
65 ::SetInput1(const TInputImage *image)
66 {
67  // The ProcessObject is not const-correct so the const_cast is required here
68  this->SetNthInput(0, const_cast<TInputImage *>(image));
69 }
70 
74 template<class TInputImage>
75 void
77 ::SetInput2(const TInputImage *image)
78 {
79  // The ProcessObject is not const-correct so the const_cast is required here
80  this->SetNthInput(1, const_cast<TInputImage *>(image));
81 }
82 
83 template<class TInputImage>
84 const TInputImage *
87 {
88  if (this->GetNumberOfInputs()<1)
89  {
90  return 0;
91  }
92  return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(0));
93 }
94 
95 template<class TInputImage>
96 const TInputImage *
99 {
100  if (this->GetNumberOfInputs()<2)
101  {
102  return 0;
103  }
104  return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(1));
105 }
106 
107 template<class TInputImage>
111 {
112  switch (output)
113  {
114  case 0:
115  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
116  break;
117  case 1:
118  case 2:
119  case 3:
120  case 4:
121  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
122  break;
123  default:
124  // might as well make an image
125  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
126  break;
127  }
128 }
129 
130 template<class TInputImage>
134 {
135  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
136 }
137 
138 template<class TInputImage>
142 {
143  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
144 }
145 
146 template<class TInputImage>
150 {
151  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
152 }
153 
154 template<class TInputImage>
158 {
159  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
160 }
161 
162 template<class TInputImage>
166 {
167  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
168 }
169 
170 template<class TInputImage>
174 {
175  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
176 }
177 
178 template<class TInputImage>
182 {
183  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
184 }
185 
186 template<class TInputImage>
190 {
191  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
192 }
193 
194 template<class TInputImage>
195 void
198 {
199  Superclass::GenerateOutputInformation();
200  if (this->GetInput())
201  {
202  this->GetOutput()->CopyInformation(this->GetInput());
203  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
204 
205  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
206  {
207  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
208  }
209  }
210 }
211 template<class TInputImage>
212 void
215 {
216  // This is commented to prevent the streaming of the whole image for the first stream strip
217  // It shall not cause any problem because the output image of this filter is not intended to be used.
218  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
219  //this->GraftOutput( image );
220  // Nothing that needs to be allocated for the remaining outputs
221 }
222 
223 template<class TInputImage>
224 void
227 {
228  int i;
229  unsigned long count;
230  unsigned long diffCount;
231  RealType squareOfDifferences, absoluteValueOfDifferences;
232 
233  int numberOfThreads = this->GetNumberOfThreads();
234 
235  PixelType minimumRef, maximumRef;
236  RealType mse;
237  RealType mae;
238  RealType psnr;
239 
240  squareOfDifferences = absoluteValueOfDifferences = itk::NumericTraits<RealType>::Zero;
241  count = 0;
242  diffCount=0;
243 
244  // Find min/max and the accumulate count and difference of squares over all threads
245  minimumRef = itk::NumericTraits<PixelType>::max();
247 
248  for (i = 0; i < numberOfThreads; ++i)
249  {
250  count += m_Count[i];
251  diffCount += m_DiffCount[i];
252  squareOfDifferences += m_SquareOfDifferences[i];
253  absoluteValueOfDifferences += m_AbsoluteValueOfDifferences[i];
254 
255  if (m_ThreadMinRef[i] < minimumRef)
256  {
257  minimumRef = m_ThreadMinRef[i];
258  }
259  if (m_ThreadMaxRef[i] > maximumRef)
260  {
261  maximumRef = m_ThreadMaxRef[i];
262  }
263  }
264 
265  // compute mse
266  mse = (count != 0) ? squareOfDifferences / static_cast<RealType>(count) : 0.;
267  // compute mse
268  mae = (count != 0) ? absoluteValueOfDifferences / static_cast<RealType>(count) : 0.;
269 
270  //compute psnr
271  psnr = (vcl_abs(mse)>0.0000000001 && (maximumRef-minimumRef)>0.0000000001) ? 10.*vcl_log10(((maximumRef-minimumRef)*(maximumRef-minimumRef))/mse):0.;
272  // Set the output
273  this->GetMSEOutput()->Set(mse);
274  this->GetMAEOutput()->Set(mae);
275  this->GetPSNROutput()->Set(psnr);
276  this->GetDiffCountOutput()->Set(static_cast<RealType>(diffCount));
277 }
278 
279 template<class TInputImage>
280 void
283 {
284  int numberOfThreads = this->GetNumberOfThreads();
285 
286  // Resize the thread temporaries
287  m_Count.SetSize(numberOfThreads);
288  m_DiffCount.SetSize(numberOfThreads);
289  m_SquareOfDifferences.SetSize(numberOfThreads);
290  m_AbsoluteValueOfDifferences.SetSize(numberOfThreads);
291 
292  m_ThreadMinRef.SetSize(numberOfThreads);
293  m_ThreadMaxRef.SetSize(numberOfThreads);
294 
295  // Initialize the temporaries
296  m_Count.Fill(itk::NumericTraits<long>::Zero);
297  m_DiffCount.Fill(itk::NumericTraits<long>::Zero);
298  m_SquareOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
299  m_AbsoluteValueOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
300  m_ThreadMinRef.Fill(itk::NumericTraits<PixelType>::max());
301  m_ThreadMaxRef.Fill(itk::NumericTraits<PixelType>::NonpositiveMin());
302 }
303 
304 template<class TInputImage>
305 void
308 {
309  if (m_PhysicalSpaceCheck)
310  Superclass::VerifyInputInformation();
311 }
312 
313 template<class TInputImage>
314 void
316 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
317  itk::ThreadIdType threadId)
318 {
322  InputImagePointer inputPtr1 = const_cast<TInputImage *>(this->GetInput(0));
323  InputImagePointer inputPtr2 = const_cast<TInputImage *>(this->GetInput(1));
325 
326  // support progress methods/callbacks
327  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
328 
329  RealType realValue1, realValue2;
330  PixelType value1, value2;
331 
332  itk::ImageRegionConstIterator<TInputImage> it1(inputPtr1, outputRegionForThread);
333  itk::ImageRegionConstIterator<TInputImage> it2(inputPtr2, outputRegionForThread);
334 
335  it1.GoToBegin();
336  it2.GoToBegin();
337  // do the work
338  while (!it1.IsAtEnd() && !it2.IsAtEnd())
339  {
340  value1 = it1.Get();
341  realValue1 = static_cast<RealType>(value1);
342 
343  value2 = it2.Get();
344  realValue2 = static_cast<RealType>(value2);
345 
346  if (value1 < m_ThreadMinRef[threadId])
347  {
348  m_ThreadMinRef[threadId] = value1;
349  }
350  if (value1 > m_ThreadMaxRef[threadId])
351  {
352  m_ThreadMaxRef[threadId] = value1;
353  }
354 
355  RealType diffVal = realValue1 - realValue2;
356  m_SquareOfDifferences[threadId] += diffVal * diffVal;
357  m_AbsoluteValueOfDifferences[threadId] += vcl_abs( diffVal );
358  if (! itk::Math::FloatAlmostEqual(realValue1, realValue2))
359  {
360  m_DiffCount[threadId]++;
361  }
362  m_Count[threadId]++;
363  ++it1;
364  ++it2;
365  progress.CompletedPixel();
366  }
367 }
368 template <class TImage>
369 void
371 ::PrintSelf(std::ostream& os, itk::Indent indent) const
372 {
373  Superclass::PrintSelf(os, indent);
374 
375  os << indent << "PSNR: " << this->GetPSNR() << std::endl;
376  os << indent << "MSE: " << this->GetMSE() << std::endl;
377  os << indent << "MAE: " << this->GetMAE() << std::endl;
378  os << indent << "Count: " << this->GetDiffCount() << std::endl;
379 }
380 } // end namespace otb
381 #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)