OTB  9.0.0
Orfeo Toolbox
otbDisparityMapMedianFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 otbDisparityMapMedianFilter_hxx
22 #define otbDisparityMapMedianFilter_hxx
23 
24 #ifdef ITK_USE_CONSOLIDATED_MORPHOLOGY
25 #else
26 
27 #include "otbDisparityMapMedianFilter.h"
28 #include "itkConstNeighborhoodIterator.h"
29 #include "itkImageRegionIterator.h"
30 #include "itkImageRegionIteratorWithIndex.h"
31 #include "itkNeighborhoodAlgorithm.h"
32 #include "itkProgressReporter.h"
33 
34 namespace otb
35 {
36 
37 template <class TInputImage, class TOutputImage, class TMask>
38 DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::DisparityMapMedianFilter()
39 {
40  this->SetNumberOfRequiredInputs(1);
41  this->SetNumberOfRequiredOutputs(4);
42  this->SetNthOutput(1, TMask::New());
43  this->SetNthOutput(2, TOutputImage::New());
44  this->SetNthOutput(3, TMask::New());
45  m_Radius.Fill(3);
46  m_IncoherenceThreshold = 1.0;
47 }
48 
49 template <class TInputImage, class TOutputImage, class TMask>
50 void DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::SetMaskInput(const TMask* inputmask)
51 {
52  // Process object is not const-correct so the const casting is required.
53  this->SetNthInput(1, const_cast<TMask*>(inputmask));
54 }
55 
56 
57 template <class TInputImage, class TOutputImage, class TMask>
58 const TMask* DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GetMaskInput()
59 {
60  if (this->GetNumberOfInputs() < 2)
61  {
62  return nullptr;
63  }
64  return static_cast<const TMask*>(this->itk::ProcessObject::GetInput(1));
65 }
66 
67 template <class TInputImage, class TOutputImage, class TMask>
68 TMask* DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GetOutputMask()
69 {
70  if (this->GetNumberOfOutputs() < 2)
71  {
72  return nullptr;
73  }
74  return static_cast<TMask*>(this->itk::ProcessObject::GetOutput(1));
75 }
76 
77 
78 template <class TInputImage, class TOutputImage, class TMask>
79 TOutputImage* DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GetOutputDisparityMap()
80 {
81  if (this->GetNumberOfOutputs() < 3)
82  {
83  return nullptr;
84  }
85  return static_cast<TOutputImage*>(this->itk::ProcessObject::GetOutput(2));
86 }
87 
88 template <class TInputImage, class TOutputImage, class TMask>
89 TMask* DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GetOutputDisparityMask()
90 {
91  if (this->GetNumberOfOutputs() < 4)
92  {
93  return nullptr;
94  }
95  return static_cast<TMask*>(this->itk::ProcessObject::GetOutput(3));
96 }
97 
98 
99 template <class TInputImage, class TOutputImage, class TMask>
100 void DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GenerateOutputInformation()
101 {
102  // Call superclass implementation
103  Superclass::GenerateOutputInformation();
104 
105  // Retrieve output pointers
106  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput());
107  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
108  TMask* outputmaskPtr = this->GetOutputMask();
109  typename Superclass::OutputImagePointer outputdisparitymapPtr = this->GetOutputDisparityMap();
110  TMask* outputdisparitymaskPtr = this->GetOutputDisparityMask();
111 
112  // Update size and spacing according to grid step
113  InputImageRegionType largestRegion = outputPtr->GetLargestPossibleRegion();
114  SizeType outputSize = largestRegion.GetSize();
115 
116  // Set largest region size
117  largestRegion.SetSize(outputSize);
118  outputPtr->SetLargestPossibleRegion(largestRegion);
119  outputmaskPtr->SetLargestPossibleRegion(largestRegion);
120  outputdisparitymapPtr->SetLargestPossibleRegion(largestRegion);
121  outputdisparitymaskPtr->SetLargestPossibleRegion(largestRegion);
122 }
123 
124 
125 template <class TInputImage, class TOutputImage, class TMask>
126 void DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GenerateInputRequestedRegion()
127 {
128  // call the superclass' implementation of this method
129  Superclass::GenerateInputRequestedRegion();
130 
131  // get pointers to the input and output
132  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput());
133  TMask* inputmaskPtr = const_cast<TMask*>(this->GetMaskInput());
134  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
135  TMask* outputmaskPtr = this->GetOutputMask();
136  typename Superclass::OutputImagePointer outputdisparitymapPtr = this->GetOutputDisparityMap();
137  TMask* outputdisparitymaskPtr = this->GetOutputDisparityMask();
138 
139  if (!inputPtr || !outputPtr || !outputmaskPtr || !outputdisparitymapPtr || !outputdisparitymaskPtr)
140  {
141  return;
142  }
143 
144  if (inputmaskPtr)
145  {
146  // check that the mask has the same size as the input image
147  if (inputmaskPtr->GetLargestPossibleRegion() != inputPtr->GetLargestPossibleRegion())
148  {
149  itkExceptionMacro(<< "Input image and mask image don't have the same size ! Input image :" << inputPtr->GetLargestPossibleRegion()
150  << "; Mask image :" << inputmaskPtr->GetLargestPossibleRegion());
151  }
152  }
153 
154  // get a copy of the input requested region (should equal the output
155  // requested region)
156  typename TInputImage::RegionType inputRequestedRegion;
157  inputRequestedRegion = inputPtr->GetRequestedRegion();
158 
159  // pad the input requested region by the operator radius
160  inputRequestedRegion.PadByRadius(m_Radius);
161 
162  // crop the input requested region at the input's largest possible region
163  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
164  {
165  inputPtr->SetRequestedRegion(inputRequestedRegion);
166  if (inputmaskPtr)
167  {
168  inputmaskPtr->SetRequestedRegion(inputRequestedRegion);
169  }
170 
171  return;
172  }
173  else
174  {
175  // Couldn't crop the region (requested region is outside the largest
176  // possible region). Throw an exception.
177 
178  // store what we tried to request (prior to trying to crop)
179  inputPtr->SetRequestedRegion(inputRequestedRegion);
180 
181  // build an exception
182  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
183  e.SetLocation(ITK_LOCATION);
184  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
185  e.SetDataObject(inputPtr);
186  throw e;
187  }
188 }
189 
190 
191 template <class TInputImage, class TOutputImage, class TMask>
192 void DisparityMapMedianFilter<TInputImage, TOutputImage, TMask>::GenerateData()
193 {
194  // Allocate outputs
195  this->AllocateOutputs();
196 
197  // Get the image pointers
198  typename OutputImageType::Pointer output = this->GetOutput();
199  typename InputImageType::ConstPointer input = this->GetInput();
200  typename MaskImageType::ConstPointer inputmaskPtr = this->GetMaskInput();
201  TMask* outputmaskPtr = this->GetOutputMask();
202  TOutputImage* outputdisparitymapPtr = this->GetOutputDisparityMap();
203  TMask* outputdisparitymaskPtr = this->GetOutputDisparityMask();
204 
205  SizeType imgSize = output->GetLargestPossibleRegion().GetSize();
206 
208  itk::ConstNeighborhoodIterator<InputImageType> InputIt(m_Radius, input, input->GetRequestedRegion());
209  itk::ConstNeighborhoodIterator<TMask> MaskInputIt;
210  if (inputmaskPtr)
211  {
212  MaskInputIt.Initialize(m_Radius, inputmaskPtr, inputmaskPtr->GetRequestedRegion());
213  }
215 
217  itk::ImageRegionIteratorWithIndex<OutputImageType> outputIt(output, output->GetRequestedRegion());
218  itk::ImageRegionIterator<TMask> outputMaskIt(outputmaskPtr, output->GetRequestedRegion());
219  itk::ImageRegionIterator<OutputImageType> outputDisparityMapIt(outputdisparitymapPtr, output->GetRequestedRegion());
220  itk::ImageRegionIterator<TMask> outputDisparityMaskIt(outputdisparitymaskPtr, output->GetRequestedRegion());
221  outputIt.GoToBegin();
222  outputMaskIt.GoToBegin();
223  outputDisparityMapIt.GoToBegin();
224  outputDisparityMaskIt.GoToBegin();
226 
227  std::vector<InputPixelType> pixels;
228  while (!outputIt.IsAtEnd() && !outputMaskIt.IsAtEnd() && !outputDisparityMapIt.IsAtEnd() && !outputDisparityMapIt.IsAtEnd())
229  {
230  if (outputIt.GetIndex()[0] >= static_cast<IndexValueType>(m_Radius[0]) &&
231  outputIt.GetIndex()[0] < static_cast<IndexValueType>(imgSize[0]) - static_cast<IndexValueType>(m_Radius[0]) &&
232  outputIt.GetIndex()[1] >= static_cast<IndexValueType>(m_Radius[1]) &&
233  outputIt.GetIndex()[1] < static_cast<IndexValueType>(imgSize[1]) - static_cast<IndexValueType>(m_Radius[1]))
234  {
235  // determine pixels in the neighborhood window whose subpixel mask is not equal to 0
236  int p = 0;
237  pixels.clear();
238  if (inputmaskPtr)
239  {
240  MaskInputIt.SetLocation(outputIt.GetIndex());
241  }
242  InputIt.SetLocation(outputIt.GetIndex());
243  for (unsigned int i = 0; i < InputIt.Size(); i++)
244  {
245  if (!inputmaskPtr || (MaskInputIt.GetPixel(i) != 0))
246  {
247  p++;
248  pixels.push_back(InputIt.GetPixel(i));
249  }
250  }
251  if (p > 0)
252  {
253  // outputMaskIt.Set(itk::NumericTraits<MaskImagePixelType>::max());
254  outputMaskIt.Set(1);
255  // get the median value
256  if ((p & 0x1) == 0)
257  {
258  const unsigned int medianPosition_low = p / 2 - 1;
259  const unsigned int medianPosition_high = p / 2;
260  const typename std::vector<InputPixelType>::iterator medianIterator_low = pixels.begin() + medianPosition_low;
261  const typename std::vector<InputPixelType>::iterator medianIterator_high = pixels.begin() + medianPosition_high;
262  std::nth_element(pixels.begin(), medianIterator_low, pixels.end());
263  std::nth_element(pixels.begin(), medianIterator_high, pixels.end());
264  outputIt.Set(static_cast<typename OutputImageType::PixelType>((*medianIterator_low + *medianIterator_high) / 2));
265  }
266  else
267  {
268  const unsigned int medianPosition = p / 2;
269  const typename std::vector<InputPixelType>::iterator medianIterator = pixels.begin() + medianPosition;
270  std::nth_element(pixels.begin(), medianIterator, pixels.end());
271  outputIt.Set(static_cast<typename OutputImageType::PixelType>(*medianIterator));
272  }
273  }
274  else
275  {
276  outputIt.Set(0.0);
277  outputMaskIt.Set(0);
278  }
279  }
280  else
281  {
282  outputIt.Set(0.0);
283  outputMaskIt.Set(0);
284  }
285 
286  outputDisparityMapIt.Set(static_cast<typename OutputImageType::PixelType>(InputIt.GetCenterPixel())); // copy the input disparity map
287 
288  if (inputmaskPtr)
289  {
290  outputDisparityMaskIt.Set(MaskInputIt.GetCenterPixel()); // copy the input disparity mask
291  }
292  else
293  {
294  outputDisparityMaskIt.Set(1); // if no mask is given, fill with 1 by default
295  }
296 
297  ++outputIt;
298  ++outputMaskIt;
299  ++outputDisparityMapIt;
300  ++outputDisparityMaskIt;
301  }
302 
303  // Remove incoherences between disparity and median//
304  // creation of the auxiliary image that store positions of incoherences between the median and the input disparity map
305  MaskImagePointerType image_aux = MaskImageType::New();
306  image_aux->SetRegions(input->GetRequestedRegion());
307  image_aux->Allocate();
308  image_aux->FillBuffer(0);
309  itk::NeighborhoodIterator<TMask> image_aux_It(m_Radius, image_aux, input->GetRequestedRegion());
310  outputIt.GoToBegin();
311  outputMaskIt.GoToBegin();
312 
313  itk::ImageRegionConstIterator<OutputImageType> MedianIt(output, output->GetRequestedRegion());
314 
315  while (!outputIt.IsAtEnd() && !outputMaskIt.IsAtEnd())
316  {
317  if (outputIt.GetIndex()[0] >= static_cast<IndexValueType>(m_Radius[0]) &&
318  outputIt.GetIndex()[0] < static_cast<IndexValueType>(imgSize[0]) - static_cast<IndexValueType>(m_Radius[0]) &&
319  outputIt.GetIndex()[1] >= static_cast<IndexValueType>(m_Radius[1]) &&
320  outputIt.GetIndex()[1] < static_cast<IndexValueType>(imgSize[1]) - static_cast<IndexValueType>(m_Radius[1]))
321  {
322  InputIt.SetLocation(outputIt.GetIndex());
323  MedianIt.SetIndex(outputIt.GetIndex());
324  outputDisparityMapIt.SetIndex(outputIt.GetIndex());
325  outputDisparityMaskIt.SetIndex(outputIt.GetIndex());
326  image_aux_It.SetLocation(outputIt.GetIndex());
327  if (inputmaskPtr)
328  {
329  MaskInputIt.SetLocation(outputIt.GetIndex());
330  }
331 
332  if ((!inputmaskPtr || (MaskInputIt.GetCenterPixel() != 0)) && std::fabs(InputIt.GetCenterPixel() - MedianIt.Get()) > m_IncoherenceThreshold)
333  {
334  outputDisparityMapIt.Set(0.0); // Remove pixel from disparity map//
335  outputDisparityMaskIt.Set(0);
336  for (unsigned int i = 0; i < image_aux_It.Size(); i++)
337  {
338  image_aux_It.SetPixel(i, 1);
339  }
340  }
341  }
342  ++outputIt;
343  ++outputMaskIt;
344  }
345 
346  // Recompute median where values had been changed
347  // we use the updated sub pixel disparity map
348  itk::ConstNeighborhoodIterator<OutputImageType> updatedDisparityMapIt(m_Radius, outputdisparitymapPtr, output->GetRequestedRegion());
349  itk::ConstNeighborhoodIterator<TMask> updatedDisparityMaskIt(m_Radius, outputdisparitymaskPtr, output->GetRequestedRegion());
350 
351  outputIt.GoToBegin();
352  outputMaskIt.GoToBegin();
353  updatedDisparityMapIt.GoToBegin();
354  updatedDisparityMaskIt.GoToBegin();
355 
356  while (!updatedDisparityMapIt.IsAtEnd() && !updatedDisparityMaskIt.IsAtEnd() && !outputIt.IsAtEnd() && !outputMaskIt.IsAtEnd())
357  {
358  if (outputIt.GetIndex()[0] >= static_cast<IndexValueType>(m_Radius[0]) &&
359  outputIt.GetIndex()[0] < static_cast<IndexValueType>(imgSize[0]) - static_cast<IndexValueType>(m_Radius[0]) &&
360  outputIt.GetIndex()[1] >= static_cast<IndexValueType>(m_Radius[1]) &&
361  outputIt.GetIndex()[1] < static_cast<IndexValueType>(imgSize[1]) - static_cast<IndexValueType>(m_Radius[1]))
362  {
363  image_aux_It.SetLocation(outputIt.GetIndex());
364  if (image_aux_It.GetCenterPixel() != 0)
365  {
366  // determine pixels in the neighborhood window whose subpixel mask is not equal to 0
367  int p = 0;
368  pixels.clear();
369  for (unsigned int i = 0; i < updatedDisparityMaskIt.Size(); i++)
370  {
371  if (updatedDisparityMaskIt.GetPixel(i) != 0)
372  {
373  p++;
374  pixels.push_back(updatedDisparityMapIt.GetPixel(i));
375  }
376  }
377  if (p > 0)
378  {
379  // outputMaskIt.Set(itk::NumericTraits<MaskImagePixelType>::max());
380  outputMaskIt.Set(1);
381  // get the median value
382  if ((p & 0x1) == 0)
383  {
384  const unsigned int medianPosition_low = p / 2 - 1;
385  const unsigned int medianPosition_high = p / 2;
386  const typename std::vector<InputPixelType>::iterator medianIterator_low = pixels.begin() + medianPosition_low;
387  const typename std::vector<InputPixelType>::iterator medianIterator_high = pixels.begin() + medianPosition_high;
388  std::nth_element(pixels.begin(), medianIterator_low, pixels.end());
389  std::nth_element(pixels.begin(), medianIterator_high, pixels.end());
390  outputIt.Set(static_cast<typename OutputImageType::PixelType>((*medianIterator_low + *medianIterator_high) / 2));
391  }
392  else
393  {
394  const unsigned int medianPosition = p / 2;
395  const typename std::vector<InputPixelType>::iterator medianIterator = pixels.begin() + medianPosition;
396  std::nth_element(pixels.begin(), medianIterator, pixels.end());
397  outputIt.Set(static_cast<typename OutputImageType::PixelType>(*medianIterator));
398  }
399  }
400  else
401  {
402  outputIt.Set(0.0);
403  outputMaskIt.Set(0);
404  }
405  }
406  }
407  ++outputIt;
408  ++outputMaskIt;
409  ++updatedDisparityMapIt;
410  ++updatedDisparityMaskIt;
411  }
412 }
413 
414 
418 template <class TInputImage, class TOutput, class TMask>
419 void DisparityMapMedianFilter<TInputImage, TOutput, TMask>::PrintSelf(std::ostream& os, itk::Indent indent) const
420 {
421  Superclass::PrintSelf(os, indent);
422  os << indent << "Radius: " << m_Radius << std::endl;
423 }
425 
426 } // end namespace otb
427 
428 #endif
429 
430 #endif
otb::DisparityMapMedianFilter::SizeType
InputImageType::SizeType SizeType
Definition: otbDisparityMapMedianFilter.h:108
otb::DisparityMapMedianFilter::GetOutputDisparityMap
TOutputImage * GetOutputDisparityMap()
Definition: otbDisparityMapMedianFilter.hxx:79
otb::DisparityMapMedianFilter::GenerateInputRequestedRegion
void GenerateInputRequestedRegion() override
Definition: otbDisparityMapMedianFilter.hxx:126
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::DisparityMapMedianFilter::InputImageRegionType
InputImageType::RegionType InputImageRegionType
Definition: otbDisparityMapMedianFilter.h:106
otb::DisparityMapMedianFilter::IndexValueType
OutputImageType::IndexValueType IndexValueType
Definition: otbDisparityMapMedianFilter.h:109
otb::DisparityMapMedianFilter::GetOutputDisparityMask
TMask * GetOutputDisparityMask()
Definition: otbDisparityMapMedianFilter.hxx:89
otb::DisparityMapMedianFilter::GenerateData
void GenerateData() override
Definition: otbDisparityMapMedianFilter.hxx:192
otb::DisparityMapMedianFilter::GetOutputMask
TMask * GetOutputMask()
Definition: otbDisparityMapMedianFilter.hxx:68
otb::DisparityMapMedianFilter::GetMaskInput
const TMask * GetMaskInput()
Definition: otbDisparityMapMedianFilter.hxx:58
otbDisparityMapMedianFilter.h
otb::DisparityMapMedianFilter::SetMaskInput
void SetMaskInput(const TMask *inputmask)
Definition: otbDisparityMapMedianFilter.hxx:50
otb::DisparityMapMedianFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbDisparityMapMedianFilter.hxx:419
otb::DisparityMapMedianFilter::GenerateOutputInformation
void GenerateOutputInformation(void) override
Definition: otbDisparityMapMedianFilter.hxx:100
otb::DisparityMapMedianFilter::DisparityMapMedianFilter
DisparityMapMedianFilter()
Definition: otbDisparityMapMedianFilter.hxx:38
otb::DisparityMapMedianFilter::MaskImagePointerType
MaskImageType::Pointer MaskImagePointerType
Definition: otbDisparityMapMedianFilter.h:104
Generated at Wed Feb 28 2024 15:21:49 for OTB with   doxygen 1.8.17