Code covered by the BSD License

# Multiresolution Gabor-like transforms

### Kunal Chaudhury (view profile)

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;

```