otbBoxAndWhiskerImageFilter.hxx

Go to the documentation of this file.

/*
 * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
 * Copyright (C) 2007-2012 Institut Mines Telecom / Telecom Bretagne
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
#ifndef otbBoxAndWhiskerImageFilter_hxx
#define otbBoxAndWhiskerImageFilter_hxx
 
#include <vector>
#include <algorithm>
 
#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkProgressReporter.h"
 
#include "otbBoxAndWhiskerImageFilter.h"
#include "otbEuclideanDistanceMetricWithMissingValue.h"
 
namespace otb
{
 
template <class TInputImage>
BoxAndWhiskerImageFilter<TInputImage>::BoxAndWhiskerImageFilter()
{
  this->SetNumberOfRequiredInputs(1);
  this->SetNumberOfRequiredOutputs(1);
  this->InPlaceOn();
  m_Beta = 1.5;
  m_Radius.Fill(1);
  m_NumberFound = 0L;
}
 
template <class TInputImage>
void BoxAndWhiskerImageFilter<TInputImage>::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
{
  const TInputImage* inputPtr  = this->GetInput();
  OutputImageType*   outputPtr = this->GetOutput();
 
  // data-set boundary faces
  typedef typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> ImageBoundaryFacesCalculatorType;
  typename ImageBoundaryFacesCalculatorType::FaceListType                                   faceList;
  ImageBoundaryFacesCalculatorType                                                          boundaryCalculator;
  faceList = boundaryCalculator(inputPtr, outputRegionForThread, this->GetRadius());
 
  // face iterator
  typename ImageBoundaryFacesCalculatorType::FaceListType::iterator faceIterator;
 
  // local iterators
  itk::ImageRegionConstIterator<InputImageType> inputIter;
  itk::ImageRegionIterator<OutputImageType>     outputIter;
 
  // support progress methods/callbacks
  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
 
  for (faceIterator = faceList.begin(); faceIterator != faceList.end(); ++faceIterator)
  {
    inputIter = itk::ImageRegionConstIterator<InputImageType>(inputPtr, *faceIterator);
    inputIter.GoToBegin();
 
    outputIter = itk::ImageRegionIterator<OutputImageType>(outputPtr, *faceIterator);
    outputIter.GoToBegin();
 
    while (!outputIter.IsAtEnd())
    {
      outputIter.Set(this->PerformBoxAndWhiskerDetection(inputIter.Get()));
 
      ++inputIter;
      ++outputIter;
 
      progress.CompletedPixel();
    }
  }
}
 
template <class TInputImage>
typename BoxAndWhiskerImageFilter<TInputImage>::PixelType BoxAndWhiskerImageFilter<TInputImage>::PerformBoxAndWhiskerDetection(const PixelType& pixel)
{
  typedef std::vector<ValueType>                                              VectorType;
  typedef otb::Statistics::EuclideanDistanceMetricWithMissingValue<PixelType> OutlierType;
 
  unsigned int       i;
  const unsigned int vectorSize            = pixel.Size();
  const unsigned int medianPosition        = vectorSize / 2;
  const unsigned int firstQuartilePosition = vectorSize / 4;
  const unsigned int thirdQuartilePosition = (3 * vectorSize) / 4;
 
  VectorType data(vectorSize);
  for (i    = 0; i < vectorSize; ++i)
    data[i] = pixel[i];
 
  std::sort(data.begin(), data.end());
 
  double box    = m_Beta * static_cast<double>(data[thirdQuartilePosition] - data[firstQuartilePosition]);
  double median = ::fabs(static_cast<double>(data[medianPosition]));
 
  PixelType outputPixel(pixel);
  for (i = 0; i < vectorSize; ++i)
  {
    if (::fabs(static_cast<double>(outputPixel[i]) - box) > median)
    {
      OutlierType::SetToMissingValue(outputPixel[i]);
      m_NumberFound++;
    }
  }
 
  return outputPixel;
}
 
template <class TInputImage>
void BoxAndWhiskerImageFilter<TInputImage>::GenerateOutputInformation()
{
  Superclass::GenerateOutputInformation();
  this->GetOutput()->SetNumberOfComponentsPerPixel(this->GetInput()->GetNumberOfComponentsPerPixel());
 
  this->GetOutput()->CopyInformation(this->GetInput());
}
 
template <class TInputImage>
void BoxAndWhiskerImageFilter<TInputImage>::AllocateOutputs()
{
  OutputImageType* output = this->GetOutput();
  output->SetNumberOfComponentsPerPixel(this->GetInput()->GetNumberOfComponentsPerPixel());
  output->SetBufferedRegion(output->GetRequestedRegion());
  output->Allocate();
}
 
} // end of namespace otb
 
#endif