2D matrix indexing for element wise operation

Hello all, I am using following code to estimate pairwise angle between for a vector. Where the angle should be between 0 to 2*pi. I am able to predict the angle. But I need your advise on implementing this two for loops. This is slowing the code. Looking forward for your attention. In the code, XYdata is a n x 2 matrix representing the Cartesian coordinates of individual point [xi, yi]. I presume this slow is because of two FOR loops.

        for i2 =1:length(XYdata_new)
            for j2 = 1:length(XYdata_new)
                angle_ij(i2,j2)= Angle360(XYdata_new(i2,1), XYdata_new(i2,2), XYdata_new(j2,1),XYdata_new(j2,2));
            end
        end

Thanks in advance.

[Stephen23]

Please show us the code for Angle360.

[Raj Raj]

function  angleout = Angle360(x1,y1,x2,y2)
if x1 == x2 && y1 ==y2
    angleout = Inf;
end
angle = atan((y2-y1)/(x2-x1));
if angle > 0 
  if (y2 < y1) % y2 > y1, x2 > x1  % dy > 0, dx > 0 %first quadrant
    angleout=pi+ angle;
  else         % 
    angleout=  angle; % 3rd quadrant
  end
elseif angle <  0
  if x2 < x1 && y2 > y1
    angleout= pi + angle; % 2nd quadrant
  elseif y2 < y1 && x2 > x1 
    angleout= 3*pi/2 + (pi/2+angle); % 4th quadrant
  elseif y2 < y1 && x2-x1 == 0
    angleout= pi - angle; % 4th quadrant    
  end
elseif angle == 0
    if x2 < x1 
        angleout = pi +angle;
    else
        angleout = angle;
    end
end

[Raj Raj]

here it is

Answer 1

[Star Strider]

I would use the atan2 or atan2d functions. They are vectorised, so you can likely give your matrix columns to them directly and then do the angle difference calculations.

[Raj Raj]

Hi Star, Can I seek an extension to your reply. I wanted pairwise angle. So, I may need a nXn matrix out of nX2 coordinates.

[Star Strider]

Having some of your data would help. One option to generate your matrix would be to use the bsxfun function.

Example:

M = randi([-5 5], 10, 2);                               % Create Data
Angles = atan2d(M(:,2), M(:,1));                        % Generate Angles (°)
AngleDifferences = bsxfun(@minus, Angles, Angles');     % Pairwise Angle Differences (°)

If you want to use radian angles, use the atan2 function instead.

[Raj Raj]

Hi Star,

Thanks for your example. Let me re-phrase my problem. I have a nX2 data points. I will select a point, lets say "i" and I want to calculate the angle between "i" and all the remaining points. A vector output is also serves the purpose. All the angles supposed to be in between 0 to 2*pi.

I had used atan2d or atan2 commands. However they don't give angles between 0-2*pi. So I can use my subroutine Angle360.

Using two for loops seems to be a costly choice.

If I consider M=[0,0; 1,1; 0,1; -1,0; 0,-1]; and if I consider M(1,:) with respect to M, I should get output as [0, 0, 90, 180, 270, 360 or 0].

Having that, I tried to rephrase my code to calculate only a vector as follows:

          for i2 =1:length(XYdata_new)
              angle_ij(i2,1)= Angle360(XYdata_new(i,1), XYdata_new(i,2), XYdata_new(i2,1),XYdata_new(i2,2));
          end

Let me know, if I can enhance any way by avoiding the for loop and using self indexing. It seems using indexing, my Angle360 is not gonna work, as Angle360 is written only for individual values rather than vectors.

[Star Strider]

Two loops definitely is a costly choice.

This updated version wraps the angles so they all go from (0,360) or (0,2*pi), depending on what you want to do. If you want degrees, use atan2d, if radians use atan2, then use the appropriate conversion line:

M = randi([-5 5], 10, 2);                                       % Create Data
Angles = atan2d(M(:,2), M(:,1));                                % Generate Angles (°)
Angles360 = Angles.*(Angles>=0) + (360+Angles).*(Angles < 0);   % Convert Degree Angles To (0,360)
% Angles2pi = Angles.*(Angles>=0) + (2*pi+Angles).*(Angles < 0);  % Convert Radian Angles To (0,2*pi)
AngleDifferences = bsxfun(@minus, Angles360, Angles360');       % Pairwise Angle Differences

[Raj Raj]

Thank you Star for your patience in extending your response.

In your code, I don't see a i^th variable (with respect to which I need to find the angles). In my recent "one FOR loop version", you can see that I am specifying XYdata_new(i,:) is the point about which I am finding the angles to the all data points. This implementation is similar to the pdist2.

[Star Strider]

In the matrix bsxfun calculates, the first row are the angle differences between the first angle an all the others (so ‘AngleDifferences(1,2)’ is ‘Angles360(1)-Angles360(2)’ and so for the rest. This is the angular equivalent of the pdist function.

If you have two different sets of angles, run ‘Angles’ and ‘Angles360’ on each set separately, producing for example ‘Angles360A’ and ‘Angles360B’. The bsxfun call then becomes:

AngleDifferences = bsxfun(@minus, Angles360A(:), Angles360B(:)');               % Pairwise Angle Differences

The ‘AngleDifferences(1,2)’ matrix element would then become for example ‘Angles360A(1)-Angles360B(2)’, and so for the rest. (Note that in this line I forced them to both be column vectors.) (Make the appropriate changes to my code if you are using radian measure.)

This the angle equivalent of pdist2, so it should do what you want.

[Raj Raj]

Thanks Star, I have marked it as answered.

[Star Strider]

My pleasure!

Thank you.

EDIT — More efficient versions of the angle conversion lines are:

Angles360 = rem(360+a, 360);
Angles2pi = rem(2*pi+a, 2*pi);