I have a matrix with size 1000×1000. Assume we're at the position (x,y) and we draw a circle with radius r around us. Afterwards we want to multiply the value of position (x,y) with the values of the position, which lie in the circle and average over all the products. This procedure should be done for every position (x,y) and with multiple values for the radius r. I have tried to write a code (and it works!), but it is really slow. Any assumptions to speed it up? Maybe by avoiding for-loops? I'm pretty new to big-data analysis, so maybe anyone can help? Thank you.
Here is the code: I am only calculating the distance from position (x,y) within a square with sidelength 2*radius. Then the points which are inside the circle are calculated. Also there is the part of matrix u which is in interest for position (x,y) and radius r (variable utmp), so you dont have to multiply a value with the whole matrix u but just with a part. The largest radius will be half of one dimension of x (or y or u). I think the problem are the huge for loops for all the 1000 rows and columns, but I can't find an alternative to them…
INPUT: radius = min_radius:deltar:max_radius; (given parameters) x and y are 1000x1000 matrices containing the x and y positions of a grid u is also a 1000x1000 matrix which contains the values at the positions (x,y) cd = zeros(1, length(radius)); tmp = 0; nancount = 0; dist = cell(size(x, 1), size(x, 2)); in = cell(length(radius),size(x, 1),size(x, 2)); rcount = 1; for dr = radius disp(sprintf('Radius: %d px',dr)); for row = 1:size(x, 1) for col = 1:size(x, 2) %%==== Calculate distances from position (x,y) to any other point ====
if rcount == 1 tic if row-dr < 1 && col-dr >= 1 && col+dr <= size(x,2) dist{row,col} = ((x(row,col) - x(1:row+dr,col-dr:col+dr)).^2 +... (y(row,col) - y(1:row+dr,col-dr:col+dr)).^2).^(1/2); % distance calc.
elseif col-dr < 1 && row-dr >= 1 && row+dr <= size(x,1) dist{row,col} = ((x(row,col) - x(row-dr:row+dr,1:col+dr)).^2 +... (y(row,col) - y(row-dr:row+dr,1:col+dr)).^2).^(1/2); % distance calc. elseif row+dr > size(x,1) && col-dr >= 1 && col+dr <= size(x,2) dist{row,col} = ((x(row,col) - x(row-dr:size(x,1),col-dr:col+dr)).^2 +... (y(row,col) - y(row-dr:size(x,1),col-dr:col+dr)).^2).^(1/2); % distance calc. elseif col+dr > size(x,2) && row-dr >= 1 && row+dr <= size(x,1) dist{row,col} = ((x(row,col) - x(row-dr:row+dr,col-dr:size(x,2))).^2 +... (y(row,col) - y(row-dr:row+dr,col-dr:size(x,2))).^2).^(1/2); % distance calc. elseif row-dr < 1 && col-dr < 1 dist{row,col} = ((x(row,col) - x(1:row+dr,1:col+dr)).^2 +... (y(row,col) - y(1:row+dr,1:col+dr)).^2).^(1/2); % distance calc. elseif row-dr < 1 && col+dr > size(x,2) dist{row,col} = ((x(row,col) - x(1:row+dr,col-dr:size(x,2))).^2 +... (y(row,col) - y(1:row+dr,col-dr:size(x,2))).^2).^(1/2); % distance calc. elseif row+dr > size(x,1) && col-dr < 1 dist{row,col} = ((x(row,col) - x(row-dr:size(x,1),1:col+dr)).^2 +... (y(row,col) - y(row-dr:size(x,1),1:col+dr)).^2).^(1/2); % distance calc. elseif row+dr > size(x,1) && col+dr > size(x,2) dist{row,col} = ((x(row,col) - x(row-dr:size(x,1),col-dr:size(x,2))).^2 +... (y(row,col) - y(row-dr:size(x,1),col-dr:size(x,2))).^2).^(1/2); % distance calc. else dist{row,col} = ((x(row,col) - x(row-dr:row+dr,col-dr:col+dr)).^2 +... (y(row,col) - y(row-dr:row+dr,col-dr:col+dr)).^2).^(1/2); % distance calc.
end eltime(col*row*rcount) = toc; end% %%Calculate inner points
tic% in{rcount,row,col} = dist{row,col}<=dr;
eltime2(col*row*rcount) = toc; end %%Calculate product
tic % Look for the interesting part of Matrix u{dt}
if row-dr < 1 && col-dr >= 1 && col+dr <= size(x,2) utmp = u{dt}(1:row+dr,col-dr:col+dr); elseif col-dr < 1 && row-dr >= 1 && row+dr <= size(x,1) utmp = u{dt}(row-dr:row+dr,1:col+dr); elseif row+dr > size(x,1) && col-dr >= 1 && col+dr <= size(x,2) utmp = u{dt}(row-dr:size(x,1),col-dr:col+dr); elseif col+dr > size(x,2) && row-dr >= 1 && row+dr <= size(x,1) utmp = u{dt}(row-dr:row+dr,col-dr:size(x,2)); elseif row-dr < 1 && col-dr < 1 utmp = u{dt}(1:row+dr,1:col+dr); elseif row-dr < 1 && col+dr > size(x,2) utmp = u{dt}(1:row+dr,col-dr:size(x,2)); elseif row+dr > size(x,1) && col-dr < 1 utmp = u{dt}(row-dr:size(x,1),1:col+dr); elseif row+dr > size(x,1) && col+dr > size(x,2) utmp = u{dt}(row-dr:size(x,1),col-dr:size(x,2)); else utmp = u{dt}(row-dr:row+dr,col-dr:col+dr); end eltime3(col*row*rcount) = toc; tic m = nanmean(u{dt}(row,col).*utmp(in{rcount,row,col})); if isnan(m), m = 0; nancount = nancount + 1; end tmp = tmp + m; eltime4(col*row*rcount) = toc; end end cd(dt, rcount) = tmp/(col*row-nancount); %%OUTPUT
tmp = 0; nancount = 0; rcount = rcount + 1; end
Best Answer