Initial commit

.gitignore

@@ -0,0 +1,12 @@
+
+unitTests/swUnitTestDistTransformL1.m
+
+Docs/create_images.m
+
+unitTests/swDistTransformL1FH.m
+
+unitTests/swUnitTestDistAngles.m
+
+unitTests/swUnitTestDistTransformNE.m
+
+unitTests/swUnitTestMinL1TVCirc.m

Auxiliary/deltaSNR.m

@@ -0,0 +1,14 @@
+function dsnr = deltaSNR( groundTruth, data, estimate, dataSpace)
+%deltaSNR Computes the signal to noise ratio improvement of a restoration
+
+switch dataSpace
+    case 'circ'
+        dist = @(x,y) distAngle(x,y);
+    case 'real'
+        dist = @(x,y) abs(x - y);
+end
+
+dsnr = 10 * log10(sum( dist(groundTruth, data).^2 ) / sum( dist(groundTruth, estimate).^2 ));
+
+end
+

Auxiliary/distAngle.m

@@ -0,0 +1,14 @@
+function d = distAngle( phi, psi )
+%distAngle Distance between two angles phi and psi
+
+if (all(abs(phi) <= pi)) && (all(abs(psi) <= pi))
+    % cheap computation for angles in the interval [-pi, pi]
+    aux = phi - psi;
+    d = min(abs([aux(:) + 2 * pi, aux(:), aux(:) - 2 * pi]), [], 2);
+    d = reshape(d, size(aux));
+else
+    % more expensive calculation for angles outside the interval [-pi, pi]
+    aux = abs(angle(exp(1i * (phi - psi))));
+    d = min(aux, 2*pi-aux);
+end
+

Auxiliary/distTransformL1.m

@@ -0,0 +1,15 @@
+function B = distTransformL1( B, R, alpha )
+%distTransformL1 Fast computation of L1 distance transform
+
+K = numel(B);
+% forward pass
+for k=2:K
+    B(k) = min( B(k), B(k-1) + alpha * (R(k) - R(k-1)));
+end
+% backward pass
+for k=K-1:-1:1
+    B(k) = min( B(k), B(k+1) + alpha * (R(k+1) - R(k)));
+end
+
+end
+

Auxiliary/randCP.m

@@ -0,0 +1,9 @@
+function r = randCP( r, lambda)
+%randCP Generates a random vector distributed according to a certain 
+% compound Poisson distribution with parameter lambda
+
+x = rand(size(r)); % uniform distr. vector
+r(x(:) <= exp(-lambda)) = 0; % compound Poisson
+
+end
+

Auxiliary/randl.m

@@ -0,0 +1,10 @@
+function y = randl( varargin )
+%randl Random normal Laplacian noise, mean = 0, variance sigma^2 = 1
+% pdf: $\frac{1}{\sqrt{2}} e^{- \sqrt{2} |x|}$
+% Same input arguments as function rand
+
+x = rand( varargin{:} ) - 0.5;
+y = -sign(x) .* log(1 - 2 * abs(x)) / sqrt(2);
+
+end
+

Auxiliary/wrapAngle.m

@@ -0,0 +1,7 @@
+function x = wrapAngle( y )
+%wrapAngle Wraps the angle y to the interval [-pi, pi]
+
+x = angle(exp(1i * y));
+
+end
+

Demos/demo_L1TV_Circ.m

@@ -0,0 +1,36 @@
+%%% Demo for denoising a circle-valued signal by the L1-TV model
+
+% create random signal (smoothed pcw constant signal)
+rng(12345) % random seed for reproducibility
+N = 2000;
+t = 2*pi;
+lambda = 20 / N;
+K = 20;
+sigma= 0.3;
+innovation = randCP((rand([N, 1])-0.5) * t, lambda );
+signalUnwrapped = cumsum(innovation);
+h = fspecial('Gaussian', [N/10, 1], 10);
+smoothed = conv(signalUnwrapped, h, 'same');
+groundTruth = wrapAngle(smoothed);
+
+% add noise
+y = wrapAngle(groundTruth + sigma* randl(size(groundTruth)));
+
+% perform restoration using L1TV_Circ
+alpha = sqrt(N)*sigma;
+x = L1TV_Circ(y, alpha);
+
+% plot the results
+figure('Color', 'w')
+subplot(1,2,1)
+plot(y, '.')
+ylim([-pi,pi])
+title('Data with values on the unit circle')
+
+subplot(1,2,2)
+plot(x, '.')
+ylim([-pi,pi])
+title(sprintf('L1TV restoration, SNR improvement: %.2f dB', deltaSNR(groundTruth, y, x, 'circ')))
+
+
+

Demos/demo_L1TV_Real.m

@@ -0,0 +1,33 @@
+%%% Demo for denoising a circle-valued signal by the L1-TV model
+
+% create random signal (smoothed pcw constant signal)
+rng(12345) % random seed for reproducibility
+N = 2000;
+lambda = 20 / N;
+K = 20;
+scale = 100; % scale of signal (change to observe constrast invariance of L1TV)
+sigma= scale* 0.3;
+innovation =  randCP(randn([N, 1]), lambda );
+signal = scale * cumsum(innovation);
+h = fspecial('Gaussian', [N/10, 1], 10);
+groundTruth = conv(signal, h, 'same');
+
+% add noise
+y = groundTruth + sigma * randl(size(groundTruth));
+
+% perform restoration using L1TV_Real
+alpha = sqrt(N)*sigma/scale;
+x = L1TV_Real(y, alpha);
+
+% plot the results
+figure('Color', 'w')
+subplot(1,2,1)
+plot(y, '.')
+ylim_set = ylim;
+title('Real valued data')
+
+subplot(1,2,2)
+plot(x, '.')
+ylim(ylim_set);
+title(sprintf('L1TV restoration, SNR improvement: %.2f dB', deltaSNR(groundTruth, y, x, 'real')))
+

L1TV_Circ.m

@@ -0,0 +1,73 @@
+function x = L1TV_Circ( y, alpha, varargin )
+%L1TV_Circ Exact solver for the univariate L1TV problem with circle valued
+%data. i.e. it computes the solution of
+% min_{x \in T^n} \alpha sum_{i=1}^{n-1} d(x_i, x_{i+1}) + sum_{i=1}^{n} w_i d(x_i, y_{i}) 
+% where T is the unit circle (torus) and d(x, y) is the shortest arc
+% distance between x and y
+%
+% Input:
+% y: n-vector of phase angles between [-pi, pi]
+% alpha: regularization parameter
+% Optional arguments
+% 'weights': pointwise weights for the data fidelity (n-vector of positve numbers)
+%
+% Reference:
+% Storath, M., Weinmann, A., & Unser, M. (2016). 
+% Exact algorithms for L^1-TV regularization of real-valued or circle-valued signals. 
+% SIAM Journal on Scientific Computing, 38(1), A614-A630.
+
+
+% parse input
+ip = inputParser;
+ip.addParameter('useDistTrans', true);
+ip.addParameter('weights', ones(size(y)));
+
+parse(ip, varargin{:});
+par = ip.Results;
+
+% initialization
+yCom = exp(1i * y); % complex number representation of y
+yAng = angle(yCom); % assure angle to be in [-pi, pi]
+uniqueValues = unique(yAng); % determine uniqe angles
+V = [uniqueValues; angle( -exp(1i * uniqueValues))]; % add antipodal points
+V = sort(V); % sort candidate values
+N = numel(yAng);
+K = numel(V);
+B = zeros(K,N); 
+pen = zeros(K,1);
+
+
+% compute tabulation
+B(:,1) = par.weights(1) .* distAngle(V, yAng(1));
+if par.useDistTrans % Using distance transforms, the complexity is O(N K)
+    % auxiliary structure for distance transforms
+    Vrep = cat(1, V-2*pi, V, V+2*pi);
+    for n=2:N
+        d = par.weights(n) .* distAngle(V, yAng(n));
+        Brep = repmat(B(:,n-1), [3,1]);
+        aux = distTransformL1(Brep, Vrep, alpha);
+        pen = min(reshape(aux, K, 3), [], 2);
+        B(:,n) = d + pen;
+    end
+else
+    for n=2:N % Naive implementation is O(N K^2)
+        d = par.weights(n) .* distAngle(V, yAng(n));
+        for k=1:K
+            pen(k) = min( B(:,n-1) +  alpha * distAngle(V, V(k)));
+        end
+        B(:,n) = d + pen;
+    end
+
+end
+
+% backtracking
+[~,l] = min(B(:,N));
+x(N,1) = V(l);
+for n=N-1:-1:1
+    [~, l] = min(B(:, n) + alpha * distAngle(V, x(n+1)));
+    x(n) = V(l);
+end
+
+
+end
+

L1TV_Real.m

@@ -0,0 +1,58 @@
+function x = L1TV_Real( y, alpha, varargin )
+%L1TV_Real Exact solver for the univariate L1TV problem with real-valued
+%data. i.e. it computes the solution of
+% min_{x \in R^n} \alpha sum_{i=1}^{n-1} |x_i - x_{i+1}| + sum_{i=1}^{n} w_i |x_i - y_{i}| 
+% where R is the real numbers 
+%
+% Input:
+% y: n-vector of real numbers
+% alpha: regularization parameter
+% Optional arguments
+% 'weights': pointwise weights for the data fidelity (n-vector of positve numbers)
+%
+% Reference:
+% Storath, M., Weinmann, A., & Unser, M. (2016). 
+% Exact algorithms for L^1-TV regularization of real-valued or circle-valued signals. 
+% SIAM Journal on Scientific Computing, 38(1), A614-A630.
+
+% parse input
+ip = inputParser;
+ip.addParameter('useDistTrans', true);
+ip.addParameter('weights', ones(size(y)));
+
+parse(ip, varargin{:});
+par = ip.Results;
+
+% initialization
+uniqueValues = unique(y); % determine unique angles
+V = sort(uniqueValues); % sort values
+N = numel(y);
+K = numel(V);
+B = zeros(K,N); % tabulation
+pen = zeros(K,1);
+
+% tabulation
+B(:,1) = par.weights(1) .* abs(V - y(1));
+for n=2:N
+    d = par.weights(n) .* abs(V - y(n));
+    if par.useDistTrans % Using distance transforms, the complexity is O(NK)
+        pen = distTransformL1(B(:,n-1), V, alpha);
+    else % The naive implemention is O(NK^2)
+        for k=1:K
+            pen(k) = min( B(:,n-1) +  alpha * abs(V - V(k)));
+        end
+    end
+    B(:,n) = d + pen;
+end
+
+% backtracking
+[~,l] = min(B(:,N));
+x(N,1) = V(l);
+for n=N-1:-1:1
+    [~, l] = min(B(:, n) + alpha * abs(V - x(n+1)));
+    x(n) = V(l);
+end
+
+
+end
+

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2020 Martin Storath
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,19 @@
+Matlab functions for estimation (denoising/reconstruction) of approximately piecewise constant signals.
+The functions are reference implementations of the method described in the paper
+- M. Storath, A. Weinmann, M. Unser. 
+  " Exact algorithms for L^1-TV regularization of real-valued or circle-valued signals."
+  SIAM Journal on Scientific Computing, 38(1), A614-A630, 2016
+
+## Estimation of real-valued signals
+Estimates a real-valued signal using the L1-TV model (exact non-iterative solver)
+
+![L1 TV Denoising of real-valued signal](/Docs/L1TV_Real.png)
+
+## Estimation of circle-valued signals
+Estimates a circle-valued signal using the L1-TV model (exact non-iterative solver)
+
+![L1 TV Denoising of circle-valued signal](/Docs/L1TV_Circ.png)
+
+## Installation and usage
+- Set the Matlab path by calling setPath.m
+- Run the demos of the Demos folder

setPath.m

@@ -0,0 +1,10 @@
+% sets the paths 
+disp('Setting Matlab path...');
+folder = fileparts(which(mfilename));
+addpath(...
+    fullfile(folder, ''),...
+    fullfile(folder, 'Auxiliary') ...
+    );
+
+% save pathdef
+savepath;