Commit 93bec743 authored by Nicolas Gasnier's avatar Nicolas Gasnier
Browse files

Merge branch 'release' into 'master'

Release

See merge request !1
parents c3d5338d 95f9cfcb
%% 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;
This diff is collapsed.
This diff is collapsed.
# 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 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 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 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
File added
This diff is collapsed.
This diff is collapsed.
clean:
@find . -name '*~' -exec rm -f {} \;
# 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.
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
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');
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);
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) []);