libDXF snippets

snippets.c

/*!
 * \file snippets.c
 */

#include <stdio.h>
#include <time.h>
#include <math.h>


typedef struct
{
        double x;
        double y;
        double z;
} XYZ;

typedef struct
{
        int x;
        int y;
} Point;


/*!
 * \brief Calculate the determinant of a 2x2 (2D) matrix.
 *
\f[
\left |
\matrix{
a b \cr
c d
}
\right | = a d - b c
\f]
 *
 * \return Determinant.
 */
static double
determinant_2D
(
        double a,
        double b,
        double c,
        double d
)
{
  return (a * d - b * c);
}


/*!
 * \brief Calculate the determinant of a 3x3 (3D) matrix.
 *
\f[
\left |
\matrix{
a b c \cr
d e f \cr
g h i
}
\right | = a e i + b f g + c d h - c e g - b d i - a f h
\f]
 *
 * \return Determinant.
 */
static double
determinant_3D
(
        double a,
        double b,
        double c,
        double d,
        double e,
        double f,
        double g,
        double h,
        double i
)
{
  return (a * e * i + b * f * g + c * d * h - c * e * g - b * d * i
  - a * f * h);
}


/*!
 * \brief Return the angle between two vertices on a plane.
 *
 * The angle is from vertex 1 to vertex 2, positive counterclockwise
 * (CCW).
 *
 * \return The angle value is in the range (\f$ -\pi \cdots \pi \f$) in radians.
 */
double
angle_2D
(
        double x1,
                /*!< X-coordinate of the first vertex (of the pair). */
        double y1,
                /*!< Y-coordinate of the first vertex (of the pair). */
        double x2,
                /*!< X-coordinate of the second vertex (of the pair). */
        double y2
                /*!< Y-coordinate of the second vertex (of the pair). */
)
{
        double dtheta;
        double theta1;
        double theta2;

        theta1 = atan2 (y1, x1);
        theta2 = atan2 (y2, x2);
        dtheta = theta2 - theta1;
        while (dtheta > M_PI)
                dtheta -= 2 * M_PI;
        while (dtheta < -M_PI)
        dtheta += 2 * M_PI;

        return(dtheta);
}


/*!
 * \brief Compute if the coordinates of a point \f$ p \f$ lies inside or
 * outside a polygon \c polygon.
 *
 * \author Paul Bourke <http://www.paulbourke.net/geometry/insidepoly/>
 * Adapted for libDXF by Bert Timmerman <[email protected]>
 *
 * A solution by Philippe Reverdy is to compute the sum of the angles
 * made between the test point and each pair of points making up the
 * polygon.\n
 * If this sum is (\f$ 2\pi \f$) then the point is an interior point,
 * if 0 then the point is an exterior point.\n
 * This also works for polygons with holes given the polygon is defined
 * with a path made up of coincident edges into and out of the hole as
 * is common practice in many CAD packages.\n
 *
 * \note For most of the "point-in-polygon" algorithms above there is a
 * pathological case if the point being queries lies exactly on a
 * vertex.\n
 * The easiest way to cope with this is to test that as a separate
 * process and make your own decision as to whether you want to consider
 * them inside or outside.
 */
int
point_inside_polygon_2D
(
        Point *polygon,
                /*!< A set of polygon vertices. */
        int n,
                /*!< The number of vertices in \c polygon. */
        Point point
                /*!< The point to be tested. */
)
{
        int i;
        double angle = 0;
        Point p1;
        Point p2;

        for (i = 0; i < n; i++)
        {
                p1.x = polygon[i].x - point.x;
                p1.y = polygon[i].y - point.y;
                p2.x = polygon[(i+1)%n].x - point.x;
                p2.y = polygon[(i+1)%n].y - point.y;
                angle += angle_2D (p1.x, p1.y, p2.x, p2.y);
        }

        if (ABS (angle) < PI)
                return (FALSE);
        else
                return (TRUE);
}


/*!
 * \brief Determine the intersection point of two line segments.
 *
 * \author Paul Bourke <http://www.paulbourke.net/geometry/lineline2d/>
 *
 * This describes the technique and algorithm for determining the
 * intersection point of two lines (or line segments) in 2 dimensions.\n
 * \n
 * \image html line_line_intersection_2D.gif
 * \n
 * The equations of the lines are:\n
\f[
  P_a = P_1 + \mu_a (P_2 - P_1)
\f]
\f[
  P_b = P_3 + \mu_b (P_4 - P_3)
\f]
 * Solving for the point where \f$ P_a = P_b \f$ gives the following two
 * equations in two unknowns (\f$ \mu_a \f$ and \f$ \mu_b \f$)\n
\f[
x_1 + \mu_a (x_2 - x_1) = x_3 + \mu_b (x_4 - x_3)
\f]
 * and \n
\f[
y_1 + \mu_a (y_2 - y_1) = y_3 + \mu_b (y_4 - y_3)
\f]
 * Solving gives the following expressions for \f$ \mu_a \f$ and
 * \f$ \mu_b \f$:\n
\f[
  \mu_a = \frac
  {(x_4-x_3)(y_1-y_3)-(y_4-y_3)(x_1-x_3)}
  {(y_4-y_3)(x_2-x_1)-(x_4-x-3)(y_2-y_1)}
\f]
\f[
  \mu_b = \frac
  {(x_2-x_1)(y_1-y_3)-(y_2-y_1)(x_1-x_3)}
  {(y_4-y_3)(x_2-x_1)-(x_4-x-3)(y_2-y_1)}
\f]
 * Substituting either of these into the corresponding equation for the
 * line gives the intersection point.\n
 * For example the intersection point \f$ (x,y) \f$ is:\n
\f[
  x = x_1 + \mu_a (x_2 - x_1)
\f]
\f[
  y = y_1 + \mu_a (y_2 - y_1) 
\f]
 *
 * \note The denominators for the equations for \f$ \mu_a \f$ and
 * \f$ \mu_b \f$ are the same.
 *
 * \note If the denominator for the equations for \f$ \mu_a \f$ and
 * \f$ \mu_b \f$ is 0 then the two lines are parallel.
 *
 * \note If the denominator and numerator for the equations for
 * \f$ \mu_a \f$ and \f$ \mu_b \f$ are 0 then the two lines are coincident.
 *
 * \note The equations apply to lines, if the intersection of line
 * segments is required then it is only necessary to test if
 * \f$ \mu_a \f$ and \f$ \mu_b \f$ lie between 0 and 1.\n
 * Whichever one lies within that range then the corresponding line
 * segment contains the intersection point.\n
 * If both lie within the range of 0 to 1 then the intersection point is
 * within both line segments.
 *
 * \return FALSE if the lines don't intersect.
 */
int
line_line_intersection_point_2D
(
        double x1, double y1,
        double x2, double y2,
        double x3, double y3,
        double x4, double y4,
        double *x, double *y
)
{
        double mua;
        double mub;
        double denom;
        double numera;
        double numerb;

        denom  = (y4-y3) * (x2-x1) - (x4-x3) * (y2-y1);
        numera = (x4-x3) * (y1-y3) - (y4-y3) * (x1-x3);
        numerb = (x2-x1) * (y1-y3) - (y2-y1) * (x1-x3);

        /* Are the line coincident? */
        if (abs (numera) < EPS && abs (numerb) < EPS && abs (denom) < EPS)
        {
                *x = (x1 + x2) / 2;
                *y = (y1 + y2) / 2;
                return (TRUE);
        }

        /* Are the line parallel */
        if (abs (denom) < EPS)
        {
                *x = 0;
                *y = 0;
                return (FALSE);
        }

        /* Is the intersection along the the segments */
        mua = numera / denom;
        mub = numerb / denom;
        if (mua < 0 || mua > 1 || mub < 0 || mub > 1)
        {
                *x = 0;
                *y = 0;
                return (FALSE);
        }

        *x = x1 + mua * (x2 - x1);
        *y = y1 + mua * (y2 - y1);
        return (TRUE);
}


/*!
 * \brief Calculate the line segment \f$ P_aP_b \f$ that is the shortest
 * route between two lines \f$ P_1P_2 \f$ and \f$ P_3P_4 \f$.
 *
 * \author Paul Bourke <http://www.paulbourke.net/geometry/lineline3d/>
 *
 * Two lines in 3 dimensions generally don't intersect at a point, they
 * may be parallel (no intersections) or they may be coincident
 * (infinite intersections) but most often only their projection onto a
 * plane intersect.\n
 * \n
 * \image html line_line_closest_line_3D.gif
 * \n
 * When they don't exactly intersect at a point they can be connected by
 * a line segment, the shortest line segment is unique and is often
 * considered to be their intersection in 3D.\n
 *
 * The following will show how to compute this shortest line segment that
 * joins two lines in 3D, it will as a bi-product identify parallel
 * lines.\n
 * In what follows a line will be defined by two points lying on it, a
 * point on line "a" defined by points \f$ P_1 \f$ and \f$ P_2 \f$ has
 * an equation:\n
\f[
  P_a = P_1 + \mu_a (P_2 - P_1)
\f]
 * similarly a point on a second line "b" defined by points \f$ P_3 \f$
 * and \f$ P_4 \f$ will be written as:\n
\f[
  P_b = P_3 + \mu_b (P_4 - P_3)
\f]
 * The values of \f$ \mu_a \f$ and \f$ \mu_b \f$ range from negative to
 * positive infinity.\n
 * The line segments between \f$ P_1P_2 \f$ and \f$ P_3P_4 \f$have their
 * corresponding \f$ \mu \f$ between 0 and 1.\n
 * There are two approaches to finding the shortest line segment between
 * lines "a" and "b".\n
 * The first is to write down the length of the line segment joining the
 * two lines and then find the minimum.\n
 * That is, minimise the following:\n
\f[
  || P_b - P_a ||^2 
\f]
 * Substituting the equations of the lines gives:\n
\f[
|| P_1 - P_3 + \mu_a (P_2 - P_1) - \mu_b (P_4 - P_3) ||^2
\f]
 * The above can then be expanded out in the (x,y,z) components.\n
 * There are conditions to be met at the minimum, the derivative with
 * respect to \f$ \mu_a \f$ and \f$ \mu_b \f$ must be zero.
 *
 * \note It is easy to convince oneself that the above function only has
 * one minima and no other minima or maxima.\n
 * These two equations can then be solved for \f$ \mu_a \f$ and
 * \f$ \mu_b \f$, the actual intersection points found by substituting
 * the values of mu into the original equations of the line.\n
 *
 * An alternative approach but one that gives the exact same equations
 * is to realise that the shortest line segment between the two lines
 * will be perpendicular to the two lines.\n
 * This allows us to write two equations for the dot product as:\n
\f[
  (P_a - P_b) \cdot (P_2 - P_1) = 0
\f]
\f[
  (P_a - P_b) \cdot (P_4 - P_3) = 0
\f]
 * Expanding these given the equation of the lines: \n
\f[
  (P_1 - P_3 + \mu_a (P_2 - P_1) - \mu_b (P_4 - P_3) ) \cdot (P_2 - P_1) = 0
\f]
\f[
  (P_1 - P_3 + \mu_a (P_2 - P_1) - \mu_b (P_4 - P_3) ) \cdot (P_4 - P_3) = 0
\f]
 * Expanding these in terms of the coordinates (x,y,z) is a nightmare
 * but the result is as follows:\n
\f[
  d_{1321} + \mu_a d_{2121} - \mu_b d_{4321} = 0
\f]
\f[
  d_{1343} + \mu_a d_{4321} - \mu_b d_{4343} = 0
\f]
 * where:\n
\f[
  dmnop = (xm - xn)(xo - xp) + (ym - yn)(yo - yp) + (zm - zn)(zo - zp)
\f]
 * Note that \f$ d_{mnop} = d_{opmn} \f$ \n
 * \n
 * Finally, solving for \f$ \mu_a \f$ gives:
\f[
  \mu_a = (d_{1343} d_{4321} - d_{1321} d_{4343}) / (d_{2121} d_{4343} - d_{4321} d_{4321})
\f]
 * and back-substituting gives \f$ \mu_b \f$ :\n
\f[
  \mu_b = (d_{1343} + \mu_a d_{4321}) / d_{4343}
\f]
 *
 * \return FALSE if no solution exists.
 */
int
line_line_closest_line_3D
(
        XYZ p1,
        XYZ p2,
        XYZ p3,
        XYZ p4,
        XYZ *pa,
        XYZ *pb,
        double *mua,
        double *mub
)
{
        XYZ p13;
        XYZ p43;
        XYZ p21;
        double d1343;
        double d4321;
        double d1321;
        double d4343;
        double d2121;
        double numer;
        double denom;

   p13.x = p1.x - p3.x;
   p13.y = p1.y - p3.y;
   p13.z = p1.z - p3.z;
   p43.x = p4.x - p3.x;
   p43.y = p4.y - p3.y;
   p43.z = p4.z - p3.z;

   if (ABS (p43.x) < EPS && ABS (p43.y) < EPS && ABS (p43.z) < EPS)
      return (FALSE);

   p21.x = p2.x - p1.x;
   p21.y = p2.y - p1.y;
   p21.z = p2.z - p1.z;

   if (ABS (p21.x) < EPS && ABS (p21.y) < EPS && ABS (p21.z) < EPS)
      return (FALSE);

   d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z;
   d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z;
   d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z;
   d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z;
   d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z;

   denom = d2121 * d4343 - d4321 * d4321;

   if (ABS (denom) < EPS)
      return (FALSE);

   numer = d1343 * d4321 - d1321 * d4343;

   *mua = numer / denom;
   *mub = (d1343 + d4321 * (*mua)) / d4343;

   pa->x = p1.x + *mua * p21.x;
   pa->y = p1.y + *mua * p21.y;
   pa->z = p1.z + *mua * p21.z;
   pb->x = p3.x + *mub * p43.x;
   pb->y = p3.y + *mub * p43.y;
   pb->z = p3.z + *mub * p43.z;

   return (TRUE);
}


/*!
 * \brief Return whether a polygon in 2D is concave or convex.
 *
 * \author Paul Bourke <http://www.paulbourke.net/geometry/clockwise/>
 *
 * The following describes a method for determining whether or not a
 * polygon has its vertices ordered clockwise or anticlockwise for both
 * convex and concave polygons.\n
 * A polygon will be assumed to be described by N vertices, ordered\n
\f[
(x_0, y_0), (x_1, y_1), (x_2, y_2), \cdots (x_{n-1}, y_{n-1})
\f]
 * A simple test of vertex ordering for convex polygons is based on
 * considerations of the cross product between adjacent edges.\n
 * If the crossproduct is positive then it rises above the plane (z axis
 * up out of the plane) and if negative then the cross product is into
 * the plane.\n
 * \n
\f[
cross product = ((x_i - x_{i-1}), (y_i - y_{i-1}))
  \cdot ((x_{i+1} - x_i), (y_{i+1} - y_i))
\f]
 * \n
\f[
= (x_i - x_{i-1})
  \cdot (y_{i+1} - y_i) - (y_i - y_{i-1})
  \cdot (x_{i+1} - x_i)
\f]
 * \n
 * A positive cross product means we have a counterclockwise polygon.\n
 * \n
 * To determine the vertex ordering for concave polygons one can use a
 * result from the calculation of polygon areas, where the area is given
 * by:\n
 * \n
\f[
  A = \frac {1} {2}
  \cdot \sum _{i = 0} ^{N - 1}
  {({x_i} \cdot {y_{i + 1}} - {x_{i + 1}} \cdot {y_i})} 
\f]
 * \n
 * If the above expression is positive then the polygon is ordered
 * counter clockwise otherwise if it is negative then the polygon
 * vertices are ordered clockwise.
 *
 * \note It is assumed that the polygon is simple (does not intersect
 * itself or have holes).
 *
 * \return 0 for incomputables eg: colinear points,\n
 * CONVEX == 1 ,\n
 * CONCAVE == -1 .
 */
int
polygon_is_convex_2D
(
        XY *p,
        int n
)
{
        int i;
        int j;
        int k;
        int flag = 0;
        double z;

        if (n < 3)
        return (0);

        for (i=0;i<n;i++)
        {
                j = (i + 1) % n;
                k = (i + 2) % n;
                z  = (p[j].x - p[i].x) * (p[k].y - p[j].y);
                z -= (p[j].y - p[i].y) * (p[k].x - p[j].x);

                if (z < 0)
                        flag |= 1;
                else if (z > 0)
                        flag |= 2;

                if (flag == 3)
                        return (CONCAVE);
        }

        if (flag != 0)
                return (CONVEX);
        else
                return (0);
}


/*!
 * \brief Calculate the area of the given closed polyline.
 *
 * \author Paul Bourke <http://www.paulbourke.net/geometry/polyarea/>
 *
 * The problem of determining the area of a polygon seems at best messy
 * but the final formula is particularly simple.\n
 * Consider a polygon made up of line segments between N vertices
 * \f$ (x_i, y_i) \f$ , i = 0 to (N - 1).\n
 * The last vertex \f$ (x_N, y_N) \f$ is assumed to be the same as the
 * first, ie: the polygon is closed.\n
 * \n
 * \image html dxf_hatch_boundary_path_polyline_area1.gif
 * \n
 * The area is given by:
 * \n
\f[
  A = \frac {1} {2}
  \cdot \sum _{i = 0} ^{N - 1}
  {({x_i} \cdot {y_{i + 1}} - {x_{i + 1}} \cdot {y_i})} 
\f]
 * \n
 *
 * \note
 * <ul>
 *   <li> For polygons with holes.\n
 *        The holes are usually defined by ordering the vertices of the
 *        enclosing polygon in the opposite direction to those of the
 *        holes.\n
 *        This algorithm still works except that the absolute value
 *        should be taken after adding the polygon area to the area of
 *        all the holes.\n
 *        That is, the holes areas will be of opposite sign to the
 *        bounding polygon area.
 *   <li> The sign of the area expression above (without the absolute
 *        value) can be used to determine the ordering of the vertices
 *        of the polygon.\n
 *        If the sign is positive then the polygon vertices are ordered
 *        counterclockwise (CCW) about the normal, otherwise clockwise.
 * </ul>
 *
 * \return The area of the polygon.
 */
double
polygon_area
(
        Point *polygon,
        int N
)
{
        int i;
        int j;
        double area = 0;

        for (i = 0; i < N; i++)
        {
                j = (i + 1) % N;
                area += polygon[i].x * polygon[j].y;
                area -= polygon[i].y * polygon[j].x;
        }

        area /= 2;
        return (area < 0 ? -area : area);
}


/*!
 * \brief Determine whether or not the line segment \f$ p_1, p_2 \f$ intersects
 * the 3 vertex facet bounded by \f$ p_a, p_b, p_c \f$.
 *
 * The labeling and naming conventions for the line segment and the
 * facet are shown in the following diagram:\n
 * \n
 * \image html line_facet1.gif
 * \n
 * The equation of the line is: \n
\f[
  p = p_1 + \mu (p_2 - p_1)
\f]
 * The equation of the plane is: \n
\f[
  a \cdot x + b \cdot y + c \cdot z + d = 0
\f]
 * The following will find the intersection point (if it exists) between
 * a line segment and a planar 3 vertex facet.\n
 * The mathematics and solution can also be used to find the
 * intersection between a plane and line, a simpler problem.\n
 * The intersection between more complex polygons can be found by first
 * triangulating them into multiple 3 vertex facets.\n
 * \n
 * The procedure will be implemented given the line segment defined by
 * its two end points and the facet bounded by its three vertices.\n
 * The solution involves the following steps:\n
 * <ol>
 *   <li> Check the line and plane are not parallel.
 *   <li> Find the intersection of the line, on which the given line
 *        segment lies, with the plane containing the facet.
 *   <li> Check that the intersection point lies along the line segment.
 *   <li> Check that the intersection point lies within the facet.
 * </ol>
 * \n
 * The intersection point \f$ P \f$ is found by substituting the
 * equation for the line \f$ P = P_1 + \mu (P_2 - P_1) \f$ into the
 * equation for the plane \f$ A \cdot x + B \cdot y + C \cdot z + D = 0 \f$.
 * \n
 * Note that the values of \f$ A,B,C \f$ are the components of the
 * normal to the plane which can be found by taking the cross product of
 * any two normalised edge vectors, for example:\n
\f[
  (A,B,C) = (P_b - P_a) cross (P_c - P_a)
\f]
 * Then \f$ D \f$ is found by substituting one vertex into the equation
 * for the plane for example:\n
\f[
  A \cdot P_{ax} + B \cdot P_{ay} + C \cdot P_{az} = -D 
\f]
 * This gives an expression for \f$ \mu \f$ from which the point of
 * intersection \f$ P \f$ can be found using the equation of the line.\n
\f[
  mu = \frac
  {( D + A \cdot P_{1x} + B \cdot P_{1y} + C \cdot P_{1z} )}
  {( A \cdot (P_{1x} - P_{2x}) + B \cdot (P_{1y} - P_{2y}) + C \cdot (P_{1z} - P_{2z}) )}
\f]
 * If the denominator above is 0 then the line is parallel to the plane
 * and they don't intersect.\n
 * For the intersection point to lie on the line segment, \f$ \mu \f$
 * must be between 0 and 1.\n
 * \n
 * Lastly, we need to check whether or not the intersection point lies
 * within the planar facet bounded by \f$ P_a, P_b, P_c \f$.\n
 * \n
 * The method used here relies on the fact that the sum of the internal
 * angles of a point on the interior of a triangle is \f$ 2\pi \f$,
 * points outside the triangular facet will have lower angle sums.\n
 * \n
 * \image html line_facet2.gif
 * \n
 * If we form the unit vectors \f$ P_{a1}, P_{a2}, P_{a3} \f$ as follows
 * (P is the point being tested to see if it is in the interior):\n
\f[
  P_{a1} = (P_a - P) / |(P_a - P)|
\f]
\f[
  P_{a2} = (P_b - P) / |(P_b - P)|
\f]
\f[
  P_{a3} = (P_c - P) / |(P_c - P)|
\f]
 * the angles are: \n
\f[
  a_1 = \arccos (P_{a1} \cdot P_{a2})
\f]
\f[
  a_2 = \arccos (P_{a2} \cdot P_{a3})
\f]
\f[
  a_3 = \arccos (P_{a3} \cdot P_{a1})
\f]
 *
 * \return Return true/false and the intersection point p.
 */
int
line_facet
(
        p1,
        p2,
        pa,
        pb,
        pc,
        p
)

XYZ p1,p2,pa,pb,pc,*p;

{
        double d;
        double a1;
        double a2;
        double a3;
        double total;
        double denom;
        double mu;
        XYZ n;
        XYZ pa1;
        XYZ pa2;
        XYZ pa3;

        /* Calculate the parameters for the plane */
        n.x = (pb.y - pa.y) * (pc.z - pa.z) - (pb.z - pa.z) * (pc.y - pa.y);
        n.y = (pb.z - pa.z) * (pc.x - pa.x) - (pb.x - pa.x) * (pc.z - pa.z);
        n.z = (pb.x - pa.x) * (pc.y - pa.y) - (pb.y - pa.y) * (pc.x - pa.x);
        Normalise (&n);
        d = - n.x * pa.x - n.y * pa.y - n.z * pa.z;

        /* Calculate the position on the line that intersects the plane */
        denom = n.x * (p2.x - p1.x) + n.y * (p2.y - p1.y) + n.z * (p2.z - p1.z);

        if (abs (denom) < EPS) /* Line and plane don't intersect */
                return (FALSE);

        mu = - (d + n.x * p1.x + n.y * p1.y + n.z * p1.z) / denom;
        p->x = p1.x + mu * (p2.x - p1.x);
        p->y = p1.y + mu * (p2.y - p1.y);
        p->z = p1.z + mu * (p2.z - p1.z);

        if (mu < 0 || mu > 1) /* Intersection not along line segment */
        return (FALSE);

        /* Determine whether or not the intersection point is bounded by pa,pb,pc */
        pa1.x = pa.x - p->x;
        pa1.y = pa.y - p->y;
        pa1.z = pa.z - p->z;
        Normalise (&pa1);
        pa2.x = pb.x - p->x;
        pa2.y = pb.y - p->y;
        pa2.z = pb.z - p->z;
        Normalise (&pa2);
        pa3.x = pc.x - p->x;
        pa3.y = pc.y - p->y;
        pa3.z = pc.z - p->z;
        Normalise (&pa3);
        a1 = pa1.x * pa2.x + pa1.y * pa2.y + pa1.z * pa2.z;
        a2 = pa2.x * pa3.x + pa2.y * pa3.y + pa2.z * pa3.z;
        a3 = pa3.x * pa1.x + pa3.y * pa1.y + pa3.z * pa1.z;
        total = (acos (a1) + acos (a2) + acos (a3)) * RTOD;

        if (abs (total - 360) > EPS)
                return (FALSE);

        return(TRUE);
}


/*!
 * \brief Calculate the intersection of a line and a sphere.
 *
 * \image html line_sphere1.gif
 *
 * The line segment is defined from \f$ p_1 \f$ to \f$ p_2 \f$ .\n
 * The sphere is of radius \f$ r \f$ and centered at \f$ p_3 \f$ .\n
 * There are potentially two points of intersection given by:\n
\f[
  p = p_1 + \mu_1 (p_2 - p_1)
\f]
\f[
  p = p_1 + \mu_2 (p_2 - p_1)
\f]
 *
 * Points \f$ P (x,y) \f$ on a line defined by two points
 * \f$ P_1 (x_1, y_1, z_1) \f$ and \f$ P_2 (x_2, y_2, z_2) \f$ is
 * described by:
\f[
  P = P_1 + \mu (P_2 - P_1)
\f]
 * or in each coordinate:
\f[
  x = x_1 + \mu (x_2 - x_1)
\f]
\f[
  y = y_1 + \mu (y_2 - y_1)
\f]
\f[
  z = z_1 + \mu (z_2 - z_1)
\f]
 * A sphere centered at \f$ P_3 (x_3, y_3, z_3) \f$ with radius
 * \f$ r \f$ is described by:\n
\f[
  (x - x_3)^2 + (y - y_3)^2 + (z - z_3)^2 = r^2
\f]
 * Substituting the equation of the line into the sphere gives a
 * quadratic equation of the form:\n
\f[
  a \mu^2 + b \mu + c = 0
\f]
 * where:\n
\f[
  a = (x_2 - x_1)^2 + (y_2 - y_1)^2 + (z_2 - z_1)^2
\f]
\f[
  b = 2[ (x_2 - x_1) (x_1 - x_3) + (y_2 - y_1) (y_1 - y_3)
    + (z_2 - z_1) (z_1 - z_3) ]
\f]
\f[
  c = x_3^2 + y_3^2 + z_3^2 + x_1^2 + y_1^2 + z_1^2
    - 2[x_3 x_1 + y_3 y_1 + z_3 z_1] - r^2
\f]
 * The solutions to this quadratic are described by:\n
\f[
  \frac
    {  -b \pm \sqrt {b^2 - 4 \cdot a \cdot c}}
    {2a}
\f]
 * The exact behaviour is determined by the expression within the square
 * root:\n
\f[
  \sqrt {b^2 - 4 \cdot a \cdot c}
\f]
 * If this is less than 0 then the line does not intersect the sphere.\n
 * \n
 * If it equals 0 then the line is a tangent to the sphere intersecting
 * it at one point, namely at:
\f[
 \mu = -b/2a
\f]
 * If it is greater then 0 the line intersects the sphere at two
 * points.\n
 * \n
 * To apply this to two dimensions, that is, the intersection of a line
 * and a circle simply remove the \f$ z \f$ component from the above
 * mathematics.
 * \n
 * <h3>Line Segment</h3>
 * \n
 * For a line segment between \f$ P_1 \f$ and \f$ P_2 \f$ there are 5
 * cases to consider.\n
 * <ol>
 *   <li> Line segment doesn't intersect and on outside of sphere, in
 *        which case both values of \f$ \mu \f$ will either be less than
 *        0 or greater than 1.
 *   <li> Line segment doesn't intersect and is inside sphere, in which
 *        case one value of u will be negative and the other greater
 *        than 1.
 *   <li> Line segment intersects at one point, in which case one value
 *        of \f$ \mu \f$ will be between 0 and 1 and the other not.
 *   <li> Line segment intersects at two points, in which case both
 *        values of \f$ \mu \f$ will be between 0 and 1.
 *   <li> Line segment is tangential to the sphere, in which case both
 *        values of u will be the same and between 0 and 1.
 * <ol>
 * \n
 * When dealing with a line segment it may be more efficient to first
 * determine whether the line actually intersects the sphere or circle.\n
 * This is achieved by noting that the closest point on the line through
 * \f$ P_1P_2 \f$ to the point \f$ P_3 \f$ is along a perpendicular from
 * \f$ P_3 \f$ to the line.\n
 * In other words if \f$ P \f$ is the closest point on the line then:\n
\f[
  (P_3 - P) \cdot (P_2 - P_1) = 0
\f]
 * \n
Substituting the equation of the line into this:\n
\f[
 [P_3 - P_1 - \mu (P_2 - P_1)] \cdot (P_2 - P_1) = 0
\f]
Solving the above for:\n
\f[
  u = \frac
    {(x_3 - x_1)(x_2 - x_1) + (y_3 - y_1)(y_2 - y_1) + (z_3 - z_1)
      (z_2 - z_1)}
    {(x_2 - x_1)(x_2 - x_1) + (y_2 - y_1)(y_2 - y_1) + (z_2 - z_1)
      (z_2 - z_1)}
\f]
 * If \f$ \mu \f$ is not between 0 and 1 then the closest point is not
 * between \f$ P_1 \f$ and \f$ P_2 \f$ .\n
 * Given \f$ \mu \f$ , the intersection point can be found, it must also
 * be less than the radius \f$ r \f$ .\n
 * If these two tests succeed then the earlier calculation of the actual
 * intersection point can be applied.
 *
 * \return FALSE if the ray doesn't intersect the sphere.
 */
int line_sphere
(
        XYZ p1,
        XYZ p2,
        XYZ p3,
        double r,
        double *mu1,
        double *mu2
)
{
        double a;
        double b;
        double c;
        double bb4ac;
        XYZ dp;

        dp.x = p2.x - p1.x;
        dp.y = p2.y - p1.y;
        dp.z = p2.z - p1.z;
        a = dp.x * dp.x + dp.y * dp.y + dp.z * dp.z;
        b = 2 * (dp.x * (p1.x - p3.x)
                + dp.y * (p1.y - p3.y) + dp.z * (p1.z - p3.z));
        c = p3.x * p3.x + p3.y * p3.y + p3.z * p3.z;
        c += p1.x * p1.x + p1.y * p1.y + p1.z * p1.z;
        c -= 2 * (p3.x * p1.x + p3.y * p1.y + p3.z * p1.z);
        c -= r * r;
        bb4ac = b * b - 4 * a * c;

        if (abs (a) < EPS || bb4ac < 0)
        {
                *mu1 = 0;
                *mu2 = 0;
                return(FALSE);
        }

        *mu1 = (-b + sqrt (bb4ac)) / (2 * a);
        *mu2 = (-b - sqrt (bb4ac)) / (2 * a);
        return(TRUE);
}


/*!
 * \brief Clip a 3 vertex facet in place.
 *
 * The 3 point facet is defined by vertices \f$ p[0],p[1],p[2] \f$ .\n
 * There must be a fourth point \f$ "p[3]" \f$ as a 4 point facet may
 * result.\n
 * The normal to the plane is n.\n
 * A point on the plane is \f$ p0 \f$ .\n
 * The side of the plane containing the normal is clipped away.\n
 * \n
 * This note along with the source code at the end clips a 3 vertex
 * facet by an arbitrary plane.\n
 * The facet is described by its 3 vertices, the clipping plane is
 * specified by its normal and one point on the plane.\n
 * The clipping side of the plane is taken to be that one containing the
 * normal, in the diagram below this is the right side of the clipping
 * plane.\n
 * The standard equation of a plane is:\n
\f[
    A \cdot x + B \cdot y + C \cdot z + D = 0
\f]
 * where \f$ (A,B,C) \f$ is the unit normal.\n
 * The value of \f$ D \f$ is determined by substituting in the known
 * point \f$(P_x, P_y, P_z) \f$ on the plane, namely
\f[
    D = - (A \cdot P_x + B \cdot P_y + C \cdot P_z)
\f]
 * For a vertex \f$ (Q_x,Q_y,Q_z) \f$ the expression
\f[
    side (Q) = A \cdot Q_x + B \cdot Q_y + C \cdot Q_z + D
\f]
 * can be used to determine which side of the plane the vertex lies on.\n
 * If it is positive the point lies on the same side as the normal, if
 * negative it lies on the other side, if zero it lies on the plane.\n
 * \n
 * After determining if an edge intersects the cutting plane it is
 * necessary to calculate the intersection point as that will become a
 * vertex of the clipped facet.\n
 * Let the edge be between two points \f$ P_0 \f$ and \f$ P_1 \f$ , the
 * equation of the points along the line segment
\f[
    P = P_0 + \mu (P_1 - P_0)
\f]
 * where \f$ \mu \f$ lies between 0 and 1.\n
 * Substituting this into the expression for the plane:\n
\f[
    A (P_{0x} + \mu (P_{1x} - P_{0x}))
    + B (P_{0y} + \mu (P_{1y} - P_{0y}))
    + C (P_{0z} + \mu (P_{1z} - P_{0z})) + D = 0
\f]
 * Solving for \f$ \mu \f$ :\n
\f[
    \mu = - side (P_0) / (side (P_1) - side (P_0)) 
\f]
 * Substituting this into the equation of the line gives the actual
 * intersection point.\n
 * \n
 * There are 8 different cases to consider, they are illustrated in the
 * table below and appear in order in the source code.\n
 * The top two cases are when all the points are on either side of the
 * clipping plane.\n
 * The three cases on the left are when there is only one vertex on the
 * clipping side of the plane, the three cases on the right occur when
 * there are two vertices on the clipping side of the plane.\n
 *
 * \image html clip_facet1.gif
 *
 * \note
 * <ul>
 *   <li> When there is only one point on the clipping side the
 *        resulting facet has 4 vertices, if the underlying database
 *        only deals with 3 vertex facets then this can be bisected
 *        using a number of methods.
 *   <li> When the facets to be clipped have more than 3 vertices then
 *        they can be subdivided first and each segment clipped
 *        individually.
 *   <li> The algorithm as presented has no divide by zero problems.
 *   <li> The source below does the clipping in place, the order in
 *        which the vertices are calculated is important for the 3 cases
 *        when there is one point on the clipping side of the plane.
 * </ul>
 *
 * \return Return the number of vertices in the clipped polygon.
*/
int
clip_facet
(
        p,
        n,
        p0
)
XYZ *p;
XYZ n;
XYZ p0;
{
   double A,B,C,D;
   double l;
   double side[3];
   XYZ q;

   /*
      Determine the equation of the plane as
      Ax + By + Cz + D = 0
   */
   l = sqrt(n.x*n.x + n.y*n.y + n.z*n.z);
   A = n.x / l;
   B = n.y / l;
   C = n.z / l;
   D = -(n.x*p0.x + n.y*p0.y + n.z*p0.z);

   /*
      Evaluate the equation of the plane for each vertex
      If side < 0 then it is on the side to be retained
      else it is to be clipped
   */
   side[0] = A*p[0].x + B*p[0].y + C*p[0].z + D;
   side[1] = A*p[1].x + B*p[1].y + C*p[1].z + D;
   side[2] = A*p[2].x + B*p[2].y + C*p[2].z + D;

   /* Are all the vertices are on the clipped side */
   if (side[0] >= 0 && side[1] >= 0 && side[2] >= 0)
      return(0);

   /* Are all the vertices on the not-clipped side */
   if (side[0] <= 0 && side[1] <= 0 && side[2] <= 0)
      return(3);

   /* Is p0 the only point on the clipped side */
   if (side[0] > 0 && side[1] < 0 && side[2] < 0) {
      q.x = p[0].x - side[0] * (p[2].x - p[0].x) / (side[2] - side[0]);
      q.y = p[0].y - side[0] * (p[2].y - p[0].y) / (side[2] - side[0]);
      q.z = p[0].z - side[0] * (p[2].z - p[0].z) / (side[2] - side[0]);
      p[3] = q;
      q.x = p[0].x - side[0] * (p[1].x - p[0].x) / (side[1] - side[0]);
      q.y = p[0].y - side[0] * (p[1].y - p[0].y) / (side[1] - side[0]);
      q.z = p[0].z - side[0] * (p[1].z - p[0].z) / (side[1] - side[0]);
      p[0] = q;
      return(4);
   }

   /* Is p1 the only point on the clipped side */
   if (side[1] > 0 && side[0] < 0 && side[2] < 0) {
      p[3] = p[2];
      q.x = p[1].x - side[1] * (p[2].x - p[1].x) / (side[2] - side[1]);
      q.y = p[1].y - side[1] * (p[2].y - p[1].y) / (side[2] - side[1]);
      q.z = p[1].z - side[1] * (p[2].z - p[1].z) / (side[2] - side[1]);
      p[2] = q;
      q.x = p[1].x - side[1] * (p[0].x - p[1].x) / (side[0] - side[1]);
      q.y = p[1].y - side[1] * (p[0].y - p[1].y) / (side[0] - side[1]);
      q.z = p[1].z - side[1] * (p[0].z - p[1].z) / (side[0] - side[1]);
      p[1] = q;
      return(4);
   }

   /* Is p2 the only point on the clipped side */
   if (side[2] > 0 && side[0] < 0 && side[1] < 0) {
      q.x = p[2].x - side[2] * (p[0].x - p[2].x) / (side[0] - side[2]);
      q.y = p[2].y - side[2] * (p[0].y - p[2].y) / (side[0] - side[2]);
      q.z = p[2].z - side[2] * (p[0].z - p[2].z) / (side[0] - side[2]);
      p[3] = q;
      q.x = p[2].x - side[2] * (p[1].x - p[2].x) / (side[1] - side[2]);
      q.y = p[2].y - side[2] * (p[1].y - p[2].y) / (side[1] - side[2]);
      q.z = p[2].z - side[2] * (p[1].z - p[2].z) / (side[1] - side[2]);
      p[2] = q;
      return(4);
   }

   /* Is p0 the only point on the not-clipped side */
   if (side[0] < 0 && side[1] > 0 && side[2] > 0) {
      q.x = p[0].x - side[0] * (p[1].x - p[0].x) / (side[1] - side[0]);
      q.y = p[0].y - side[0] * (p[1].y - p[0].y) / (side[1] - side[0]);
      q.z = p[0].z - side[0] * (p[1].z - p[0].z) / (side[1] - side[0]);
      p[1] = q;
      q.x = p[0].x - side[0] * (p[2].x - p[0].x) / (side[2] - side[0]);
      q.y = p[0].y - side[0] * (p[2].y - p[0].y) / (side[2] - side[0]);
      q.z = p[0].z - side[0] * (p[2].z - p[0].z) / (side[2] - side[0]);
      p[2] = q;
      return(3);
   }

   /* Is p1 the only point on the not-clipped side */
   if (side[1] < 0 && side[0] > 0 && side[2] > 0) {
      q.x = p[1].x - side[1] * (p[0].x - p[1].x) / (side[0] - side[1]);
      q.y = p[1].y - side[1] * (p[0].y - p[1].y) / (side[0] - side[1]);
      q.z = p[1].z - side[1] * (p[0].z - p[1].z) / (side[0] - side[1]);
      p[0] = q;
      q.x = p[1].x - side[1] * (p[2].x - p[1].x) / (side[2] - side[1]);
      q.y = p[1].y - side[1] * (p[2].y - p[1].y) / (side[2] - side[1]);
      q.z = p[1].z - side[1] * (p[2].z - p[1].z) / (side[2] - side[1]);
      p[2] = q;
      return(3);
   }

   /* Is p2 the only point on the not-clipped side */
   if (side[2] < 0 && side[0] > 0 && side[1] > 0) {
      q.x = p[2].x - side[2] * (p[1].x - p[2].x) / (side[1] - side[2]);
      q.y = p[2].y - side[2] * (p[1].y - p[2].y) / (side[1] - side[2]);
      q.z = p[2].z - side[2] * (p[1].z - p[2].z) / (side[1] - side[2]);
      p[1] = q;
      q.x = p[2].x - side[2] * (p[0].x - p[2].x) / (side[0] - side[2]);
      q.y = p[2].y - side[2] * (p[0].y - p[2].y) / (side[0] - side[2]);
      q.z = p[2].z - side[2] * (p[0].z - p[2].z) / (side[0] - side[2]);
      p[0] = q;
      return(3);
   }

   /* Shouldn't get here */
   return(-1);
}


/*!
 * \brief Function for parsing an int while reading a DxfHeader
 *
 * \return \c EXIT_SUCCESS when done, or \c EXIT_FAILURE when an error
 * occurred.
 */
static int
dxf_read_header_parse_int
(
        DxfFile *fp,  /*!< DXF file handler.\n */
        const char *header_var,
        int *value,
        int acad_version_number,
        int valid_acad_version
)
{
#if DEBUG
        DXF_DEBUG_BEGIN
#endif
        char temp_string[255];
        int n;

        if (strcmp (temp_string, header_var)
          && acad_version_number >= valid_acad_version)
        {
                if (fscanf (fp->fp, "%i\n%i\n", &n, value) >= 1)
                        return (EXIT_SUCCESS);
                else
                        return (EXIT_FAILURE);
        }
#if DEBUG
        DXF_DEBUG_END
#endif
        return (EXIT_SUCCESS);
}


/*!
 * \brief Function to generate and return random numbers.
 */
int *
get_random ()
{
        static int r[10];
        int i;

        /* set the seed */
        srand ((unsigned) time (NULL));
        for (i = 0; i < 10; ++i)
        {
                r[i] = rand ();
                printf ("%d\n", r[i] );
        }
        return r;
}
 
/*!
 * \brief Test function to call above defined  get_random()function.
 */
int
test_get_random ()
{
        /* a pointer to an int */
        int *p;
        int i;

        p = get_random ();
        for (i = 0; i < 10; i++)
        {
                printf ("*(p + [%d]) : %d\n", i, *(p + i));
        }
        return 0;
}


double
magnitude (XYZ *Point1, XYZ *Point2)
{
        XYZ Vector;

        Vector.X = Point2->X - Point1->X;
        Vector.Y = Point2->Y - Point1->Y;
        Vector.Z = Point2->Z - Point1->Z;

        return (double) sqrt (Vector.X * Vector.X + Vector.Y * Vector.Y + Vector.Z * Vector.Z);
}


int
distance_point_line (XYZ *Point, XYZ *LineStart, XYZ *LineEnd, float *Distance)
{
        double LineMag;
        double U;
        XYZ Intersection;
 
        LineMag = magnitude (LineEnd, LineStart);
        U = (((Point->X - LineStart->X) * (LineEnd->X - LineStart->X))
          + ((Point->Y - LineStart->Y) * (LineEnd->Y - LineStart->Y))
          + ((Point->Z - LineStart->Z) * (LineEnd->Z - LineStart->Z)))
          / (LineMag * LineMag);
        if (U < 0.0f || U > 1.0f)
                return 0; // closest point does not fall within the line segment
        Intersection.X = LineStart->X + U * (LineEnd->X - LineStart->X);
        Intersection.Y = LineStart->Y + U * (LineEnd->Y - LineStart->Y);
        Intersection.Z = LineStart->Z + U * (LineEnd->Z - LineStart->Z);
        *Distance = magnitude (Point, &Intersection);
        return 1;
}


void
test_distance_point_line (void)
{
        XYZ LineStart, LineEnd, Point;
        double Distance;

        LineStart.X = 50.0;
        LineStart.Y = 80.0;
        LineStart.Z = 300.0;
        LineEnd.X = 50.0;
        LineEnd.Y = -800.0;
        LineEnd.Z = 1000.0;
        Point.X = 20.0;
        Point.Y = 1000.0;
        Point.Z = 400.0;
        if (distance_point_line (&Point, &LineStart, &LineEnd, &Distance))
                printf ("closest point falls within line segment, distance = %f\n", Distance);
        else
                printf ("closest point does not fall within line segment\n");
        LineStart.X = 0.0;
        LineStart.Y = 0.0;
        LineStart.Z = 50.0;
        LineEnd.X = 0.0;
        LineEnd.Y = 0.0;
        LineEnd.Z = -50.0;
        Point.X = 10.0;
        Point.Y = 50.0;
        Point.Z = 10.0;
        if (distance_point_line (&Point, &LineStart, &LineEnd, &Distance))
                printf ("closest point falls within line segment, distance = %f\n", Distance);
        else
                printf ("closest point does not fall within line segment\n");
}


/* EOF */