diff --git a/Code/Fusion/otbLmvmPanSharpeningFusionImageFilter.h b/Code/Fusion/otbLmvmPanSharpeningFusionImageFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..822647d89a390b09ad2fd61bbf4dc79a8f4c64b0 --- /dev/null +++ b/Code/Fusion/otbLmvmPanSharpeningFusionImageFilter.h @@ -0,0 +1,237 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + Some parts of this code are derived from ITK. See ITKCopyright.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 __otbLmvmPanSharpeningFusionImageFilter_h +#define __otbLmvmPanSharpeningFusionImageFilter_h + +#include "itkImageToImageFilter.h" +#include "otbConvolutionImageFilter.h" +#include "itkTernaryFunctorImageFilter.h" +#include "otbPerBandVectorImageFilter.h" +#include "itkNoiseImageFilter.h" +#include "otbVectorImage.h" + +#include "itkProgressAccumulator.h" + + +namespace otb +{ +/** + * \class LmvmPanSharpeningFusionImageFilter + * \brief This class performs a Local Mean and Variance Matching (LMVM) Pan sharpening operation + * + * Given a Pan image and the corresponding Xs image (oversampled to have the + * same number of pixels), this filter realizes an Local Mean and Variance + * Matching Pan sharpening operation as described in: + * + * de Bethune, S. and Muller, F. and Donnay, J.P. 1998. Fusion of multispectral and panchromatic images + * by local mean and variance matching filtering techniques. Fusion of Earth Data, 28-30 August 1998. + * + * \ingroup Streamed + * \ingroup Multithreaded + * \ingroup Fusion + * + * + **/ + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision = float> +class ITK_EXPORT LmvmPanSharpeningFusionImageFilter : + public itk::ImageToImageFilter<TXsImageType, TOutputImageType> +{ +public: + /** Standard class typedefs */ + typedef LmvmPanSharpeningFusionImageFilter Self; + typedef itk::ImageToImageFilter<TXsImageType, TOutputImageType> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Internal image type used as Pan smoothing and local standard deviation filter output */ + typedef otb::Image<TInternalPrecision, + TPanImageType::ImageDimension> InternalImageType; + + /** Single band Xs image type */ + typedef otb::Image<TXsImageType::PixelType, + TXsImageType::ImageDimension> XsBandImageType; + + /** Internal image type used as Xs smoothing and local standard deviation filter output */ + typedef otb::VectorImage<TInternalPrecision, + TXsImageType::ImageDimension> InternalVectorImageType; + + /** Typedef for the radius of the smoothing and local standard deviation filter */ + typedef typename itk::Array<TInternalPrecision> ArrayType; + + /** Method for creation through object factory */ + itkNewMacro(Self); + + /** Run-time type information */ + itkTypeMacro(LmvmPanSharpeningFusionImageFilter, + itk::ImageToImageFilter); + + /** Define the radius type for the smoothing and local standard deviation operation */ + typedef typename InternalImageType::SizeType RadiusType; + + /** Set the smoothing and local standard deviation filter radius */ + itkGetMacro(Radius, RadiusType); + itkSetMacro(Radius, RadiusType); + + /** Set the kernel used for the smoothing and local standard deviation filter */ + itkSetMacro(Filter, ArrayType); + itkGetConstReferenceMacro(Filter, ArrayType); + + virtual void SetPanInput(const TPanImageType * image); + const TPanImageType * GetPanInput(void) const; + + virtual void SetXsInput(const TXsImageType * path); + const TXsImageType * GetXsInput(void) const; + +protected: + /** Constructor */ + LmvmPanSharpeningFusionImageFilter(); + + /** Destructor */ + virtual ~LmvmPanSharpeningFusionImageFilter() {}; + + /** Call to generate data, wiring composite internal minipipeline */ + void GenerateData(); + + /** PrintSelf method */ + void PrintSelf(std::ostream& os, itk::Indent indent) const; + +private: + LmvmPanSharpeningFusionImageFilter(Self &); // intentionally not implemented + void operator =(const Self&); // intentionally not implemented + + /** \class FusionFunctor + * This functor applies the LMVM + * operation. It is intended for internal use only. + */ + class FusionFunctor + { + public: + // Implement the fusion as a six arguments operator + typename TOutputImageType::PixelType operator()(const typename TXsImageType::PixelType& xsPixel, + const TInternalPrecision& smoothPanchroPixel, + const typename InternalVectorImageType::PixelType& smoothXsPixel, + const typename InternalVectorImageType::PixelType& stdXsPixel, + const TInternalPrecision& stdPanchroPixel, + const typename TPanImageType::PixelType& sharpPanchroPixel) const + { + // Build output pixel + typename TOutputImageType::PixelType output(xsPixel.Size()); + + TInternalPrecision scale = 1.; + + if(vcl_abs(stdPanchroPixel) > 1e-10) + { + scale = 1.0/stdPanchroPixel; + } + + // Perform fusion for each band with appropriate casting + for(unsigned int i = 0; i < xsPixel.Size(); ++i) + { + output[i] = static_cast<typename TOutputImageType::InternalPixelType>( + (((sharpPanchroPixel - (smoothPanchroPixel * stdXsPixel[i])) * scale) + smoothXsPixel[i])); + } + // Returns the output pixel + return output; + } + }; + + /** + * Typedef of the TernaryFunctorImageFilter applying the fusion functor to + * p, p_smooth, p_std, xs_smooth, xs_std and xs. + */ + typedef itk::TernaryFunctorImageFilter<TXsImageType, + InternalImageType, + InternalVectorImageType, + InternalVectorImageType, + InternalImageType, + TPanImageType, + TOutputImageType, + FusionFunctor> FusionFilterType; + + /** Pointer to the fusion filter */ + typename FusionFilterType::Pointer m_FusionFilter; + + /** Typedef of the convolution filter performing Pan smoothing */ + typedef otb::ConvolutionImageFilter + <TPanImageType, + InternalImageType, + itk::ZeroFluxNeumannBoundaryCondition<TPanImageType>, + TInternalPrecision> PanConvolutionFilterType; + + /** Pointer to the internal Pan convolution filter */ + typename PanConvolutionFilterType::Pointer m_PanConvolutionFilter; + + /** Typedef of the convolution filter performing Xs smoothing */ + typedef otb::ConvolutionImageFilter + <XsBandImageType, + InternalImageType, + itk::ZeroFluxNeumannBoundaryCondition<TPanImageType>, + TInternalPrecision> XsConvolutionFilterType; + + /** Pointer to the internal Xs convolution filter */ + typename XsConvolutionFilterType::Pointer m_XsConvolutionFilter; + + /** Typedef of the Pan local standard deviation filter*/ + typedef itk::NoiseImageFilter + <TPanImageType, InternalImageType> PanNoiseFilterType; + + /** Pointer to the Pan local standard deviation filter */ + typename PanNoiseFilterType::Pointer m_PanNoiseFilter; + + /** Typedef of the Xs local standard deviation filter*/ + typedef itk::NoiseImageFilter + <XsBandImageType, InternalImageType> XsNoiseFilterType; + + /** Pointer to the Xs local standard deviation filter */ + typename XsNoiseFilterType::Pointer m_XsNoiseFilter; + + /** Typedef of a helper filter to apply to a vector image*/ + typedef otb::PerBandVectorImageFilter + < TXsImageType, InternalVectorImageType, XsConvolutionFilterType> XsVectorConvolutionFilterType; + + /** Pointer to the helper vector image filter */ + typename XsVectorConvolutionFilterType::Pointer m_XsVectorConvolutionFilter; + + /** Typedef of a helper filter to apply to a vector image*/ + typedef otb::PerBandVectorImageFilter + < TXsImageType, InternalVectorImageType, XsNoiseFilterType> XsVectorNoiseFilterType; + + /** Pointer to the helper vector image filter */ + typename XsVectorNoiseFilterType::Pointer m_XsVectorNoiseFilter; + + /** Radius used for the smoothing filter */ + RadiusType m_Radius; + + /** Kernel used for the smoothing filter */ + ArrayType m_Filter; + + /** The internal progress accumulator */ + typename itk::ProgressAccumulator::Pointer m_ProgressAccumulator; +}; + +} // end namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbLmvmPanSharpeningFusionImageFilter.txx" +#endif + +#endif diff --git a/Code/Fusion/otbLmvmPanSharpeningFusionImageFilter.txx b/Code/Fusion/otbLmvmPanSharpeningFusionImageFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..73b3983f8666bb3e32299389b52354a8aad54987 --- /dev/null +++ b/Code/Fusion/otbLmvmPanSharpeningFusionImageFilter.txx @@ -0,0 +1,181 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + Some parts of this code are derived from ITK. See ITKCopyright.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 __otbLmvmPanSharpeningFusionImageFilter_txx +#define __otbLmvmPanSharpeningFusionImageFilter_txx + +#include "otbLmvmPanSharpeningFusionImageFilter.h" + +namespace otb +{ +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::LmvmPanSharpeningFusionImageFilter() +{ + // Fix number of required inputs + this->SetNumberOfRequiredInputs(2); + + // Instantiate convolution filters + m_PanConvolutionFilter = PanConvolutionFilterType::New(); + m_PanConvolutionFilter->NormalizeFilterOn(); + m_XsConvolutionFilter = XsConvolutionFilterType::New(); + m_XsConvolutionFilter->NormalizeFilterOn(); + m_PanNoiseFilter = PanNoiseFilterType::New(); + m_XsNoiseFilter = XsNoiseFilterType::New(); + m_XsVectorConvolutionFilter = XsVectorConvolutionFilterType::New(); + m_XsVectorNoiseFilter = XsVectorNoiseFilterType::New(); + + + // Set-up default parameters + m_Radius.Fill(3); + m_Filter.SetSize(7 * 7); + m_Filter.Fill(1); + + // Instantiate fusion filter + m_FusionFilter = FusionFilterType::New(); + + // Set-up progress reporting + m_ProgressAccumulator = itk::ProgressAccumulator::New(); + m_ProgressAccumulator->SetMiniPipelineFilter(this); + m_ProgressAccumulator->RegisterInternalFilter(m_PanConvolutionFilter, 0.2); + m_ProgressAccumulator->RegisterInternalFilter(m_PanNoiseFilter, 0.2); + m_ProgressAccumulator->RegisterInternalFilter(m_XsVectorConvolutionFilter, 0.2); + m_ProgressAccumulator->RegisterInternalFilter(m_XsVectorNoiseFilter, 0.2); + m_ProgressAccumulator->RegisterInternalFilter(m_FusionFilter, 0.2); +} + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +void +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::SetPanInput(const TPanImageType *image) +{ + // We have 2 inputs: an image and a vector image + + // Process object is not const-correct so the const_cast is required here + this->itk::ProcessObject::SetNthInput(1, + const_cast<TPanImageType*>(image)); + this->Modified(); +} + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +const TPanImageType * +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::GetPanInput(void) const +{ + if (this->GetNumberOfInputs() < 2) + { + return 0; + } + + return static_cast<const TPanImageType *> + (this->itk::ProcessObject::GetInput(1)); +} + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +void +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::SetXsInput(const TXsImageType *image) +{ + // We have 2 inputs: an image and a vector image + + // Process object is not const-correct so the const_cast is required here + this->itk::ProcessObject::SetNthInput(0, + const_cast<TXsImageType*>(image)); + this->Modified(); +} + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +const TXsImageType * +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::GetXsInput(void) const +{ + if (this->GetNumberOfInputs() < 1) + { + return 0; + } + + return static_cast<const TXsImageType *> + (this->itk::ProcessObject::GetInput(0)); +} + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +void +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::GenerateData() +{ + //Check if size is correct + typename TPanImageType::SizeType sizePan; + typename TXsImageType::SizeType sizeXs; + sizePan = this->GetPanInput()->GetLargestPossibleRegion().GetSize(); + sizeXs = this->GetXsInput()->GetLargestPossibleRegion().GetSize(); + if ((sizePan[0] != sizeXs[0]) || (sizePan[1] != sizeXs[1])) + { + itkExceptionMacro(<< "LmvmPanSharpeningFusionImageFilter: Wrong Pan/Xs size"); + } + + //Process the fusion + m_PanConvolutionFilter->SetInput(this->GetPanInput()); + m_PanConvolutionFilter->SetRadius(this->m_Radius); + m_PanConvolutionFilter->SetFilter(this->m_Filter); + m_PanNoiseFilter->SetRadius(this->m_Radius); + m_PanNoiseFilter->SetInput(this->GetPanInput()); + m_XsConvolutionFilter->SetRadius(this->m_Radius); + m_XsConvolutionFilter->SetFilter(this->m_Filter); + m_XsVectorConvolutionFilter->SetInput(this->GetXsInput()); + m_XsVectorConvolutionFilter->SetFilter(m_XsConvolutionFilter); + m_XsNoiseFilter->SetRadius(this->m_Radius); + m_XsVectorNoiseFilter->SetInput(this->GetXsInput()); + m_XsVectorNoiseFilter->SetFilter(m_XsNoiseFilter); + + + m_FusionFilter->SetInput1(this->GetXsInput()); + m_FusionFilter->SetInput2(m_PanConvolutionFilter->GetOutput()); + m_FusionFilter->SetInput3(m_XsVectorConvolutionFilter->GetOutput()); + m_FusionFilter->SetInput4(m_XsVectorNoiseFilter->GetOutput()); + m_FusionFilter->SetInput5(m_PanNoiseFilter->GetOutput()); + m_FusionFilter->SetInput6(this->GetPanInput()); + + // Wire composite filter + m_FusionFilter->GraftOutput(this->GetOutput()); + m_FusionFilter->Update(); + this->GraftOutput(m_FusionFilter->GetOutput()); +} + +template <class TPanImageType, class TXsImageType, class TOutputImageType, class TInternalPrecision> +void +LmvmPanSharpeningFusionImageFilter +<TPanImageType, TXsImageType, TOutputImageType, TInternalPrecision> +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os + << indent << "Radius:" << this->m_Radius + << std::endl; +} + +} // end namespace otb + +#endif