#include "inside.h"


#define	M_PI		3.14159265358979323846	/* pi */
#define	M_PI_2		1.57079632679489661923	/* pi/2 */




/***************************************************************************\
*                                                                           *
*   Module Name:        angle                                               *
*                                                                           *
*   Programmer(s):      Unknown                                             *
*                                                                           *
*   Date Written:       Unknown                                             *
*                                                                           *
*   Modified:           Jan C. Depner, ported to C.                         *
*                                                                           *
*   Date Modified:      August 1992                                         *
*                                                                           *
*   Purpose:            Computes the angle between the three geographic     *
*                       points specified.                                   *
*                                                                           *
*   Arguments:          x0                  -   x position of point 0       *
*                       y0                  -   y position of point 0       *
*                       x1                  -   x position of point 1       *
*                       y1                  -   y position of point 1       *
*                       x2                  -   x position of point 2       *
*                       y2                  -   y position of point 2       *
*                                                                           *
*   Return Value:       NV_FLOAT64          -   angle                       *
*                                                                           *
*   Method:             Computes the cross and dot products of the two      *
*                       vectors.  The cross product divided by the dot      *
*                       product gives the tangent of the angle.  Then the   *
*                       arc tangent is taken to get the angle.              *
*                                                                           *
\***************************************************************************/

static NV_FLOAT64 angle (NV_FLOAT64 x0, NV_FLOAT64 y0, NV_FLOAT64 x1, NV_FLOAT64 y1, 
                         NV_FLOAT64 x2, NV_FLOAT64 y2)
{
  NV_FLOAT64          xx1, yy1, xx2, yy2, cross, dot, local_angle;

    xx1 = x1 - x0;
    yy1 = y1 - y0;
    xx2 = x2 - x0;
    yy2 = y2 - y0;
    cross = xx1 * yy2 - xx2 * yy1;
    dot = xx1 * xx2 + yy1 * yy2;

    if (fabs (dot) <= 1.0e-30)
    {
        local_angle = M_PI_2;
    }
    else
    {
        local_angle = atan (fabs (cross / dot));
        if (dot < 0.0) local_angle = M_PI - local_angle;
    }

    if (cross < 0.0) local_angle = -local_angle;

    return (local_angle);
}



/***************************************************************************\
*                                                                           *
*   Module Name:        inside                                              *
*                                                                           *
*   Programmer(s):      Unknown                                             *
*                                                                           *
*   Date Written:       Unknown                                             *
*                                                                           *
*   Modified:           Jan C. Depner, ported to C.                         *
*                                                                           *
*   Date Modified:      August 1992                                         *
*                                                                           *
*   Purpose:            Checks a point to see if it falls within the        *
*                       specified polygon.                                  *
*                                                                           *
*   Arguments:          ax                  -   x points                    *
*                       ay                  -   y points                    *
*                       count               -   number of points            *
*                       x                   -   x position of point         *
*                       y                   -   y position of point         *
*                                                                           *
*   Return Value:       NV_INT32            -   1 if inside, 0 if not       *
*                                                                           *
*   Method:             Sums all of the angles between the point and each   *
*                       pair of polygon vertices.  If the sum is greater    *
*                       than PI, the point is inside the polygon.           *
*                                                                           *
*   Restrictions:       The vertices of the polygon must be in order around *
*                       the polygon.                                        *
*                                                                           *
\***************************************************************************/

NV_INT32 inside (NV_FLOAT64 *ax, NV_FLOAT64 *ay, NV_INT32 count, NV_FLOAT64 x, NV_FLOAT64 y)
{
    NV_FLOAT64          angle_sum;
    NV_INT32            i;


    /*  There have to be at least three points in the polygon.          */

    if (count > 2)
    {
        angle_sum = 0.0;

        for (i = 1 ; i < count ; i++)
        {
            angle_sum += angle (x, y, ax[i], ay[i], ax[i - 1], ay[i - 1]);
        }
        
        angle_sum += angle (x, y, ax[0], ay[0], ax[count - 1], ay[count - 1]);

        if (fabs (angle_sum) > M_PI) return (1);
    }
    
    return (0);
}