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.

analysis(inImg, J, alpha, tau)
function  [w, lowpass, L1, L2, L3, L4]  =  analysis(inImg, J, alpha, tau)
%
% First, project the image into the four (tensor) approximation space.
% Then a J-level decomposition is performed yielding a total of 4 x 3 = 12
% wavelet subbands.
% The subands are approapiately combined to yield six complex wavelet
% subbands corresponding to the six directionl Gabor wavelets.
%
% pre-filtering
[im1, L1] = prefilter(inImg, alpha, tau, 1);
[im2, L2] = prefilter(inImg, alpha, tau, 2);
[im3, L3] = prefilter(inImg, alpha, tau, 3);
[im4, L4] = prefilter(inImg, alpha, tau, 4);


% complex wavelet coefficients (images)
w = cell(J,1);

% analysis filters for the four channels
shift = 1/2;
M     = length(inImg);
[H1, G1] = filters(M, alpha, tau, 0);
[H2, G2] = filters(M, alpha, tau + shift, 0);

H1 = conj(H1);
G1 = conj(G1);
H2 = conj(H2);
G2 = conj(G2);

for depth = 1 : J,
    
    % channel 1
    for i = 1 : M
        X = fft(im1(i, : ),M);
        Y = H1.*X;
        Z = G1.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im1(i, : ) = [ifft(Y,M/2)  z];
    end
    
    for i = 1 : M
        X = fft(im1( : ,i).');
        Y = H1.*X;
        Z = G1.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im1( : ,i) = [ifft(Y,M/2)  z];
    end
    
    
    % channel 2
    for i = 1 : M
        X = fft(im2(i, : ),M);
        Y = H2.*X;
        Z = G2.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im2(i, : ) = [ifft(Y,M/2)  z];
    end
    
    for i = 1 : M
        X = fft(im2( : ,i).');
        Y = H2.*X;
        Z = G2.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im2( : ,i) = [ifft(Y,M/2)  z];
    end
    
 
    % channel 3
    for i = 1 : M
        X = fft(im3(i,:),M);
        Y = H1.*X;
        Z = G1.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im3(i, : ) = [ ifft(Y,M/2) z ];
    end
    for i = 1 : M
        X = fft(im3( : ,i).');
        Y = H2.*X;
        Z = G2.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im3( : ,i) = [ ifft(Y,M/2) z ];
    end
    
    
    % channel 4
    for i = 1 : M
        X = fft(im4(i, : ),M);
        Y = H2.*X;
        Z = G2.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im4(i, : )  =  [ifft(Y,M/2)  z];
    end
    for i = 1 : M
        X = fft(im4(:,i).');
        Y = H1.*X;
        Z = G1.*X;
        Y = 1/2 * (Y(1 : M/2) + Y(M/2 + (1 : M/2)));
        Z = 1/2 * (Z(1 : M/2) + Z(M/2 + (1 : M/2)));
        z = ifft(Z,M/2);
        im4( : ,i) = [ifft(Y,M/2)  z];
    end
    
    
    % complex wavelet coefficients 
    w{depth}(:, :, 1) = im1(1 : M/2,M/2 + (1:M/2)) + ...
        sqrt(-1) * im4(1 : M/2, M/2 + (1:M/2));
    w{depth}(:, :, 2) = im3(1 : M/2,M/2 + (1:M/2)) + ...
        sqrt(-1) * im2(1 : M/2, M/2 + (1:M/2));
    w{depth}(:, :, 3) = im1(M/2 + (1:M/2),1 : M/2) + ...
        sqrt(-1) * im3(M/2 + (1 : M/2),1 : M/2);
    w{depth}(:, :, 4) = im4(M/2 + (1:M/2),1 : M/2) + ...
        sqrt(-1) * im2(M/2 + (1:M/2),1 : M/2);
    w{depth}(:, :, 5) = (im1(M/2 + (1:M/2),M/2 + (1:M/2)) ...
        - im2(M/2 + (1:M/2),M/2 + (1:M/2)))/sqrt(2) ...
        + sqrt(-1) * (im3(M/2 + (1:M/2),M/2 + (1:M/2)) ...
        + im4(M/2 + (1:M/2),M/2 + (1:M/2)))/sqrt(2);
    w{depth}(:, :, 6) = (im1(M/2 + (1:M/2),M/2 + (1:M/2)) ...
        + im2(M/2 + (1:M/2),M/2 + (1:M/2)))/sqrt(2) ...
        + sqrt(-1) * (im3(M/2 + (1:M/2),M/2 + (1:M/2)) ...
        - im4(M/2 + (1:M/2),M/2 + (1:M/2)))/sqrt(2);
    
    
    im1 = im1(1 : M/2,1 : M/2);
    im2 = im2(1 : M/2,1 : M/2);
    im3 = im3(1 : M/2,1 : M/2);
    im4 = im4(1 : M/2,1 : M/2);
    
    M = M/2;
    H1 = H1(1 : 2 : length(H1));
    G1 = G1(1 : 2 : length(G1));
    H2 = H2(1 : 2 : length(H2));
    G2 = G2(1 : 2 : length(G2));
    
end


% lowpass subbands
[p, q]  = size(im1);
lowpass = zeros(p, q, 4);
lowpass(:,:,1) = im1;
lowpass(:,:,2) = im2;
lowpass(:,:,3) = im3;
lowpass(:,:,4) = im4;

Contact us