Commit a52635d8 authored by Cédric Traizet's avatar Cédric Traizet

Merge branch 'local_rx_detector' into 'develop'

Local Rx detector

See merge request !353
parents 09e69e47 5ff7ff95
......@@ -28,7 +28,13 @@ otb_create_application(
SOURCES otbVertexComponentAnalysis.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
otb_create_application(
NAME LocalRxDetection
SOURCES otbLocalRxDetection.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
otb_create_application(
NAME EndmemberNumberEstimation
SOURCES otbEndmemberNumberEstimation.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* 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.
*/
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbLocalRxDetectorFilter.h"
#include "otbFunctorImageFilter.h"
namespace otb
{
namespace Wrapper
{
class LocalRxDetection : public Application
{
public:
/** Standard class typedefs. */
typedef LocalRxDetection Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(LocalRxDetection, otb::Application);
/** Image typedefs */
typedef DoubleVectorImageType VectorImageType;
typedef DoubleImageType ImageType;
/** Filter typedefs */
typedef otb::LocalRxDetectorFilter<VectorImageType, ImageType> LocalRxDetectorFilterType; //TODO remove this
private:
void DoInit() override
{
SetName("LocalRxDetection");
SetDescription("Performs local Rx score computation on an hyperspectral image.");
// Documentation
SetDocName("Local Rx Detection");
SetDocLongDescription("Performs local Rx score computation on an input "
"hyperspectral image. For each hyperspectral pixel, the Rx score is "
"computed using statistics computed on a dual neighborhood. The dual "
"neighborhood is composed of all pixel that are in between two radiuses "
"around the center pixel. This score can then be used to detect "
"anomalies in the image, this can be done for example by thresholding "
"the result of this application with the BandMath application.");
SetDocLimitations("None");
SetDocAuthors("OTB-Team");
SetDocSeeAlso("BandMath");
AddDocTag(Tags::Hyperspectral);
AddParameter(ParameterType_InputImage, "in", "Input Image");
SetParameterDescription("in","Input hyperspectral data cube");
AddParameter(ParameterType_OutputImage, "out", "Output Image");
SetParameterDescription("out","Output Rx score image");
MandatoryOff("out");
AddParameter(ParameterType_Int, "ir", "Internal radius");
SetParameterDescription("ir", "Internal radius in pixel");
SetDefaultParameterInt("ir", 1);
AddParameter(ParameterType_Int, "er", "External radius");
SetParameterDescription("er","External radius in pixel");
SetDefaultParameterInt("er", 5);
AddRAMParameter();
// Doc example parameter settings
SetDocExampleParameterValue("in", "cupriteSubHsi.tif");
SetDocExampleParameterValue("out", "LocalRxScore.tif");
SetDocExampleParameterValue("ir", "1");
SetDocExampleParameterValue("er", "5");
SetOfficialDocLink();
}
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
}
void DoExecute() override
{
auto inputImage = GetParameterDoubleVectorImage("in");
inputImage->UpdateOutputInformation();
// The localRxDetectionFilter can be replaced by a functorImageFilter using the appropriate
// functor. However using functorImageFilter with neighborhood is buggy (see issue #1802). Still,
// the functor has been implemented and localRxDetectionFilter will be deprecated when the
// bug is corrected.
#if 1 // Using localRxDetectionFilter
auto localRxDetectionFilter = LocalRxDetectorFilterType::New();
localRxDetectionFilter->SetInput(inputImage);
unsigned int externalRadius = GetParameterInt("er");
unsigned int internalRadius = GetParameterInt("ir");
localRxDetectionFilter->SetInternalRadius(internalRadius);
localRxDetectionFilter->SetExternalRadius(externalRadius);
#else // Using a functorImageFilter
Functor::LocalRxDetectionFunctor<double> detectorFunctor;
detectorFunctor.SetInternalRadius(GetParameterInt("ir"), GetParameterInt("ir"));
auto localRxDetectionFilter = otb::NewFunctorFilter
(detectorFunctor ,{{GetParameterInt("er"),GetParameterInt("er")}});
localRxDetectionFilter->SetInputs(inputImage);
#endif
SetParameterOutputImage("out", localRxDetectionFilter->GetOutput());
RegisterPipeline();
}
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::LocalRxDetection)
......@@ -26,6 +26,8 @@ otb_module(OTBAppHyperspectral
OTBApplicationEngine
OTBEndmembersExtraction
OTBUnmixing
OTBAnomalyDetection
OTBFunctor
TEST_DEPENDS
OTBTestKernel
OTBCommandLine
......
......@@ -73,4 +73,16 @@ otb_test_application(NAME apTvHyEndmemberNumberEstimation_elm
VALID --compare-ascii ${EPSILON_7}
${BASELINE_FILES}/aptTvHyEndmemberNumberEstimation.txt
${TEMP}/aptTvHyEndmemberNumberEstimation_elm.txt
)
\ No newline at end of file
)
#----------- LocalRxDetection TESTS ------------------------
otb_test_application(NAME apTvHyLocalRxDetection
APP LocalRxDetection
OPTIONS -in ${OTB_DATA_ROOT}/Input/Hyperspectral/synthetic/hsi_cube.tif?&bands=1:10
-out ${TEMP}/apTvHyLocalRxDetection.tif double
-ir 1
-er 3
VALID --compare-image ${EPSILON_9}
${BASELINE}/apTvHyLocalRxDetection.tif
${TEMP}/apTvHyLocalRxDetection.tif
)
......@@ -114,6 +114,118 @@ private:
};
/** \class LocalRxDetectionFunctor
* \brief This functor computes a local Rx score on an input neighborhood. Pixel of the neighborhood
* inside the internal radius are not considered during the computation of local statistics.
*
* \ingroup ImageFilters
*
* \ingroup OTBAnomalyDetection
*/
namespace Functor
{
template<typename TInput, typename TOutput =TInput>
class LocalRxDetectionFunctor
{
public:
/** typedef */
typedef typename itk::Neighborhood<itk::VariableLengthVector<TInput>>::OffsetType OffsetType;
typedef typename itk::VariableLengthVector<TInput> VectorMeasurementType;
typedef itk::Statistics::ListSample<VectorMeasurementType> ListSampleType;
typedef itk::Statistics::CovarianceSampleFilter<ListSampleType> CovarianceCalculatorType;
typedef typename CovarianceCalculatorType::MeasurementVectorRealType MeasurementVectorRealType;
typedef typename CovarianceCalculatorType::MatrixType MatrixType;
private:
// Internal radius along the X axis
unsigned int m_InternalRadiusX;
// Internal radius along the Y axis
unsigned int m_InternalRadiusY;
public:
LocalRxDetectionFunctor() : m_InternalRadiusX(1), m_InternalRadiusY(1) {};
void SetInternalRadius(const unsigned int internalRadiusX, const unsigned int internalRadiusY)
{
m_InternalRadiusX = internalRadiusX;
m_InternalRadiusY = internalRadiusY;
};
int GetInternalRadiusX()
{
return m_InternalRadiusX;
};
int GetInternalRadiusY()
{
return m_InternalRadiusY;
};
auto operator()(const itk::Neighborhood<itk::VariableLengthVector<TInput>> & in) const
{
// Create a list sample with the pixels of the neighborhood located between
// the two radius.
typename ListSampleType::Pointer listSample = ListSampleType::New();
// The pixel on whih we will compute the Rx score, we load it now to get the input vector size.
auto centerPixel = in.GetCenterValue();
listSample->SetMeasurementVectorSize(centerPixel.Size());
OffsetType off;
// Cache radiuses attributes for threading performances
const int internalRadiusX = m_InternalRadiusX;
const int internalRadiusY = m_InternalRadiusY;
auto externalRadius = in.GetRadius();
for (int y = -externalRadius[1]; y <= static_cast<int>(externalRadius[1]); y++)
{
off[1] = y;
for (int x = -externalRadius[0]; x <= static_cast<int>(externalRadius[0]); x++)
{
off[0] = x;
if ((abs(x) > internalRadiusX) || (abs(y) > internalRadiusY))
{
listSample->PushBack(in[off]);
}
}
}
// Compute mean & inverse covariance matrix
typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
covarianceCalculator->SetInput(listSample);
covarianceCalculator->Update();
MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
{
meanVec.SetElement(i, meanVector.GetElement(i));
}
MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();
typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
{
centeredTestPixMat.put(i, 0, (centerPixel.GetElement(i) - meanVector.GetElement(i)));
}
// Rx score computation
typename MatrixType::InternalMatrixType rxValue
= centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;
return static_cast<TOutput> (rxValue.get(0, 0));
}
};
} // end namespace functor
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
......
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* 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 otbLocalRxDetectorNonThreadFilter_h
#define otbLocalRxDetectorNonThreadFilter_h
#include "itkImageToImageFilter.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkImageRegionIterator.h"
#include "itkListSample.h"
#include "itkCovarianceSampleFilter.h"
namespace otb
{
/** \class otbLocalRxDetectorNonThreadFilter
* \brief Local-RX detector algorithm with multichannel VectorImage data as input
*
*
* \ingroup ImageFilters
*
* \ingroup OTBAnomalyDetection
*/
template <class TInputImage, class TOutputImage>
class ITK_EXPORT LocalRxDetectorNonThreadFilter:
public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:
/** Standard class typedefs. */
typedef LocalRxDetectorNonThreadFilter Self;
typedef itk::ImageToImageFilter< TInputImage, TOutputImage > Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Type macro */
itkNewMacro(Self);
/** Creation through object factory macro */
itkTypeMacro(LocalRxDetectorNonThreadFilter, ImageToImageFilter);
/** typedef related to input and output images */
typedef TInputImage InputImageType;
typedef typename InputImageType::Pointer InputPointerType;
typedef typename InputImageType::ConstPointer InputConstPointerType;
typedef typename InputImageType::IndexType InputIndexType;
typedef typename InputImageType::SizeType InputSizeType;
typedef TOutputImage OutputImageType;
typedef typename OutputImageType::Pointer OutputPointerType;
typedef typename OutputImageType::IndexType OutputIndexType;
typedef typename OutputImageType::OffsetType OutputOffsetType;
typedef typename OutputImageType::SizeType OutputSizeType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
/** typedef related to iterators */
typedef itk::ConstShapedNeighborhoodIterator<InputImageType> ConstShapedNeighborhoodIteratorType;
typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> VectorFaceCalculatorType;
typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<OutputImageType> FaceCalculatorType;
typedef itk::ImageRegionIterator<OutputImageType> ImageRegionIteratorType;
/** typedef related to statistics */
typedef typename InputImageType::PixelType VectorMeasurementType;
typedef itk::Statistics::ListSample<VectorMeasurementType> ListSampleType;
typedef itk::Statistics::CovarianceSampleFilter<ListSampleType> CovarianceCalculatorType;
typedef typename CovarianceCalculatorType::MeasurementVectorRealType MeasurementVectorRealType;
typedef typename CovarianceCalculatorType::MatrixType MatrixType;
/** Getter and Setter */
itkSetMacro(InternalRadius, int);
itkGetMacro(InternalRadius, int);
itkSetMacro(ExternalRadius, int);
itkGetMacro(ExternalRadius, int);
/** Main computation method */
void GenerateInputRequestedRegion() override;
void GenerateData() override;
// virtual void BeforeThreadedGenerateData();
// virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId);
protected:
LocalRxDetectorNonThreadFilter();
~LocalRxDetectorNonThreadFilter() override {}
void PrintSelf(std::ostream& os, itk::Indent indent) const override;
private:
LocalRxDetectorNonThreadFilter(const Self&) = delete;
void operator=(const Self&) = delete;
int m_InternalRadius;
int m_ExternalRadius;
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbLocalRxDetectorNonThreadFilter.hxx"
#endif
#endif
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* 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 otbLocalRxDetectorNonThreadFilter_hxx
#define otbLocalRxDetectorNonThreadFilter_hxx
#include "otbLocalRxDetectorNonThreadFilter.h"
namespace otb
{
/**
*
*/
template <class TInputImage, class TOutputImage>
LocalRxDetectorNonThreadFilter<TInputImage, TOutputImage>
::LocalRxDetectorNonThreadFilter()
{
this->m_ExternalRadius = 0;
this->m_InternalRadius = 0;
}
/**
*
*/
template <class TInputImage, class TOutputImage>
void
LocalRxDetectorNonThreadFilter<TInputImage, TOutputImage>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "Internal Radius: " << m_InternalRadius << std::endl;
os << indent << "External Radius: " << m_ExternalRadius << std::endl;
}
/**
*
*/
template <class TInputImage, class TOutputImage>
void
LocalRxDetectorNonThreadFilter<TInputImage, TOutputImage>
::GenerateData()
{
this->AllocateOutputs();
// Get the input and output pointers
InputPointerType inputPtr =
const_cast< InputImageType * >( this->GetInput());
OutputPointerType outputPtr = this->GetOutput();
inputPtr->Update();
// Fill the output buffer with black pixels
outputPtr->FillBuffer(0);
// Iterator on input region
typename ConstShapedNeighborhoodIteratorType::RadiusType radius;
radius.Fill(m_ExternalRadius);
VectorFaceCalculatorType vectorFaceCalculator;
typename VectorFaceCalculatorType::FaceListType vectorFaceList;
typename VectorFaceCalculatorType::FaceListType::iterator vectorFit;
vectorFaceList = vectorFaceCalculator(inputPtr, inputPtr->GetRequestedRegion(), radius);
vectorFit = vectorFaceList.begin(); // Only the first face is used
ConstShapedNeighborhoodIteratorType inputIt(radius, inputPtr, *vectorFit);
// Neighborhood Configuration
typename ConstShapedNeighborhoodIteratorType::OffsetType off;
for (int y = -m_ExternalRadius; y <= m_ExternalRadius; y++)
{
off[1] = y;
for (int x = -m_ExternalRadius; x <= m_ExternalRadius; x++)
{
off[0] = x;
if ((abs(x) > m_InternalRadius) || (abs(y) > m_InternalRadius))
{
inputIt.ActivateOffset(off);
}
}
}
// Iterator on output region
FaceCalculatorType faceCalculator;
typename FaceCalculatorType::FaceListType faceList;
typename FaceCalculatorType::FaceListType::iterator fit;
faceList = faceCalculator(outputPtr, inputPtr->GetRequestedRegion(), radius);
fit = faceList.begin(); // Only the first face is used
ImageRegionIteratorType outputIt(outputPtr, *fit);
// Run Input Image
for (inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
{
// Create ListSample
typename ListSampleType::Pointer listSample = ListSampleType::New();
listSample->SetMeasurementVectorSize(inputPtr->GetNumberOfComponentsPerPixel());
// Run neighborhood
typename ConstShapedNeighborhoodIteratorType::ConstIterator ci;
for (ci = inputIt.Begin(); !ci.IsAtEnd(); ++ci)
{
// Pushback element in listSample
//std::cout << "pixel of shaped iteror : " << ci.Get() << std::endl;
listSample->PushBack(ci.Get());
}
// Compute mean & covariance matrix
typename CovarianceCalculatorType::Pointer covarianceCalculator = CovarianceCalculatorType::New();
covarianceCalculator->SetInput(listSample);
covarianceCalculator->Update();
MeasurementVectorRealType meanVector = covarianceCalculator->GetMean();
MatrixType covarianceMatrix = covarianceCalculator->GetCovarianceMatrix();
// Compute RX value
typename MatrixType::InternalMatrixType invCovMat = covarianceMatrix.GetInverse();
VectorMeasurementType testPixVec;
testPixVec = inputPtr->GetPixel(inputIt.GetIndex());
VectorMeasurementType meanVec(meanVector.GetNumberOfElements());
for(unsigned int i = 0; i < meanVector.GetNumberOfElements(); ++i)
{
meanVec.SetElement(i, meanVector.GetElement(i));
}
typename MatrixType::InternalMatrixType centeredTestPixMat(meanVector.GetNumberOfElements(), 1);
for (unsigned int i = 0; i < centeredTestPixMat.rows(); ++i)
{
centeredTestPixMat.put(i, 0, (testPixVec.GetElement(i) - meanVector.GetElement(i)));
}
typename MatrixType::InternalMatrixType rxValue = centeredTestPixMat.transpose() * invCovMat * centeredTestPixMat;
outputIt.Set(rxValue.get(0, 0));
}
}
/**
*
*/
template <class TInputImage, class TOutputImage>
void
LocalRxDetectorNonThreadFilter<TInputImage, TOutputImage>
::GenerateInputRequestedRegion()
{
// Call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// Get pointers to the input and output
InputPointerType inputPtr =
const_cast< InputImageType * >( this->GetInput());
OutputPointerType outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
// Get a copy of the input requested region (should equal the output
// requested region)
typename TInputImage::RegionType inputRequestedRegion;
inputRequestedRegion = inputPtr->GetRequestedRegion();
// Pad the input requested region by the operator radius
inputRequestedRegion.PadByRadius( m_ExternalRadius );
// Crop the input requested region at the input's largest possible region
if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
{
inputPtr->SetRequestedRegion( inputRequestedRegion );
return;
}
else
{
// Couldn't crop the region (requested region is outside the largest
// possible region). Throw an exception.
// Store what we tried to request (prior to trying to crop)
inputPtr->SetRequestedRegion( inputRequestedRegion );
// Build an exception
itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
e.SetLocation(ITK_LOCATION);
e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
e.SetDataObject(inputPtr);
throw e;
}
}
} // end namespace otb
#endif
......@@ -22,7 +22,6 @@ otb_module_test()
set(OTBAnomalyDetectionTests
otbAnomalyDetectionTestDriver.cxx
otbLocalRxDetectorRoiTest.cxx
otbLocalRxDetectorTest.cxx
)
......@@ -32,3 +31,13 @@ otb_module_target_label(otbAnomalyDetectionTestDriver)
# Tests Declaration
otb_add_test(NAME hyTvLocalRxDetectorFilter COMMAND otbAnomalyDetectionTestDriver
--compare-image ${EPSILON_9}
${BASELINE}/hyTvLocalRxDetectorFilter.tif
${TEMP}/hyTvLocalRxDetectorFilter.tif
LocalRXDetectorTest
${INPUTDATA}/cupriteSubHsi.tif
${TEMP}/hyTvLocalRxDetectorFilter.tif
3
1
)
\ No newline at end of file
......@@ -22,6 +22,5 @@
void RegisterTests()
{
REGISTER_TEST(LocalRXDetectorROITest);
REGISTER_TEST(LocalRXDetectorTest);
}
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* 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