Commit 89359fb0 authored by Julien Michel's avatar Julien Michel
Browse files

Merge branch 'refactor-sinclairimagefilter' into 'develop'

Refactor filters in Polarimetry module to use new FunctorFilter

See merge request !318
parents b6432bd9 36dd35bf
......@@ -21,20 +21,17 @@
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbReciprocalHAlphaImageFilter.h"
#include "otbReciprocalBarnesDecompImageFilter.h"
#include "otbReciprocalHuynenDecompImageFilter.h"
#include "otbReciprocalPauliDecompImageFilter.h"
#include "otbReciprocalHAlphaImageFilter.h"
#include "otbSinclairReciprocalImageFilter.h"
#include "otbSinclairToReciprocalCoherencyMatrixFunctor.h"
#include "otbPerBandVectorImageFilter.h"
#include "itkMeanImageFilter.h"
// #include "otbNRIBandImagesToOneNComplexBandsImage.h"
#include "otbImageListToVectorImageFilter.h"
#include "otbImageList.h"
#include "otbSinclairToReciprocalCoherencyMatrixImageFilter.h"
namespace otb
{
......@@ -50,28 +47,22 @@ public:
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef otb::Functor::SinclairToReciprocalCoherencyMatrixFunctor<ComplexDoubleImageType::PixelType,
ComplexDoubleImageType::PixelType,
ComplexDoubleImageType::PixelType,
ComplexDoubleVectorImageType::PixelType> FunctorType;
typedef SinclairReciprocalImageFilter<ComplexDoubleImageType,
ComplexDoubleImageType,
ComplexDoubleImageType,
ComplexDoubleVectorImageType,
FunctorType > SRFilterType;
using SRFilterType = otb::SinclairToReciprocalCoherencyMatrixImageFilter<ComplexDoubleImageType, ComplexDoubleVectorImageType>;
typedef itk::MeanImageFilter<ComplexDoubleImageType, ComplexDoubleImageType> MeanFilterType;
typedef otb::PerBandVectorImageFilter<ComplexDoubleVectorImageType, ComplexDoubleVectorImageType, MeanFilterType> PerBandMeanFilterType;
//typedef otb::NRIBandImagesToOneNComplexBandsImage<DoubleVectorImageType, ComplexDoubleVectorImageType> NRITOOneCFilterType;
typedef otb::ImageList<ComplexDoubleImageType> ImageListType;
typedef ImageListToVectorImageFilter<ImageListType, ComplexDoubleVectorImageType > ListConcatenerFilterType;
typedef otb::ReciprocalHAlphaImageFilter<ComplexDoubleVectorImageType, ComplexDoubleVectorImageType> HAFilterType;
typedef otb::ReciprocalBarnesDecompImageFilter<ComplexDoubleVectorImageType, ComplexDoubleVectorImageType> BarnesFilterType;
typedef otb::ReciprocalHuynenDecompImageFilter<ComplexDoubleVectorImageType, ComplexDoubleVectorImageType> HuynenFilterType;
......@@ -91,17 +82,17 @@ private:
// Documentation
SetDocName("SARDecompositions");
SetDocLongDescription("From one-band complex images (HH, HV, VH, VV), returns the selected decomposition.\n \n"
"All the decompositions implemented are intended for the mono-static case (transmitter and receiver are co-located).\n"
"There are two kinds of decomposition: coherent ones and incoherent ones.\n"
"In the coherent case, only the Pauli decomposition is available.\n"
"In the incoherent case, there the decompositions available: Huynen, Barnes, and H-alpha-A.\n"
"User must provide three one-band complex images HH, HV or VH, and VV (mono-static case <=> HV = VH).\n"
"Incoherent decompositions consist in averaging 3x3 complex coherency/covariance matrices; the user must provide the size of the averaging window, thanks to the parameter inco.kernelsize."
);
SetDocLongDescription(
"From one-band complex images (HH, HV, VH, VV), returns the selected decomposition.\n \n"
"All the decompositions implemented are intended for the mono-static case (transmitter and receiver are co-located).\n"
"There are two kinds of decomposition: coherent ones and incoherent ones.\n"
"In the coherent case, only the Pauli decomposition is available.\n"
"In the incoherent case, there the decompositions available: Huynen, Barnes, and H-alpha-A.\n"
"User must provide three one-band complex images HH, HV or VH, and VV (mono-static case <=> HV = VH).\n"
"Incoherent decompositions consist in averaging 3x3 complex coherency/covariance matrices; the user must provide the size of the averaging window, "
"thanks to the parameter inco.kernelsize.");
SetDocLimitations("Some decompositions output real images, while this application outputs complex images for general purpose.\n"
"Users should pay attention to extract the real part of the results provided by this application.\n");
SetDocAuthors("OTB-Team");
......@@ -111,21 +102,21 @@ private:
AddParameter(ParameterType_InputImage, "inhh", "Input Image");
SetParameterDescription("inhh", "Input image (HH)");
AddParameter(ParameterType_InputImage, "inhv", "Input Image");
SetParameterDescription("inhv", "Input image (HV)");
MandatoryOff("inhv");
AddParameter(ParameterType_InputImage, "invh", "Input Image");
SetParameterDescription("invh", "Input image (VH)");
MandatoryOff("invh");
AddParameter(ParameterType_InputImage, "invv", "Input Image");
SetParameterDescription("invv", "Input image (VV)");
AddParameter(ParameterType_OutputImage, "out", "Output Image");
SetParameterDescription("out", "Output image");
AddParameter(ParameterType_Choice, "decomp", "Decompositions");
AddChoice("decomp.haa","H-alpha-A incoherent decomposition");
SetParameterDescription("decomp.haa","H-alpha-A incoherent decomposition");
......@@ -135,10 +126,9 @@ private:
SetParameterDescription("decomp.huynen","Huynen incoherent decomposition");
AddChoice("decomp.pauli","Pauli coherent decomposition");
SetParameterDescription("decomp.pauli","Pauli coherent decomposition");
AddParameter(ParameterType_Group,"inco","Incoherent decompositions");
SetParameterDescription("inco","This group allows setting parameters related to the incoherent decompositions.");
AddParameter(ParameterType_Int, "inco.kernelsize", "Kernel size for spatial incoherent averaging");
SetParameterDescription("inco.kernelsize", "Minute (0-59)");
SetMinimumParameterIntValue("inco.kernelsize", 1);
......@@ -167,83 +157,83 @@ private:
void DoExecute() override
{
bool inhv = HasUserValue("inhv");
bool inhv = HasUserValue("inhv");
bool invh = HasUserValue("invh");
if ( (!inhv) && (!invh) )
if ( (!inhv) && (!invh) )
otbAppLogFATAL( << "Parameter inhv or invh not set. Please provide a HV or a VH complex image.");
m_SRFilter = SRFilterType::New();
m_HAFilter = HAFilterType::New();
m_MeanFilter = PerBandMeanFilterType::New();
m_HAFilter = HAFilterType::New();
m_MeanFilter = PerBandMeanFilterType::New();
MeanFilterType::InputSizeType radius;
m_BarnesFilter = BarnesFilterType::New();
m_HuynenFilter = HuynenFilterType::New();
m_PauliFilter = PauliFilterType::New();
m_Concatener = ListConcatenerFilterType::New();
m_ImageList = ImageListType::New();
m_ImageList = ImageListType::New();
switch (GetParameterInt("decomp"))
{
case 0: // H-alpha-A
if (inhv)
m_SRFilter->SetInputHV_VH(GetParameterComplexDoubleImage("inhv"));
else if (invh)
m_SRFilter->SetInputHV_VH(GetParameterComplexDoubleImage("invh"));
m_SRFilter->SetInputHH(GetParameterComplexDoubleImage("inhh"));
m_SRFilter->SetInputVV(GetParameterComplexDoubleImage("invv"));
radius.Fill( GetParameterInt("inco.kernelsize") );
m_MeanFilter->GetFilter()->SetRadius(radius);
m_MeanFilter->SetInput(m_SRFilter->GetOutput());
m_HAFilter->SetInput(m_MeanFilter->GetOutput());
SetParameterOutputImage("out", m_HAFilter->GetOutput() );
break;
case 1: // Barnes
if (inhv)
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hv_or_vh>(GetParameterComplexDoubleImage("inhv"));
else if (invh)
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hv_or_vh>(GetParameterComplexDoubleImage("invh"));
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hh>(GetParameterComplexDoubleImage("inhh"));
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::vv>(GetParameterComplexDoubleImage("invv"));
radius.Fill(GetParameterInt("inco.kernelsize"));
m_MeanFilter->GetFilter()->SetRadius(radius);
m_MeanFilter->SetInput(m_SRFilter->GetOutput());
m_HAFilter->SetVariadicInput<0>(m_MeanFilter->GetOutput());
SetParameterOutputImage("out", m_HAFilter->GetOutput());
break;
case 1: // Barnes
if (inhv)
m_SRFilter->SetInputHV_VH(GetParameterComplexDoubleImage("inhv"));
else if (invh)
m_SRFilter->SetInputHV_VH(GetParameterComplexDoubleImage("invh"));
m_SRFilter->SetInputHH(GetParameterComplexDoubleImage("inhh"));
m_SRFilter->SetInputVV(GetParameterComplexDoubleImage("invv"));
radius.Fill( GetParameterInt("inco.kernelsize") );
m_MeanFilter->GetFilter()->SetRadius(radius);
m_MeanFilter->SetInput(m_SRFilter->GetOutput());
m_BarnesFilter->SetInput(m_MeanFilter->GetOutput());
SetParameterOutputImage("out", m_BarnesFilter->GetOutput() );
break;
case 2: // Huynen
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hv_or_vh>(GetParameterComplexDoubleImage("inhv"));
else if (invh)
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hv_or_vh>(GetParameterComplexDoubleImage("invh"));
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hh>(GetParameterComplexDoubleImage("inhh"));
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::vv>(GetParameterComplexDoubleImage("invv"));
radius.Fill(GetParameterInt("inco.kernelsize"));
m_MeanFilter->GetFilter()->SetRadius(radius);
m_MeanFilter->SetInput(m_SRFilter->GetOutput());
m_BarnesFilter->SetVariadicInput<0>(m_MeanFilter->GetOutput());
SetParameterOutputImage("out", m_BarnesFilter->GetOutput());
break;
case 2: // Huynen
if (inhv)
m_SRFilter->SetInputHV_VH(GetParameterComplexDoubleImage("inhv"));
else if (invh)
m_SRFilter->SetInputHV_VH(GetParameterComplexDoubleImage("invh"));
m_SRFilter->SetInputHH(GetParameterComplexDoubleImage("inhh"));
m_SRFilter->SetInputVV(GetParameterComplexDoubleImage("invv"));
radius.Fill( GetParameterInt("inco.kernelsize") );
m_MeanFilter->GetFilter()->SetRadius(radius);
m_MeanFilter->SetInput(m_SRFilter->GetOutput());
m_HuynenFilter->SetInput(m_MeanFilter->GetOutput());
SetParameterOutputImage("out", m_HuynenFilter->GetOutput() );
break;
case 3: // Pauli
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hv_or_vh>(GetParameterComplexDoubleImage("inhv"));
else if (invh)
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hv_or_vh>(GetParameterComplexDoubleImage("invh"));
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::hh>(GetParameterComplexDoubleImage("inhh"));
m_SRFilter->SetVariadicNamedInput<polarimetry_tags::vv>(GetParameterComplexDoubleImage("invv"));
radius.Fill(GetParameterInt("inco.kernelsize"));
m_MeanFilter->GetFilter()->SetRadius(radius);
m_MeanFilter->SetInput(m_SRFilter->GetOutput());
m_HuynenFilter->SetVariadicInput<0>(m_MeanFilter->GetOutput());
SetParameterOutputImage("out", m_HuynenFilter->GetOutput());
break;
case 3: // Pauli
m_ImageList->PushBack(GetParameterComplexDoubleImage("inhh"));
......@@ -253,18 +243,16 @@ private:
m_ImageList->PushBack(GetParameterComplexDoubleImage("invh"));
m_ImageList->PushBack(GetParameterComplexDoubleImage("invv"));
m_Concatener->SetInput( m_ImageList );
m_PauliFilter->SetInput(m_Concatener->GetOutput());
SetParameterOutputImage("out", m_PauliFilter->GetOutput() );
break;
m_Concatener->SetInput(m_ImageList);
m_PauliFilter->SetVariadicInput<0>(m_Concatener->GetOutput());
SetParameterOutputImage("out", m_PauliFilter->GetOutput() );
break;
}
}
//MCPSFilterType::Pointer m_MCPSFilter;
SRFilterType::Pointer m_SRFilter;
HAFilterType::Pointer m_HAFilter;
BarnesFilterType::Pointer m_BarnesFilter;
......@@ -273,8 +261,7 @@ private:
PerBandMeanFilterType::Pointer m_MeanFilter;
ListConcatenerFilterType::Pointer m_Concatener;
ImageListType::Pointer m_ImageList;
};
};
} //end namespace Wrapper
} //end namespace otb
......
......@@ -22,12 +22,13 @@
#ifndef otbMuellerToPolarisationDegreeAndPowerImageFilter_h
#define otbMuellerToPolarisationDegreeAndPowerImageFilter_h
#include "otbUnaryFunctorImageFilter.h"
#include "itkNumericTraits.h"
#include "itkMatrix.h"
#include "itkVector.h"
#include "otbMath.h"
#include "otbFunctorImageFilter.h"
namespace otb
{
......@@ -67,6 +68,8 @@ namespace Functor {
* - channel #2 : \f$ DegP_{min} \f$
* - channel #3 : \f$ DegP_{max} \f$
*
* Use otb::MuellerToPolarisationDegreeAndPowerImageFilter to apply
*
* \ingroup Functor
* \ingroup SARPolarimetry
*
......@@ -84,8 +87,8 @@ public:
typedef itk::Matrix<double, 4, 4> MuellerMatrixType;
typedef itk::Vector<double, 4> StokesVectorType;
inline TOutput operator()( const TInput & Mueller ) const
{
inline void operator()(TOutput& result, const TInput& Mueller) const
{
double P;
double deg_pol;
double tau;
......@@ -98,9 +101,6 @@ public:
double l_PolarisationDegreeMin(itk::NumericTraits<double>::max());
double l_PolarisationDegreeMax(itk::NumericTraits<double>::min());
TOutput result;
result.SetSize(m_NumberOfComponentsPerPixel);
MuellerMatrixType muellerMatrix;
muellerMatrix[0][0] = Mueller[0];
muellerMatrix[0][1] = Mueller[1];
......@@ -172,70 +172,35 @@ public:
result[1] = l_PowerMax;
result[2] = l_PolarisationDegreeMin;
result[3] = l_PolarisationDegreeMax;
return result;
}
unsigned int GetOutputSize()
{
return m_NumberOfComponentsPerPixel;
constexpr size_t OutputSize(...) const
{
// Size of the result
return 4;
}
/** Constructor */
MuellerToPolarisationDegreeAndPowerFunctor() : m_NumberOfComponentsPerPixel(4), m_Epsilon(1e-6), m_PI_90(2*CONST_PI_180) {}
/** Destructor */
virtual ~MuellerToPolarisationDegreeAndPowerFunctor() {}
private:
unsigned int m_NumberOfComponentsPerPixel;
const double m_Epsilon;
const double m_PI_90;
static constexpr double m_Epsilon = 1e-6;
static constexpr double m_PI_90 = 2 * CONST_PI_180;
};
}
} // namespace Functor
/** \class otbMuellerToPolarisationDegreeAndPowerImageFilter
* \brief Compute the polarization degree and power (4 channels : Power min and max, Polarization degree min and max)
* from the Mueller image (16 real channels)
* For more details, please refer to the class MuellerToPolarisationDegreeAndPowerFunctor.
*
* \ingroup SARPolarimetry
* \sa MuellerToPolarisationDegreeAndPowerFunctor
/**
* \typedef MuellerToPolarisationDegreeAndPowerImageFilter
* \brief Applies otb::Functor::MuellerToPolarisationDegreeAndPowerFunctor
* \sa otb::Functor::MuellerToPolarisationDegreeAndPowerFunctor
*
* Set inputs with:
* \code
* SetVariadicInput<0>(inputPtr);
* \endcode
*
* \ingroup OTBPolarimetry
*/
template <class TInputImage, class TOutputImage>
class ITK_EXPORT MuellerToPolarisationDegreeAndPowerImageFilter :
public UnaryFunctorImageFilter<TInputImage, TOutputImage, Functor::MuellerToPolarisationDegreeAndPowerFunctor<
typename TInputImage::PixelType, typename TOutputImage::PixelType> >
{
public:
/** Standard class typedefs. */
typedef MuellerToPolarisationDegreeAndPowerImageFilter Self;
typedef typename Functor::MuellerToPolarisationDegreeAndPowerFunctor<
typename TInputImage::PixelType, typename TOutputImage::PixelType> FunctionType;
typedef UnaryFunctorImageFilter<TInputImage, TOutputImage, FunctionType> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
template <typename TInputImage, typename TOutputImage>
using MuellerToPolarisationDegreeAndPowerImageFilter =
FunctorImageFilter<Functor::MuellerToPolarisationDegreeAndPowerFunctor<typename TInputImage::PixelType, typename TOutputImage::PixelType>>;
/** Runtime information support. */
itkTypeMacro(MuellerToPolarisationDegreeAndPowerImageFilter, UnaryFunctorImageFilter);
protected:
MuellerToPolarisationDegreeAndPowerImageFilter() {}
~MuellerToPolarisationDegreeAndPowerImageFilter() override {}
private:
MuellerToPolarisationDegreeAndPowerImageFilter(const Self&) = delete;
void operator=(const Self&) = delete;
};
} // end namespace otb
......
......@@ -22,7 +22,7 @@
#ifndef otbMuellerToReciprocalCovarianceImageFilter_h
#define otbMuellerToReciprocalCovarianceImageFilter_h
#include "otbUnaryFunctorImageFilter.h"
#include "otbFunctorImageFilter.h"
namespace otb
{
......@@ -54,6 +54,8 @@ namespace Functor {
* The output image has 6 channels : the diagonal and the upper element of the reciprocal matrix.
* Element are stored from left to right, line by line.
*
* Use otb::MuellerToReciprocalCovarianceImageFilter to apply
*
* \ingroup Functor
* \ingroup SARPolarimetry
*
......@@ -71,12 +73,8 @@ public:
typedef std::complex<double> ComplexType;
typedef typename TOutput::ValueType OutputValueType;
inline TOutput operator()( const TInput & Mueller ) const
inline void operator()(TOutput& result, const TInput& Mueller) const
{
TOutput result;
result.SetSize(m_NumberOfComponentsPerPixel);
// Keep the upper diagonal of the matrix
const double M11 = static_cast<double>(Mueller[0]);
const double M12 = static_cast<double>(Mueller[1]);
......@@ -89,7 +87,7 @@ public:
const double M34 = static_cast<double>(Mueller[11]);
const double M44 = static_cast<double>(Mueller[15]);
const ComplexType A(0.5*(M11+M22+2*M12));
const ComplexType B(0.5*std::sqrt(2.0)*(M13+M23), 0.5*std::sqrt(2.0)*(M14+M24));
const ComplexType C(-0.5*(M33+M44), -M34);
......@@ -103,67 +101,32 @@ public:
result[3] = static_cast<OutputValueType>( E );
result[4] = static_cast<OutputValueType>( F );
result[5] = static_cast<OutputValueType>( I );
return result;
}
unsigned int GetOutputSize()
{
return m_NumberOfComponentsPerPixel;
}
/** Constructor */
MuellerToReciprocalCovarianceFunctor() : m_NumberOfComponentsPerPixel(6) {}
/** Destructor */
virtual ~MuellerToReciprocalCovarianceFunctor() {}
private:
unsigned int m_NumberOfComponentsPerPixel;
constexpr size_t OutputSize(...) const
{
// Size of the result
return 6;
}
};
}
} // namespace Functor
/** \class otbMuellerToReciprocalCovarianceImageFilter
* \brief Compute the MLC image
* from the Mueller image (16 real channels)
* For more details, lease refer to the class MuellerToReciprocalCovarianceFunctor.
*
* \ingroup SARPolarimetry
* \sa MuellerToReciprocalCovarianceFunctor
/**
* \typedef MuellerToReciprocalCovarianceImageFilter
* \brief Applies otb::Functor::MuellerToReciprocalCovarianceFunctor
* \sa otb::Functor::MuellerToReciprocalCovarianceFunctor
*
* Set inputs with:
* \code
* SetVariadicInput<0>(inputPtr);
* \endcode
*
* \ingroup OTBPolarimetry
*/
template <class TInputImage, class TOutputImage>
class ITK_EXPORT MuellerToReciprocalCovarianceImageFilter :
public UnaryFunctorImageFilter<TInputImage, TOutputImage, Functor::MuellerToReciprocalCovarianceFunctor<
typename TInputImage::PixelType, typename TOutputImage::PixelType> >
{
public:
/** Standard class typedefs. */
typedef MuellerToReciprocalCovarianceImageFilter Self;
typedef Functor::MuellerToReciprocalCovarianceFunctor<
typename TInputImage::PixelType, typename TOutputImage::PixelType> FunctionType;
typedef UnaryFunctorImageFilter<TInputImage, TOutputImage, FunctionType> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
template <typename TInputImage, typename TOutputImage>
using MuellerToReciprocalCovarianceImageFilter =
FunctorImageFilter<Functor::MuellerToReciprocalCovarianceFunctor<typename TInputImage::PixelType, typename TOutputImage::PixelType>>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Runtime information support. */
itkTypeMacro(MuellerToReciprocalCovarianceImageFilter, UnaryFunctorImageFilter);
protected:
MuellerToReciprocalCovarianceImageFilter() {}
~MuellerToReciprocalCovarianceImageFilter() override {}
private:
MuellerToReciprocalCovarianceImageFilter(const Self&) = delete;
void operator=(const Self&) = delete;
};
} // end namespace otb
......
/*
* Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/