How do I replace my for loops using matrix manipulation?

I want a more efficient way to do a window search between two images. I currently have two images; a small clip, and a larger image the clip should be a part of. I want to compare the small clip to every possible subset of the larger image to build a map of mutual information values. I can do this with a for loop as shown below.

    dim_Large = size(Large); 
    dim_Clip = size(Clip);
    for j = 1:dim_Large(1)-dim_Clip(1)
        for k = 1:dim_Large(2)-dim_Clip(2)
            map(j,k) = ...
                ent(Large(j:j+dim_Clip(1)-1,k:k+dim_Clip(2)-1),Clip);
        end 
    end

The 'ent' function is code to calculate the mutual information between two matrices of the same size. This does what I need it to do, but is very slow.

I'm sure there is a better way to write this for MatLab, something like

    map = ent( Large(1:dim_Large(1)-dim_Clip(1)-1,1:dim_Large(2)-dim_Clip(2)-1), Clip);

but I can't get it to work.

Is there a more efficient way to write this code, preferably one that uses no for loops?

[Rik]

This depends entirely on the contents of the ent function. Without it, it is impossible to answer this question.

Answer 1

[Jan]

The speed of your code depends on two points:

1. The speed of the function ent()
2. The number of times this function is called

There is no magic way to reduce one or both.

Use the profiler to find the bottleneck of the code. Then start with improving the corresponding code lines.

But I assume, a real speed-up needs a vectorization, e.g.

ent(Large(:, 1:dim_Large(2)-dim_Clip(2)-1), Clip)

Now the complete columns of Large are used and one loop can vanish. If this is possible depends on the contents of ent(), and it is not guaranteed, that it is faster.

[Joshua Knicely]

Thank you for the quick response and information, Jan. The code is from the MATLAB file exchange https://www.mathworks.com/matlabcentral/fileexchange/41714-fast-mutual-information--joint-entropy--and-joint-histogram-calculation-for-n-d-images?focused=3800593&tab=function

I tried your above vectorization, but it seems the code wasn't built with that kind of behavior in mind and it resulted in an error.

    Matrix dimensions must agree. 
    Error in ent (line 18)
    xx = xx + yy;

I'm making some minor edits to the code to speed it up, but I think I'm going to have to deal with the long computation time.

[Jan]

If you only need the first output of ent(), omit the expensive calculation of the 2 others. By the way: x.^-1 is much more expensive than 1./x.

Please use tic/toc to find if this improves the speed:

function jhist = ent2(J, K)
dimen = 256;
x = numel(J);
t = 1:x;
xx = J(:)+1;
yy = dimen*K(:);
xx = sort(xx + yy);
yy(1:x-1) = xx(2:x);
zz = yy - xx;
zz(x) = 1;
zz = t(zz ~=0);
yy = xx(zz);
t = numel(zz);
zz(2:t) = zz(2:t)-zz(1:t-1);  % Better: diff(zz)
jhist = zeros(dimen);
jhist(yy) = zz / x;
end

[Joshua Knicely]

I need the output MI from ent( ). I removed the other two outputs and commented out the portion that calculates the joint entropy to speed it up a bit.

Your suggested change of x.^-1 to 1./x sped up the code. It went from 0.624 seconds per 100 runs to 0.569 seconds per 100 runs, about a 10% increase in speed, so thanks for that.

Change 1 [xx = xx+yy; xx = sort(xx); -> xx = sort(xx+yy)] actually resulted in slowing down the code very slightly, which seems odd. It went from 0.556 seconds per 100 runs to 0.560 seconds per 100 runs. Less than 1%, so I imagine it doesn't matter.

Change 2 [zz(2:t)-zz(1:t-1) -> diff(zz)] increased the speed by a tiny margin (<1%).

Change 3 [ xx = zz/x; jhist(yy) = xx; -> jhist(yy) = zz/x;] increased run speed by about 2%.

Thank you very much for the extremely specific ways in which to speed up this code.