Release

Example.m

+%% Author: Nicolas Gasnier
+%% Copyright 2020 CS GROUP - France, Centre national d'études spatiales, Télécom Paris 
+%% nicolas.gasnier@telecom-paris.fr
+
+%% Adapted from the original RABASAR code written by Charles Deledalle (https://bitbucket.org/cdeledalle/rabasar)
+
+
+
+%% Example of RABASAR modified to use a geometric mean superimage on a stack of 20 Sentinel-1 SLC from St Nazaire (France)
+clear all
+close all
+set(0, 'DefaultAxesFontSize',30)
+addpath(genpath('.'))
+
+disp('Please note that this implementation uses DDID denoiser and not BM2D for copyright reasons. The results presented in the article have been obtained with BM3D denoiser. See README.md for more information')
+% Load SAR time series
+% v is a crop of a S1 SLC time series spatially undersampled by a factor 2
+% to limit noise correlation
+
+verbose = false; % Set to true to display more results
+
+
+load('Example_S1_20_images_undersampled.mat')
+
+v = v.*v; % Amplitude to intensity
+
+
+t = 11;
+L = 1;
+
+% Rabasar
+h = robustwaitbar(0);
+tic
+[u_hat_t, u_hat_si, L_hat_si, u_hat_dsi, L_hat_dsi] = ...
+    rabasar_geom(v, L, t, ...
+            'waitbar', @(p) robustwaitbar(p, h));
+toc
+close(h);
+u_hat_t_g = u_hat_t;
+u_hat_si_g = u_hat_si;
+L_hat_si_g = L_hat_si;
+u_hat_dsi_g = u_hat_dsi;
+L_hat_dsi_g = L_hat_dsi;
+
+%% Visualization
+min1 = 0.1;
+max1 = 10;
+fancyfigure;
+h1 = plotimagesar(u_hat_si_g, 'alpha', 2/3);
+title(sprintf('Noisy Geometric mean'));
+fancyfigure;
+plotimagesar(u_hat_dsi_g, 'rangeof', h1);
+title(sprintf('Denoised geometric mean'));
+fancyfigure;
+plotimagesar(v(:,:,t), 'rangeof', h1);
+title(sprintf('Noisy image at date: t=%d', t));
+fancyfigure;
+plotimagesar(v(:,:,17), 'rangeof', h1);
+title(sprintf('Noisy image at date: t=%d', 17));
+fancyfigure;
+plotimagesar(u_hat_t, 'rangeof', h1);
+title(sprintf('Rabasar result (with geometric mean) at date: t=%d', t));
+
+if verbose == true
+    fancyfigure;
+    imagesc(v(:,:,t)./u_hat_t,[min1,max1]);
+    title(sprintf('Ratio noisy/denoised (with geometric mean) at date: t=%d', t));
+    colormap gray
+    fancyfigure;
+    histogram(v(:,:,t)./u_hat_t,200,'BinLimits',[0.1,10]);
+    title(sprintf('Ratio noisy/denoised (with geometric mean) at date: t=%d', t));
+    linkaxes;
+    mean(mean(v(:,:,t)./u_hat_t))
+end
+
+%%
+
+tic
+
+clear u_hat_t u_hat_si L_hat_si u_hat_dsi L_hat_dsi h
+h = robustwaitbar(0);
+
+[u_hat_t, u_hat_si, L_hat_si, u_hat_dsi, L_hat_dsi] = ...
+    rabasar_arith(v, L, t, ...
+            'waitbar', @(p) robustwaitbar(p, h));
+toc
+close(h);
+
+u_hat_t_a = u_hat_t;
+u_hat_si_a = u_hat_si;
+L_hat_si_a = L_hat_si;
+u_hat_dsi_a = u_hat_dsi;
+L_hat_dsi_a = L_hat_dsi;
+
+%% Visualization
+min1 = 0.1;
+max1 = 10;
+fancyfigure;
+plotimagesar(u_hat_si_a, 'rangeof', h1);
+title(sprintf('Arithmetic mean'));
+fancyfigure;
+plotimagesar(u_hat_dsi_a, 'rangeof', h1);
+title(sprintf('Denoised arithmetic mean'));
+fancyfigure;
+plotimagesar(u_hat_t, 'rangeof', h1);
+title(sprintf('Rabasar result (with arithmetic mean) at date: t=%d', t));
+
+if verbose == true
+    fancyfigure;
+    imagesc(v(:,:,t)./u_hat_t,[min1,max1]);
+    title(sprintf('Ratio noisy/denoised (with arithmetic mean) at date: t=%d', t));
+    colormap gray
+    fancyfigure;
+    histogram(v(:,:,t)./u_hat_t,200,'BinLimits',[0.1,10]);
+    title(sprintf('Ratio noisy/denoised (with arithmetic mean) at date: t=%d', t));
+    mean(mean(v(:,:,t)./u_hat_t))
+end
+%%
+
+
+    min1 = 0.1;
+    max1 = 7;
+    
+    fancyfigure;
+    imagesc(u_hat_dsi_a./u_hat_dsi_g,[min1,max1]);
+    title(sprintf('Ratio denoised arithmetic mean / denoised arithmetic mean'));
+    colormap gray
+    
+    fancyfigure;
+    imagesc(u_hat_si_a./u_hat_si_g,[min1,max1]);
+    title(sprintf('Ratio noisy arithmetic mean / noisy geometric mean'));
+    colormap gray
+    
+
+
+if verbose == true
+    
+    
+    min1 = 0.3;
+    max1 = 3;
+    
+    fancyfigure;
+    imagesc(u_hat_t_a./u_hat_t_g,[min1,max1]);
+    title(sprintf('Ratio arith/geom results at date: t=%d', t));
+    colormap gray
+    
+    fancyfigure;
+    imagesc(u_hat_t_g./u_hat_t_a,[min1,max1]);
+    title(sprintf('Ratio geom/arith results at date: t=%d', t));
+    colormap gray
+    
+
+    
+    fancyfigure;
+    histogram(u_hat_t_a./u_hat_t_g,200,'BinLimits',[0.1,10]);
+    title(sprintf('Ratio arith/geom at date: t=%d', t));
+    colormap gray
+    
+    fancyfigure;
+    histogram(u_hat_dsi_a./u_hat_dsi_g,200,'BinLimits',[0.1,10]);
+    title(sprintf('Ratio arith/geom dsi'));
+    colormap gray
+end
+    linkaxes;
+    

README.md

 # Geometric_Mean_Denoising
 
-Matlab implementation associated with the article "On the use and denoising of the temporal geometric mean for SAR time series", submitted to IEEE Geoscience and Remote sensing Letters by N.Gasnier, L.Denis and F.Tupin.
+Matlab implementation associated with the article "On the use and denoising of the temporal geometric mean for SAR time series", published in IEEE Geoscience and Remote sensing Letters by N.Gasnier, L.Denis and F.Tupin ([https://ieeexplore.ieee.org/document/9339951](url))
 
 
-The code is still under review for licence compliance and cannot be released yet.
+Author: Nicolas Gasnier
+contact : [nicolas.gasnier@telecom-paris.fr](mailto:nicolas.gasnier@telecom-paris.fr)
 
-If you want to be notified when it is released, please send a mail to [nicolas.gasnier@telecom-paris.fr](mailto:nicolas.gasnier@telecom-paris.fr)
+**Usage**
+
+The main results of the article can be reproduced by running the `example.m` script.
+
+The _t_-th image of a SAR intensity images stack `v` with an equivalent number of looks `L` can be denoised by using the following command in matlab:
+
+`h = robustwaitbar(0);
+[u_hat_t, u_hat_si, L_hat_si, u_hat_dsi, L_hat_dsi] = rabasar_geom(v, L, t,'waitbar', @(p) robustwaitbar(p, h));`
+
+It will return `u_hat_t` the denoised image at date `t`, `u_hat_si` the geometric mean of the time series, `L_hat_si` the ENL of the geometric mean, `u_hat_dsi` the denoised geometric mean (superimage) and its estimated ENL `L_hat_dsi`.
+
+
+
+**Using the BM3D denoiser**
+
+Please note that because of copyright issues, we could not to release our code along with the BM3D implementation used for our experiments. A DDID denoiser is used instead and its results are differents. To use a BM3D denoiser, you can copy the content of the denoisers folder that can be found in the original mulog repository:[https://bitbucket.org/charles_deledalle/mulog/src/master/](url)) in mulog/denoisers and change the choice of denoiser by uncommenting line 90 and commenting line 91 in both rabasar/rabasar_arith.m and rabasar/rabasar_geom.m
+
+
+
+**Copyright and license**
+
+Copyright (C) Copyright 2020 CS GROUP - France, Centre national d'études spatiales, Télécom Paris.
+
+The Geometric Mean Denoising software as a whole is distributed under the CeCILL licence, version 2.0. A copy of this license is available in the LICENSE file. 
+
+All code snippets provided in the documentation, tutorials and examples (available in the Readme and in the example script) are intended to show how to use Geometric Mean Denoising and are not, strictly speaking, part of Geometric Mean Denoising project. In order to facilitate their reuse, they are released under the “Zero-Clause BSD” license (aka BSD 0-Clause or 0BSD). This is a permissive license that allows you to copy and paste this source code without preserving copyright notices.

check/check_approx_inv_psi.m

+%% Check the approximation of the invert of psi used for ENL estimation
+
+clear all
+close all
+
+siglist = linspace(.01, 60, 1000);
+stdlog = @(L) sqrt(psi(1, L));
+for k = 1:length(siglist)
+    sig = siglist(k);
+    ENLn(k) = exp(fminbnd(@(logL) (stdlog(exp(logL)) - sig).^2, log(0.001), log(200000)));
+
+    k2 = sig^2 / 4;
+    ENLa(k) = ...
+        (1.+3.5232*k2+1.57472*k2.^2+4.8288*10^(-2)*k2.^3) ./ ...
+        (-1.705*10^(-5)+4.004*k2+5.9856*k2.^2+8.6208*10^(-1)*k2.^3);
+end
+
+fancyfigure;
+loglog(siglist, ENLn);
+hold on
+plot(siglist, ENLa, '--');
+legend('Exact', 'Approx');

check/check_mulog_vs_molog.m

+%% Check that MoLoG=MuLoG
+
+clear all
+close all
+
+init();
+
+%% Setting
+L = 2;
+denoiser = @bm3d;
+
+%% Gaussian noise simulation setting
+x = loadimage('data/cameraman.png').^2;
+[m, n] = size(x);
+y = x .* mean((randn(m, n, L).^2 + randn(m, n, L).^2) / 2, 3);
+
+%% Run MuLoG with embedded BM3D Gaussian denoiser
+disp('Run MuLoG with embedded BM3D Gaussian denoiser');
+tic;
+h = robustwaitbar(0);
+xhat = mulog(y, L, denoiser, ...
+            'waitbar', @(p) robustwaitbar(p, h));
+close(h);
+disp(sprintf('  Elapsed time %.2f s', toc));
+
+%% Run MoLoG with embedded BM3D Gaussian denoiser
+disp('Run MoLoG with embedded BM3D Gaussian denoiser');
+tic;
+h = robustwaitbar(0);
+xhat2 = molog(y, L, denoiser, ...
+              'waitbar', @(p) robustwaitbar(p, h));
+close(h);
+disp(sprintf('  Elapsed time %.2f s', toc));
+
+%% Display results
+f = fancyfigure;
+subplot(1, 3, 1);
+h = plotimagesar(y, 'range', [0 255]);
+title(sprintf('Noisy image (psnr %.2f)', perfs(y, x, L, 'psnr')));
+subplot(1, 3, 2);
+plotimagesar(xhat, 'rangeof', h);
+title(sprintf('MuLoG+BM3D (psnr %.2f)', perfs(xhat, x, L, 'psnr')));
+subplot(1, 3, 3);
+plotimagesar(xhat2, 'rangeof', h);
+title(sprintf('MoLoG-BM3D (psnr %.2f)', perfs(xhat2, x, L, 'psnr')));
+linkaxes

check/check_rulog_vs_molog.m

+%% Check that MuLoG=RuLog when Lden is large
+
+clear all
+close all
+
+init();
+
+%% Setting
+L = 2;
+denoiser = @bm3d;
+
+%% Gaussian noise simulation setting
+x = loadimage('data/cameraman.png').^2;
+[m, n] = size(x);
+y = x .* mean((randn(m, n, L).^2 + randn(m, n, L).^2) / 2, 3);
+
+%% Run MoLoG with embedded BM3D Gaussian denoiser
+disp('Run MoLoG with embedded BM3D Gaussian denoiser');
+tic;
+h = robustwaitbar(0);
+xhat = molog(y, L, denoiser, ...
+            'waitbar', @(p) robustwaitbar(p, h));
+close(h);
+disp(sprintf('  Elapsed time %.2f s', toc));
+
+%% Run RuLoG with embedded BM3D Gaussian denoiser
+disp('Run RuLoG with embedded BM3D Gaussian denoiser');
+tic;
+h = robustwaitbar(0);
+xhat2 = rulog(y, L, 1000, denoiser, ...
+              'waitbar', @(p) robustwaitbar(p, h));
+close(h);
+disp(sprintf('  Elapsed time %.2f s', toc));
+
+%% Display results
+f = fancyfigure;
+subplot(1, 3, 1);
+h = plotimagesar(y, 'range', [0 255]);
+title(sprintf('Noisy image (psnr %.2f)', perfs(y, x, L, 'psnr')));
+subplot(1, 3, 2);
+plotimagesar(xhat, 'rangeof', h);
+title(sprintf('MoLoG+BM3D (psnr %.2f)', perfs(xhat, x, L, 'psnr')));
+subplot(1, 3, 3);
+plotimagesar(xhat2, 'rangeof', h);
+title(sprintf('RuLoG-BM3D (psnr %.2f)', perfs(xhat2, x, L, 'psnr')));
+linkaxes

mulog/.gitignore

+*~
+#*

mulog/Makefile

+clean:
+	@find . -name '*~' -exec rm -f {} \;

mulog/README.md

+# MuLoG
+
+This is a Matlab script implementing the MuLoG algorithm for speckle reduction
+in intensity and mutivariate SAR data. For more details about the approach
+please refer to: https://www.charles-deledalle.fr/pages/mulog
+
+
+## Installation
+
+Just clone the project.
+
+## Usage
+
+Have a look at the `example_*.m` files.
+
+## License
+
+This software is a computer program whose purpose is to provide a suite of tools
+to manipulate SAR images.
+
+This software is governed by the CeCILL license under French law and
+abiding by the rules of distribution of free software. You can use,
+modify and/ or redistribute the software under the terms of the CeCILL
+license as circulated by CEA, CNRS and INRIA at the following URL
+"http://www.cecill.info".
+
+As a counterpart to the access to the source code and rights to copy,
+modify and redistribute granted by the license, users are provided only
+with a limited warranty and the software's author, the holder of the
+economic rights, and the successive licensors have only limited
+liability.
+
+In this respect, the user's attention is drawn to the risks associated
+with loading, using, modifying and/or developing or reproducing the
+software by the user in light of its specific status of free software,
+that may mean that it is complicated to manipulate, and that also
+therefore means that it is reserved for developers and experienced
+professionals having in-depth computer knowledge. Users are therefore
+encouraged to load and test the software's suitability as regards their
+requirements in conditions enabling the security of their systems and/or
+data to be ensured and, more generally, to use and operate it in the
+same conditions as regards security.
+
+The fact that you are presently reading this means that you have had
+knowledge of the CeCILL license and that you accept its terms.
+

mulog/addpathrec.m

+function addpathrec(path)
+%% Add all directories and subdirectories recursively to matlab
+%  pathdefs variable
+%
+%
+% License
+%
+% This software is governed by the CeCILL license under French law and
+% abiding by the rules of distribution of free software. You can use,
+% modify and/ or redistribute the software under the terms of the CeCILL
+% license as circulated by CEA, CNRS and INRIA at the following URL
+% "http://www.cecill.info".
+%
+% As a counterpart to the access to the source code and rights to copy,
+% modify and redistribute granted by the license, users are provided only
+% with a limited warranty and the software's author, the holder of the
+% economic rights, and the successive licensors have only limited
+% liability.
+%
+% In this respect, the user's attention is drawn to the risks associated
+% with loading, using, modifying and/or developing or reproducing the
+% software by the user in light of its specific status of free software,
+% that may mean that it is complicated to manipulate, and that also
+% therefore means that it is reserved for developers and experienced
+% professionals having in-depth computer knowledge. Users are therefore
+% encouraged to load and test the software's suitability as regards their
+% requirements in conditions enabling the security of their systems and/or
+% data to be ensured and, more generally, to use and operate it in the
+% same conditions as regards security.
+%
+% The fact that you are presently reading this means that you have had
+% knowledge of the CeCILL license and that you accept its terms.
+%
+% Copyright 2017 Charles Deledalle
+% Email charles-alban.deledalle@math.u-bordeaux.fr
+
+
+
+list = dir(path);
+for k = 1:length(list)
+    if ~strcmp(list(k).name(1), '.') && list(k).isdir
+        addpath([path '/' list(k).name]);
+        addpathrec([path '/' list(k).name]);
+    end
+end

mulog/comparison_mulog_with_bm3d_and_tv.m

+clear all
+close all
+
+addpathrec('.');
+
+%% Load a polsar image
+C = loadsarimg('data/tsukaba_pisar.mat');
+
+%% Setting
+L            = 2;       % input number of looks
+R            = 1;       % input rank
+
+%% Run MuLoG
+h = robustwaitbar(0);
+
+%% ... with TV
+disp('Run MuLoG with embedded TV');
+tic;
+Shat_tv = ...
+    mulog(C, L, @tv, ...
+          'waitbar', @(p) robustwaitbar(p/2, h), ...
+          'R', R);
+disp(sprintf('  Elapsed time %.2f s', toc));
+
+%% ... with BM3D
+disp('Run MuLoG with embedded BM3D');
+tic;
+Shat_bm3d = ...
+    mulog(C, L, @bm3d, ...
+          'waitbar', @(p) robustwaitbar(p/2 + 1/2, h), ...
+          'R', R);
+disp(sprintf('  Elapsed time %.2f s', toc));
+
+close(h);
+
+%% Display results
+close all
+f = fancyfigure;
+subplot(1, 3, 1);
+h = plotimagesar(C, 'beta', 3, 'alpha', 0.7, 'adjust', 'm2s');
+title('Input SAR image');
+subplot(1, 3, 2);
+plotimagesar(Shat_tv, 'rangeof', h);
+title('MuLoG with TV');
+subplot(1, 3, 3);
+plotimagesar(Shat_bm3d, 'rangeof', h);
+title('MuLoG with BM3D');
+linkaxes;
+
+%% Save figure
+%savesubfig(f, '/tmp/result');

mulog/denoisers/bm3d.m

+function x = bm3d(y, sig, varargin)
+%% Wraper to the BM3D implementation of the authors
+%
+% Input/Ouput:
+%
+%    Y          the observation (size MxNxK)
+%
+%    X          the solution (size MxNxK)
+%
+%    SIG        the standard deviation of the noise
+%
+%
+% Reference
+%
+%    K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian,
+%    Image Denoising by Sparse 3D Transform-Domain Collaborative Filtering,
+%    IEEE Transactions on Image Processing, vol. 16, no. 8, August, 2007.
+%
+%
+% License
+%
+% This software is governed by the CeCILL license under French law and
+% abiding by the rules of distribution of free software. You can use,
+% modify and/ or redistribute the software under the terms of the CeCILL
+% license as circulated by CEA, CNRS and INRIA at the following URL
+% "http://www.cecill.info".
+%
+% As a counterpart to the access to the source code and rights to copy,
+% modify and redistribute granted by the license, users are provided only
+% with a limited warranty and the software's author, the holder of the
+% economic rights, and the successive licensors have only limited
+% liability.
+%
+% In this respect, the user's attention is drawn to the risks associated
+% with loading, using, modifying and/or developing or reproducing the
+% software by the user in light of its specific status of free software,
+% that may mean that it is complicated to manipulate, and that also
+% therefore means that it is reserved for developers and experienced
+% professionals having in-depth computer knowledge. Users are therefore
+% encouraged to load and test the software's suitability as regards their
+% requirements in conditions enabling the security of their systems and/or
+% data to be ensured and, more generally, to use and operate it in the
+% same conditions as regards security.
+%
+% The fact that you are presently reading this means that you have had
+% knowledge of the CeCILL license and that you accept its terms.
+%
+% Copyright 2017 Charles Deledalle
+% Email charles-alban.deledalle@math.u-bordeaux.fr
+
+
+
+options   = makeoptions(varargin{:});
+cbwaitbar = getoptions(options, 'waitbar', @(dummy) []);
+
+x = zeros(size(y));
+m = min(y(:));
+M = max(y(:));
+K = size(y, 3);
+if isscalar(sig)
+    sig = sig * ones(K ,1);
+end
+parfor k = 1:K
+    tau = 255 * sig(k) / (M - m);
+    if tau > 40
+        ratio = 40 / tau;
+    else
+        ratio = 1;
+    end
+    [~, x(:,:,k)] = BM3D(1, ratio * (y(:,:,k) - m) / (M - m), ...
+                         ratio * tau, 'np', 0);
+    x(:, :, k) = (M - m) * x(:, :, k) / ratio + m;
+    cbwaitbar((K-k+1) / K);
+end
+x(isnan(x)) = 0;
+cbwaitbar(1);

mulog/denoisers/ddid.m

+function x = ddid(y, sigma, varargin)
+%% Wraper to the DDID implementation of the authors
+%
+% Input/Ouput:
+%
+%    Y          the observation (size MxNxK)
+%
+%    X          the solution (size MxNxK)
+%
+%    SIG        the standard deviation of the noise
+%
+%
+% Reference
+%
+%    Knaus, Claude, and Matthias Zwicker,
+%    Dual-domain image denoising." Image Processing (ICIP),
+%    20th IEEE International Conference on. IEEE, 2013.
+%
+%
+% License
+%
+% This software is governed by the CeCILL license under French law and
+% abiding by the rules of distribution of free software. You can use,
+% modify and/ or redistribute the software under the terms of the CeCILL
+% license as circulated by CEA, CNRS and INRIA at the following URL
+% "http://www.cecill.info".
+%
+% As a counterpart to the access to the source code and rights to copy,
+% modify and redistribute granted by the license, users are provided only
+% with a limited warranty and the software's author, the holder of the
+% economic rights, and the successive licensors have only limited
+% liability.
+%
+% In this respect, the user's attention is drawn to the risks associated
+% with loading, using, modifying and/or developing or reproducing the
+% software by the user in light of its specific status of free software,
+% that may mean that it is complicated to manipulate, and that also
+% therefore means that it is reserved for developers and experienced
+% professionals having in-depth computer knowledge. Users are therefore
+% encouraged to load and test the software's suitability as regards their
+% requirements in conditions enabling the security of their systems and/or
+% data to be ensured and, more generally, to use and operate it in the
+% same conditions as regards security.
+%
+% The fact that you are presently reading this means that you have had
+% knowledge of the CeCILL license and that you accept its terms.
+%
+% Copyright 2017 Charles Deledalle
+% Email charles-alban.deledalle@math.u-bordeaux.fr
+
+
+
+options   = makeoptions(varargin{:});
+cbwaitbar = getoptions(options, 'waitbar', @(dummy) []);
+[m, n, K] = siz