Zip file

Convex Unwrapping Matlab Demonstration files (zip, 511kB)

Description

The script cvx_unwrapping.m (developped under MATLAB® 7.14) in the zip file above shows the behavior of the phase unwrapping algorithm described in (Gonzalez Gonzalez and Jacques 2014) whose abtract is:

The 2-D phase unwrapping problem aims at retrieving a “phase” image from its modulo $2\pi$ observations. Many applications, such as interferometry or synthetic aperture radar imaging, are concerned by this problem since they proceed by recording complex or modulated data from which a “wrapped” phase is extracted. Although 1-D phase unwrapping is trivial, a challenge remains in higher dimensions to overcome two common problems: noise and discontinuities in the true phase image. In contrast to state-of-the-art techniques, this work aims at simultaneously unwrap and denoise the phase image. We propose a robust convex optimization approach that enforces data fidelity constraints expressed in the corrupted phase derivative domain while promoting a sparse phase prior. The resulting optimization problem is solved by the Chambolle-Pock primal-dual scheme. We show that under different observation noise levels, our approach compares favorably to those that perform the unwrapping and denoising in two separate steps.

The test of this script is performed on a synthetic phase image defined in a $256 \times 256$ pixel grid (see below).

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 (Gonzalez Gonzalez and Jacques 2014) allows solving the following minimization problem:

$$ [x^*,v^*] = \arg\min_{x,v} ||\Psi^\top x ||_1\quad {\rm s.t.}\quad || v ||_2 < \epsilon_1,\ || q - \nabla(x + v)||_1 < \epsilon_2,\ x(1) = 0 $$

with $y$ the observation (the wrapped phase), and $q(y) = {\rm wrap}(\nabla(y))$

This code relies on the following toolbox and references:

“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.
The PUMA algorithm provided in 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.
Rice Wavelet Toolbox available at 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},
}

Demo file

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

[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

Convex Unwrapping Matlab Demonstration files

by Adriana Gonzalez Gonzalez & Laurent Jacques