demo_videodenoising.m

A set of C++ and Matlab routines implementing the surfacelet transform surfacelet的一个非常好用的工具箱

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%	SurfBox-MATLAB (c)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	Yue Lu and Minh N. Do
%%
%%	Department of Electrical and Computer Engineering
%%	Coordinated Science Laboratory
%%	University of Illinois at Urbana-Champaign
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%	demo_DenoisingVideo.m
%%	
%%	First created: 04-22-06
%%	Last modified: 04-23-06
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ');
disp(' ');
disp('In this demo, we apply the 3-D surfacelet transform to remove additive');
disp('white Gaussian noise added to a video sequence. This is achieved by simply');
disp('applying hard-thresholding on the surfacelet coefficients.');
disp(' ');
disp('Note: This demo is memory intensive, and might need 1.5 GB of memory.');
disp('A quick workaround is to start Matlab without the Java virtual machine,');
disp('i.e., use Matlab -nojvm');
disp(' ');
disp(['Make sure the Matlab directory of SurfBox is on the search path.']);
disp(' ');

%% We first load a video sequence "Mobile and Calendar"
disp(' ');
disp('Step 1: Load a video sequence "Mobile and Calendar"');
load mobile; %% An array X of size 192 * 192 * 192 have been loaded
r = input('Press <enter> to play the original sequence ...');
PlayImageSequence(X);

%% We add Gaussian noise to the video sequence
disp(' ');
disp('Step 2: Add white Gaussian noise to the sequence.');
sigma = 40;
Xn = double(X) + sigma * randn(size(X));
clear X
r = input('Press <enter> to play the noisy sequence ...');
PlayImageSequence(uint8(Xn));

%% Surfacelet Denoising

%% Some parameters
Pyr_mode = 1.5; %% Note: better performance can be achieved if we use Pyr_mode = 1,
%% but that requires even more memory.
bo = 8;
Level_64 = [-1 3 3; 3 -1 3; 3 3 -1];
Level_16 = [-1 2 2; 2 -1 2; 2 2 -1];
Level_4 =  [-1 1 1; 1 -1 1; 1 1 -1];
Lev_array = {Level_64, Level_64, Level_16, Level_4};

%% Note: different surfacelets have slightly different L^2 different norm
%% values. For best performance, we need to fist estimate these norm
%% values. The estimation can then be saved and used by future runs of the
%% denoising routine.

%% Note: this step may take a while.
disp(' ');
disp('Step 3: Estimate the L^2 norm of surfacelets. This may take a while.');
disp('Note: In real applicaitons, we can skip this step by using pre-computed estimates.');
r = input('Press <enter> to continue ...');
h = waitbar(0, 'Estimating the L^2 norm of surfacelets. Please wait ...');
drawnow;
niter = 30;
for n = 1 : niter
    Noise = randn(size(Xn));
    
    [Y, Recinfo] = Surfdec(Noise, Pyr_mode, Lev_array, 'ritf', 'bo', bo);
    clear Noise;
    if n == 1
        Ep = Y;
        for j = 1 : length(Y) - 1
            for dd = 1 : length(Y{j})
                for sd = 1 : length(Y{j}{dd})
                    Ep{j}{dd}{sd} = Ep{j}{dd}{sd} .^ 2;
                end
            end
        end
        
    else
        for j = 1 : length(Y) - 1
            for dd = 1 : length(Y{j})
                for sd = 1 : length(Y{j}{dd})
                    Ep{j}{dd}{sd} = Ep{j}{dd}{sd} + Y{j}{dd}{sd} .^ 2;
                end
            end
        end
    end
    
    clear Y;
    
    waitbar(n / niter, h);
    drawnow;
end

for j = 1 : length(Ep) - 1
    for dd = 1 : length(Ep{j})
        for sd = 1 : length(Ep{j}{dd})
            Ep{j}{dd}{sd} = sqrt(Ep{j}{dd}{sd} ./ (niter - 1));
        end
    end
end

close(h);
drawnow;

%% Take the surfacelet transform on the noisy sequence
disp(' ');
disp('Step 4: Apply surfacelet transform on the noisy sequence and hard-threshold the coefficients');
r = input('Press <enter> to continue ...');

h = waitbar(0, 'Please wait ...');
drawnow;
c = 1;
for m = 1 : 2
    if (m == 2)
        %% Optional step: apply the same algorithm on a shifted volume to
        %% improve denoising performance.
        Xn = circshift(Xn, [1 1 1]);
    end
    [Y, Recinfo] = Surfdec(Xn, Pyr_mode, Lev_array, 'ritf', 'bo', bo);
    
    %% applying hard-thresholding
    for n = 1 : length(Y) - 1
        thresh = 3 * sigma + (n==1) * sigma;
        for dd = 1 : length(Y{n})
            for sd = 1 : length(Y{n}{dd})
                Y{n}{dd}{sd} = Y{n}{dd}{sd} .* (abs(Y{n}{dd}{sd}) > thresh * Ep{n}{dd}{sd});
            end
        end
    end
    
    waitbar(c / 4, h);
    c = c + 1;
    drawnow;

    %% surfacelet reconstruction
    if (m == 1)
        Xd1 = Surfrec(Y, Recinfo);
    else
        Xd2 = Surfrec(Y, Recinfo);
        Xd2 = circshift(Xd2, [-1 -1 -1]);
    end
        
    waitbar(c / 4, h);
    c = c + 1;
    drawnow;
end

close(h);
drawnow;

clear Ep;

%% Take the avarage of two different copies of the denoised signal
Xd = (Xd1 + Xd2) / 2.0;
clear Xd1 Xd2;

disp(' ');
disp('Done.');
r = input('Press <enter> to show the denoised sequence ...');
PlayImageSequence(uint8(Xd(:,:, 10 : end - 10)));

%%	This software is provided "as-is", without any express or implied
%%	warranty. In no event will the authors be held liable for any 
%%	damages arising from the use of this software.