Matlab : Vectorize technics in 3 dimensions matrix

Question

I actually vectorizing one of my code and I have some issues.

This is my initial code:

CoordVorBd = random(N+1,3)
CoordCP = random(N,3)
v = random(1,3)
for i = 1 : N
    for j = 1 : N
            ri1j = (-CoordVorBd (i,:) + CoordCP(j,:));

            vij(i,j,:) = cross(v,ri1j))/(norm(ri1j)
     end     

end

I have start to vectorize that creating some matrix that contains 3*1 Vectors. My size of matrix is N*N*3.

CoordVorBd1(1:N,:) = CoordVorBd(2:N+1,:);
CoordCP_x= CoordCP(:,1);
CoordCP_y= CoordCP(:,2);
CoordCP_z= CoordCP(:,3);
CoordVorBd_x = CoordVorBd([1:N],1);
CoordVorBd_y = CoordVorBd([1:N],2);
CoordVorBd_z = CoordVorBd([1:N],3);
CoordVorBd1_x = CoordVorBd1(:,1);
CoordVorBd1_y = CoordVorBd1(:,2);
CoordVorBd1_z = CoordVorBd1(:,3);
[X,Y] = meshgrid (1:N);

ri1j_x = (-CoordVorBd_x(X) + CoordCP_x(Y));
ri1j_y = (-CoordVorBd_y(X) + CoordCP_y(Y));
ri1j_z = (-CoordVorBd_z(X) + CoordCP_z(Y));

ri1jmat(:,:,1) = ri1j_x(:,:);
ri1jmat(:,:,2) = ri1j_y(:,:);
ri1jmat(:,:,3) = ri1j_z(:,:);
vmat(:,:,1) = ones(N)*v(1);
vmat(:,:,2) = ones(N)*v(2);
vmat(:,:,3) = ones(N)*v(3);

This code works but is heavy in terms of variable creation. I did'nt achieve to apply the vectorization to all the matrix in one time.

The formule like

ri1jmat(X,Y,1:3) = (-CoordVorBd (X,:) + CoordCP(Y,:));

doesn't work... If someone have some ideas to have something cleaner.

At this point I have a N*N*3 matrix ri1jmat with all my vectors.

I want to compute the N*N rij1norm matrix that is the norm of the vectors

rij1norm(i,j) = norm(ri1jmat(i,j,1:3))

to be able to vectorize the vij matrix.

vij(:,:,1:3) = (cross(vmat(:,:,1:3),ri1jmat(:,:,1:3))/(ri1jmatnorm(:,:));

The cross product works.

I tried numbers of method without achieve to have this rij1norm matrix without doing a double loop.

If someone have some tricks, thanks by advance.

score 0 · Answer 1 · answered Mar 06 '17 at 19:53

Here's a vectorized version. Note your original loop didn't include the last column of CoordVorBd, so if that was intentional you need to remove it from the below code as well. I assumed it was a mistake.

CoordVorBd = rand(N+1,3);
CoordCP = rand(N,3);
v = rand(1,3);

repCoordVor=kron(CoordVorBd', ones(1,size(CoordCP,1)))'; %based on http://stackoverflow.com/questions/16266804/matlab-repeat-every-column-sequentially-n-times
repCoordCP=repmat(CoordCP, size(CoordVorBd,1),1); %repeat matrix
V2=-repCoordVor + repCoordCP; %your ri1j 
nrm123=sqrt(sum(V2.^2,2)); %vectorized norm for each row
vij_unformatted=cat(3,(v(:,2).*V2(:,3) - V2(:,2).*v(:,3))./nrm123,(v(:,3).*V2(:,1) - V2(:,3).*v(:,1))./nrm123,(v(:,1).*V2(:,2) - V2(:,1).*v(:,2))./nrm123); % cross product, expanded, and each term divided by norm, could use bsxfun(@rdivide,cr123,nrm123) instead, if cr123 is same without divisions
vij=permute(reshape( vij_unformatted,N,N+1,3),[2,1,3]); %reformat to match your vij

score 0 · Answer 2 · answered May 11 '17 at 06:31

Here is another way to do it using arrayfun

% Define a meshgrid of indices to run over
[I, J] = meshgrid(1:N, 1:(N+1));
% Calculate ril for each index
rilj = arrayfun(@(x, y) -CoordVorBd (y,:) + CoordCP(x,:), I, J, 'UniformOutput', false);
%Calculate vij for each point
temp_vij1 = arrayfun(@(x, y) cross(v, rilj{x, y}) / norm(rilj{x, y}), J, I, 'UniformOutput', false);
%Reshape the matrix into desired format
temp_vij2 = cell2mat(temp_vij1);
vij = cat(3, temp_vij2(:, 1:3:end), temp_vij2(:, 2:3:end), temp_vij2(:, 3:3:end));

Matlab : Vectorize technics in 3 dimensions matrix

2 Answers2