OTB  5.0.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 "otbMacro.h"
25 
26 namespace otb
27 {
28 
29 template<class TInputImage>
31 ::PersistentCompareImageFilter() : m_SquareOfDifferences(1), m_AbsoluteValueOfDifferences(1),
32  m_ThreadMinRef(1), m_ThreadMaxRef(1), m_Count(1), m_PhysicalSpaceCheck(true)
33 {
34  this->SetNumberOfRequiredInputs( 2 );
35  // first output is a copy of the image, DataObject created by
36  // superclass
37 
38  // allocate the data objects for the outputs which are
39  // just decorators around real types
40  for (int i = 1; i < 4; ++i)
41  {
42  typename RealObjectType::Pointer output
43  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
45  }
46 
47  this->GetPSNROutput()->Set(itk::NumericTraits<RealType>::max());
48  this->GetMSEOutput()->Set(itk::NumericTraits<RealType>::max());
49  this->GetMAEOutput()->Set(itk::NumericTraits<RealType>::max());
50 
51  this->Reset();
52 }
53 
57 template<class TInputImage>
58 void
60 ::SetInput1(const TInputImage *image)
61 {
62  // The ProcessObject is not const-correct so the const_cast is required here
63  this->SetNthInput(0, const_cast<TInputImage *>(image));
64 }
65 
69 template<class TInputImage>
70 void
72 ::SetInput2(const TInputImage *image)
73 {
74  // The ProcessObject is not const-correct so the const_cast is required here
75  this->SetNthInput(1, const_cast<TInputImage *>(image));
76 }
77 
78 template<class TInputImage>
79 const TInputImage *
82 {
83  if (this->GetNumberOfInputs()<1)
84  {
85  return 0;
86  }
87  return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(0));
88 }
89 
90 template<class TInputImage>
91 const TInputImage *
94 {
95  if (this->GetNumberOfInputs()<2)
96  {
97  return 0;
98  }
99  return static_cast<const TInputImage *>(this->itk::ProcessObject::GetInput(1));
100 }
101 
102 template<class TInputImage>
106 {
107  switch (output)
108  {
109  case 0:
110  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
111  break;
112  case 1:
113  case 2:
114  case 3:
115  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
116  break;
117  default:
118  // might as well make an image
119  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
120  break;
121  }
122 }
123 
124 template<class TInputImage>
128 {
129  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
130 }
131 
132 template<class TInputImage>
136 {
137  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
138 }
139 
140 template<class TInputImage>
144 {
145  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
146 }
147 
148 template<class TInputImage>
152 {
153  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
154 }
155 
156 template<class TInputImage>
160 {
161  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
162 }
163 
164 template<class TInputImage>
168 {
169  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
170 }
171 
172 template<class TInputImage>
173 void
176 {
177  Superclass::GenerateOutputInformation();
178  if (this->GetInput())
179  {
180  this->GetOutput()->CopyInformation(this->GetInput());
181  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
182 
183  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
184  {
185  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
186  }
187  }
188 }
189 template<class TInputImage>
190 void
193 {
194  // This is commented to prevent the streaming of the whole image for the first stream strip
195  // It shall not cause any problem because the output image of this filter is not intended to be used.
196  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
197  //this->GraftOutput( image );
198  // Nothing that needs to be allocated for the remaining outputs
199 }
200 
201 template<class TInputImage>
202 void
205 {
206  int i;
207  long count;
208  RealType squareOfDifferences, absoluteValueOfDifferences;
209 
210  int numberOfThreads = this->GetNumberOfThreads();
211 
212  PixelType minimumRef, maximumRef;
213  RealType mse;
214  RealType mae;
215  RealType psnr;
216 
217  squareOfDifferences = absoluteValueOfDifferences = itk::NumericTraits<RealType>::Zero;
218  count = 0;
219 
220  // Find min/max and the accumulate count and difference of squares over all threads
221  minimumRef = itk::NumericTraits<PixelType>::max();
223 
224  for (i = 0; i < numberOfThreads; ++i)
225  {
226  count += m_Count[i];
227  squareOfDifferences += m_SquareOfDifferences[i];
228  absoluteValueOfDifferences += m_AbsoluteValueOfDifferences[i];
229 
230  if (m_ThreadMinRef[i] < minimumRef)
231  {
232  minimumRef = m_ThreadMinRef[i];
233  }
234  if (m_ThreadMaxRef[i] > maximumRef)
235  {
236  maximumRef = m_ThreadMaxRef[i];
237  }
238  }
239 
240  // compute mse
241  mse = (count != 0) ? squareOfDifferences / static_cast<RealType>(count) : 0.;
242  // compute mse
243  mae = (count != 0) ? absoluteValueOfDifferences / static_cast<RealType>(count) : 0.;
244 
245  //compute psnr
246  psnr = (vcl_abs(mse)>0.0000000001 && (maximumRef-minimumRef)>0.0000000001) ? 10.*vcl_log10(((maximumRef-minimumRef)*(maximumRef-minimumRef))/mse):0.;
247  // Set the output
248  this->GetMSEOutput()->Set(mse);
249  this->GetMAEOutput()->Set(mae);
250  this->GetPSNROutput()->Set(psnr);
251 }
252 
253 template<class TInputImage>
254 void
257 {
258  int numberOfThreads = this->GetNumberOfThreads();
259 
260  // Resize the thread temporaries
261  m_Count.SetSize(numberOfThreads);
262  m_SquareOfDifferences.SetSize(numberOfThreads);
263  m_AbsoluteValueOfDifferences.SetSize(numberOfThreads);
264 
265  m_ThreadMinRef.SetSize(numberOfThreads);
266  m_ThreadMaxRef.SetSize(numberOfThreads);
267 
268  // Initialize the temporaries
269  m_Count.Fill(itk::NumericTraits<long>::Zero);
270  m_SquareOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
271  m_AbsoluteValueOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
272  m_ThreadMinRef.Fill(itk::NumericTraits<PixelType>::max());
273  m_ThreadMaxRef.Fill(itk::NumericTraits<PixelType>::NonpositiveMin());
274 }
275 
276 template<class TInputImage>
277 void
280 {
281  if (m_PhysicalSpaceCheck)
282  Superclass::VerifyInputInformation();
283 }
284 
285 template<class TInputImage>
286 void
288 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
289  itk::ThreadIdType threadId)
290 {
294  InputImagePointer inputPtr1 = const_cast<TInputImage *>(this->GetInput(0));
295  InputImagePointer inputPtr2 = const_cast<TInputImage *>(this->GetInput(1));
297 
298  // support progress methods/callbacks
299  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
300 
301  RealType realValue1, realValue2;
302  PixelType value1, value2;
303 
304  itk::ImageRegionConstIterator<TInputImage> it1(inputPtr1, outputRegionForThread);
305  itk::ImageRegionConstIterator<TInputImage> it2(inputPtr2, outputRegionForThread);
306 
307  it1.GoToBegin();
308  it2.GoToBegin();
309  // do the work
310  while (!it1.IsAtEnd() && !it2.IsAtEnd())
311  {
312  value1 = it1.Get();
313  realValue1 = static_cast<RealType>(value1);
314 
315  value2 = it2.Get();
316  realValue2 = static_cast<RealType>(value2);
317 
318  if (value1 < m_ThreadMinRef[threadId])
319  {
320  m_ThreadMinRef[threadId] = value1;
321  }
322  if (value1 > m_ThreadMaxRef[threadId])
323  {
324  m_ThreadMaxRef[threadId] = value1;
325  }
326 
327  m_SquareOfDifferences[threadId] += ( realValue1 - realValue2 ) * ( realValue1 - realValue2 );
328  m_AbsoluteValueOfDifferences[threadId] += vcl_abs( realValue1 - realValue2 );
329  m_Count[threadId]++;
330  ++it1;
331  ++it2;
332  progress.CompletedPixel();
333  }
334 }
335 template <class TImage>
336 void
338 ::PrintSelf(std::ostream& os, itk::Indent indent) const
339 {
340  Superclass::PrintSelf(os, indent);
341 
342  os << indent << "PSNR: " << this->GetPSNR() << std::endl;
343  os << indent << "MSE: " << this->GetMSE() << std::endl;
344  os << indent << "MAE: " << this->GetMAE() << std::endl;
345 }
346 } // end namespace otb
347 #endif
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx)
ObjectType * GetPointer() const
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId)
static T max(const T &)
DataObject * GetInput(const DataObjectIdentifierType &key)
void PrintSelf(std::ostream &os, itk::Indent indent) const
itk::NumericTraits< PixelType >::RealType RealType
static T NonpositiveMin()
bool IsAtEnd(void) const
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
unsigned int ThreadIdType
DataObject * GetOutput(const DataObjectIdentifierType &key)