What operations/functions can I use to speed up these nested for loops?

I tried to solve my problem very sequentially. I know that I am doing this incorrectly (inefficiently), and I need a push in the right direction. If looking at this makes you cry, I apologize and promise to do better with the knowledge bestowed (hopefully) upon me after receiving help.

in the snippet below, source_r and source_c are the rows and colums of a pixel which I deem to be a 'heat source'. from each source is a normal curve radiating in all directions. A pixel some distance away is affected by this source and has a value corresponding to the normal curve.

for c = 1:100
    for r = 1:100
        for i = 1:length(source_r)
            dist = sqrt(((r - source_r(i))^2 +(c - source_c(i))^2));
            normvalue = 1-1/(sig*sqrt(2*pi()))*e^(-(sqrt(dist)-mean)^2/(2*sig^2));
            Z(c,r) = Z(c,r) + normvalue;
        end
    end
end

Works fine for 10,000 pixels, (100x100) but eventually I need to do it on 2 million. Attached below is an example of what I want to end up with.

What kind of functions should I use to simplify my loops?

[Aaron Liao]

So @DGM I think bwdist does what I want, but I'd like it for each source point as well... and then I'd like to sum up all of the values from all of the source points after. I've played around with bwdist with adjacent 0's and I'm not sure if its want...

[DGM]

Yeah. I noticed that. That's why i got rid of the comment. What values for sig and mean are you using?

[Aaron Liao]

Thought it would be useful to @ you just to have some discussion. At the moment, those are arbitrary as I really don't have a good idea of how the heat actually spreads in the sample that I'm trying to model. For this sig = 20, and mean = 10.

Answer 1

[DGM]

Try this:

A = im2double(rgb2gray(imread('eyetest.png')));

avg = 10;
sig = 20;

% create gaussian filter
r = sig*5;
xx = -r:r;
yy = xx.';
R = sqrt(xx.^2 + yy.^2);
k = 1/(sig*sqrt(2*pi));
f = 1-k*exp(-(sqrt(R)-avg).^2/(2*sig^2));
% f = sum(f(:)); % do you want to normalize?

% apply filter
Z = imfilter(A,f);
imshow(Z,[])

That should reduce exec time by about 95-98%

You can use immse() to compare the results to yours (should be near zero), though bear in mind that you'll need to fix

Z(c,r) = Z(c,r) + normvalue;

to Z(r,c) so that your image isn't transposed.

[Aaron Liao]

So I guess what I unknowingly described is a gaussian filter. I have a lot to learn. Thank you for your help!

Thank you also for the correction on the (row,column) mistake.

[DGM]

Sometimes it's easy to get lost in an idea and make it complicated.

[Aaron Liao]

Hey so turns out, what I was looking for was not the filter function. I believe that what I was describing was actually convolution with the mask being a normal curve. This makes a little more sense when considering the final value of each pixel given its surrounding values added as a function of the distance from the sources around it.

Just thought I'd share. Thanks for your help again!

[DGM]

By default, imfilter() does correlation, but it can also do convolution. In this context, the only difference between the two is a 180 degree rotation of the gaussian. Since the gaussian is rotationally symmetric, the result should be the same. IIRC, both conv2 and imfilter default to doing zero-padding for edge treatment, but imfilter has other options. For the most part, you can think of imfilter as a wrapper for conv2 with a bunch of extra convenience options.

You're right though. Each point is the integral of a product of two curves -- convolution should come to mind.