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 = imread('tumor-1.jpg');
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 :
Okay, I looked into the problem a bit deeper. 


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. 
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.
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.
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 = imread('tumor-1.jpg');
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(:) &gt; 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(:) &gt; 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(:) &gt; 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 &lt; 5*eps))
all(all(timgD_LH - imgD_LH &lt; 5*eps))
all(all(timgD_HL - imgD_HL &lt; 5*eps))
all(all(timgD_HH - imgD_HH &lt; 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');