Commit a9db248c authored by Gregoire Mercier's avatar Gregoire Mercier

ADD: Box and Whisker filters

parent e4c5deae
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Copyright (c) Institut Telecom ; Telecom bretagne. All rights reserved.
See ITCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbBoxAndWhiskerImageFilter__h
#define __otbBoxAndWhiskerImageFilter__h
#include "itkMacro.h"
#include "itkUnaryFunctorImageFilter.h"
namespace otb {
/** \class BoxAndWhiskerImageFilter
* \brief This class performs the detection of outlier with the Box and Whisker technique
*
* This is appropriated to perform simple outlier detection over vector data.
* The input image has to contain vector pixel through a VectorImage type.
* When an outlier is detected, it is treated to as an missing value through the function
* \code{otb::EuclideanDistanceWithMissingValue::SetToMissingValue()}.
*
* A component is considered to as a missing value when
* \f$ \| x_i - \beta \left( x_{3/4} - x_{1/4} \right) \| > \| x_{1/2} \| \f$
* where \f$ x_{1/4}, x_{1/2}, x_{3/4} \f$ stand for the first, second (median) and
* thrid quantile values.
*
* \ingroup Streamed
* \sa EuclideanDistanceWithMissingValue
*/
template < class TInputImage >
class ITK_EXPORT BoxAndWhiskerImageFilter
: public itk::InPlaceImageFilter< TInputImage >
{
public :
/** Standard class typedefs. */
typedef BoxAndWhiskerImageFilter Self;
typedef typename itk::InPlaceImageFilter< TInputImage > Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(BoxAndWhiskerImageFilter, InPlaceImageFilter);
/** Template parameters typedefs */
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImagePointer OutputImagePointer;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
typedef typename Superclass::OutputImagePixelType OutputImagePixelType;
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::InputImagePointer InputImagePointer;
typedef typename Superclass::InputImageConstPointer InputImageConstPointer;
typedef typename Superclass::InputImageRegionType InputImageRegionType;
typedef typename Superclass::InputImagePixelType InputImagePixelType;
typedef typename InputImageType::PixelType PixelType;
typedef typename PixelType::ValueType ValueType;
typedef typename InputImageType::SizeType SizeType;
typedef typename InputImageType::RegionType RegionType;
/** Dimension */
itkStaticConstMacro(InputImageDimension, unsigned int, InputImageType::ImageDimension);
itkStaticConstMacro(OutputImageDimension, unsigned int, OutputImageType::ImageDimension);
/** Get the Radius (all to 1) */
itkGetConstReferenceMacro(Radius, SizeType);
/** Fix the \em whisker */
itkGetConstMacro(Beta,double);
itkSetMacro(Beta,double);
itkGetConstMacro(NumberFound,unsigned int);
protected:
BoxAndWhiskerImageFilter ();
virtual ~BoxAndWhiskerImageFilter () {}
/** Main computation method implemented as a multithreaded filter */
virtual void ThreadedGenerateData ( const OutputImageRegionType& outputRegionForThread, int threadId );
virtual void GenerateOutputInformation();
virtual void AllocateOutputs();
/** Perform the outlier detection */
PixelType PerformBoxAndWhiskerDetection ( const PixelType & pixel );
private:
BoxAndWhiskerImageFilter ( const Self & );
void operator= ( const Self & ); // not implemented
SizeType m_Radius;
double m_Beta;
long int m_NumberFound;
}; // end of class BoxAndWhiskerImageFilter
} // end of namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbBoxAndWhiskerImageFilter.txx"
#endif
#endif
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Copyright (c) Institut Telecom ; Telecom bretagne. All rights reserved.
See ITCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbBoxAndWhiskerImageFilter__txx
#define __otbBoxAndWhiskerImageFilter__txx
#include <vector>
#include <algorithm>
#include "itkConstNeighborhoodIterator.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkProgressReporter.h"
#include "otbBoxAndWhiskerImageFilter.h"
#include "otbEuclideanDistanceWithMissingValue.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,
int 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() );
/**
* Process each of the boundary faces.
* These are N-d regions which border the edge of the buffer.
*/
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::EuclideanDistanceWithMissingValue< 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() );
this->GetOutput()->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
}
template < class TInputImage >
void
BoxAndWhiskerImageFilter < TInputImage >
::AllocateOutputs ()
{
OutputImageType * output = this->GetOutput();
output->SetNumberOfComponentsPerPixel(
this->GetInput()->GetNumberOfComponentsPerPixel() );
output->Allocate();
}
} // end of namespace otb
#endif
......@@ -9,11 +9,11 @@
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
See ITCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
......@@ -21,84 +21,76 @@
#ifndef __otbEuclideanDistanceWithMissingValue_h
#define __otbEuclideanDistanceWithMissingValue_h
#include "itkEuclideanDistance.h"
#include "otbEuclideanDistanceWithMissingValuePow2.h"
namespace otb
{
namespace otb {
namespace Statistics
{
namespace Statistics {
/** \class EuclideanDistanceWithMissingValue
* \brief Euclidean distance function facing missing value.
*
* This class is derived from EuclideanDistance class. In addition
* to the initial Evaluate method, the class does not perform calculation
* when a component does not contain any data.
*
* This class is derived from EuclideanDistanceWithMissingValuePow2 class that handle the missing value
* functonnalities. Here, the square root is included in the evaluation...
*
* The class can be templated over any container that holds data elements, as
* for template of EuclideanDistance.
* for template of EuclideanDistance.
*
* The only restriction is that elemnts have to support \code NaN \endcode.
* The only restriction is that elemnts have to support '\code{NaN}'.
*
* \sa EuclideanDistance
* \sa EuclideanDistanceWithMissingValuePow2
*/
template< class TVector >
class ITK_EXPORT EuclideanDistanceWithMissingValue :
public itk::Statistics::EuclideanDistance< TVector >
public otb::Statistics::EuclideanDistanceWithMissingValuePow2< TVector >
{
public:
/** Standard "Self" typedef. */
typedef EuclideanDistanceWithMissingValue Self;
typedef itk::Statistics::EuclideanDistance< TVector > Superclass;
typedef itk::SmartPointer< Self > Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef typename Superclass::MeasurementVectorSizeType MeasurementVectorSizeType;
public:
/** Standard "Self" typedef. */
typedef EuclideanDistanceWithMissingValue Self;
typedef otb::Statistics::EuclideanDistanceWithMissingValuePow2< TVector >
Superclass;
typedef itk::SmartPointer< Self > Pointer ;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Run-time type information (and related methods). */
itkTypeMacro(EuclideanDistanceWithMissingValue, EuclideanDistance);
typedef typename Superclass::MeasurementVectorSizeType
MeasurementVectorSizeType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(EuclideanDistanceWithMissingValue, EuclideanDistance);
/** Type of the component of a vector */
typedef typename TVector::ValueType ValueType;
/** Method for creation through the object factory. */
itkNewMacro(Self) ;
/** Gets the distance between the origin and x */
double Evaluate(const TVector &x) const;
/** Type of the component of a vector */
typedef typename TVector::ValueType ValueType ;
/** Gets the distance between x1 and x2 */
double Evaluate(const TVector &x1, const TVector &x2) const;
/** Gets the distance between the origin and x */
double Evaluate(const TVector &x) const {
return ::vcl_sqrt( Superclass::Evaluate( x ) ); }
/** Gets the cooridnate distance between a and b. NOTE: a and b
* should be type of component */
double Evaluate(const ValueType &a, const ValueType &b) const;
/** Gets the distance between x1 and x2 */
double Evaluate(const TVector &x1, const TVector &x2) const {
return ::vcl_sqrt( Superclass::Evaluate( x1, x2 ) ); }
/** Returns true if the distance between x and the origin is less
* than radius */
bool IsWithinRange(const TVector &x, const double radius) const
{
return Superclass::IsWithinRange( x, radius );
}
/** Gets the cooridnate distance between a and b. NOTE: a and b
* should be type of component */
double Evaluate(const ValueType &a, const ValueType &b) const {
return ::vcl_sqrt( Superclass::Evaluate( a, b ) ); }
/** Check if a value is NaN or not */
static bool IsMissingValue ( const ValueType & v);
/** Returns true if the distance between x and the origin is less
* than radius */
bool IsWithinRange(const TVector &x, const double radius) const {
return Superclass::IsWithinRange( x, radius ); }
/** Set a value to Nan */
static void SetToMissingValue ( ValueType & v );
protected:
EuclideanDistanceWithMissingValue() {}
virtual ~EuclideanDistanceWithMissingValue() {}
}; // end of class
protected:
EuclideanDistanceWithMissingValue() {}
virtual ~EuclideanDistanceWithMissingValue() {}
} ; // end of class
} // end namespace statistics
} // end namespace otb
#ifndef ITK_MANUAL_INSTANTIATION
#include "otbEuclideanDistanceWithMissingValue.txx"
#endif
#endif
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
See ITCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbEuclideanDistanceWithMissingValuePow2_h
#define __otbEuclideanDistanceWithMissingValuePow2_h
#include <itkEuclideanDistance.h>
namespace otb {
namespace Statistics {
/** \class EuclideanDistanceWithMissingValuePow2
* \brief Euclidean comparison distance function facing missing value. The square root is not performed in this class.
*
* This class is derived from EuclideanDistance class. In addition
* to the initial Evaluate method, the class does not perform calculation
* when a component does not contain any data. The 'no data' case is perform through the \code{SetToMissingValue()}
* which performs a \code{Nan} affectation.
*
* The class can be templated over any container that holds data elements, as
* for template of EuclideanDistance.
*
* The only restriction is that elemnts have to support '\code{NaN}'.
*
* \sa EuclideanDistance
* \sa EuclideanDistanceWithMissingValue
*/
template< class TVector >
class ITK_EXPORT EuclideanDistanceWithMissingValuePow2 :
public itk::Statistics::EuclideanDistance< TVector >
{
public:
/** Standard "Self" typedef. */
typedef EuclideanDistanceWithMissingValuePow2 Self;
typedef itk::Statistics::EuclideanDistance< TVector > Superclass;
typedef itk::SmartPointer< Self > Pointer ;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef typename Superclass::MeasurementVectorSizeType MeasurementVectorSizeType;
/** Run-time type information (and related methods). */
itkTypeMacro(EuclideanDistanceWithMissingValuePow2, EuclideanDistance);
/** Method for creation through the object factory. */
itkNewMacro(Self) ;
/** Type of the component of a vector */
typedef typename TVector::ValueType ValueType ;
/** Gets the distance between the origin and x */
double Evaluate(const TVector &x) const ;
/** Gets the distance between x1 and x2 */
double Evaluate(const TVector &x1, const TVector &x2) const ;
/** Gets the cooridnate distance between a and b. NOTE: a and b
* should be type of component */
double Evaluate(const ValueType &a, const ValueType &b) const ;
/** Returns true if the distance between x and the origin is less
* than radius */
bool IsWithinRange(const TVector &x, const double radius) const {
return Superclass::IsWithinRange( x, radius ); }
/** Check if a value is NaN or not */
static bool IsMissingValue ( const ValueType & v) ;
/** Set a value to Nan */
static void SetToMissingValue ( ValueType & v );
protected:
EuclideanDistanceWithMissingValuePow2() {}
virtual ~EuclideanDistanceWithMissingValuePow2() {}
} ; // end of class
} // end namespace statistics
} // end namespace otb
#ifndef ITK_MANUAL_INSTANTIATION
#include "otbEuclideanDistanceWithMissingValuePow2.txx"
#endif
#endif
......@@ -9,107 +9,134 @@
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
See ITCopyright.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbEuclideanDistanceWithMissingValue_txx
#define __otbEuclideanDistanceWithMissingValue_txx
#ifndef __otbEuclideanDistanceWithMissingValuePow2_txx
#define __otbEuclideanDistanceWithMissingValuePow2_txx
#include "itkNumericTraits.h"
#include "vnl/vnl_math.h"
namespace otb
{
namespace otb {
namespace Statistics
{
namespace Statistics {
template< class TVector >
inline double
EuclideanDistanceWithMissingValue< TVector >
EuclideanDistanceWithMissingValuePow2< TVector >
::Evaluate(const TVector &x1, const TVector &x2) const
{
if ( itk::MeasurementVectorTraits::GetLength( x1 ) !=
itk::MeasurementVectorTraits::GetLength( x2 ) )
if( itk::MeasurementVectorTraits::GetLength( x1 ) !=
itk::MeasurementVectorTraits::GetLength( x2 ) )
{
itkExceptionMacro( << "Vector lengths must be equal." );
}
double temp, distance = itk::NumericTraits< double >::Zero;
double temp, distance = itk::NumericTraits< double >::Zero ;
for (unsigned int i = 0; i < x1.Size(); ++i )
for(unsigned int i = 0 ; i < x1.Size(); i++ )
{
if ( !IsMissingValue( x1[i] ) && !IsMissingValue( x2[i] ) )
{
temp = x1[i] - x2[i];
distance += temp * temp;
temp = x1[i] - x2[i] ;
distance += temp * temp ;
}
}
return ::vcl_sqrt(distance);
return distance ;
}
template< class TVector >
inline double
EuclideanDistanceWithMissingValue< TVector >
EuclideanDistanceWithMissingValuePow2< TVector >
::Evaluate(const TVector &x) const
{
MeasurementVectorSizeType
measurementVectorSize = this->GetMeasurementVectorSize();
if (measurementVectorSize == 0)
MeasurementVectorSizeType
measurementVectorSize = this->GetMeasurementVectorSize();
if(measurementVectorSize == 0)
{
itkExceptionMacro( << "Please set the MeasurementVectorSize first" );
}
itk::MeasurementVectorTraits::Assert( this->m_Origin, measurementVectorSize,
"EuclideanDistance::Evaluate Origin and input vector have different lengths");
itk::MeasurementVectorTraits::Assert( this->m_Origin, measurementVectorSize,
"EuclideanDistance::Evaluate Origin and input vector have different lengths");
double temp, distance = itk::NumericTraits< double >::Zero;
double temp, distance = itk::NumericTraits< double >::Zero ;
for (unsigned int i = 0; i < measurementVectorSize; ++i )
for(unsigned int i = 0 ; i < measurementVectorSize ; i++ )
{
if ( !IsMissingValue( this->GetOrigin()[i] ) && !IsMissingValue( x[i] ) )
{
temp = this->GetOrigin()[i] - x[i];
distance += temp * temp;
temp = this->GetOrigin()[i] - x[i] ;
distance += temp * temp ;
}
}
return ::vcl_sqrt(distance);
return distance ;