Code covered by the BSD License  

Highlights from
Multiresolution Gabor-like transforms

image thumbnail

Multiresolution Gabor-like transforms

by

 

11 May 2012 (Updated )

Matlab implementation of the multiresolution Gabor filters in 1 and 2 dimensions.

generateGaborWavelets.m
function   generateGaborWavelets (  ~  )
%
% generates and plots the real components of the four (from a total of six) 
% Gabor-likewavelets. 
%
fprintf('Synthesizing the four  directional wavelets...\n');

% set parameters
len   =  1024;
alpha =  12;
tau   =  0;
shift =  1/2;
J     =  6;
M     =  len / 2^J;

% get filters
[GG0, GG1] = filters(len, alpha, tau, 0);
[HH0, HH1] = filters(len, alpha, tau - shift, 0);

% set up the impulse (Kronecker delta) images
Yl_A = zeros(M, M);
LH_A = zeros(M,M);
HL_A = zeros(M,M);
HH_A = zeros(M,M);

LH_A(M/2, M/2) = 1;

Yl_B = zeros(M,M);
LH_B = zeros(M,M);
HL_B = zeros(M,M);
HH_B = zeros(M,M);

HL_B(M/2, M/2) = 1;


Yl_C = zeros(M,M);
LH_C = zeros(M,M);
HL_C = zeros(M,M);
HH_C = zeros(M,M);

HH_C(M/2, M/2) = 1;


Yl_D = zeros(M,M);
LH_D = zeros(M,M);
HL_D = zeros(M,M);
HH_D = zeros(M,M);

HH_D(M/2, M/2) = 1;


% synthesis
for depth = J : -1 : 1
    
    G0 = GG0(1 : 2^(depth-1) : length(GG0));
    G1 = GG1(1 : 2^(depth-1) : length(GG1));
    H0 = HH0(1 : 2^(depth-1) : length(HH0));
    H1 = HH1(1 : 2^(depth-1) : length(HH1));
    
    
    % channel 1
    R_A=[Yl_A  LH_A;  HL_A HH_A];
    
    for i = 1 : 2*M
        Y  = fft(reshape(R_A(1:M,i),1,M),M);
        v  = reshape(R_A((M + 1) : 2*M, i),1,M);
        V  = fft(v,M);
        Y0 = G0(1:M).*Y + G1(1:M).*V;
        Y1 = G0(M + (1:M)).*Y + G1(M + (1:M)).*V;
        Y  = [Y0  Y1];
        R_A(:,i) = reshape(real(ifft(Y,2*M)),2*M,1);
    end
    
    for i = 1 : 2*M
        Y  = fft(R_A(i,1:M),M);
        v  = R_A(i,(M + 1) : 2*M);
        V  = fft(v,M);
        Y0 = G0(1:M).*Y + G1(1:M).*V;
        Y1 = G0(M + (1:M)).*Y + G1(M + (1:M)).*V;
        Y  = [Y0  Y1];
        R_A(i,:) = real(ifft(Y,2*M));
    end
    
    
    
    % channel 2
    R_B=[Yl_B  LH_B; HL_B HH_B];
    
    for i = 1: 2*M
        Y  = fft(reshape(R_B(1:M,i),1,M),M);
        v  = reshape(R_B((M + 1) : 2*M,i),1,M);
        V  = fft(v,M);
        Y0 = G0(1:M).*Y + G1(1:M).*V;
        Y1 = G0(M + (1:M)).*Y + G1(M + (1:M)).*V;
        Y  = [Y0 Y1];
        R_B(:,i) = reshape(real(ifft(Y,2*M)),2*M,1);
    end
    
    for i = 1 : 2*M
        Y  = fft(R_B(i,1:M),M);
        v  = R_B(i,(M+1):2*M);
        V  = fft(v,M);
        Y0 = G0(1:M).*Y+G1(1:M).*V;
        Y1 = G0(M+(1:M)).*Y+G1(M+(1:M)).*V;
        Y  = [Y0 Y1];
        R_B(i, : ) = real(ifft(Y,2*M));
    end
    
    
    % channel 3
    R_C=[Yl_C  LH_C; HL_C HH_C];
    
    
    for i = 1 : 2*M
        Y  = fft(reshape(R_C(1:M,i),1,M),M);
        v  = reshape(R_C((M + 1) : 2*M,i),1,M);
        V  = fft(v,M);
        Y0 = H0(1:M).*Y + H1(1:M).*V;
        Y1 = H0(M + (1:M)).*Y + H1(M + (1:M)).*V;
        Y  = [Y0 Y1];
        R_C(:,i) = reshape(real(ifft(Y,2*M)),2*M,1);
    end
    
    for i = 1 : 2*M
        Y  = fft(R_C(i,1 : M),M);
        v  = R_C(i,(M + 1) : 2*M);
        V  = fft(v,M);
        Y0 = H0(1:M).*Y + H1(1:M).*V;
        Y1 = H0(M + (1:M)).*Y + H1(M + (1:M)).*V;
        Y  = [Y0  Y1];
        R_C(i,:) = real(ifft(Y,2*M));
    end
    
    
    % channel 4
    R_D=[Yl_D  LH_D; HL_D HH_D];
    
    
    for i = 1 : 2*M
        Y  = fft(reshape(R_D(1:M,i),1,M),M);
        v  = reshape(R_D((M + 1) : 2*M,i),1,M);
        V  = fft(v,M);
        Y0 = G0(1:M).*Y + G1(1:M).*V;
        Y1 = G0(M + (1:M)).*Y + G1(M + (1:M)).*V;
        Y  = [Y0 Y1];
        R_D(:,i) = reshape(real(ifft(Y,2*M)),2*M,1);
    end
    
    
    for i = 1 : 2*M
        Y  = fft(R_D(i,1:M),M);
        v  = R_D(i,(M + 1) : 2*M);
        V  = fft(v,M);
        Y0 = G0(1:M).*Y + G1(1:M).*V;
        Y1 = G0(M + (1:M)).*Y + G1(M + (1:M)).*V;
        Y  = [Y0  Y1];
        R_D(i,:) = real(ifft(Y,2*M));
    end
    
    
    M    = 2*M;
    Yl_A = R_A;
    
    LH_A = zeros(M,M);
    HL_A = zeros(M,M);
    HH_A = zeros(M,M);
    
    Yl_B = R_B;
    
    LH_B = zeros(M,M);
    HL_B = zeros(M,M);
    HH_B = zeros(M,M);
    
    
    Yl_C = R_C;
    
    LH_C = zeros(M,M);
    HL_C = zeros(M,M);
    HH_C = zeros(M,M);
    
    Yl_D = R_D;
    
    LH_D = zeros(M,M);
    HL_D = zeros(M,M);
    HH_D = zeros(M,M);
    
end

% four directional wavelets
D1 = R_A;
D2 = R_B;
D3 = (R_C+R_D) / sqrt(2);
D4 = (R_C-R_D) / sqrt(2);


% display
close all force; clc;
figure,
subplot(2,2,1), imagesc(D1); colormap gray; axis image; axis off;
title('Horizontal wavelet', 'fontsize', 15);
subplot(2,2,2),imagesc(D2); colormap gray; axis image;axis off;
title('Vertical wavelet', 'fontsize', 15);
subplot(2,2,3),imagesc(D3); colormap gray; axis image;axis off;
title('Diagonal wavelet 1', 'fontsize', 15);
subplot(2,2,4),imagesc(D4); colormap gray; axis image;axis off;
title('Diagonal wavelet 2', 'fontsize', 15);





Contact us