# Wavelet decomposition of images

I am working on wavelet decomposition and found a formula for it and wrote it by myself for computation.

The formula is:

$$I_{LH}(i,j) = \sum_{p=1}^{N_L} \sum_{q=1}^{N_H}\ L(p)\ H(q)\ I(i+p, j+q)$$

This is LH filtered and I do it for LL, LH, HL and HH so i get 4 images.

Can someone check the for loop in it and help me improve the computation time of it?

Elapsed time is 10.040462 seconds. I want to reduce time as much as possible as I need to work on hundreds of images.

%   Wavelet Decomposition
%Decomposition wavelet 'coif1' filters
[Lo_D,Hi_D] = wfilters('coif1','d');
img=  im2double(img);
figure, imshow(img), title('original image');
[nrows ncols] = size(img);
imgg = zeros(nrows+6,ncols+6);
imgg(1:nrows,1:ncols) = img(:,:);
figure, imshow(imgg);

% wavelet decomposition.
imgD_LL = zeros(size(img));
imgD_LH = zeros(size(img));
imgD_HL = zeros(size(img));
imgD_HH = zeros(size(img));

% Deconstruction
tic
for i = 1: nrows
for j = 1:ncols
for p = 1:numel(Lo_D)
for q = 1:numel(Hi_D)
tempqD_LL(q) = Lo_D(p)*Lo_D(q)*imgg(i+p,j+q);
tempqD_LH(q) = Lo_D(p)*Hi_D(q)*imgg(i+p,j+q);
tempqD_HL(q) = Hi_D(p)*Lo_D(q)*imgg(i+p,j+q);
tempqD_HH(q) = Hi_D(p)*Hi_D(q)*imgg(i+p,j+q);
end
if (length(tempqD_LL) == q)
tempQD_LL(p) = sum(tempqD_LL);
tempQD_LH(p) = sum(tempqD_LH);
tempQD_HL(p) = sum(tempqD_HL);
tempQD_HH(p) = sum(tempqD_HH);
end

end
if (length(tempQD_LL) == p )
imgD_LL(i,j) = sum(tempQD_LL);
imgD_LH(i,j) = sum(tempQD_LH);
imgD_HL(i,j) = sum(tempQD_HL);
imgD_HH(i,j) = sum(tempQD_HH);
end
end
end
toc
% Plots Deconstruction..
figure,
subplot(2,2,1), imshow(imgD_LL, [min(min(imgD_LL)) max(max(imgD_LL))]), title('Img Deconstruction - LL filter');
subplot(2,2,2), imshow(imgD_LH, [min(min(imgD_LH)) max(max(imgD_LH))]); title('Img Deconstruction - LH filter');
subplot(2,2,3), imshow(imgD_HL, [min(min(imgD_HL)) max(max(imgD_HL))]); title('Img Deconstruction - HL filter');
subplot(2,2,4), imshow(imgD_HH, [min(min(imgD_HH)) max(max(imgD_HH))]); title('Img Deconstruction - HH filter');


The image i used for above is : • I know nothing about wavelet decomposition, but does the wavedec2 function not do what you want? – 200_success Sep 21 '16 at 16:58
• It does similar thing but not the one i want.there are various once but nothing like the one i wrote above,,,All i need to do is improve computation time of the for loop..which is quite time taking – Gopi Sep 21 '16 at 17:54
• The computation of Lo_D(p)*Lo_D(q), Lo_D(p)*Hi_D(q), etc.. can be done outside the 4 loops right? Since it only depends on p and q, which are the length of Lo_D and Hi_D. So you are doing multiple multiplications already numerous times. So just create 4 new variables which holds the result of Lo_D(p)*Lo_D(q), Lo_D(p)*Hi_D(q), etc... In addition it is better to already allocate space for the variables tempqD_LL, tempqD_LH, etc... now the variables are changing in size dynamically. – WG- Sep 22 '16 at 8:24

Okay, I looked into the problem a bit deeper.

1. Currently, with your use case, you are doing a lot of multiplications which end up in zero anyway. Most of the image is namely black. You can easily skip all these multiplications by checking if the image block you are looping over is bigger than zero.
2. The two loops basically do a matrix multiplication. You are multiplying an image block with two vector values over which you loop. You can create a matrix which has all the possible combinations of Lo_D(p)*Lo_D(q) by taking the Kronecker product of Lo_D. Then you can simply do a matrix multiplication. and sum all elements over the matrix. However, it seems that this optimization is actually slower then the first one in this case. It might be the case that it is faster in other cases.
3. You are calculating each index of an array each time and afterwards you are summing it, this can be done in one line as well. The performance is however does not seem to increase or decrease.
4. Finally, you can do one minor optimization the for loop where you loop over Lo_D, you can also there check if an entire row is equal to zero. This is a small optimization but might give you something. Doing this another time in the next loop is useless.

Full code:

clear all;
close all;
clc

[Lo_D,Hi_D] = wfilters('coif1','d');

img = im2double(img);
figure;
imshow(img);
title('original image');

[nrows ncols] = size(img);
imgg = zeros(nrows+6,ncols+6);
imgg(1:nrows,1:ncols) = img(:,:);
figure;
imshow(imgg);

% Wavelet decomposition.
imgD_LL = zeros(size(img));
imgD_LH = zeros(size(img));
imgD_HL = zeros(size(img));
imgD_HH = zeros(size(img));

timgD_LL = zeros(size(img));
timgD_LH = zeros(size(img));
timgD_HL = zeros(size(img));
timgD_HH = zeros(size(img));

%% Optimization #1
tic
for i = 1:nrows
for j = 1:ncols

imgBlock = imgg((i+1):(i+numel(Lo_D)),(j+1):(j+numel(Hi_D)));

if any(imgBlock(:) > 0)

for p = 1:numel(Lo_D)
for q = 1:numel(Hi_D)
tempqD_LL(q) = Lo_D(p)*Lo_D(q)*imgg(i+p,j+q);
tempqD_LH(q) = Lo_D(p)*Hi_D(q)*imgg(i+p,j+q);
tempqD_HL(q) = Hi_D(p)*Lo_D(q)*imgg(i+p,j+q);
tempqD_HH(q) = Hi_D(p)*Hi_D(q)*imgg(i+p,j+q);
end
if (length(tempqD_LL) == q)
tempQD_LL(p) = sum(tempqD_LL);
tempQD_LH(p) = sum(tempqD_LH);
tempQD_HL(p) = sum(tempqD_HL);
tempQD_HH(p) = sum(tempqD_HH);
end

end
if (length(tempQD_LL) == p)
timgD_LL(i,j) = sum(tempQD_LL);
timgD_LH(i,j) = sum(tempQD_LH);
timgD_HL(i,j) = sum(tempQD_HL);
timgD_HH(i,j) = sum(tempQD_HH);
end
end
end
end
toc

%% Optimization #2
% % Kronecker product
% LL = kron(Lo_D,Lo_D);
% LH = kron(Hi_D,Lo_D);
% HL = kron(Lo_D,Hi_D);
% HH = kron(Hi_D,Hi_D);
%
% % Reshape to matrix
% LL = reshape(LL,numel(Lo_D)*[1 1]);
% LH = reshape(LH,numel(Lo_D)*[1 1]);
% HL = reshape(HL,numel(Lo_D)*[1 1]);
% HH = reshape(HH,numel(Lo_D)*[1 1]);
%
% tic
% for i = 1:nrows
%     for j = 1:ncols
%
%         imgBlock = imgg((i+1):(i+numel(Lo_D)),(j+1):(j+numel(Hi_D)));
%
%         if any(imgBlock(:) > 0)
%             timgD_LL(i,j) = sum(sum(LL.*imgBlock));
%             timgD_LH(i,j) = sum(sum(LH.*imgBlock));
%             timgD_HL(i,j) = sum(sum(HL.*imgBlock));
%             timgD_HH(i,j) = sum(sum(HH.*imgBlock));
%         else
%             timgD_LL(i,j) = 0;
%             timgD_LH(i,j) = 0;
%             timgD_HL(i,j) = 0;
%             timgD_HH(i,j) = 0;
%         end
%     end
% end
% toc
%
%% Optimization #3
% tic
% for i = 1:nrows
%     for j = 1:ncols
%         tLL = 0;
%         tLH = 0;
%         tHL = 0;
%         tHH = 0;
%
%         imgBlock = imgg((i+1):(i+numel(Lo_D)),(j+1):(j+numel(Hi_D)));
%
%         if any(imgBlock(:) > 0)
%             for p = 1:numel(Lo_D)
%                 for q = 1:numel(Hi_D)
%                     tLL = Lo_D(p)*Lo_D(q)*imgg(i+p,j+q) + tLL;
%                     tLH = Lo_D(p)*Hi_D(q)*imgg(i+p,j+q) + tLH;
%                     tHL = Hi_D(p)*Lo_D(q)*imgg(i+p,j+q) + tHL;
%                     tHH = Hi_D(p)*Hi_D(q)*imgg(i+p,j+q) + tHH;
%                 end
%             end
%         end
%
%         timgD_LL(i,j) = tLL;
%         timgD_LH(i,j) = tLH;
%         timgD_HL(i,j) = tHL;
%         timgD_HH(i,j) = tHH;
%     end
% end
% toc
%
%% Original implementation
tic
for i = 1:nrows
for j = 1:ncols
for p = 1:numel(Lo_D)
for q = 1:numel(Hi_D)
tempqD_LL(q) = Lo_D(p)*Lo_D(q)*imgg(i+p,j+q);
tempqD_LH(q) = Lo_D(p)*Hi_D(q)*imgg(i+p,j+q);
tempqD_HL(q) = Hi_D(p)*Lo_D(q)*imgg(i+p,j+q);
tempqD_HH(q) = Hi_D(p)*Hi_D(q)*imgg(i+p,j+q);
end
if (length(tempqD_LL) == q)
tempQD_LL(p) = sum(tempqD_LL);
tempQD_LH(p) = sum(tempqD_LH);
tempQD_HL(p) = sum(tempqD_HL);
tempQD_HH(p) = sum(tempqD_HH);
end

end
if (length(tempQD_LL) == p)
imgD_LL(i,j) = sum(tempQD_LL);
imgD_LH(i,j) = sum(tempQD_LH);
imgD_HL(i,j) = sum(tempQD_HL);
imgD_HH(i,j) = sum(tempQD_HH);
end
end
end
toc

%%
% Check if they are the same, calculate difference and check if it is
% smaller then 5 machine epsilon because of numerical differences
all(all(timgD_LL - imgD_LL < 5*eps))
all(all(timgD_LH - imgD_LH < 5*eps))
all(all(timgD_HL - imgD_HL < 5*eps))
all(all(timgD_HH - imgD_HH < 5*eps))

% Plots Deconstruction..
figure;
subplot(2,2,1), imshow(timgD_LL, [min(min(timgD_LL)) max(max(timgD_LL))]), title('Img Deconstruction - LL filter');
subplot(2,2,2), imshow(timgD_LH, [min(min(timgD_LH)) max(max(timgD_LH))]); title('Img Deconstruction - LH filter');
subplot(2,2,3), imshow(timgD_HL, [min(min(timgD_HL)) max(max(timgD_HL))]); title('Img Deconstruction - HL filter');
subplot(2,2,4), imshow(timgD_HH, [min(min(timgD_HH)) max(max(timgD_HH))]); title('Img Deconstruction - HH filter');
figure;
subplot(2,2,1), imshow(imgD_LL, [min(min(imgD_LL)) max(max(imgD_LL))]), title('Img Deconstruction - LL filter');
subplot(2,2,2), imshow(imgD_LH, [min(min(imgD_LH)) max(max(imgD_LH))]); title('Img Deconstruction - LH filter');
subplot(2,2,3), imshow(imgD_HL, [min(min(imgD_HL)) max(max(imgD_HL))]); title('Img Deconstruction - HL filter');
subplot(2,2,4), imshow(imgD_HH, [min(min(imgD_HH)) max(max(imgD_HH))]); title('Img Deconstruction - HH filter');

• It was really useful and time saving the optimization #1 you gave worked perfectly other two didn't work as the same and images are not identical... I tried it for a whole image which has more pixels which are not zero and the computation time was improved... Orignal Code: Elapsed time is 9.765361 seconds. and Optimization code #1: Elapsed time is 6.419653 seconds. and for image like i posted with question it was Elapsed time is 0.714627 seconds... Thank You.. – Gopi Sep 22 '16 at 15:24
• If you have time can you explain the usefulness of this code you wrote: 'about smaller than 5 machine epsilon' I am thinking you are checking all pixel values in both and check if they are close enough all your codes agree with this but when i see ' isequal(timgD_LL,imgD_LL)' for optimization #2 and #3 its not same..for #1 its perfect... can you explain this..will it make a difference if i use #2 or #3 ,becuase Elapsed time is 1.598046 seconds. and Elapsed time is 1.860081 seconds. for #2 and #3 repectively which is really huge improvement compared to orignal and opt#1.. – Gopi Sep 22 '16 at 15:49
• The machine epsilon is the upper bound on the relative error due to rounding in floating point arithmetic. In my case eps (you can type this in Matlab) has the value of 2.2204e-16 ≈ 2^(−52). Because the calculation is different the numerical error will propagate differently. Therefor the very digits at the end of the number will be different and you cannot compare the two images by isequal(timgD_LL,imgD_LL). Because of this I subtract them from each other and check that the difference between the two is small enough, then the two images are essentially equal to each other. – WG- Sep 23 '16 at 8:49
• Personally I would use #2 because in my opinion it is cleaner code and also more clear what is happening. But I can imagine that for people who do not do a lot of matrix calculations it can be harder to understand. In addition, I am still a bit confused why the matrix calculation seems to be a bit slower then the implementation in which you have the two additional for loops. Normally matrix calculations are faster. – WG- Sep 23 '16 at 8:53
• I like the Kronecker product method too..Thanks for your explanation – Gopi Sep 23 '16 at 15:42