MATLAB - get angle between arrays of points

Question

For the sake of illumination analysis, based on this document, I am trying to determine three things for an array of lights and a series of points on a solid surface:

(Image key: big blue points are lights with illumination direction shown, small points are the points on my surface)

• 1) The distances between each of the lights and each of the points,
• 2) the angles between the direction each light is facing and the normal vectors of all of the points:

Note in this image I have replicated the normal vector and moved it to more clearly show the angle.

• 3) the angles between the direction each light is facing, and the vector from that light to all of the points on the solid:

Originally I had nested for loops iterating through all of the lights and points on the solid, but am now doing my best to do it in true MATLAB style with matrices:

I have found the distances between all the points with the pdist2 function, but have not managed to find a similar method to find the angles between the lights and all the points, nor the lights and the normal vectors of the points. I would prefer to do this with matrix methods rather than with iteration as I have been using.

Considering I have data set out, where each column of Lmat has my x,y,z position vectors of my lights; Dmat gives x,y,z directions of each light, thus the combination of each row from both of these matrices fully define the light and the direction it is facing. Similarly, Omega and nmat do the same for the points on the surface.

I am fairly sure that to get angles I want to do something along the lines of:

distMatrix = pdist2(Omega, Lmat);

LmatNew = zeros(numPoints, numLights, 3);
DmatNew = zeros(numPoints, numLights, 3);
OmegaNew = zeros(numPoints, numLights, 3);
nmatNew = zeros(numPoints, numLights, 3);
for i = 1:numLights
    LmatNew(:,i,1) = Lmat(i,1);
    LmatNew(:,i,2) = Lmat(i,2);
    LmatNew(:,i,3) = Lmat(i,3);

    DmatNew(:,i,1) = Dmat(i,1);
    DmatNew(:,i,2) = Dmat(i,2);
    DmatNew(:,i,3) = Dmat(i,3);
end

for j = 1:numPoints
    OmegaNew(j,:,1) = Omega(j,1);
    OmegaNew(j,:,2) = Omega(j,2);
    OmegaNew(j,:,3) = Omega(j,3);

    DmatNew(:,i,1) = Dmat(i,1);
    DmatNew(:,i,2) = Dmat(i,2);
    DmatNew(:,i,3) = Dmat(i,3);
end


angleMatrix = -dot(LmatNew-OmegaNew, DmatNew, 3);
angleMatrix = atand(angleMatrix);

angleMatrix = angleMatrix.*(angleMatrix > 0);

But I am getting conceptually stuck trying to get my head around what to do after my dot product.

Am I on the right track? Is there an inbuilt angle equivalent of pdist2 that I am overlooking?

Thanks all for your help, and sorry for the paint images!

Context: This image shows my lights (big blue points), the directions the lights are facing (little black traces), and my model.

Answer 1

According to MathWorks, there is no built-in function to calculate the angle between vectors. However, you can use trigonometry to calculate the angles.

Inputs

Since you unfortunately didn't explain your input data in great detail, I'm going to assume that you have a matrix Lmat containing a location vector of a light source in each row and a matrix Dmat containing the directional vectors for the light sources, both of size n×3, where n is the number of light sources in your scene.

The matrices Omega and Nmat supposedly are of size m×3 and contain the location vectors and normal vectors of all m surface points. The desired result are the angles between all light direction vectors and surface normal vectors, of which there are n⋅m, and the angles between the light direction vectors and the vectors connecting the light to each point on the surface, of which there are n⋅m as well.

To get results for all combinations of light sources and surface points, the input matrices have to be repeated vertically:

Lmat  = repmat(Lmat,  size(Omega,1), 1);
Dmat  = repmat(Dmat,  size(Omega,1), 1);
Omega = repmat(Omega, size(Lmat,1),  1);
Nmat  = repmat(Nmat,  size(Lmat,1),  1);

Using the inner product / dot product

The definition of the inner product of two vectors is

where θ is the angle between the two vectors. Reordering the equation yields

You can therefore calculate the angles between your directional vectors Dmat and your normal vectors Nmat like this:

normProd = sqrt(sum(Dmat.^2,2)).*sqrt(sum(Nmat.^2,2));
anglesInDegrees = acos(dot(Dmat.',Nmat.')' ./ normProd) * 180 / pi;

To calculate the angles between the light-to-point vectors and the directional vectors, just replace Nmat with Omega - Lmat.

Using the vector product / cross product

It has been mentioned that the above method will have problems with accuracy for very small (θ ≈ 0°) or very large (θ ≈ 180°) angles. The suggested solution is calculating the angles using the cross product and the inner product.

The norm of the vector product of two vectors is

You can combine this with the above definition of the inner product to get

which can obviously be reordered to this:

The corresponding MATLAB code looks like this:

normCross = sqrt(sum(cross(Dmat,Nmat,2).^2,2));
anglesInDegrees = atan2(normCross,dot(Dmat.',Nmat.')') * 180/pi;