Edges with wavelets


This is an example taken from my signal processing notes. If we apply k times certain wavelet transform to an image NxN with N a power of 2, then the part not including details, the approximation part, concentrates in the N/2kxN/2k left upper corner. If we set it to zero and invert the process then we only preserve the details to some level and we get a primitive edge detector.
If the number of iterations is large we enter into a coarse scale. In this case a curious metallic effect is obtained.

To make the images more appealing it is considered the negative of the absolute value of the result. In this way the background color corresponds to the maximal value and gives white.

The octave code is included below.


armour

Original image


ed_d4_2

D4 with k=2


ed_d4_3

D4 with k=3


ed_d6_3

D6 with k=3


ed_d8_7

D8 with k=7


ed_d12_6

D12 with k=6


ed_d4_6

D4 with k=6


ed_d6_4

D6 with k=4




The code

This is the Octave code for the images:


pkg load image

clear all

name = './armour_head.png'

% original image
ima = imread( name );
ima = im2double(ima);
N = size(ima,1)


% D4
fd4 = [ 1-sqrt(3), sqrt(3)-3, 3+sqrt(3), -1-sqrt(3)]/4/sqrt(2);

% D6
fd6 = [0.47046721, 1.14111692 , 0.650365, -0.19093442, -0.12083221, 0.0498175 ];
fd6 = -fliplr(fd6)/sqrt(2).*(-1).^[1:size(fd6,2)];

% D8
fd8 = [0.32580343 , 1.01094572 , 0.89220014, -0.03957503 , -0.26450717 , 0.0436163 , 0.0465036, -0.01498699];
fd8 = -fliplr(fd8)/sqrt(2).*(-1).^[1:size(fd8,2)];

% D12
fd12 = [0.15774243 , 0.69950381 ,1.06226376 , 0.44583132, -0.31998660, -0.18351806 , 0.13788809 , 0.03892321 , -0.04466375 , 7.83251152e-4, 6.75606236e-3, -1.52353381e-3 ];
fd12 = -fliplr(fd12)/sqrt(2).*(-1).^[1:size(fd12,2)];

% Haar
haar = [1,-1]/sqrt(2);



%%%%%%%%%%%%%%%%%%%%%%%
% D4 k=2
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 2;
f = fd4;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = im2double(mat2gray( -abs(wav_im_i(temp,f,k)) ));
imwrite(temp,'./ed_d4_2.png');


%%%%%%%%%%%%%%%%%%%%%%%
% D4 k=3
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 3;
f = fd4;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = im2double(mat2gray( -abs(wav_im_i(temp,f,k)) ));
imwrite(temp,'./ed_d4_3.png');


%%%%%%%%%%%%%%%%%%%%%%%
% D6 k=3
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 3;
f = fd6;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = im2double(mat2gray( -abs(wav_im_i(temp,f,k)) ));
imwrite(temp,'./ed_d6_3.png');



%%%%%%%%%%%%%%%%%%%%%%%
% D8 k=7 (metallic, maximal)
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 7;
f = fd8;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = im2double(mat2gray( -abs(wav_im_i(temp,f,k)) ));
imwrite(temp,'./ed_d8_7.png');

  
%%%%%%%%%%%%%%%%%%%%%%%
% D12 k=6 (metallic, maximal)
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 6;
f = fd12;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = im2double(mat2gray( -abs(wav_im_i(temp,f,k)) ));
imwrite(temp,'./ed_d12_6.png');


  
%%%%%%%%%%%%%%%%%%%%%%%
% D4 k=6
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 6;
f = fd4;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = wav_im_i(temp,f,k);
imwrite(temp,'./ed_d4_6.png');

  
%%%%%%%%%%%%%%%%%%%%%%%
% D6 k=4
%%%%%%%%%%%%%%%%%%%%%%%

% parameters
k = 4;
f = fd6;

% algorithm
temp = wav_im_t(ima,f,k);
temp(1:N/2^k,1:N/2^k) = zeros(N/2^k);
temp = im2double(mat2gray( -abs(wav_im_i(temp,f,k)) ));
imwrite(temp,'./ed_d6_4.png');

  

It calls to three functions:
File wav_im_i.m File wav_im_t.m File wav_t_n.m
% k iterations of the inverse wavelet transform associated to a filter f


function y = wav_im_i(ima,f,k)

N = size(ima,1);

Ntemp = N/2^(k-1);
Atemp = ima;

for j = 1:k
  T_N = wav_t_n( f, Ntemp);
  Atemp(1:Ntemp, 1:Ntemp) = T_N' * Atemp(1:Ntemp, 1:Ntemp) *T_N;
  Ntemp = 2*Ntemp;
endfor

y = Atemp;

end
% k iterations of the wavelet transform associated to a filter f


function y = wav_im_t(ima,f,k)

N = size(ima,1);

Ntemp = N;
Atemp = ima;

for j = 1:k
  T_N = wav_t_n( f, Ntemp);
  Atemp(1:Ntemp, 1:Ntemp) = T_N * Atemp(1:Ntemp, 1:Ntemp) *T_N';
  Ntemp = Ntemp/2;
endfor

y = Atemp;

end
% Wavelet transform: f=filter nx=dimension (power of two)

function y = wav_t_n(f,nx)
  
  Lf = size(f,2);
  
  ro = zeros(1, nx);
  ro(1:Lf) = f;
  
  F = ro;
  for j= 2:nx/2
    F = [F; circshift(ro, [0,2*j-2])];
  endfor
  
  
  f = fliplr(f).*(-1).^[1:Lf];
  ro(1:Lf) = f;
  
  Fp = ro;
  for j= 2:nx/2
    Fp = [Fp; circshift(ro, [0,2*j-2])];
  end

  
  y = [Fp;F];
  
end