Commit 62529436 authored by Guillaume Pasero's avatar Guillaume Pasero

Merge branch 'fast_nlmeans_filter' into 'develop'

Fast nlmeans filter

See merge request !644
parents 0d4fae35 b47f546e
Pipeline #3208 passed with stages
in 21 minutes and 6 seconds
......@@ -27,3 +27,8 @@ otb_create_application(
NAME ContrastEnhancement
SOURCES otbContrastEnhancement.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
otb_create_application(
NAME FastNLMeans
SOURCES otbFastNLMeans.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
/*
* Copyright (C) 2005-2018 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 "otbFastNLMeansImageFilter.h"
namespace otb
{
namespace Wrapper
{
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
{
// 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)
......@@ -30,6 +30,7 @@ otb_module(OTBAppFiltering
OTBStatistics
OTBStreaming
OTBFunctor
OTBSmoothing
TEST_DEPENDS
OTBTestKernel
......
......@@ -101,6 +101,19 @@ otb_test_application(NAME apTvUtContrastTest_lum
VALID --compare-image ${NOTOL}
${BASELINE}/apTvUtContrastTest_lum.tif
${TEMP}/apTvUtContrastTest_lum.tif)
#----------- NL Means TESTS ----------------
otb_test_application(NAME nlMeansTest_base
APP FastNLMeans
OPTIONS -in ${INPUTDATA}/GomaAvant.tif
-out ${TEMP}/GomaAvant_NLMeans.tif
-patchradius 2
-searchradius 11
-thresh 30
VALID --compare-image ${EPSILON_7}
${BASELINE}/GomaAvant_NLMeans.tif
${TEMP}/GomaAvant_NLMeans.tif)
......
/*
* Copyright (C) 2005-2018 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 otbFastNLMeansImageFilter_h
#define otbFastNLMeansImageFilter_h
#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
*/
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;
}
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;
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbFastNLMeansImageFilter.hxx"
#endif
#endif
......@@ -25,6 +25,7 @@ otbSmoothingTestDriver.cxx
otbMeanShiftSmoothingImageFilter.cxx
otbMeanShiftSmoothingImageFilterSpatialStability.cxx
otbMeanShiftSmoothingImageFilterThreading.cxx
otbFastNLMeansImageFilter.cxx
)
add_executable(otbSmoothingTestDriver ${OTBSmoothingTests})
......@@ -148,3 +149,10 @@ otb_add_test(NAME bfTvMeanShiftSmoothingImageFilterThreadingNonOpt COMMAND otbSm
${TEMP}/bfMeanShiftSmoothingImageFilterMultiThread.tif
4 50 0
)
otb_add_test(NAME fastNLMeansImageFilter COMMAND otbSmoothingTestDriver
otbFastNLMeansImageFilter
${INPUTDATA}/GomaAvant.tif
${TEMP}/GomaAvant_FastNLMeansFilter.tif
2 11 30
)
/*
* 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 "itkMacro.h"
#include <iostream>
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otbFastNLMeansImageFilter.h"
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;
typedef otb::Image<InputPixelType> InputImageType;
typedef otb::Image<OutputPixelType> OutputImageType;
typedef otb::ImageFileReader<InputImageType> ReaderType;
typedef otb::ImageFileWriter<OutputImageType> WriterType;
typedef otb::NLMeansFilter<InputImageType, OutputImageType> FilterType;
FilterType::Pointer filterNLMeans = FilterType::New();
filterNLMeans->SetHalfWindowSize(HalfPatchSize);
filterNLMeans->SetHalfSearchSize(HalfSearchSize);
filterNLMeans->SetCutOffDistance(Thresh);
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
reader->SetFileName(inputFilename);
writer->SetFileName(outputFilename);
filterNLMeans->SetInput(reader->GetOutput());
writer->SetInput(filterNLMeans->GetOutput());
writer->Update();
return EXIT_SUCCESS;
}
......@@ -25,4 +25,5 @@ void RegisterTests()
REGISTER_TEST(otbMeanShiftSmoothingImageFilter);
REGISTER_TEST(otbMeanShiftSmoothingImageFilterSpatialStability);
REGISTER_TEST(otbMeanShiftSmoothingImageFilterThreading);
REGISTER_TEST(otbFastNLMeansImageFilter);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment