How to calculate on quick way angle (the closest to one of eight offered values) between point and org(0,0) without math functions(arctan)?

Question

How to calculate on quick way angle (the closest to one of eight offered values) between point and org(0,0) without math functions(arctan) ? I have split xy coordinate system in 8 segments ( 0, 45, 90, 135, 180, 225, 270, 315), I need to find angle between point and org( not be accurate just the closest value of eight above) without heavy math functions. I can find like

std::pair<float,float> point;
float angle =arctan(point.second/point.first);
int index =static_cast<int>( (angle+22.5)/45);

and then read from array by index. ( I have added 22.5 degree because to [-22.5, 22.5)=>0, [22.5,67.5)=>45,[67.5,112.5)=>90... ) Is there quicker way, any idea(time of execution is very important)?

Answer 1

This solution is based on maximizing cos(theta-theta_i) where theta_i is 0,45,90,...,315 degrees. I used the angle-difference identity cos(a-b)=cos(a)cos(b)-sin(a)sin(b) and worked with x=r*cos(a) and y=r*sin(a) to simplify the different expressions, so no actual trig functions have to be used:

int getIndex_simple(float x, float y){

  float rCosThetaMinusThetaI[8];
  rCosThetaMinusThetaI[0]=x;
  rCosThetaMinusThetaI[1]=(x+y)*M_SQRT1_2;
  rCosThetaMinusThetaI[2]=y;
  rCosThetaMinusThetaI[3]=(y-x)*M_SQRT1_2;
  rCosThetaMinusThetaI[4]=-x;
  rCosThetaMinusThetaI[5]=(-x-y)*M_SQRT1_2;
  rCosThetaMinusThetaI[6]=-y;
  rCosThetaMinusThetaI[7]=(x-y)*M_SQRT1_2;

  int best_i=0;
  float best_rCosThetaMinusThetaI=rCosThetaMinusThetaI[0];
  for(int i=1; i<8; i++)
    if( rCosThetaMinusThetaI[i]>best_rCosThetaMinusThetaI ){
      best_i=i;
      best_rCosThetaMinusThetaI=rCosThetaMinusThetaI[i];
    }

  return best_i;

}

There are a few simple things you can do to speed this up. For example rCosThetaMinusThetaI[i] == -rCosThetaMinusThetaI[i-4] for i>=4, so you only need four variables. Obviously you don't have to use an array, I just did this to make it easy to read. Also, since rCosThetaMinusThetaI[0]==x and rCosThetaMinusThetaI[2]==y, you only really need two variables. So we can rewrite the above as just a series of eight if statements:

int getIndex_optimized(float x, float y){

  float cosTheta1 = (x+y)*M_SQRT1_2;
  float cosTheta3 = (y-x)*M_SQRT1_2;

  int best_i=0;
  float best_cosTheta=x;
  if( best_cosTheta < cosTheta1 ){ best_i=1; best_cosTheta=cosTheta1; }
  if( best_cosTheta < y         ){ best_i=2; best_cosTheta=y; }
  if( best_cosTheta < cosTheta3 ){ best_i=3; best_cosTheta=cosTheta3; }
  if( best_cosTheta < -x        ){ best_i=4; best_cosTheta=-x; }
  if( best_cosTheta < -cosTheta1){ best_i=5; best_cosTheta=-cosTheta1; }
  if( best_cosTheta < -y        ){ best_i=6; best_cosTheta=-y; }
  if( best_cosTheta < -cosTheta3){ best_i=7; best_cosTheta=-cosTheta3; }

  return best_i;

}

These two functions are C99 code; you can compile them with -std=c99. You need to include math.h to get the constant M_SQRT1_2=0.7071...

Answer 2

Given a point (a,b) with a<b, a>=0 and b>=0, is it closer to 45 degrees or 0 degrees?

Well, tan(theta) = opposite/adjacent, and in the range of 0 degrees to 45 degrees tan is monotonically increasing.

tan(22.5 degrees) =~ 207107/500000. So if a/b > 207107/500000, the closest angle is 45 degrees. If a/b < 207107/500000, the closest angle is 0 degrees. We can even do this without floating point mathematics by saying 500000*a < 207107*b.

For arbitrary points (a,b), we can figure out what quadrant it is in via the signs on a and b. We can rotate (by negation) the problem into the positive-positive quadrant, then invert that rotation on the resulting angle (which is a really simple map).

For arbitrary (a,b) in the positive-positive quadrant, if a>b just reverse a and b, solve as above, and the "closer to 0 degrees" corresponds to "closer to 90 degrees".

Some of the above is overly branchy, but you should be able to turn these branches into integer ops and end with an array access.

Now, note that on some systems, trig function intrinsics can be blazingly fast, much faster than a pile of branchless integer ops and an array lookup. Your first step should be to see if you can replace your arctan with a faster arctan.

bool neg_a = a<0;
bool neg_b = b<0;
a *= (1-2*neg_a);
b *= (1-2*neg_b);
bool near_0 = 500000*a<207107*b; // a/b < 207107/500000
bool near_90 = 207107*a>500000*b; // a/b > 500000/207107
bool near_45 = !near_0 & !near_90;
// 3 CW 2     1
//   -+ | ++
//  2-4 | 0-2 CCW
//4 ----+---- 0
//CCW-- | +- CW
//  4-6 | 6-8
// 5    6     7

// 0 1 or 2
int index = near_45 + 2*near_90;
// negating a or b reverses angle
index *= (1-2*neg_a);
index *= (1-2*neg_b);
// base is 4 if a is negative:
index += 4*(neg_a);
// base is 8 if b is negative, and a is not negative:
index += 8*(neg_b&!neg_a);
index &= 7;

return index;

which is pretty ridiculous, but branch-free. Also not debugged.

Answer 3

Of the 8 rays you are using, one set of four are used when the point is closer (or same distance) to the x-axis than the y-axis. The other four are used otherwise.
Of those 4, the right two rays are used when x is positive (or zero). The other two are used otherwise.
Of those 2, the upper ray is used when y is positive. Otherwise the lower one is used.
So use those simple tests to form an index in a lookup table.

float ClosestAngle(std::pair<float,float> point) {
    float angle[8] = { 247.5 , 112.5f , 292.5 , 67.5f , 202.5 , 157.5 , 337.5 , 22.5f };

    int closerToXAxisThanYAxis = (abs(point.first) <= abs(point.second)) ? 4 : 0;
    int xIsPositive = (point.first >= 0) ? 2 : 0;
    int yIsPositive = (point.second >= 0) ? 1 : 0;

    return angle[closerToXAxisThanYAxis + xIsPositive + yIsPositive];
}

Answer 4

You can calculate which octant the point is in by comparing the gradient of its ray with the gradients of the octant borders, at 22.5, 67.5 degrees, etc. So:

static const float borders[] = { tan(-3 * PI / 8), tan(-PI / 8), tan(PI / 8), tan(3 * PI / 8) };
static const int angles[] = { 270, 315, 0, 45, 90, 135, 180, 225 };
float gradient = y / x;
int i = std::distance(borders, std::lower_bound(std::begin(borders), std::end(borders), gradient)) + (x < 0 ? 4 : 0);
int angle = angles[i];