xylimeng/warp

Wavelets with adaptive recursive partitioning

WARP

Wavelets with adaptive recursive partitioning applied to image reconstruction (both 2D and 3D). Matlab users can directly use the code for real data analysis. The C++ source code is also available, which is portable toR and Matlab.

Instructions for Matlab users

• Install the c++ library armadillo(http://arma.sourceforge.net) to your computer;
• Let path be the directory containing the header file of armadillo;
• Compile the c++ code for matlab by running the code mex_this(path).

Insturctions for R users

• Please refer to https://github.com/dcbdan/WARP for the R package.

2D example

The following code is available in the file Demo.m. It usually takes 20sec if ran in a Macbook Pro with usual specifications, which increases by 1min if the step size in cycling spinning (step) is changed to 5.

rng(0);
load('lena_data.mat');
obs_true = double(lena) ./ 255; % ground truth
MSE = @(x) mean(mean((x - obs_true).^2, 1));

sigma = 0.2;
obs_raw = obs_true + randn(size(obs_true, 1)) .* sigma; % 2D observation
obs = obs_raw(:); % vectorize obs
dimension = size(obs_raw)'; % obtain size information as a column vector

tic
% default method to select parameters in the model
hyper0 = hyper_default(obs(:), dimension);
t1 = toc;
% Bayesian model averaging without cycle spinning
BMA_no_cs = treeFit(obs, dimension, hyper0, 0);
BMA_no_cs = reshape(BMA_no_cs, dimension'); 
t2 = toc;
% BMA with cycle spinning;
step = 1; % increase `step` will lead to better performance (smaller MSEs) but takes longer time. 
BMA_cs = treeFit(obs, dimension, hyper0, step);
BMA_cs = reshape(BMA_cs, dimension'); 
t3 = toc;

figure; 
subplot(2, 2, 1)
imshow(obs_true); title('True'); 
subplot(2, 2, 2)
imshow(obs_raw); 
title(sprintf('Noisy observation \n MSE = %.4f', MSE(obs_raw))); 
subplot(2, 2, 3)
imshow(BMA_no_cs); 
title(sprintf('WARPed Haar without CS \n MSE = %.4f | Time = %.1fs', MSE(BMA_no_cs), t2)); 
subplot(2, 2, 4)
imshow(BMA_cs); 
title(sprintf('WARPed Haar with %d shifts \n MSE = %.4f | Time = %.1fs', (2 * step + 1)^2, MSE(BMA_cs), t3 - t2 + t1));

Although we do not recommend it, if you'd like to use a full optimization to select hyperparameters rather than the default method implemented in hyper_default, use the following code to obtain hyper0.

idx = (1:2:numel(obs));
a = (obs(idx) - obs(idx + 1)) ./ sqrt(2);
sigma_hat = mad(a(:), 1) * 1.4826;
x0 = [0, 0, 0, 0, 0]';
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
f = @(x) -treeLikelihood(obs,[x; log(sigma_hat)]);
[x, fval] = fminunc(f,x0, options);
hyper0 = [x; log(sigma_hat)];