How can I vectorise the next for - loop?

Hello, I have the following for loop and I would like to know if there is a way to vectorise it:

JM = 91;
stepVec = 10;
JN = 100; 

x = zeros(JM,JM);
y = zeros(JM,JM);
z = zeros(JM,JM);
v = zeros(1,3,JN);
kk=1;
for i=1:stepVec:JM
 for ii=1:stepVec:JM
    
     theta=ii*pi/180;
     phi=i*pi/180;
     x(i,ii) = r*sin(theta)*sin(phi);
     y(i,ii) = r*sin(theta)*cos(phi);
     z(i,ii) = r*cos(theta);
    
    v(:,:,kk) = [x(i,ii) y(i,ii) z(i,ii)];
    kk=kk+1;
 end
end

Thank you for your answers, I would like to know if I can vectorise the next two too, although I don't think it's possible.

Rotat = zeros(4,4,length(t));
RotatVector = cell(1, JN);
F = struct('cdata',cell(1,length(t)),'colormap',cell(1,length(t)));
Faxis = cell(1,JN);
for j=1:JN
 for i=1:length(t)
     RR = makehgtform('axisrotate',[v(1,1,j) v(1,2,j) v(1,3,j)],i*beta); 
     set(hgtc,'Matrix',RR);
     Rotat(:,:,i) = RR;
     RotatVector{1,j} = Rotat;
          
    
     width = 300;
     height = 225;
     F(i) = getframe(gcf,[140 109.5 width height]); 

     Faxis{j} = F;
     drawnow;
 end
end

% --------------------------------------------------------------------------

grayFrame = zeros(282,375,length(F));
grayFrameNEW = cell(1,JN);
meanGrayLevels = cell (1,length(F));
meanGrayLevelsCELL = cell(1,length(F),JN);
for j=1:JN
 for i=1:length(F)    
     % convert the image to a frame
     grayFrame(:,:,i) = rgb2gray(Faxis{j}(i).cdata);  
     grayFrameNEW{j}=grayFrame;
     % Calculate the mean gray level
     meanGrayLevels{i} = mean2(grayFrameNEW{1,j}(:,:,i));
     meanGrayLevelsCELL(:,:,j) = meanGrayLevels;
 end
end

[Walter Roberson]

You cannot vectorize getframe(), so you will not be able to avoid looping for that.

However, your line

Faxis{j} = F;

should be after the "for i" loop.

[Raúl Alonso Merino]

Thank you, I have corrected that, about the third loop is there something possible?

Answer 1

[Walter Roberson]

for i=1:stepVec:JM
 for ii=1:stepVec:JM
    
     theta=ii*pi/180;
     phi=i*pi/180;
     x(i,ii) = r*sin(theta)*sin(phi);
     y(i,ii) = r*sin(theta)*cos(phi);
     z(i,ii) = r*cos(theta);
    
    v(:,:,kk) = [x(i,ii) y(i,ii) z(i,ii)];
    kk=kk+1;
 end
end

can be replaced with

angles = (1:stepVec:JM) * pi/180;
sinang = sin(angles);
cosang = cos(angles);
nang = length(angles);
rsin = r * sinang;
rcos = r * cosang;
k = 1;
for i = 1 : nang
    x(i, :) = rsin(i) * sinang;
    y(i, :) = rsin(i) * cosang;
    z(i, :) = rcos;
    v(1, :, k:k+nang-1) = [x y z];
    k = k + nang;
end

and this can be further vectorized:

angles = (1:stepVec:JM) * pi/180;
nang = length(angles);
[sTheta, sPhi] = ndgrid(sin(angles));
[cTheta, cPhi] = ndgrid(cos(angles));
v = r .* [sTheta(:) .* sPhi(:), sTheta(:) .* cPhi(:), cTheta(:)];  %(nang^2) x 3 x 1
v = permute(v, [3 1 2]);     %1 x 3 x (nang^2)

[Raúl Alonso Merino]

Thank you for your answer, my problem using the last vectorization now is that I get this error:

Index in position 3 exceeds array bounds (must not exceed 3).

Error in PruebaDosAngulos (line 89)
     RR = makehgtform('axisrotate',[v(1,1,j) v(1,2,j) v(1,3,j)],i*beta);  %Rotation matrix

And it is because the second loop I posted, the one with the getframe, uses the variable v. How can I solve this error so the next for loop uses the factorization you propose?

Thank you.

[Walter Roberson]

Loop version:

r = 20;
JM = 91;
stepVec = 10;
JN = 100; 

angles = (1:stepVec:JM) * pi/180;
sinang = sin(angles);
cosang = cos(angles);
nang = length(angles);

x = zeros(nang, nang);
y = zeros(nang, nang);
z = zeros(nang, nang);
v = zeros(1, 3, nang*nang);

rsin = r * sinang;
rcos = r * cosang;
k = 1;
for i = 1 : nang
    x(i, :) = rsin(i) * sinang;
    y(i, :) = rsin(i) * cosang;
    z(i, :) = rcos;
    v(1, :, k:k+nang-1) = [x(i,:); y(i,:); z(i,:)];
    k = k + nang;
end

vectorized version:

r = 20;
JM = 91;
stepVec = 10;
JN = 100;

angles = (1:stepVec:JM) * pi/180;
nang = length(angles);
[sTheta, sPhi] = ndgrid(sin(angles));
[cTheta, cPhi] = ndgrid(cos(angles));
v = r .* [sTheta(:) .* sPhi(:), sTheta(:) .* cPhi(:), cTheta(:)];  %(nang^2) x 3 x 1
v = permute(v, [3 2 1]);     %1 x 3 x (nang^2)

Note: I do not understand why you want a scalar first dimension for v. It makes things awkward.

[Raúl Alonso Merino]

Thank you for your answer for the first loop, I have the v with those dimensions to have everything in the code more organised although I know I could propbably make it quite better but I'm not good at programming and I do not know how to improve my code because right now is Running out of memory in my computer.

Note: I leave the question marked without an answer for now to see if someone knows how to improve the second and third loop. Thank you.

[Walter Roberson]

You are running out of memory because you are saving a copy of every frame in memory.

If your purpose is to create a series of movies, then after the for i loop, instead of doing

Faxis{j} = F;

use VideoWriter or equivalent to save F to file instead of saving it in Faxis .

If you are creating one long movie out of all of these, then you would open the movie before starting for j and after the for i you would append what you just captured.

Or instead of saving everything, just calculate the gray levels as you go along and do not save anything you do not need to save.

I cannot see any good reason to organize v as 3D. You are not using the first dimension for any purpose that I can see.

[Raúl Alonso Merino]

I was saving some variables to see them later but you are right I do not need to save every frame, I thought I needed to save them so I can compute later the gray levels. How can I go along as you say without saving for example the frames?

Thank you again so much, you are helping me a lot because I was really stuck.

And also how to use later the v vector saving it as 3D? I have changed so many times my code I have forgotten how I had some parts.

[Walter Roberson]

You can initialize

meanGrayLevels = zeros(length(t), JN);

Then as you capture a frame into F (probably without index),

meanGrayLevels(i, j) = mean2( rgb2gray(F.cdata) );

[Raúl Alonso Merino]

Thanks that worked perfectly, it took it 54 hours to compute but it was because I choose a JN = 129600 as I have two angles of 360 degrees each and I wanted to rotate for all the degrees. I will have to think now how to reduce time on it but thanks.

I'll give your answer as solution as you were really helpful. Thank you a lot.