cvx_unwrapping

REQUIREMENTS
INITIALIZATION
GROUND TRUTH (GT)
ADDING GAUSSIAN NOISE
WRAPPED DATA
COMPUTING THE NOISE BOUND
CVX UNWRAPPING

% Robust Phase Unwrapping by Convex Optimization
%
% This script shows the behavior of the phase unwrapping algorithm
% 'cvx_unwrapping.m' described in [1].
%
% The test is performed on a synthetic phase image defined in a 256x256
% pixel grid. This image is simulated by a 2-D Gaussian function of height
% 0.9*rho*pi (rho is a parameter to modify the height of the Gaussian), and
% standard deviations of 40 px horizontally and 25 px vertically. It then
% adds normal distributed noise and compute the observations using the
% wrapping function.
%
% The unwrapping algorithm in 'cvx_unwrapping.m' (see [1]) allows solving
% the following minimization problem:
%
% [x*,v*] = argmin_[x,v] || Psi^T x ||_1
%                   s.t. || v ||_2 < epsilon1
%                        || q - nabla(x + v)||_1 < epsilon2
%                        x(1) = 0
% with y the observation (wrapped phase) and q(y) = wrap(grad(y))
%
% References:
%
% [1] A. Gonzalez and L. Jacques. ''Robust Phase Unwrapping by Convex
% Optimization,'' in IEEE International Conference on Image Processing
% (ICIP14), Paris, France. October 27-30, 2014.
%
% [2] SPARCO Toolbox. E. v. Berg, M. P. Friedlander, G. Hennenfent, F.
% Herrmann, R. Saab, and O. Yilmaz, ''Sparco: A testing framework for
% sparse reconstruction,'' Dept. Computer Science, University of British
% Columbia, Vancouver, Tech. Rep. TR-2007-20, October 2007.
%
% [3] J. Bioucas-Dias and G. Valad˜ao, ''Phase unwrapping via graph-cuts,''
% IEEE Transactions on Image Processing, vol. 16, no. 3, pp. 698–709, 2007.
%
% [4] Rice Wavelet Toolbox.
% http://dsp.rice.edu/software/rice-wavelet-toolbox
%
% BibTeX Citation when using this code:
% @INPROCEEDINGS{7025343,
% author={Gonzalez, A. and Jacques, L.},
% booktitle={Image Processing (ICIP), 2014 IEEE International Conference on},
% title={Robust phase unwrapping by convex optimization},
% year={2014},
% pages={1713-1717},
% doi={10.1109/ICIP.2014.7025343},
% month={Oct},}
%
% Copyright 2015 Adriana Gonzalez
% PhD Student at Universite catholique de Louvain (UCL), Belgium

REQUIREMENTS

% The Rice Wavelet Toolbox (RWT) must be installed and must be located in
% the following path
addpath(genpath(strcat(pwd,'\rwt')))

% The PUMA algorithm must be located in the following path
addpath(genpath(strcat(pwd,'\PUMA')))

INITIALIZATION

% Intitialization of the algorithm ****************************************
init = 'puma_denoised';
            % 'y': initialized using the observed data
            % 'zeros': initialized using a vector of zeros
            % 'puma': initialized using the output of PUMA
            % 'puma_denoised': initialized using the denoised PUMA

% Computation of the noise bound ******************************************
type_epsilon1 = 'Oracle';
            % 'GT':     computing the noise bound using the GT of the noise
            %           variance
            % 'Oracle': computing the noise bound using the variance
            %           estimated by the robust median estimator

% Computation of the error bound ******************************************
type_epsilon2 = 'Oracle';
            % 'GT': computing the error bound from the GT image
            % 'Puma': computing the error bound from the PUMA output

% Parameters to adjust ****************************************************

% Parameter that allows analyzing the behavior of the algorithm with
% respect to the amount of ''wraps''. rho = 1 provides no ''wraps''
rho = 10;

% Boolean variable to decide whether to show the images or not
show_im = 1;

% Standard deviation of the Gaussian noise to add to the phase
std_dev = 0.35;

% Data size
n = 256;
N = n^2;

% Gradient Operator
opG = opGrad_AG(n);

% Wavelet Operator
w.name = 'daubechies';  % Mother wavelet family: 'daubechies' or 'haar'
w.filter = 14;          % Number of taps of the filter
w.level = 7;            % Number of decomposition levels
w.type = 1;             % 1:DWT, 2:UDWT
opW = opTranspose(opWavelet_AG(n,n,w.name,w.filter,w.level,'min',w.type));

GROUND TRUTH (GT)

[xx,yy] = meshgrid(-n/2+1:n/2,-n/2+1:n/2); sx = 40; sy = 25;
im_noiseless = 0.9*pi*(exp(-xx.^2/(2*sx^2)-yy.^2/(2*sy^2)));
clear xx yy sx sy top

im_noiseless = reshape(im_noiseless*rho,N,1);
im_noiseless_aux = im_noiseless - mean(im_noiseless);

if show_im
    figure; colormap(pink(256)); imagesc(reshape(im_noiseless,n,n));
    colorbar; axis image;
    title('Continuous phase image (x)')
    xlabel('Pixels'); ylabel('Pixels')
end

ADDING GAUSSIAN NOISE

noise = std_dev*randn(N,1);
im = im_noiseless + noise;
im_aux = im - mean(im);

MSNR = 20*log10(norm(im_noiseless)/norm(noise));
fprintf('MSNR = %3.5fdB\n', MSNR);

if show_im
    figure; colormap(pink(256)); imagesc(reshape(im,n,n));
    colorbar; axis image;
    title(sprintf('Noisy continuous phase image (x+n) - MSNR = %2.2fdB',MSNR))
    xlabel('Pixels'); ylabel('Pixels')
end

MSNR = 24.93145dB

WRAPPED DATA

% Observations
y = mod(im+pi,2*pi)-pi;

% Indirect observations
q = mod(opG(y(:),1)+pi,2*pi) - pi;

if show_im

    figure;colormap(pink(256));imagesc(reshape(y,n,n));
    colorbar; axis image;
    title('Measurement - Wrapped phase - y = (x+n) mod 2\pi')
    xlabel('Pixels'); ylabel('Pixels')

end

COMPUTING THE NOISE BOUND

switch type_epsilon1

    case 'GT'

        gx = opG(im_noiseless,1);
        epsilon1 = std_dev*sqrt(N + 2*sqrt(N));

    case 'Oracle'

        % Computation of the variance of the observation noise

        Wn = 0.4;  % Normalized cut-off frequency to filter the noise from
                   % the noisy signal inside the robust median estimator

        QX = reshape(q(1:N),n,n);
        QY = reshape(q(N+1:2*N),n,n);

        aux_X = 0; aux_Y = 0;
        for ii=1:n
            aux_X = aux_X + sqrt(median_estimator((n-1),QX(ii,:),Wn));
            aux_Y = aux_Y + sqrt(median_estimator((n-1),QY(ii,:),Wn));
        end
        sigma_X = aux_X/n;
        sigma_Y = aux_Y/n;

        epsilon1 = (sigma_X + sigma_Y)*sqrt(N + 2*sqrt(N));

end

fprintf('|| n ||_2 = %e\n', epsilon1);

|| n ||_2 = 9.027190e+01

CVX UNWRAPPING

% Algorithm parameters ****************************************************

param.iter = [];
param.ThMode = 3; % 1 => Stability, 2 => Residuals, 3 => Both
param.Th = [1.5e-4 5e-1 5e-1]; % [Stability Residuals]
param.max_iter = 1e4;
param.res_norm = 2;
param.t = 1e3;
param.print = 0;

% Compute the output of PUMA if it is needed for initialization or for
% computing the error bound ***********************************************

if strcmp(type_epsilon2,'Oracle')&&(strcmp(init,'puma_denoised')||strcmp(init,'puma_denoised'))

    % PUMA output
    [unwph,~,~] = puma_ho(reshape(y,n,n),2,'verbose','no'); xp = unwph(:);

    % Denoising
    fun = @(lambda) abs(norm(xp-opW(SoftTh(opW(xp,1),lambda),2))-epsilon1);
    lambda_min = fminsearch(fun,epsilon1);
    xp_denoised = opW(SoftTh(opW(xp,1),lambda_min),2);

end

% Computing the error bound ***********************************************

switch type_epsilon2
    case 'GT'
        epsilon2 = norm(q - opG(im,1),1);
    case 'Oracle'
        epsilon2 = norm(q - opG(xp,1),1);
end

fprintf('|| q - nabla(x+n) ||_1 = %e\n', epsilon2);

% Ground Truth ************************************************************

gt = [im_noiseless; noise];

% Initialization **********************************************************

switch init
    case 'zeros'
        x0 = zeros(2*N,1);
    case 'y'
        x0 = [y; zeros(N,1)];
    case 'puma_denoised'
        x0 = [xp_denoised; zeros(N,1)];
    case 'puma'
        x0 = [xp; zeros(N,1)];
end

% Algorithm ***************************************************************

[est,iter,T,SNR,p,d,condition,time,l1norm_x,l1norm_gt] = ...
    cvx_unwrapping(q,epsilon1,epsilon2,w,gt,x0,param);

x = est(1:N);

% Quality *****************************************************************

x_aux = x - mean(x);

SNRx = 20*log10(norm(im_noiseless_aux)/norm(im_noiseless_aux - x_aux));
RMSEx = std(im_noiseless_aux-x_aux);

fprintf('\ncvx unwrapping: RSNR = %2.4f dB\n',SNRx)

% Results *****************************************************************

if show_im

    figure; colormap(pink(256)); imagesc(reshape(x,n,n));
    colorbar; axis image;
    title('Unwrapped Phase')
    xlabel('Pixels'); ylabel('Pixels')

end

|| q - nabla(x+n) ||_1 = 1.125100e-12

|| W^T x ||_1 = 2.476665e+03
|| W^T x* ||_1 = 2.837763e+03

epsilon = 9.027190e+01
|| v* ||_2 = 8.814280e+01

|| q - grad(x*+v*) ||_1 = 2.919999e+02

cvx unwrapping: RSNR = 40.8506 dB

Contents

REQUIREMENTS

INITIALIZATION

GROUND TRUTH (GT)

ADDING GAUSSIAN NOISE

WRAPPED DATA

COMPUTING THE NOISE BOUND

CVX UNWRAPPING