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;