Vectorizing a For-Loop

7 Ansichten (letzte 30 Tage)
Pyroholic
Pyroholic am 5 Jul. 2015
Kommentiert: Pyroholic am 6 Jul. 2015
So I have the following code:
[x y] = meshgrid(1:size(A,1), 1:size(A,2));
for i=1:n
for j=1:m
disk_indices = sqrt((x-disks_center_coordinates(1,j)).^2+...
(y-disks_center_coordinates(2,j)).^2) <= disks_radii(j);
new_matrix(i,j) = sum(sum((A(disk_indices)))) ;
end
end
Here, assume that disks_center_coordinates, and disks_radii are available somewhere.
As you can see, what I am trying to do here is to define circular regions (disks), sum all the numbers in the matrix inside those disks, and store them in the new_matrix.
Of course, this code consumes a lot of time. When I tested the code on 2000 disks, it took 7 minutes for the code to do the job.
I know that a vectorized implementation will save a lot of time, but I need some guidance on how to vectorize such code. (what I mean by vectorization is the use of vectors instead of for loops)
I will be glad to clarify any ambiguities.

Antworten (2)

Jan
Jan am 5 Jul. 2015
The disk_indices do not depend on the outer loop over i. This is strange, but if this is intented you can move the expensive calculation before the "for j" loop.
sqrt is very expensive. Avoid it by squaring both sides:
((x-disks_center_coordinates(1,j)).^2 + ...
(y-disks_center_coordinates(2,j)).^2) <= disks_radii(j)^2
You blow up the vectors 1:size(A,1), 1:size(A,2) at first by meshgrid. But it would be much cheaper to keep the vectors instead of creating two matrices.
x = 1:size(A,1);
y = 1:size(A,2);
Now bsxfun allows to calculate the (squared !) distances much cheaper. I cannot suggest an explicit solution, because I'm convinced that the posted code is buggy due to the missing dependency to i .
Because disk_indices is a logical vector, A(disk_indices) is a vector also. Then one sum command is sufficient.
  1 Kommentar
Pyroholic
Pyroholic am 6 Jul. 2015
Hello, Jan. Thank you for your time.
Regarding the non-dependency issue, I totally missed that. Thank you for pointing it out. I also got rid of the sqrt as per your suggestion. I also took the suggestion of Geoff regarding the counter being 'i' and 'j' and changed it to 'b' and 'a'.This is what I have now:
[x y] = meshgrid(1:size(A,1), 1:size(A,2));
for a=1:m
disk_indices = (x-disks_center_coordinates(1,a)).^2+...
(y-disks_center_coordinates(2,a)).^2 <= disks_radii(a)^2;
for b=1:n
HSIZE = 2 * ceil(2 * SIGMA(b-1)) + 1 ;
gaussian_blur_mask = fspecial('gaussian',HSIZE,SIGMA(b-1)) ;
A = imfilter(A,gaussian_blur_mask) ;
new_matrix(b,a) = sum(A(disk_indices)) ;
end
end
Regarding the bsxfun implementation, it cannot be done in this case as the dimensions of the input matrices do not match. This is why I am having trouble vectorizing the loops: the unmatching dimensions.

Melden Sie sich an, um zu kommentieren.


Geoff Hayes
Geoff Hayes am 5 Jul. 2015
Pyroholic - unless there is a typo in your code, it would appear that the dependence in the inner for loop on the index variable i (from the outer for loop) is only for updating new_matrix. So, for example, new_matrix(1,j) is identical to new_matrix(2,j) for all j. Is this the case? Because if so, then you could reduce your code to
new_matrix = zeros(n,m); % pre-size the matrix
for k=1:m
disk_indices = sqrt((x-disks_center_coordinates(1,k)).^2+...
(y-disks_center_coordinates(2,k)).^2) <= disks_radii(k);
new_matrix(:,k) = sum(sum((A(disk_indices)))) ;
end
Note the use of k instead of i or j as MATLAB uses these last two as representations of the imaginary number. Also we pre-size the new_matrix to avoid the resizing (growth) of the matrix on each iteration of the for loop.
The next step is to try to simplify the above (single) loop.
  1 Kommentar
Pyroholic
Pyroholic am 6 Jul. 2015
Hello, Geoff. Thank you for taking the time.
I took your suggestion regarding the counting variables being 'i' and 'j'.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by