Vectorization over different array sizes

Question

Cyrill Slezak 2022-3-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1668039-vectorization-over-different-array-sizes

回答： Suman Sahu 2023-2-13

I'm trying to linearize (i.e. speed up) the following looped version of my code

for x=-range:1:range
   for y=-range:1:range
      for z=-range:1:range
        canvas(index(1)+x,index(2)+y,index(3)+z)+=const*exp(-sum(((delta *[x,y,z]*RotationMatix)./[a,a,b]).^2));
      end  for
   end for
end for

The element vise evaluation within the sum works well but I'm trying to also vectorize the remaining loops. I would appreciate any help and if you could point out why my aproach (below) fails I'd appreciate it!

[x y z] = ndgrid(-range:1:range,-range:1:range, -range:1:range);
canvas(index(1)-range:index(1)+range,index(2)-range:index(2)+range,index(3)-range:index(3)+range)+= ...
        const*exp(-sum(((delta*[x(:) y(:) z(:)]*RotationMatix)./[a,a,b]).^2));

Thank you!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Suman Sahu 2023-2-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1668039-vectorization-over-different-array-sizes#answer_1170240

在 MATLAB Online 中打开

The issue with your approach is that you are trying to use matrix addition between two arrays with different shapes. The first loop results in a 3D matrix of size (2 * range + 1)^3, whereas the second approach results in a 1D vector of size (2 * range + 1). You cannot add these two arrays together because they have different shapes.

To solve this problem, you can try to reshape the 1D vector into a 3D matrix of the same size as the output of the first loop before adding it to the canvas matrix. The following code should work:

[x, y, z] = ndgrid(-range:1:range, -range:1:range, -range:1:range);
result = const * exp(-sum(((delta * [x(:), y(:), z(:)] * RotationMatix) ./ [a, a, b]).^2)); 
%reshape so the result can be added to the canvas
result = reshape(result, size(x));
canvas(index(1)-range:index(1)+range, index(2)-range:index(2)+range, index(3)-range:index(3)+range) += result;