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