54 views (last 30 days)

I start with a 3 x 46 array where each 3 x1 column in the array is viewed as a position vector of a point, let's say the first vector in the array is p_1, the next one is p_2.

I also have a 3 x 100 array, which I am viewing as the position vectors of 100 points (where these points are always of a lower 'radius' than the points in the 3 x 46 array). I basically need to create a matrix which looks as follows, where each distance is the row vector obtained by vector subtraction of the two vectors:

[Distance from p_1 to first vector in the 3 x 100 array Distance from p_1 to second vector etc ]

[Distance from p_2 to first vector in the 3 x 100 array Distance from p_2 to second vector etc ]

[ Distance from p_3 ..... ]

and so on, until all the positions in the 3 x 46 array are covered, let me know if my intention is not clear.

Guillaume
on 8 Aug 2019

Edited: Guillaume
on 8 Aug 2019

This can be achieved easily without a loop:

%demo data

m = 46

n = 100

A = rand(3,m)

B = rand(3,n)

distance = sqrt(sum((permute(A, [2, 3, 1]) - permute(B, [3, 2, 1])).^2, 3));

The result is a mxn matrix. It basically moves the (X, Y, Z) of each point into the 3rd dimension for both matrix, and move the 2nd dimension of the 1st matrix into the 1st. So points in A are along the rows, and points in B are along the column (and coordinates across the 3rd for both). We now have compatible array sizes, so can create a 3D matrix of (Xa-Xb, Ya-Yb, Za-Zb) for all the points. Square that, sum across (X, Y, Z) so it reduces to a 2D matrix, and take the square root for the distance.

edit: Note that if you want to store the vectors between the points as suggested by Jon instead of the distance (norm of the vectors), then it's simply:

vectors = permute(A, [2, 3, 1]) - permute(B, [3, 2, 1])

vectors(p1, p2, :) is the vector between point A(:, p1) and point B(:, p2).

Guillaume
on 12 Aug 2019

arrayfun is just a for loop in disguise. It just saves you on having to compute the bounds yourself (more on that later). With either method, if speed is not an issue, I'd create a cell array rather than index into the destination. The code is a lot easier to read this way.

So, I'd do:

matrixbuilder = @(rij, normal) [gradletSlip(rij, normal), thermletSlip(rij, normal);

gradletJump(rij, normal), thermletJump(rij, normal)];

[colp, colsing] = ndgrid(1:size(colp, 2), 1:size(colsing, 2))

distance = arrayfun(@(colp, colsing) matrixbuilder((p(:, colp) - sing(:, colsing))', p(:, colp)' / norm(p(:, colp))), colp, colsing, 'UniformOutput', false)

distmat = cell2mat(distance);

Sign in to comment.

Jon
on 8 Aug 2019

Edited: Jon
on 8 Aug 2019

So from your description you actually want to find and store the vectors that go from each point in the first matix to each point in the second matrix. Here is a way you could do that

m = 46

n = 100

A = rand(3,m)

B = rand(3,n)

% preallocate 3d matrix to hold results

C = zeros(m,3,n);

% loop through columns of A calculating "distance" vectors and storing

% them

for k = 1:m

C(k,:,:) = B - A(:,k);

end

% to look, for example, at "distance vectors" from the 3rd column of A to all of the

% columns of B you can use

disp(squeeze(C(3,:,:)))

Note, the squeeze function gets rid of extraneous single dimension in multidimensional arrays, which would otherwise make it difficult to see the result

If instead you just want the scalar Euclidean distance from each point in A to each point in B, you could do the following

D = zeros(m,n);

for k = 1:m

D(k,:) = vecnorm(B - A(:,k));

end

If you wanted both you could put both calculations in the same loop.

It may be possible to vectorize this further and eliminate the loop, but it is not immediately obvious to me how to do this. Maybe someone else has a suggestion.

Jon
on 8 Aug 2019

I don't know what your application is so I can't give you any advice on whether you want to define the "distance" vectors as P-S or S-P. I don't think it changes the main logic that I outlined. If I'm understanding correctly it just changes some algebraic signs.

You could probably work out how to permute the 4d array that I have obtained above to get the 3n x 3m array that you want. Otherwise, here is how you can store it this way in the first place. You should check the details, to make sure it is right, but this should give you enough of an idea to go on.

m = 4

n = 5

A = rand(3,m)

B = rand(3,n)

% with 3d result for each pairing

% preallocate 2d array to hold results

J = zeros(n*3,m*3);

for i = 1:m

r = B - A(:,i);

rnorm = vecnorm(r);

for j = 1:n

startRow = (j-1)*3 + 1; % starting row where this 3x3 result is stored

startCol = (i-1)*3 + 1; % starting column where this 3x3 result is stored

J(startRow:startRow+2,startCol:startCol+2) = ...

1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/rnorm(j)+(r(:,j)'*r(:,j))/rnorm(j)^3);

end

end

Guillaume
on 9 Aug 2019

Working in N-D is no more harder than working in 2-D. And in your case it would make it much easier to find the matrix related to a pair of point. You can either have the indices of the points as the first 2 indices, so:

distmat = squeeze(J(p1, p2, :, :)); %3x3 distance matrix between A(:, p1) and B(:, p2)

or as the last 2 indices:

distmat = J(:, :, p1, p2); %3x3 distance matrix between A(:, p1) and B(:, p2)

Storing the same in a 2D array, it's a lot less clear:

distmat = J((p1-1)*3 + (1:3), (p2-1)*3 + (1:3));

Another option, as I suggested is to store the matrices in a 2D cell array, in which case:

distmat = J{p1, p2};

The disadvantage of cell arrays is the memory and speed overhead. For an easy way to do this, see comment to my answer.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.