diff --git a/Anima/filtering/rician_to_gaussian_noise/CMakeLists.txt b/Anima/filtering/rician_to_gaussian_noise/CMakeLists.txt
new file mode 100644
index 000000000..27eb120b4
--- /dev/null
+++ b/Anima/filtering/rician_to_gaussian_noise/CMakeLists.txt
@@ -0,0 +1,37 @@
+if(BUILD_TOOLS AND USE_ARB)
+
+project(animaRicianToGaussianNoise)
+
+## #############################################################################
+## List Sources
+## #############################################################################
+
+list_source_files(${PROJECT_NAME}
+  ${CMAKE_CURRENT_SOURCE_DIR}
+  )
+
+## #############################################################################
+## add executable
+## #############################################################################
+
+add_executable(${PROJECT_NAME}
+  ${${PROJECT_NAME}_CFILES}
+  )
+
+
+## #############################################################################
+## Link
+## #############################################################################
+
+target_link_libraries(${PROJECT_NAME}
+  ${ITKIO_LIBRARIES}
+  AnimaSpecialFunctions
+  )
+
+## #############################################################################
+## install
+## #############################################################################
+
+set_exe_install_rules(${PROJECT_NAME})
+
+endif()
diff --git a/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.h b/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.h
new file mode 100644
index 000000000..a968590d1
--- /dev/null
+++ b/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.h
@@ -0,0 +1,104 @@
+#pragma once
+
+#include <animaMaskedImageToImageFilter.h>
+
+#include <boost/math/distributions/normal.hpp>
+#include <boost/math/distributions/rayleigh.hpp>
+
+namespace anima
+{
+template <class ComponentType,unsigned int ImageDimension>
+class RicianToGaussianImageFilter :
+public anima::MaskedImageToImageFilter<itk::Image<ComponentType,ImageDimension>,itk::Image<ComponentType,ImageDimension> >
+{
+public:
+    /** Standard class typedefs. */
+    typedef RicianToGaussianImageFilter Self;
+    
+    typedef itk::Image<ComponentType,ImageDimension> InputImageType;
+    typedef typename InputImageType::PixelType InputPixelType;
+    
+    typedef itk::Image<ComponentType,ImageDimension> OutputImageType;
+    typedef typename OutputImageType::PixelType OutputPixelType;
+    
+    typedef anima::MaskedImageToImageFilter<InputImageType, OutputImageType> Superclass;
+    
+    typedef itk::SmartPointer<Self> Pointer;
+    typedef itk::SmartPointer<const Self>  ConstPointer;
+    
+    /** Method for creation through the object factory. */
+    itkNewMacro(Self)
+    
+    /** Run-time type information (and related methods) */
+    itkTypeMacro(RicianToGaussianImageFilter, MaskedImageToImageFilter)
+    
+    /** Superclass typedefs. */
+    typedef typename Superclass::InputImageRegionType InputImageRegionType;
+    typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
+    typedef typename Superclass::MaskImageType MaskImageType;
+    typedef typename MaskImageType::SizeType SizeType;
+    
+    typedef typename InputImageType::Pointer InputPointerType;
+    typedef typename OutputImageType::Pointer OutputPointerType;
+    
+    itkSetMacro(MaximumNumberOfIterations, unsigned int)
+    itkSetMacro(Epsilon, double)
+    itkSetMacro(Sigma, double)
+    itkSetMacro(MeanImage, InputPointerType)
+    itkSetMacro(VarianceImage, InputPointerType)
+    
+    itkSetMacro(Scale, double)
+    itkGetConstMacro(Scale, double)
+    
+    itkSetMacro(Alpha, double)
+    itkGetConstMacro(Alpha, double)
+    
+    OutputPointerType GetLocationImage() {return this->GetOutput(0);}
+    OutputPointerType GetScaleImage() {return this->GetOutput(1);}
+    OutputPointerType GetGaussianImage() {return this->GetOutput(2);}
+    
+protected:
+    RicianToGaussianImageFilter()
+    : Superclass()
+    {
+        m_MaximumNumberOfIterations = 100;
+        m_Epsilon = 1.0e-8;
+        m_Sigma = 1.0;
+        m_Scale = 0.0;
+        m_Alpha = 0.005;
+        m_ThreadScaleSamples.clear();
+        m_NeighborWeights.clear();
+        m_Radius.Fill(0);
+        m_MeanImage = NULL;
+        m_VarianceImage = NULL;
+        
+        this->SetNumberOfRequiredOutputs(3);
+        this->SetNthOutput(0, this->MakeOutput(0));
+        this->SetNthOutput(1, this->MakeOutput(1));
+        this->SetNthOutput(2, this->MakeOutput(2));
+    }
+    
+    virtual ~RicianToGaussianImageFilter() {}
+    
+    void BeforeThreadedGenerateData(void) ITK_OVERRIDE;
+    void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+                              itk::ThreadIdType threadId) ITK_OVERRIDE;
+    void AfterThreadedGenerateData(void) ITK_OVERRIDE;
+    
+private:
+    ITK_DISALLOW_COPY_AND_ASSIGN(RicianToGaussianImageFilter);
+    
+    unsigned int m_MaximumNumberOfIterations;
+    double m_Epsilon, m_Sigma, m_Scale, m_Alpha;
+    std::vector<std::vector<double> > m_ThreadScaleSamples;
+    SizeType m_Radius;
+    std::vector<double> m_NeighborWeights;
+    InputPointerType m_MeanImage, m_VarianceImage;
+    
+    boost::math::normal_distribution<> m_NormalDistribution;
+    boost::math::rayleigh_distribution<> m_RayleighDistribution;
+};
+    
+} // end of namespace anima
+
+#include "animaRicianToGaussianImageFilter.hxx"
diff --git a/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.hxx b/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.hxx
new file mode 100644
index 000000000..49ff4d5c3
--- /dev/null
+++ b/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianImageFilter.hxx
@@ -0,0 +1,247 @@
+#pragma once
+#include "animaRicianToGaussianImageFilter.h"
+#include <animaChiDistribution.h>
+#include <animaVectorOperations.h>
+
+#include <itkConstNeighborhoodIterator.h>
+#include <itkImageRegionConstIterator.h>
+#include <itkImageRegionIterator.h>
+#include <itkGaussianOperator.h>
+
+namespace anima
+{
+
+template <class ComponentType,unsigned int ImageDimension>
+void
+RicianToGaussianImageFilter<ComponentType,ImageDimension>
+::BeforeThreadedGenerateData(void)
+{
+    // Compute spatial weights of neighbors beforehand
+    typename InputImageType::SpacingType spacing = this->GetInput()->GetSpacing();
+    typedef itk::GaussianOperator<double,ImageDimension> GaussianOperatorType;
+    std::vector<GaussianOperatorType> gaussianKernels(ImageDimension);
+    
+    for (unsigned int i = 0;i < ImageDimension;++i)
+    {
+        unsigned int reverse_i = ImageDimension - i - 1;
+        double stddev = m_Sigma / spacing[i];
+        gaussianKernels[reverse_i].SetDirection(i);
+        gaussianKernels[reverse_i].SetVariance(stddev * stddev);
+        gaussianKernels[reverse_i].SetMaximumError(1e-3);
+        gaussianKernels[reverse_i].CreateDirectional();
+        gaussianKernels[reverse_i].ScaleCoefficients(1.0e4);
+        m_Radius[i] = gaussianKernels[reverse_i].GetRadius(i);
+    }
+    
+    m_NeighborWeights.clear();
+    for (unsigned int i = 0;i < gaussianKernels[0].Size();++i)
+    {
+        for (unsigned int j = 0;j < gaussianKernels[1].Size();++j)
+        {
+            for (unsigned int k = 0;k < gaussianKernels[2].Size();++k)
+            {
+                double weight = gaussianKernels[0][i] * gaussianKernels[1][j] * gaussianKernels[2][k];
+                m_NeighborWeights.push_back(weight);
+            }
+        }
+    }
+    
+    // Initialize thread containers for global scale estimation
+    unsigned int numThreads = this->GetNumberOfThreads();
+    m_ThreadScaleSamples.resize(numThreads);
+    this->Superclass::BeforeThreadedGenerateData();
+}
+
+template <class ComponentType,unsigned int ImageDimension>
+void
+RicianToGaussianImageFilter<ComponentType,ImageDimension>
+::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
+                       itk::ThreadIdType threadId)
+{
+    typedef itk::ConstNeighborhoodIterator<InputImageType> InputIteratorType;
+    typedef itk::ConstNeighborhoodIterator<MaskImageType> MaskIteratorType;
+    typedef itk::ImageRegionConstIterator<InputImageType> InputSimpleIteratorType;
+    typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
+    
+    InputIteratorType inputItr(m_Radius, this->GetInput(), outputRegionForThread);
+    unsigned int neighborhoodSize = inputItr.Size();
+    
+    MaskIteratorType maskItr(m_Radius, this->GetComputationMask(), outputRegionForThread);
+    
+    InputSimpleIteratorType meanItr, varItr;
+    if (m_MeanImage)
+        meanItr = InputSimpleIteratorType(m_MeanImage, outputRegionForThread);
+    if (m_VarianceImage)
+        varItr = InputSimpleIteratorType(m_VarianceImage, outputRegionForThread);
+    
+    OutputIteratorType locationItr(this->GetOutput(0), outputRegionForThread);
+    OutputIteratorType scaleItr(this->GetOutput(1), outputRegionForThread);
+    OutputIteratorType signalItr(this->GetOutput(2), outputRegionForThread);
+    
+    std::vector<double> samples, weights;
+    bool isInBounds;
+    typename InputImageType::IndexType currentIndex, neighborIndex;
+    typename InputImageType::PointType currentPoint, neighborPoint;
+    
+    while (!maskItr.IsAtEnd())
+    {
+        // Discard voxels outside of the brain
+        if (maskItr.GetCenterPixel() == 0)
+        {
+            locationItr.Set(0);
+            scaleItr.Set(0);
+            signalItr.Set(0);
+            
+            ++inputItr;
+            ++maskItr;
+            if (m_MeanImage)
+                ++meanItr;
+            if (m_VarianceImage)
+                ++varItr;
+            ++locationItr;
+            ++scaleItr;
+            ++signalItr;
+            
+            continue;
+        }
+        
+        // Get signal at current central voxel
+        double inputSignal = inputItr.GetCenterPixel();
+        
+        // Rice-corrupted signals should all be positive
+        if (inputSignal <= 0)
+        {
+            locationItr.Set(0);
+            scaleItr.Set(0);
+            signalItr.Set(0);
+            
+            ++inputItr;
+            ++maskItr;
+            if (m_MeanImage)
+                ++meanItr;
+            if (m_VarianceImage)
+                ++varItr;
+            ++locationItr;
+            ++scaleItr;
+            ++signalItr;
+            
+            continue;
+        }
+        
+        // Estimation of location and scale
+        double location = 0;
+        double scale = m_Scale;
+        
+        if (m_MeanImage && m_VarianceImage)
+        {
+            // If mean and variance images are available, use them instead of neighborhood
+            double sigmaValue = std::sqrt(varItr.Get());
+            double rValue = meanItr.Get() / sigmaValue;
+            double k1Value = 0;
+            double thetaValue = anima::FixedPointFinder(rValue, 1, k1Value);
+            k1Value = anima::KummerFunction(-thetaValue * thetaValue / 2.0, -0.5, 1);
+            
+            scale = sigmaValue / std::sqrt(anima::XiFunction(thetaValue * thetaValue, 1, k1Value, M_PI / 2.0));
+            location = thetaValue * scale;
+        }
+        else
+        {
+            // Use neighbors to create samples
+            samples.clear();
+            weights.clear();
+            
+            for (unsigned int i = 0; i < neighborhoodSize; ++i)
+            {
+                double tmpVal = static_cast<double>(inputItr.GetPixel(i, isInBounds));
+                
+                if (isInBounds && !std::isnan(tmpVal) && std::isfinite(tmpVal))
+                {
+                    if (maskItr.GetPixel(i) != maskItr.GetCenterPixel())
+                        continue;
+                    
+                    double weight = m_NeighborWeights[i];
+                    
+                    if (weight < m_Epsilon)
+                        continue;
+                    
+                    samples.push_back(tmpVal);
+                    weights.push_back(weight);
+                }
+            }
+            
+            if (samples.size() == 1)
+                location = inputSignal;
+            else
+                anima::GetRiceParameters(samples, weights, location, scale);
+        }
+        
+        // Transform Rice signal in Gaussian signal
+        double outputSignal = inputSignal;
+        double snrValue = location / scale;
+        
+        if (!std::isfinite(snrValue) || std::isnan(snrValue))
+        {
+            std::cout << snrValue << " " << location << " " << scale << std::endl;
+            itkExceptionMacro("Estimated SNR is invalid");
+        }
+        
+        if (snrValue <= 0)
+            outputSignal = 0.0;
+        else if (snrValue <= m_Epsilon)
+        {
+            double unifSignal = boost::math::cdf(m_RayleighDistribution, inputSignal / scale);
+            
+            if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha)
+                unifSignal = boost::math::cdf(m_RayleighDistribution, snrValue);
+            
+            outputSignal = scale * boost::math::quantile(m_NormalDistribution, unifSignal);
+        }
+        else if (snrValue <= 600) // if SNR if > 600 keep signal as is, else...
+        {
+            double unifSignal = anima::GetRiceCDF(inputSignal, location, scale);
+            
+            if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha)
+                unifSignal = anima::GetRiceCDF(location, location, scale);
+            
+            outputSignal = location + scale * boost::math::quantile(m_NormalDistribution, unifSignal);
+        }
+        
+        m_ThreadScaleSamples[threadId].push_back(scale);
+        
+        locationItr.Set(static_cast<OutputPixelType>(location));
+        scaleItr.Set(static_cast<OutputPixelType>(scale));
+        signalItr.Set(static_cast<OutputPixelType>(outputSignal));
+        
+        ++inputItr;
+        ++maskItr;
+        if (m_MeanImage)
+            ++meanItr;
+        if (m_VarianceImage)
+            ++varItr;
+        ++locationItr;
+        ++scaleItr;
+        ++signalItr;
+    }
+}
+
+template <class ComponentType,unsigned int ImageDimension>
+void
+RicianToGaussianImageFilter<ComponentType,ImageDimension>
+::AfterThreadedGenerateData(void)
+{
+    if (m_Scale == 0)
+    {
+        std::vector<double> scaleSamples;
+        for (unsigned int i = 0;i < this->GetNumberOfThreads();++i)
+            scaleSamples.insert(scaleSamples.end(), m_ThreadScaleSamples[i].begin(), m_ThreadScaleSamples[i].end());
+        
+        m_Scale = anima::GetMedian(scaleSamples);
+    }
+    
+    m_ThreadScaleSamples.clear();
+    m_NeighborWeights.clear();
+    
+    this->Superclass::AfterThreadedGenerateData();
+}
+    
+} // end of namespace anima
diff --git a/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianNoise.cxx b/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianNoise.cxx
new file mode 100644
index 000000000..f955a49dc
--- /dev/null
+++ b/Anima/filtering/rician_to_gaussian_noise/animaRicianToGaussianNoise.cxx
@@ -0,0 +1,324 @@
+#include <animaRicianToGaussianImageFilter.h>
+#include <animaReadWriteFunctions.h>
+#include <animaRetrieveImageTypeMacros.h>
+
+#include <tclap/CmdLine.h>
+
+//Update progression of the process
+void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData)
+{
+    itk::ProcessObject * processObject = (itk::ProcessObject*) caller;
+    std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"<<std::flush;
+}
+
+struct arguments
+{
+    std::string input, output, mask, mean, variance;
+    double epsilon, sigma;
+    unsigned int maxiters, nthreads;
+};
+
+template <class ComponentType>
+void
+concatenate(const std::vector<typename itk::Image<ComponentType,3>::Pointer> &inputData, const std::string &output)
+{
+    const unsigned int InputDim = 3;
+    typedef itk::Image<ComponentType,InputDim> InputImageType;
+    typedef itk::Image<ComponentType,4> OutputImageType;
+
+    typename OutputImageType::Pointer baseImg;
+    
+    typename OutputImageType::RegionType finalRegion;
+    typename OutputImageType::SizeType finalSize;
+    typename OutputImageType::IndexType finalIndex;
+    typename OutputImageType::PointType finalOrigin;
+    typename OutputImageType::SpacingType finalSpacing;
+    typename OutputImageType::DirectionType finalDirection;
+
+    for (unsigned int d = 0;d < InputDim;++d)
+    {
+        finalIndex[d] = inputData[0]->GetLargestPossibleRegion().GetIndex()[d];
+        finalSize[d] = inputData[0]->GetLargestPossibleRegion().GetSize()[d];
+        finalOrigin[d] = inputData[0]->GetOrigin()[d];
+        finalSpacing[d] = inputData[0]->GetSpacing()[d];
+        for(unsigned int i = 0;i < InputDim;++i)
+            finalDirection[d][i] = inputData[0]->GetDirection()[d][i];
+    }
+    finalIndex[InputDim] = 0;
+    finalSize[InputDim] = 1;
+    finalOrigin[InputDim] = 0;
+    finalSpacing[InputDim] = 1;
+
+    for(unsigned int i = 0; i < InputDim; ++i)
+    {
+        finalDirection[InputDim][i] = 0;
+        finalDirection[i][InputDim] = 0;
+    }
+    finalDirection[InputDim][InputDim] = 1;
+
+    finalRegion.SetIndex(finalIndex);
+    finalRegion.SetSize(finalSize);
+
+    baseImg = OutputImageType::New();
+    baseImg->Initialize();
+    baseImg->SetRegions(finalRegion);
+    baseImg->SetOrigin(finalOrigin);
+    baseImg->SetSpacing(finalSpacing);
+    baseImg->SetDirection(finalDirection);
+    baseImg->SetNumberOfComponentsPerPixel(1);
+    baseImg->Allocate();
+
+    typedef itk::ImageRegionIterator <OutputImageType> FillIteratorType;
+    FillIteratorType fillItr(baseImg, baseImg->GetLargestPossibleRegion());
+    typedef itk::ImageRegionConstIterator<InputImageType> SourceIteratorType;
+    SourceIteratorType srcItr(inputData[0], inputData[0]->GetLargestPossibleRegion());
+    while(!srcItr.IsAtEnd())
+    {
+        fillItr.Set(srcItr.Get());
+        ++fillItr; ++srcItr;
+    }
+
+    // allocate output
+    typename OutputImageType::Pointer outputImg;
+
+    for (unsigned int d = 0; d < InputDim; ++d)
+    {
+        finalIndex[d] = baseImg->GetLargestPossibleRegion().GetIndex()[d];
+        finalSize[d] = baseImg->GetLargestPossibleRegion().GetSize()[d];
+        finalOrigin[d] = baseImg->GetOrigin()[d];
+        finalSpacing[d] = baseImg->GetSpacing()[d];
+
+        for(unsigned int i = 0;i < InputDim;++i)
+            finalDirection[d][i] = baseImg->GetDirection()[d][i];
+    }
+
+    for(unsigned int i = 0;i < InputDim + 1;++i)
+    {
+        finalDirection[InputDim][i] = baseImg->GetDirection()[InputDim][i];
+        finalDirection[i][InputDim] = baseImg->GetDirection()[i][InputDim];
+    }
+
+    finalIndex[InputDim] = baseImg->GetLargestPossibleRegion().GetIndex()[InputDim];
+    finalSize[InputDim] = baseImg->GetLargestPossibleRegion().GetSize()[InputDim] + inputData.size() - 1;
+    finalOrigin[InputDim] = baseImg->GetOrigin()[InputDim];
+    finalSpacing[InputDim] = baseImg->GetSpacing()[InputDim];
+
+    finalRegion.SetIndex(finalIndex);
+    finalRegion.SetSize(finalSize);
+
+    outputImg = OutputImageType::New();
+    outputImg->Initialize();
+    outputImg->SetRegions(finalRegion);
+    outputImg->SetOrigin(finalOrigin);
+    outputImg->SetSpacing(finalSpacing);
+    outputImg->SetDirection(finalDirection);
+    outputImg->SetNumberOfComponentsPerPixel(1);
+    outputImg->Allocate();
+
+    // fill with the base:
+    typedef itk::ImageRegionIterator <OutputImageType> FillIteratorType;
+    fillItr = FillIteratorType(outputImg, outputImg->GetLargestPossibleRegion());
+    typedef itk::ImageRegionConstIterator<OutputImageType> BaseIteratorType;
+    BaseIteratorType baseItr(baseImg, baseImg->GetLargestPossibleRegion());
+    while(!baseItr.IsAtEnd())
+    {
+        fillItr.Set(baseItr.Get());
+        ++fillItr; ++baseItr;
+    }
+
+    //fill with the inputs:
+    for(unsigned int i = 1; i < inputData.size();++i)
+    {
+        typedef itk::ImageRegionConstIterator<InputImageType> InputIteratorType;
+        InputIteratorType inputItr(inputData[i], inputData[i]->GetLargestPossibleRegion());
+        while(!inputItr.IsAtEnd())
+        {
+            fillItr.Set(inputItr.Get());
+            ++fillItr; ++inputItr;
+        }
+    }
+
+    // write res
+    anima::writeImage<OutputImageType>(output, outputImg);
+}
+
+template <class ComponentType, unsigned int ImageDimension>
+typename itk::Image<ComponentType,ImageDimension>::Pointer
+performTransformation(const typename itk::Image<ComponentType,ImageDimension>::Pointer inputImage,
+                      const typename itk::Image<ComponentType,ImageDimension>::Pointer meanImage,
+                      const typename itk::Image<ComponentType,ImageDimension>::Pointer varImage,
+                      const arguments &args)
+{
+    typedef anima::RicianToGaussianImageFilter<ComponentType,ImageDimension> RiceToGaussianImageFilterType;
+    typedef typename RiceToGaussianImageFilterType::InputImageType InputImageType;
+    typedef typename RiceToGaussianImageFilterType::OutputImageType OutputImageType;
+    typedef typename RiceToGaussianImageFilterType::MaskImageType MaskImageType;
+    
+    typename MaskImageType::Pointer maskImage;
+    if (args.mask != "")
+        maskImage = anima::readImage<MaskImageType>(args.mask);
+    
+    itk::CStyleCommand::Pointer callback;
+    typename RiceToGaussianImageFilterType::Pointer mainFilter;
+    
+    // Estimate scale parameter using SNR fixed point strategy
+    mainFilter = RiceToGaussianImageFilterType::New();
+    mainFilter->SetInput(inputImage);
+    mainFilter->SetMaximumNumberOfIterations(args.maxiters);
+    mainFilter->SetEpsilon(args.epsilon);
+    mainFilter->SetSigma(args.sigma);
+    mainFilter->SetNumberOfThreads(args.nthreads);
+    
+    if (maskImage)
+        mainFilter->SetComputationMask(maskImage);
+
+    if (meanImage)
+        mainFilter->SetMeanImage(meanImage);
+
+    if (varImage)
+        mainFilter->SetVarianceImage(varImage);
+    
+    callback = itk::CStyleCommand::New();
+    callback->SetCallback(eventCallback);
+    mainFilter->AddObserver(itk::ProgressEvent(), callback);
+    
+    std::cout << "First run for estimating the scale parameter..." << std::endl;
+
+    mainFilter->Update();
+
+    double scale = mainFilter->GetScale();
+    
+    std::cout << "\nDone. Estimated scale parameter: " << scale << std::endl;
+    
+    // Now, transform Rician noise to Gaussian noise using fixed scale parameter
+    mainFilter = RiceToGaussianImageFilterType::New();
+    mainFilter->SetInput(inputImage);
+    mainFilter->SetMaximumNumberOfIterations(args.maxiters);
+    mainFilter->SetEpsilon(args.epsilon);
+    mainFilter->SetSigma(args.sigma);
+    mainFilter->SetScale(scale);
+    mainFilter->SetNumberOfThreads(args.nthreads);
+    
+    if (maskImage)
+        mainFilter->SetComputationMask(maskImage);
+    
+    callback = itk::CStyleCommand::New();
+    callback->SetCallback(eventCallback);
+    mainFilter->AddObserver(itk::ProgressEvent(), callback);
+    
+    mainFilter->Update();
+    
+    std::cout << "\n" << std::endl;
+
+    return mainFilter->GetGaussianImage();
+}
+
+template <class ComponentType, unsigned int ImageDimension>
+void
+transformNoise(itk::ImageIOBase::Pointer imageIO, const arguments &args)
+{
+    if (ImageDimension == 4)
+    {
+        typedef itk::Image<ComponentType,4> Image4DType;
+        typedef itk::Image<ComponentType,3> Image3DType;
+        typename Image4DType::Pointer input = anima::readImage<Image4DType>(args.input);
+        std::vector <typename Image3DType::Pointer> inputData = anima::getImagesFromHigherDimensionImage<Image4DType,Image3DType>(input);
+        unsigned int nbImages = inputData.size();
+        std::vector <typename Image3DType::Pointer> meanData(nbImages), varData(nbImages);
+        if (args.mean != "")
+            meanData = anima::getImagesFromHigherDimensionImage<Image4DType,Image3DType>(anima::readImage<Image4DType>(args.mean));
+        if (args.variance != "")
+            varData = anima::getImagesFromHigherDimensionImage<Image4DType,Image3DType>(anima::readImage<Image4DType>(args.variance));
+        
+        for (unsigned int i = 0;i < nbImages;++i)
+            inputData[i] = performTransformation<ComponentType,3>(inputData[i], meanData[i], varData[i], args);
+        concatenate<ComponentType>(inputData,args.output);
+        return;
+    }
+
+    typedef itk::Image<ComponentType,ImageDimension> InputImageType;
+    typename InputImageType::Pointer inputImage = anima::readImage<InputImageType>(args.input);
+    typename InputImageType::Pointer meanImage, varImage;
+
+    if (args.mean != "")
+        meanImage = anima::readImage<InputImageType>(args.mean);
+
+    if (args.variance != "")
+        varImage = anima::readImage<InputImageType>(args.variance);
+
+    inputImage = performTransformation<ComponentType,ImageDimension>(inputImage,meanImage,varImage,args);
+    
+    anima::writeImage<InputImageType>(args.output, inputImage);
+}
+
+template <class ComponentType>
+void
+retrieveNbDimensions(itk::ImageIOBase::Pointer imageIO, const arguments &args)
+{
+    ANIMA_RETRIEVE_NUMBER_OF_DIMENSIONS(imageIO, ComponentType, transformNoise, imageIO, args)
+}
+
+int main(int ac, const char** av)
+{
+    
+    TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS Team", ' ',ANIMA_VERSION);
+    
+    TCLAP::ValueArg<std::string> inputArg("i", "input", "Input Rice-corrupted image", true, "", "Input Rice-corrupted image", cmd);
+    TCLAP::ValueArg<std::string> outputArg("o", "output", "Output Gaussian-corrupted image", true, "", "Output Gaussian-corrupted image", cmd);
+    
+    TCLAP::ValueArg<std::string> maskArg("m", "mask", "Segmentation mask", false, "", "Segmentation mask", cmd);
+    TCLAP::ValueArg<std::string> meanArg("", "mean", "Mean image", false, "", "Mean image", cmd);
+    TCLAP::ValueArg<std::string> varArg("", "variance", "Variance image", false, "", "Variance image", cmd);
+    
+    TCLAP::ValueArg<double> epsArg("E", "epsilon", "Minimal absolute value difference betweem fixed point iterations (default: 1e-8)", false, 1.0e-8, "Minimal absolute value difference betweem fixed point iterations", cmd);
+    TCLAP::ValueArg<double> sigmaArg("S", "sigma", "Gaussian standard deviation for defining neighbor weights (default: 1.0)", false, 1.0, "Gaussian standard deviation for defining neighbor weights", cmd);
+    
+    TCLAP::ValueArg<unsigned int> maxiterArg("I", "maxiter", "Maximum number of iterations (default: 100)", false, 100, "Maximum number of iterations", cmd);
+    TCLAP::ValueArg<unsigned int> nbpArg("T", "nbp", "Number of threads to run on -> default : automatically determine", false, itk::MultiThreader::GetGlobalDefaultNumberOfThreads(), "Number of threads", cmd);
+    
+    try
+    {
+        cmd.parse(ac,av);
+    }
+    catch (TCLAP::ArgException& e)
+    {
+        std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    // Retrieve image info
+    itk::ImageIOBase::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(inputArg.getValue().c_str(),
+                                                                           itk::ImageIOFactory::ReadMode);
+    if (!imageIO)
+    {
+        std::cerr << "Itk could not find suitable IO factory for the input" << std::endl;
+        return EXIT_FAILURE;
+    }
+    
+    // Now that we found the appropriate ImageIO class, ask it to read the meta data from the image file.
+    imageIO->SetFileName(inputArg.getValue());
+    imageIO->ReadImageInformation();
+    
+    arguments args;
+    args.input = inputArg.getValue();
+    args.output = outputArg.getValue();
+    args.mask = maskArg.getValue();
+    args.mean = meanArg.getValue();
+    args.variance = varArg.getValue();
+    args.epsilon = epsArg.getValue();
+    args.sigma = sigmaArg.getValue();
+    args.maxiters = maxiterArg.getValue();
+    args.nthreads = nbpArg.getValue();
+    
+    try
+    {
+        ANIMA_RETRIEVE_COMPONENT_TYPE(imageIO, retrieveNbDimensions, imageIO, args);
+    }
+    catch (itk::ExceptionObject &err)
+    {
+        std::cerr << err << std::endl;
+        return EXIT_FAILURE;
+    }
+    
+    return EXIT_SUCCESS;
+}