OTB  7.2.0
Orfeo Toolbox
otbStreamingCompareImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2020 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_hxx
22 #define otbStreamingCompareImageFilter_hxx
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  : m_SquareOfDifferences(1), m_AbsoluteValueOfDifferences(1), m_ThreadMinRef(1), m_ThreadMaxRef(1), m_Count(1), m_DiffCount(1), m_PhysicalSpaceCheck(true)
36 {
37  this->SetNumberOfRequiredInputs(2);
38  // first output is a copy of the image, DataObject created by
39  // superclass
40 
41  // allocate the data objects for the outputs which are
42  // just decorators around real types
43  for (int i = 1; i < 5; ++i)
44  {
45  typename RealObjectType::Pointer output = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
46  this->itk::ProcessObject::SetNthOutput(i, output.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  this->GetDiffCountOutput()->Set(itk::NumericTraits<RealType>::Zero);
53 
54  this->Reset();
55 }
56 
60 template <class TInputImage>
62 {
63  // The ProcessObject is not const-correct so the const_cast is required here
64  this->SetNthInput(0, const_cast<TInputImage*>(image));
65 }
66 
70 template <class TInputImage>
72 {
73  // The ProcessObject is not const-correct so the const_cast is required here
74  this->SetNthInput(1, const_cast<TInputImage*>(image));
75 }
76 
77 template <class TInputImage>
79 {
80  if (this->GetNumberOfInputs() < 1)
81  {
82  return 0;
83  }
84  return static_cast<const TInputImage*>(this->itk::ProcessObject::GetInput(0));
85 }
86 
87 template <class TInputImage>
89 {
90  if (this->GetNumberOfInputs() < 2)
91  {
92  return 0;
93  }
94  return static_cast<const TInputImage*>(this->itk::ProcessObject::GetInput(1));
95 }
96 
97 template <class TInputImage>
99 {
100  switch (output)
101  {
102  case 0:
103  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
104  break;
105  case 1:
106  case 2:
107  case 3:
108  case 4:
109  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
110  break;
111  default:
112  // might as well make an image
113  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
114  break;
115  }
116 }
117 
118 template <class TInputImage>
120 {
121  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
122 }
123 
124 template <class TInputImage>
126 {
127  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(1));
128 }
129 
130 template <class TInputImage>
132 {
133  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
134 }
135 
136 template <class TInputImage>
138 {
139  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(2));
140 }
141 
142 template <class TInputImage>
144 {
145  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
146 }
147 
148 template <class TInputImage>
150 {
151  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
152 }
153 
154 template <class TInputImage>
156 {
157  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
158 }
159 
160 template <class TInputImage>
162 {
163  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
164 }
165 
166 template <class TInputImage>
168 {
169  Superclass::GenerateOutputInformation();
170  if (this->GetInput())
171  {
172  this->GetOutput()->CopyInformation(this->GetInput());
173  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
174 
175  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
176  {
177  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
178  }
179  }
180 }
181 template <class TInputImage>
183 {
184  // This is commented to prevent the streaming of the whole image for the first stream strip
185  // It shall not cause any problem because the output image of this filter is not intended to be used.
186  // InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
187  // this->GraftOutput( image );
188  // Nothing that needs to be allocated for the remaining outputs
189 }
190 
191 template <class TInputImage>
193 {
194  int i;
195  unsigned long count;
196  unsigned long diffCount;
197  RealType squareOfDifferences, absoluteValueOfDifferences;
198 
199  int numberOfThreads = this->GetNumberOfThreads();
200 
201  PixelType minimumRef, maximumRef;
202  RealType mse;
203  RealType mae;
204  RealType psnr;
205 
206  squareOfDifferences = absoluteValueOfDifferences = itk::NumericTraits<RealType>::Zero;
207  count = 0;
208  diffCount = 0;
209 
210  // Find min/max and the accumulate count and difference of squares over all threads
211  minimumRef = itk::NumericTraits<PixelType>::max();
212  maximumRef = itk::NumericTraits<PixelType>::NonpositiveMin();
213 
214  for (i = 0; i < numberOfThreads; ++i)
215  {
216  count += m_Count[i];
217  diffCount += m_DiffCount[i];
218  squareOfDifferences += m_SquareOfDifferences[i];
219  absoluteValueOfDifferences += m_AbsoluteValueOfDifferences[i];
220 
221  if (m_ThreadMinRef[i] < minimumRef)
222  {
223  minimumRef = m_ThreadMinRef[i];
224  }
225  if (m_ThreadMaxRef[i] > maximumRef)
226  {
227  maximumRef = m_ThreadMaxRef[i];
228  }
229  }
230 
231  // compute mse
232  mse = (count != 0) ? squareOfDifferences / static_cast<RealType>(count) : 0.;
233  // compute mse
234  mae = (count != 0) ? absoluteValueOfDifferences / static_cast<RealType>(count) : 0.;
235 
236  // compute psnr
237  psnr = (std::abs(mse) > 0.0000000001 && (maximumRef - minimumRef) > 0.0000000001)
238  ? 10. * std::log10(((maximumRef - minimumRef) * (maximumRef - minimumRef)) / mse)
239  : 0.;
240  // Set the output
241  this->GetMSEOutput()->Set(mse);
242  this->GetMAEOutput()->Set(mae);
243  this->GetPSNROutput()->Set(psnr);
244  this->GetDiffCountOutput()->Set(static_cast<RealType>(diffCount));
245 }
246 
247 template <class TInputImage>
249 {
250  int numberOfThreads = this->GetNumberOfThreads();
251 
252  // Resize the thread temporaries
253  m_Count.SetSize(numberOfThreads);
254  m_DiffCount.SetSize(numberOfThreads);
255  m_SquareOfDifferences.SetSize(numberOfThreads);
256  m_AbsoluteValueOfDifferences.SetSize(numberOfThreads);
257 
258  m_ThreadMinRef.SetSize(numberOfThreads);
259  m_ThreadMaxRef.SetSize(numberOfThreads);
260 
261  // Initialize the temporaries
262  m_Count.Fill(itk::NumericTraits<long>::Zero);
263  m_DiffCount.Fill(itk::NumericTraits<long>::Zero);
264  m_SquareOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
265  m_AbsoluteValueOfDifferences.Fill(itk::NumericTraits<RealType>::Zero);
266  m_ThreadMinRef.Fill(itk::NumericTraits<PixelType>::max());
267  m_ThreadMaxRef.Fill(itk::NumericTraits<PixelType>::NonpositiveMin());
268 }
269 
270 template <class TInputImage>
272 {
274  Superclass::VerifyInputInformation();
275 }
276 
277 template <class TInputImage>
278 void PersistentCompareImageFilter<TInputImage>::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
279 {
283  InputImagePointer inputPtr1 = const_cast<TInputImage*>(this->GetInput(0));
284  InputImagePointer inputPtr2 = const_cast<TInputImage*>(this->GetInput(1));
286 
287  // support progress methods/callbacks
288  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
289 
290  RealType realValue1, realValue2;
291  PixelType value1, value2;
292 
293  itk::ImageRegionConstIterator<TInputImage> it1(inputPtr1, outputRegionForThread);
294  itk::ImageRegionConstIterator<TInputImage> it2(inputPtr2, outputRegionForThread);
295 
296  it1.GoToBegin();
297  it2.GoToBegin();
298  // do the work
299  while (!it1.IsAtEnd() && !it2.IsAtEnd())
300  {
301  value1 = it1.Get();
302  realValue1 = static_cast<RealType>(value1);
303 
304  value2 = it2.Get();
305  realValue2 = static_cast<RealType>(value2);
306 
307  if (value1 < m_ThreadMinRef[threadId])
308  {
309  m_ThreadMinRef[threadId] = value1;
310  }
311  if (value1 > m_ThreadMaxRef[threadId])
312  {
313  m_ThreadMaxRef[threadId] = value1;
314  }
315 
316  RealType diffVal = realValue1 - realValue2;
317  m_SquareOfDifferences[threadId] += diffVal * diffVal;
318  m_AbsoluteValueOfDifferences[threadId] += std::abs(diffVal);
319  if (!itk::Math::FloatAlmostEqual(realValue1, realValue2))
320  {
321  m_DiffCount[threadId]++;
322  }
323  m_Count[threadId]++;
324  ++it1;
325  ++it2;
326  progress.CompletedPixel();
327  }
328 }
329 template <class TImage>
330 void PersistentCompareImageFilter<TImage>::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  os << indent << "Count: " << this->GetDiffCount() << std::endl;
338 }
339 } // end namespace otb
340 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const override
itk::SimpleDataObjectDecorator< RealType > RealObjectType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
itk::NumericTraits< PixelType >::RealType RealType
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType