Orfeo Toolbox  4.0
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"
24 #include "itkNumericTraits.h"
25 #include "itkProgressReporter.h"
26 #include "otbMacro.h"
27 
28 namespace otb
29 {
30 
31 template<class TInputImage>
33 ::PersistentCompareImageFilter() : m_SquareOfDifferences(1), m_AbsoluteValueOfDifferences(1),
34  m_ThreadMinRef(1), m_ThreadMaxRef(1), m_Count(1)
35 {
36  this->SetNumberOfRequiredInputs( 2 );
37  // first output is a copy of the image, DataObject created by
38  // superclass
39 
40  // allocate the data objects for the outputs which are
41  // just decorators around real types
42  for (int i = 1; i < 4; ++i)
43  {
44  typename RealObjectType::Pointer output
45  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
47  }
48 
49  this->GetPSNROutput()->Set(itk::NumericTraits<RealType>::max());
50  this->GetMSEOutput()->Set(itk::NumericTraits<RealType>::max());
51  this->GetMAEOutput()->Set(itk::NumericTraits<RealType>::max());
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>
107 ::MakeOutput(unsigned int output)
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  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
118  break;
119  default:
120  // might as well make an image
121  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
122  break;
123  }
124 }
125 
126 template<class TInputImage>
130 {
131  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
132 }
133 
134 template<class TInputImage>
138 {
139  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
140 }
141 
142 template<class TInputImage>
146 {
147  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
148 }
149 
150 template<class TInputImage>
154 {
155  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
156 }
157 
158 template<class TInputImage>
162 {
163  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
164 }
165 
166 template<class TInputImage>
170 {
171  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
172 }
173 
174 template<class TInputImage>
175 void
178 {
179  Superclass::GenerateOutputInformation();
180  if (this->GetInput())
181  {
182  this->GetOutput()->CopyInformation(this->GetInput());
183  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
184 
185  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
186  {
187  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
188  }
189  }
190 }
191 template<class TInputImage>
192 void
195 {
196  // This is commented to prevent the streaming of the whole image for the first stream strip
197  // It shall not cause any problem because the output image of this filter is not intended to be used.
198  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
199  //this->GraftOutput( image );
200  // Nothing that needs to be allocated for the remaining outputs
201 }
202 
203 template<class TInputImage>
204 void
207 {
208  int i;
209  long count;
210  RealType squareOfDifferences, absoluteValueOfDifferences;
211 
212  int numberOfThreads = this->GetNumberOfThreads();
213 
214  PixelType minimumRef, maximumRef;
215  RealType mse;
216  RealType mae;
217  RealType psnr;
218 
219  squareOfDifferences = absoluteValueOfDifferences = itk::NumericTraits<RealType>::Zero;
220  count = 0;
221 
222  // Find min/max and the accumulate count and difference of squares over all threads
223  minimumRef = itk::NumericTraits<PixelType>::max();
225 
226  for (i = 0; i < numberOfThreads; ++i)
227  {
228  count += m_Count[i];
229  squareOfDifferences += m_SquareOfDifferences[i];
230  absoluteValueOfDifferences += m_AbsoluteValueOfDifferences[i];
231 
232  if (m_ThreadMinRef[i] < minimumRef)
233  {
234  minimumRef = m_ThreadMinRef[i];
235  }
236  if (m_ThreadMaxRef[i] > maximumRef)
237  {
238  maximumRef = m_ThreadMaxRef[i];
239  }
240  }
241 
242  // compute mse
243  mse = (count != 0) ? squareOfDifferences / static_cast<RealType>(count) : 0.;
244  // compute mse
245  mae = (count != 0) ? absoluteValueOfDifferences / static_cast<RealType>(count) : 0.;
246 
247  //compute psnr
248  psnr = (vcl_abs(mse)>0.0000000001 && (maximumRef-minimumRef)>0.0000000001) ? 10.*vcl_log10(((maximumRef-minimumRef)*(maximumRef-minimumRef))/mse):0.;
249  // Set the output
250  this->GetMSEOutput()->Set(mse);
251  this->GetMAEOutput()->Set(mae);
252  this->GetPSNROutput()->Set(psnr);
253 }
254 
255 template<class TInputImage>
256 void
259 {
260  int numberOfThreads = this->GetNumberOfThreads();
261 
262  // Resize the thread temporaries
263  m_Count.SetSize(numberOfThreads);
264  m_SquareOfDifferences.SetSize(numberOfThreads);
265  m_AbsoluteValueOfDifferences.SetSize(numberOfThreads);
266 
267  m_ThreadMinRef.SetSize(numberOfThreads);
268  m_ThreadMaxRef.SetSize(numberOfThreads);
269 
270  // Initialize the temporaries
271  m_Count.Fill(itk::NumericTraits<long>::Zero);
272  m_SquareOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
273  m_AbsoluteValueOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
274  m_ThreadMinRef.Fill(itk::NumericTraits<PixelType>::max());
275  m_ThreadMaxRef.Fill(itk::NumericTraits<PixelType>::NonpositiveMin());
276 }
277 
278 template<class TInputImage>
279 void
281 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
282  itk::ThreadIdType threadId)
283 {
287  InputImagePointer inputPtr1 = const_cast<TInputImage *>(this->GetInput(0));
288  InputImagePointer inputPtr2 = const_cast<TInputImage *>(this->GetInput(1));
289 
290  // support progress methods/callbacks
291  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
292 
293  RealType realValue1, realValue2;
294  PixelType value1, value2;
295 
296  itk::ImageRegionConstIterator<TInputImage> it1(inputPtr1, outputRegionForThread);
297  itk::ImageRegionConstIterator<TInputImage> it2(inputPtr2, outputRegionForThread);
298 
299  it1.GoToBegin();
300  it2.GoToBegin();
301  // do the work
302  while (!it1.IsAtEnd() && !it2.IsAtEnd())
303  {
304  value1 = it1.Get();
305  realValue1 = static_cast<RealType>(value1);
306 
307  value2 = it2.Get();
308  realValue2 = static_cast<RealType>(value2);
309 
310  if (value1 < m_ThreadMinRef[threadId])
311  {
312  m_ThreadMinRef[threadId] = value1;
313  }
314  if (value1 > m_ThreadMaxRef[threadId])
315  {
316  m_ThreadMaxRef[threadId] = value1;
317  }
318 
319  m_SquareOfDifferences[threadId] += ( realValue1 - realValue2 ) * ( realValue1 - realValue2 );
320  m_AbsoluteValueOfDifferences[threadId] += vcl_abs( realValue1 - realValue2 );
321  m_Count[threadId]++;
322  ++it1;
323  ++it2;
324  progress.CompletedPixel();
325  }
326 }
327 template <class TImage>
328 void
330 ::PrintSelf(std::ostream& os, itk::Indent indent) const
331 {
332  Superclass::PrintSelf(os, indent);
333 
334  os << indent << "PSNR: " << this->GetPSNR() << std::endl;
335  os << indent << "MSE: " << this->GetMSE() << std::endl;
336  os << indent << "MAE: " << this->GetMAE() << std::endl;
337 }
338 } // end namespace otb
339 #endif

Generated at Sat Mar 8 2014 16:19:32 for Orfeo Toolbox with doxygen 1.8.3.1