diff --git a/Modules/Applications/AppFiltering/app/otbFastNLMeans.cxx b/Modules/Applications/AppFiltering/app/otbFastNLMeans.cxx index 2463dee25c365721d462d8a71539361dcfd032fa..3ed27ed028bfe43d0e1c27e626c9ab6d155125a3 100644 --- a/Modules/Applications/AppFiltering/app/otbFastNLMeans.cxx +++ b/Modules/Applications/AppFiltering/app/otbFastNLMeans.cxx @@ -28,105 +28,110 @@ namespace otb namespace Wrapper { - class FastNLMeans : public Application +class FastNLMeans : public Application +{ +public: + typedef FastNLMeans Self; + typedef Application Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + itkNewMacro(Self); + itkTypeMacro(FastNLMeans, otb::Wrapper::Application); + + // Define image types + typedef float PixelType; + typedef FloatImageType ImageType; + + // Define filter + typedef NLMeansFilter<ImageType, ImageType> NLMeansFilterType; + +private: + void DoInit() override + { + SetName("FastNLMeans"); + SetDescription("Apply NL Means filter to an image."); + + SetDocLongDescription("Implementation is an approximation of NL Means, which is faster."); + + // Optional descriptors + SetDocLimitations( + "This filter relies on integral images. Overflow may happen though the risk is limited " + " by OTB mechanism which process data by chunks."); + SetDocAuthors("OTB-Team"); + SetDocSeeAlso(" "); + AddDocTag(Tags::Filter); + + // Parameter declarations + AddParameter(ParameterType_InputImage, "in", "Input image"); + SetParameterDescription("in", "Input image to denoise"); + + AddParameter(ParameterType_OutputImage, "out", "Output Image"); + SetParameterDescription("out", "Output image."); + + AddParameter(ParameterType_Int, "patchradius", "Patch radius (patch is a square)"); + SetParameterDescription("patchradius", "Full patch will have a size of 2*patchradius +1."); + SetDefaultParameterInt("patchradius", 2); + SetMinimumParameterIntValue("patchradius", 0); + MandatoryOff("patchradius"); + + AddParameter(ParameterType_Int, "searchradius", "Search window radius (search window is a square)"); + SetParameterDescription("searchradius", "Search window is used to find similar patches. Its size will be 2*searchradius+1."); + SetDefaultParameterInt("searchradius", 7); + SetMinimumParameterIntValue("searchradius", 0); + MandatoryOff("searchradius"); + + AddParameter(ParameterType_Float, "sig", "Standard deviation in image"); + SetParameterDescription("sig", + "Noise standard deviation estimated in image. This parameter is used to correct for the expected difference between two patches. " + "This filter works fine without using this tuning."); + SetDefaultParameterFloat("sig", 0); + SetMinimumParameterFloatValue("sig", 0); + MandatoryOff("sig"); + + AddParameter(ParameterType_Float, "thresh", "Similarity threshold"); + SetParameterDescription("thresh", + "Factor influencing similarity score of two patches. The higher the threshold, the more permissive the filter. It is common to set " + "this threshold slightly below the standard deviation (for Gaussian noise), at about 0.8*sigma."); + SetDefaultParameterFloat("thresh", 1.0); + SetMinimumParameterFloatValue("thresh", 0.); + MandatoryOff("thresh"); + + AddRAMParameter(); + + SetDocExampleParameterValue("in", "GomaAvant.tif"); + SetDocExampleParameterValue("out", "denoisedImage_NLMeans.tif"); + } + + void DoUpdateParameters() override + { + // Nothing to do here : all parameters are independent + } + + void DoExecute() override { - public: - typedef FastNLMeans Self; - typedef Application Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - itkNewMacro(Self); - itkTypeMacro(FastNLMeans, otb::Wrapper::Application); - - // Define image types - typedef float PixelType; - typedef FloatImageType ImageType; - - // Define filter - typedef NLMeansFilter<ImageType, ImageType> NLMeansFilterType; - - private: - void DoInit() override - { - SetName("FastNLMeans"); - SetDescription("Apply NL Means filter to an image."); - - SetDocLongDescription("Implementation is an approximation of NL Means, which is faster."); - - //Optional descriptors - SetDocLimitations("This filter relies on integral images. Overflow may happen though the risk is limited " - " by OTB mechanism which process data by chunks."); - SetDocAuthors("OTB-Team"); - SetDocSeeAlso(" "); - AddDocTag(Tags::Filter); - - //Parameter declarations - AddParameter(ParameterType_InputImage, "in", "Input image"); - SetParameterDescription("in", "Input image to denoise"); - - AddParameter(ParameterType_OutputImage, "out", "Output Image"); - SetParameterDescription("out","Output image."); - - AddParameter(ParameterType_Int, "patchradius", "Patch radius (patch is a square)"); - SetParameterDescription("patchradius", "Full patch will have a size of 2*patchradius +1."); - SetDefaultParameterInt("patchradius", 2); - SetMinimumParameterIntValue("patchradius", 0); - MandatoryOff("patchradius"); - - AddParameter(ParameterType_Int, "searchradius", "Search window radius (search window is a square)"); - SetParameterDescription("searchradius", "Search window is used to find similar patches. Its size will be 2*searchradius+1."); - SetDefaultParameterInt("searchradius", 7); - SetMinimumParameterIntValue("searchradius", 0); - MandatoryOff("searchradius"); - - AddParameter(ParameterType_Float, "sig", "Standard deviation in image"); - SetParameterDescription("sig", "Noise standard deviation estimated in image. This parameter is used to correct for the expected difference between two patches. This filter works fine without using this tuning."); - SetDefaultParameterFloat("sig", 0); - SetMinimumParameterFloatValue("sig", 0); - MandatoryOff("sig"); - - AddParameter(ParameterType_Float, "thresh", "Similarity threshold"); - SetParameterDescription("thresh", "Factor influencing similarity score of two patches. The higher the threshold, the more permissive the filter. It is common to set this threshold slightly below the standard deviation (for Gaussian noise), at about 0.8*sigma."); - SetDefaultParameterFloat("thresh", 1.0); - SetMinimumParameterFloatValue("thresh", 0.); - MandatoryOff("thresh"); - - AddRAMParameter(); - - SetDocExampleParameterValue("in", "GomaAvant.tif"); - SetDocExampleParameterValue("out","denoisedImage_NLMeans.tif"); - } - - void DoUpdateParameters() override - { - // Nothing to do here : all parameters are independent - } - - void DoExecute() override - { - // Get the input image - ImageType::Pointer imIn = this->GetParameterFloatImage("in"); - float sigma = this->GetParameterFloat("sig"); - float cutoffDistance = this->GetParameterFloat("thresh"); - int halfPatchSize = this->GetParameterInt("patchradius"); - int halfSearchSize = this->GetParameterInt("searchradius"); - NLMeansFilterType::Pointer nlMeansFilter = NLMeansFilterType::New(); - nlMeansFilter->SetInput(imIn); - nlMeansFilter->SetSigma(sigma); - nlMeansFilter->SetHalfWindowSize(halfPatchSize); - nlMeansFilter->SetHalfSearchSize(halfSearchSize); - nlMeansFilter->SetCutOffDistance(cutoffDistance); - - m_FilterRef = nlMeansFilter; - SetParameterOutputImage("out", nlMeansFilter->GetOutput()); - RegisterPipeline(); - } - - itk::LightObject::Pointer m_FilterRef; - }; // end class - -} -} + // Get the input image + ImageType::Pointer imIn = this->GetParameterFloatImage("in"); + float sigma = this->GetParameterFloat("sig"); + float cutoffDistance = this->GetParameterFloat("thresh"); + int halfPatchSize = this->GetParameterInt("patchradius"); + int halfSearchSize = this->GetParameterInt("searchradius"); + NLMeansFilterType::Pointer nlMeansFilter = NLMeansFilterType::New(); + nlMeansFilter->SetInput(imIn); + nlMeansFilter->SetSigma(sigma); + nlMeansFilter->SetHalfWindowSize(halfPatchSize); + nlMeansFilter->SetHalfSearchSize(halfSearchSize); + nlMeansFilter->SetCutOffDistance(cutoffDistance); + + m_FilterRef = nlMeansFilter; + SetParameterOutputImage("out", nlMeansFilter->GetOutput()); + RegisterPipeline(); + } + + itk::LightObject::Pointer m_FilterRef; +}; // end class + +} // namespace Wrapper +} // namespace otb OTB_APPLICATION_EXPORT(otb::Wrapper::FastNLMeans) diff --git a/Modules/Filtering/Smoothing/include/otbFastNLMeansImageFilter.h b/Modules/Filtering/Smoothing/include/otbFastNLMeansImageFilter.h index e6ac95fb585352d6cbc899d98496d926c3cd956b..a3c33a68112ab9784d1be27087826ffb6ef4e432 100644 --- a/Modules/Filtering/Smoothing/include/otbFastNLMeansImageFilter.h +++ b/Modules/Filtering/Smoothing/include/otbFastNLMeansImageFilter.h @@ -23,144 +23,138 @@ #include "itkImageToImageFilter.h" -namespace otb{ - - /** \class NLMeansFilter - * \brief This class implements a fast approximated version of NL Means denoising - * algorithm. This implementation is based on code in scikit module skimage. - * This fast version of the NL Means algorithm - * has been described in the following papers : - * - * J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, - * Fast nonlocal filtering applied to electron cryomicroscopy, - * in 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, - * 2008, pp. 1331-1334. - * - * Jacques Froment. - * Parameter-Free Fast Pixelwise Non-Local Means Denoising. - * Image Processing On Line, 2014, vol. 4, p. 300-326. - * - * \ingroup OTBSmoothing - */ +namespace otb +{ + +/** \class NLMeansFilter + * \brief This class implements a fast approximated version of NL Means denoising + * algorithm. This implementation is based on code in scikit module skimage. + * This fast version of the NL Means algorithm + * has been described in the following papers : + * + * J. Darbon, A. Cunha, T.F. Chan, S. Osher, and G.J. Jensen, + * Fast nonlocal filtering applied to electron cryomicroscopy, + * in 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, + * 2008, pp. 1331-1334. + * + * Jacques Froment. + * Parameter-Free Fast Pixelwise Non-Local Means Denoising. + * Image Processing On Line, 2014, vol. 4, p. 300-326. + * + * \ingroup OTBSmoothing + */ + +template <class TInputImage, class TOutputImage> +class ITK_EXPORT NLMeansFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage> +{ +public: + /** Standard class typedefs */ + typedef NLMeansFilter Self; + typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Typedef to image types */ + typedef TInputImage InImageType; + typedef TOutputImage OutImageType; + + /** Typedef to describe the image pointer types. */ + typedef typename InImageType::Pointer InImagePointerType; + typedef typename InImageType::ConstPointer InImageConstPointerType; + typedef typename InImageType::RegionType InRegionType; + typedef typename InImageType::IndexType InIndexType; + typedef typename InImageType::SizeType InSizeType; + typedef typename InImageType::OffsetType InOffsetType; + typedef typename OutImageType::Pointer OutImagePointerType; + typedef typename OutImageType::RegionType OutRegionType; + typedef typename OutImageType::PixelType OutPixelType; + typedef typename OutImageType::SizeType OutSizeType; + typedef typename OutImageType::IndexType OutIndexType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + itkTypeMacro(NLMeansFilter, ImageToImageFilter); + + void SetSigma(const float sigma) + { + m_Var = 2.0 * sigma * sigma; + } + + void SetHalfWindowSize(const unsigned int hws) + { + m_HalfPatchSize.Fill(hws); + // Update NormalizeDistance + m_NormalizeDistance = (2 * hws + 1) * (2 * hws + 1) * m_CutoffDistance * m_CutoffDistance; + } - template <class TInputImage, class TOutputImage> - class ITK_EXPORT NLMeansFilter : - public itk::ImageToImageFilter<TInputImage, TOutputImage> + void SetHalfSearchSize(const unsigned int hss) + { + m_HalfSearchSize.Fill(hss); + } + void SetCutOffDistance(const float thresh) { - public: - /** Standard class typedefs */ - typedef NLMeansFilter Self; - typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - /** Typedef to image types */ - typedef TInputImage InImageType; - typedef TOutputImage OutImageType; - - /** Typedef to describe the image pointer types. */ - typedef typename InImageType::Pointer InImagePointerType; - typedef typename InImageType::ConstPointer InImageConstPointerType; - typedef typename InImageType::RegionType InRegionType; - typedef typename InImageType::IndexType InIndexType; - typedef typename InImageType::SizeType InSizeType; - typedef typename InImageType::OffsetType InOffsetType; - typedef typename OutImageType::Pointer OutImagePointerType; - typedef typename OutImageType::RegionType OutRegionType; - typedef typename OutImageType::PixelType OutPixelType; - typedef typename OutImageType::SizeType OutSizeType; - typedef typename OutImageType::IndexType OutIndexType; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - itkTypeMacro(NLMeansFilter, ImageToImageFilter); - - void SetSigma(const float sigma) - { - m_Var = 2.0 * sigma * sigma; - } - - void SetHalfWindowSize(const unsigned int hws) - { - m_HalfPatchSize.Fill(hws); - // Update NormalizeDistance - m_NormalizeDistance = (2*hws+1)*(2*hws+1) * - m_CutoffDistance*m_CutoffDistance; - } - - void SetHalfSearchSize(const unsigned int hss) - { - m_HalfSearchSize.Fill(hss); - } - void SetCutOffDistance(const float thresh) - { - m_CutoffDistance = thresh; - // Update NormalizeDistance - m_NormalizeDistance = (2*m_HalfPatchSize[0]+1)* - (2*m_HalfPatchSize[1]+1) * - m_CutoffDistance*m_CutoffDistance; - } - - protected: - /** Constructor */ - NLMeansFilter(); - /** Destructor */ - ~NLMeansFilter() override = default; - - void ThreadedGenerateData(const OutRegionType& outputRegionForThread, - itk::ThreadIdType itkNotUsed(threadId)) override; - - void GenerateInputRequestedRegion() override; - - /** Compute the requested input region, given an output region. - * If the input requested region is outside the largest input region, a mirror padding - * is necessary. The returned tuple is composed of the following paramters : - * * input requested region (always lie inside the largest input region) - * * top rows, left cols, bottom rows, right cols : numbers of rows/cols to add with a mirror padding - * * boolean : if true, a mirror padding (in at least one direction) has to be computed - */ - std::tuple<InRegionType, int, int, int, int, bool> OutputRegionToInputRegion - (const OutRegionType& outputRegion) const; - - /** PrintSelf method */ - void PrintSelf(std::ostream& os, itk::Indent indent) const override; - - private: - NLMeansFilter(const Self&) = delete; //purposely not implemented - NLMeansFilter& operator=(const Self&) = delete; //purposely not implemented - - /** For a given shift in rows and cols, this function computes - * the squared difference between the image and its shifted version. - * Results are added to form an integral image. - */ - void ComputeIntegralImage( - const std::vector<double> & dataInput, /**< input data stored in a vector */ - std::vector<double> &imIntegral, /**< output parameter. Contains the integral image of squared difference */ - const OutIndexType shift, /**< Shift (dcol, drow) to apply to compute the difference */ - const InSizeType sizeIntegral, /**< Integral image size */ - const InSizeType sizeInput /**< input data image size */ - ) const; - - /** This function computes the squared euclidean distance - * between a patch and its shifted version. - * Computation relies on the integral image obtained before. - */ - OutPixelType ComputeDistance(const unsigned int row, /**< Upper left corner row coordinate of patch*/ - const unsigned int col, /**< Upper left corner col coordinate of patch*/ - const std::vector<double>& imIntegral, /**< Integral image of squared difference*/ - const unsigned int nbCols /**< Integral image number of columns */ - ) const; - - // Define class attributes - InSizeType m_HalfSearchSize; - InSizeType m_HalfPatchSize; - float m_Var; - float m_CutoffDistance; - float m_NormalizeDistance; //cutoff**2 * windowSize**2 - - static const int m_ROW=1; - static const int m_COL=0; - }; + m_CutoffDistance = thresh; + // Update NormalizeDistance + m_NormalizeDistance = (2 * m_HalfPatchSize[0] + 1) * (2 * m_HalfPatchSize[1] + 1) * m_CutoffDistance * m_CutoffDistance; + } + +protected: + /** Constructor */ + NLMeansFilter(); + /** Destructor */ + ~NLMeansFilter() override = default; + + void ThreadedGenerateData(const OutRegionType& outputRegionForThread, itk::ThreadIdType itkNotUsed(threadId)) override; + + void GenerateInputRequestedRegion() override; + + /** Compute the requested input region, given an output region. + * If the input requested region is outside the largest input region, a mirror padding + * is necessary. The returned tuple is composed of the following paramters : + * * input requested region (always lie inside the largest input region) + * * top rows, left cols, bottom rows, right cols : numbers of rows/cols to add with a mirror padding + * * boolean : if true, a mirror padding (in at least one direction) has to be computed + */ + std::tuple<InRegionType, int, int, int, int, bool> OutputRegionToInputRegion(const OutRegionType& outputRegion) const; + + /** PrintSelf method */ + void PrintSelf(std::ostream& os, itk::Indent indent) const override; + +private: + NLMeansFilter(const Self&) = delete; // purposely not implemented + NLMeansFilter& operator=(const Self&) = delete; // purposely not implemented + + /** For a given shift in rows and cols, this function computes + * the squared difference between the image and its shifted version. + * Results are added to form an integral image. + */ + void ComputeIntegralImage(const std::vector<double>& dataInput, /**< input data stored in a vector */ + std::vector<double>& imIntegral, /**< output parameter. Contains the integral image of squared difference */ + const OutIndexType shift, /**< Shift (dcol, drow) to apply to compute the difference */ + const InSizeType sizeIntegral, /**< Integral image size */ + const InSizeType sizeInput /**< input data image size */ + ) const; + + /** This function computes the squared euclidean distance + * between a patch and its shifted version. + * Computation relies on the integral image obtained before. + */ + OutPixelType ComputeDistance(const unsigned int row, /**< Upper left corner row coordinate of patch*/ + const unsigned int col, /**< Upper left corner col coordinate of patch*/ + const std::vector<double>& imIntegral, /**< Integral image of squared difference*/ + const unsigned int nbCols /**< Integral image number of columns */ + ) const; + + // Define class attributes + InSizeType m_HalfSearchSize; + InSizeType m_HalfPatchSize; + float m_Var; + float m_CutoffDistance; + float m_NormalizeDistance; // cutoff**2 * windowSize**2 + + static const int m_ROW = 1; + static const int m_COL = 0; +}; } // end namespace otb #ifndef OTB_MANUAL_INSTANTIATION diff --git a/Modules/Filtering/Smoothing/test/otbFastNLMeansImageFilter.cxx b/Modules/Filtering/Smoothing/test/otbFastNLMeansImageFilter.cxx index f88c1f91fac3a4c541bdc663a705a4a61a6ebf26..659b76528f7f3fd1b503477f7aa1e3f583058727 100644 --- a/Modules/Filtering/Smoothing/test/otbFastNLMeansImageFilter.cxx +++ b/Modules/Filtering/Smoothing/test/otbFastNLMeansImageFilter.cxx @@ -32,12 +32,12 @@ int otbFastNLMeansImageFilter(int itkNotUsed(argc), char* argv[]) const char* inputFilename = argv[1]; const char* outputFilename = argv[2]; - int HalfPatchSize = atoi(argv[3]); - int HalfSearchSize = atoi(argv[4]); - float Thresh = atof(argv[5]); - - typedef float InputPixelType; - typedef float OutputPixelType; + int HalfPatchSize = atoi(argv[3]); + int HalfSearchSize = atoi(argv[4]); + float Thresh = atof(argv[5]); + + typedef float InputPixelType; + typedef float OutputPixelType; typedef otb::Image<InputPixelType> InputImageType; typedef otb::Image<OutputPixelType> OutputImageType;