long auto_arc_tan( long v_dx, long v_dy )
{
// V1b
// (C) Jon P, 1998.
//
// ** INTEGER VERSION **
//
// Returns an angle from 0 to 4095 (which equates to: 0 to nearly 360 degrees)
// Angles are measured clockwise from vertical.
//
// This function will find the angle of orientation of a line specified
// by parameters dx and dy. Imagine this line as the hypotenuse of a
// right angled triangle, dx and dy would be the lengths of the other
// two sides. Although dx and dy aren't strictly just lengths -- they
// need a sign as well.... decide which end of your line is the origin
// and calculate dx and dy as follows:-
// o Measure dx from the origin towards the right.
// (If the line extends left of the origin then dx should be negative.)
// o Measure dy from the origin upwards.
// (If the line extends below the origin then dy should be negative.)
// In other words this function assumes that you're using a co-ordinate
// system where Y values increase in the upward direction and X values
// increase towards the right.
//
// "The 1 o'clock Example".
// For the most straight-forward example imagine a triangle where
// the hypotenuse was the hour hand of a clock pointing to 1 o'clock and
// the second side was a vertical line from the centre of the clock face
// upwards, and the third side was a horizontal line from the tip
// of the hour hand projecting to the left intersecting with the
// vertical side at right angles. The origin is the centre of the clock
// face. This function would give the answer of the equivelent of 30
// degrees (approx. 341 since the results are 0 to 4095 and not in degrees)
//
// This function works by emulating an 'arc tangent' (the opposite
// of a tangent) where a ratio of dx against dy is created and from that
// an angle is returned. This ratio is called the gradient (for which
// the variable name in this code is 'grad').
//
// Arc tangent is accurate and easy to use to calculate angles from 0 to
// 45 degrees (gradients 0.0 to 1.0) but arc tangent becomes increasingly
// problematic to use when the gradient becomes larger...and for 90 degrees
// the ratio is infinity which computers cannot simply handle.
// And to calculate angles approaching 90 degrees dy has to be measured
// with very high accuracy to get good results.
//
// ... so that is one reason why this function only emulates arc tangent
// for the ratios 0.0 to 1.0. The other reason is that this region of the
// arc tangent curve is much easier to emulate. Actually, for extremly
// rough requirements this part of the arc tangent curve could even be
// replaced with a straight line... but this function approximates the
// curve.
//
// "The 2 o'clock Example".
// However this function will calculate any angle from 0.0 to 360.0 degrees
// (0 to 4095 inclusive). So how can it do this when it only emulates
// arc tangents that result in 0 to 45 degrees? This is done by dividing
// up the full 360 arc into 8 regions of 45 degrees each. Up until #5# all
// gradients and angles are for 0 to 45 degrees. To demonstrate this
// imagine the hour hand of the clock was pointing to 2 o'clock (instead
// of 1 o'clock)... this would result in 60 degrees -- an angle greater
// than 45 degrees and hence a gradient greater than 1.0. However if we
// were to do dx/dy (instead of dy/dx) the calculation would result in the
// same gradient & angle as our first example: 30 degrees. The trick is then
// to simply take 30 degrees away from 90 degrees... the result is the
// required answer of 60 degrees (#6# is the line of code that would
// perform this subtraction for this example). All we've done is measured
// the angle from 90 degrees anti-clockwise (instead of the usual way by
// measuring from zero degrees clockwise). When an angle is beyond 90
// degrees the same process applies yet an additional amount 90, 180 or
// 270 degrees has to be added to get the end result (this is done in the
// code after #5#).
//
// The first process in dividing the calculateion into 8 regions is in the
// 'if' statement marked #1#. This makes sure the following code deals
// only with a gradient between 0.0 and 1.0 (which is represented in
// this code as 0 to 4096). It does this by deciding wheather to do
// dy/dx (as would be for the 1 o'clock example) or dx/dy (as it would be
// for the 2 o'clock example). By the time execution has reached #2# or
// #3# 'grad' holds a number between 0 and 4096 (representing a gradient
// between 0.0 and 1.0).
//
// The line just after #2# and #3# sets a flag called 'half'. If half
// is set to 1 this will indicate to the code following #5# that your line
// is in the first half of a quadrant. In the 1 o'clock example
// the #2# branch is executed, the gradient is calculated by dy/dx,
// and 'half' is set to 1.
//
// The 4 lines of code just after the #4# comment is the formula that
// emulates the arc tangent curve between 0.0 and 1.0.
// There is little that can be explained for this since I found the formula
// via trial and error by trying various curves and adjusting their
// parameters until I found one that would approximate the arc tangent
// curve. A look-up table could instead be used here and would be slightly
// quicker yet use more memory. This section of code calculates 'ang' which
// is an angle between 0 and 45 degrees (an actual value of 0 to 512).
//
// The code after #5# reconfigures the 0 to 45 degree angle stored in 'ang'
// to be the full angle 0 to 360 required. as described above.
//
// A similar function to this one is supplied as a standard part of C
// and is called atan2() (or something similar in other computer languages).
// However that function is usually a floating-point function and I needed
// an integer function because the platform I was working on didn't
// have a floating-point processor... so that's why I wrote this function.
long grad, ang, abs_dx, abs_dy, half;
abs_dx = v_dx;
if ( abs_dx < 0 ) abs_dx = 0 - abs_dx;
abs_dy = v_dy;
if ( abs_dy < 0 ) abs_dy = 0 - abs_dy;
// ... make some unsigned versions of dx and dy.
// find the gradient between 0 and one (or 0 and 4096 in this impelmentation)....
if ( abs_dy <= abs_dx ) // #1#
{
if ( abs_dx == 0 )
{ //exact angle can be returned... also this prevents divide by 0 error...
if ( v_dy >= 0 )
return 0;
else
return 2048;
}
grad = abs_dy;
grad <<= 12; // This line multiplies by 4096 making 'grad' a fraction in 4096ths.
grad /= abs_dx; // 'grad' is 0 to 4096 (for gradients 0.0 to 1.0)
// #2#
half = 1;
}
else
{
if ( abs_dy == 0 )
{ //exact angle can be returned... also this prevents divide by 0 error...
if ( v_dx >= 0 )
return 1024;
else
return 3072;
}
grad = abs_dx;
grad <<= 12; // This line multiplies by 4096 making 'grad' a fraction in 4096ths.
grad /= abs_dy; // 'grad' is 0 to 4096 (for gradients 0.0 to 1.0)
// #3#
half = 0;
}
// #4# This is the formula that emulates arc tangent for gradients
// 0.0 to 1.0 ('grad' = 0 to 4096) ...
ang = ( grad * 2409 ) >> 12;
ang *= 4096 - ang;
ang += grad * 3097;
ang >>= 15;
// #5# Find which quadrant the line is in by examining the
// signs of dx and dy. Note that this function assumes Y ordinates
// increase in the upwards direction, and X ordinates increase
// in the right-hand direction.
//
if ( v_dx >= 0 )
{ // handle 0 to 180 degrees...
if ( v_dy >= 0 )
{ // handle 0 to 90 degrees...
if ( half ) ang = 1024 - ang; // #6#
}
else
{ // handle over 90 degrees to 180 degrees...
if ( half )
ang += 1024;
else
ang = 2048 - ang;
}
}
else // handle over 180 degrees to nearly 360 degrees...
{
if ( v_dy >= 0 )
{ // handle 270 degrees to nearly 360 degrees...
if ( half )
ang += 3072;
else
ang = 4096 - ang;
}
else
{ // handle over 180 degrees to nearly 270 degrees...
if ( half )
ang = 3072 - ang;
else
ang += 2048;
}
}
return ang;
}
I developed a new more acurate (but much slower) arc tangent emulation formula based on the curve "n * ( 1.0 - n )".
The code for it is listed below and can be substituted for the section marked #4# above.
ang = grad * grad;
ang = ( grad * -9378 + ang ) / -5282 *
( 4096 - ( ( grad * -921 + ang ) / 3175 ) );
ang /= 25885;
ang += grad;
ang >>= 3;