Computational Geometry

A tour of the CompG package

Adriano Pascoletti

This notebook gives a thorough description of the package CompG.m which provides a rather complete set of basic operations arising from problems of computational geometry  in the plane.
The basic entities defined by the package are points, lines, segments, paths, polygons, and circular arcs, as well as angles, signed areas, affine transformations, cartesian and parametric equations of lines and axes. 
Mathematica's Sort and Ordering functions are generalized to polar and lexicographic sorting and ordering. 
The point location problems considered are those involving lines, halflines, segments, convex polygons, simple polygons, and circles. 
The package can compute the intersection of any pair of straight objects determining both the intersection set (empty, single point, segment,...) and the robustness of the intersection set w.r. to perturbations (strong and weak intersections).
The last two sections of the tour present two nontrivial applications: a nonoptimal  (O((n^2))) triangulation technique for simple polygons based on the existence of diagonals, and an optimal (O(n)) triangulation of  y-monotone polygons. 

Mathematica version 5 is required.

Revised June 2007

Preliminaries

References

[1] T.H. Cormen, C.E. Leiserson and R.L. Rivest, Introduction to Algorithms, The MIT Press, 1991

[2] M. de Berg, M. van Kreveld, M. Overmars, O. Scwarzkopf, Computational geometry - Algorithms and Applications, Springer, 1997

[3] J. O’Rourke, Computational Geometry in C, Cambridge University Press, 1994

[4] A. Pascoletti, Basic Plane Computational Geometry, MathSource, 2000, (http://library.wolfram.com/infocenter/MathSource/805/)

[5] F.P. Preparata and M.I. Shamos,  Computational Geometry: An Introduction, Springer Verlag, 1985

Load the package

If the package CompG.m doesn't reside in a directory searched by Mathematica  one must append to the Mathematica's variable $Path the directory containing CompG.m

AppendTo[$Path, DirectoryName[ ... ]];

and load the package

<< CompG`

The complete list of symbols defined by the package is easily got:

?CompG`*

Points, segments, halflines, lines

Points

Points are represented by vectors with two numeric or nonnumeric components. Two points coincide if their difference equals the zero vector (see Equal), so two points whose difference is {0, 0}, {0. , 0}, {0, 0. }, or  {0. , 0. }  are considered as coincident.

PointQ[list] yields True if list is a two-component vector.

NPointQ[list] yields True if list is a vector with two numeric components.

OriginQ[list] yields True if list is a vector with two components drawn from the set {0, 0. }.

A line is represented by an ordered pair of distinct points and it is naturally oriented  from the first to the second point.

An ordered pair of distinct points also defines a half-line emanating from the first and passing through the second point, and a (nondegenerate) directed  segment.

A segment is degenerate if the difference of its endpoints is the null vector. The empty list {} is considered as a degenerate segment.

SegmentQ[s] gives True if s is a list consisting of two points.

NSegmentQ[s] gives True if s is a list consisting of two points with numeric components.

HalfLineQ[s], LineQ[s] give True if s is a list consisting of two distinct points.

NHalfLineQ[s], NLineQ[s] give True if s is a list consisting of two distinct points with numeric components.

Test

Testing points

Clear[a, b]; 
    PointQ /@ {{2, 3}, {a, b}, {2, a}}
NPointQ /@ {{2, 3}, {a, b}, {2, a}}

The following points are the four different representations of the origin

OriginQ /@ {{0, 0}, {0. , 0. }, {0, 0. }, {0. , 0}}

Testing pairs of points

In the following list of segments s only the first one has numeric endpoints:

Clear[a, b, c, d]; 
    s = {{{2, 3}, {3, 4}}, {{a, b}, {c, d}}, {{2, 3}, {c, d}}, {{a, b}, {a, 
            b}}}; 
{SegmentQ /@ slist,  NSegmentQ /@ slist}

The pair {{2,3},{3,4}} defines a line, a halfline, and a (nondegenerate) segment; the second and the third pairs of points of s define nonnumeric lines, halflines and segments. The last pair does not determine a valid line or halfline, but it defines a degenerate nonnumeric segment:

{LineQ /@ s, NLineQ /@ s}
{HalfLineQ /@ s, NHalfLineQ /@ s}
{SegmentQ /@ s, NSegmentQ /@ s}

Paths

Description

A path is a list  {(Subscript[p, 1]), …, (Subscript[p, n])} of three or more points. A path is open if  OriginQ[(Subscript[p, n]) - (Subscript[p, 1])] yields False,  otherwise it is closed. A numeric path is a path defined by a list of numeric points.
A path is degenerate if {(Subscript[p, i]), (Subscript[p, i + 1])} is a degenerate segment for some index i (with the obvious identification (Subscript[p, n + 1]) = (Subscript[p, 1])). 
An open path can be closed by appending to its first point, and a closed one can be opened by removing its last element as long as it equals the first point .

PathQ[p] yields True if p is a list of three or more points.

NPathQ[p]  yields True if p is a list of three or more numeric points.

OpenPathQ[p] yields True if p is a path with distinct endpoints.

ClosedPathQ[p] yields True if p is a path with coincident endpoints.

NondegeneratePathQ[p] yields True if the path  p defines a sequence of nondegenerate segments.

DegeneratePathQ[p] yields True if the path  p contains a degenerate segment.

ClosePath[p] appends the first endpoint of a path p  to p if p is open.

OpenPath[p] removes the last point of a closed path until it becomes different from the first one. Opening degenerate paths may lead to non-paths, i.e. lists with less than three points.

Test

Here p defines a nondegenerate nonnumeric path

p = {{x, y}, {2, 3}, {0, 5}}; {PathQ[p], NondegeneratePathQ[p], 
      NPathQ[p]}

Open paths can be closed with ClosePath and closed paths can be opened with OpenPath

p = {{1, 2}, {2, 3}, {0, 5}}; {NPathQ[p], OpenPathQ[p], 
      NondegeneratePathQ[p]}

Take the open path p, close  and then open it

ClosePath[p]
OpenPath[%]

Opening degenerate paths may lead to non-paths, i.e. lists with less than three points.

OpenPath[{{1, 2}, {1, 1}, {1, 2. }}]

Unit vectors, angles, areas

Norm and unit vector

The unit vector (w.r.to th euclidean norm) associated to a nonzero point p is p/(Sqrt[(Subscript[p, 1]^2) + (Subscript[p, 2]^2)]).

VectorNorm[p] yields the euclidean norm of point p

UnitVector[p] yields the unit vector defined by a nonzero point  p

Test

Clear[a, b]; 
    VectorNorm /@ {{a, b}, {4, -3}}
UnitVector /@ {{a, b}, {4, -3}}

Note that the Mathematica built-in function Norm handles complex vectors.

Wedge product

The wedge product p⋀q of two points {(Subscript[p, x]), (Subscript[p, y])} and {(Subscript[q, x]), (Subscript[q, y])} is (Subscript[p, x]) (Subscript[q, y]) - (Subscript[p, y]) (Subscript[q, x]). It equals the signed double area of the triangle  {{0, 0}, p, q}.

WedgeProduct[p, q] yields the wedge product of two points p and q.

Test wedge product

Fig.  shows the values of the wedge products {1, 1}⋀p as p varies in the following list

s = {{1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 
          0}, {-1, -1}, {0, -1}, {1, -1}}; WedgeProduct[{1, 
            1}, #] & /@ s

Block[{ray = Line[{{0, 0}, #}] &, 
        pts = {{1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 
              0}, {-1, -1}, {0, -1}, {1, -1}}}, 
      Show[Graphics[{ray /@ pts, 
            MapThread[
              Text[#1, 1.1  #2] &, {WedgeProduct[{1, 1}, #] & /@ pts, 
                pts}], RGBColor[1, 0, 0], ray[{1, 1}]}]]];

Figure

Signed (double) area

The signed double area of a polygon with vertices {(Subscript[p, 1]), …, (Subscript[p, n])} (n ≥ 3) is given by

(Subscript[p, 1])⋀(Subscript[p, 2]) + 
      (Subscript[p, 2])⋀(Subscript[p, 3]) + ⋯ + 
      (Subscript[p, n - 1])⋀(Subscript[p, n]).

DoubleArea[p] computes the signed double area of a list of three or more points.

DoubleArea3[{(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3])}] efficiently computes the signed double area of a triangle.

Test DoubleArea

DoubleArea takes a list of three or more points and yields the double area of the polygon with those vertices.

{DoubleArea[#], DoubleArea[Reverse@#]} &[{{2, 1}, {2, 2}, {1, 
          2}, {0, 0}}]

Compare the performances of DoubleArea and DoubleArea3 on 50000 random triangles

Module[{d, t1, t2, r1, r2}, 
      d = Table[{{Random[], Random[]}, {Random[], Random[]}, {Random[], 
              Random[]}}, {50000}]; {t1, r1} = 
        Timing[DoubleArea /@ d]; {t2, r2} = Timing[DoubleArea3 /@ d]; {t1, 
        t2, MatchQ[Chop[r1 - r2], {(0 | 0. ) .. }]}]

Turns

An ordered triple of points uvw with u, w ≠ v  has a left turn at v if (v - u)⋀(w - v) > 0, a right turn if (v - u)⋀(w - v) < 0, and no turn if u, v, w are collinear.

If u = v or w = v the path uvw is degenerate and NoTurnQ[{u, v, w}] yields True.

LeftTurnQ[{u, v, w}] yields True if moving from u to w a left turn occurs at v.

RightTurnQ[{u, v, w}] yields True  if moving from u to w a right turn occurs at v..

NoTurnQ[{u, v, w}] yields True if  u, v, w are collinear.

Test turns

Going from the first to the last point of p a left turn occurs at the central point:

p = {{1, 1}, {-1, 1}, {-1, -1}}; LeftTurnQ[p]

Reversing the order a left turn becomes a right one:

RightTurnQ[Reverse@p]

A no turn occurs with collinear points:

paths with nondegenerate segments

{NoTurnQ[{{1, 1}, {0, 0}, {-1. , -1}}], 
      NoTurnQ[{{1, 1}, {0, 0}, {2, 2}}], 
      NoTurnQ[{{1, 1}, {0, 0}, {1, 1}}]}

path with degenerate segments

{NoTurnQ[{{1, 1}, {0, 0}, {0, 0}}], NoTurnQ[{{0, 0}, {0, 0}, {1, 1}}], 
      NoTurnQ[{{1, 1}, {1, 1}, {1, 1}}]}

Angles

Angles are measured in radians,  counterclockwise from the positive horizontal axis. The function Angle can be called with one, two, or three arguments.

If (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]) are points with numeric components

Angle[(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3])] computes the angle formed by an ordered triple of points s.t. (Subscript[p, 1]), (Subscript[p, 3]) ≠ (Subscript[p, 2]).

Angle[(Subscript[p, 1]), (Subscript[p, 2])] with (Subscript[p, 1]), (Subscript[p, 2]) ≠ {0, 0} yields Angle[(Subscript[p, 1]), {0, 0}, (Subscript[p, 2])].

Angle[(Subscript[p, 1])] with (Subscript[p, 1]) ≠ {0, 0} yields Angle[{1, 0}, {0, 0}, (Subscript[p, 1])], i.e. the polar angle of (Subscript[p, 1]).

One argument (polar angle):

{Angle[{(1/2), ((Sqrt[3])/2)}], Angle[{0, 2}]}

Two arguments:

{Angle[{1, 1}, {1, 0}], Angle[{1, 1}, {0, 1}]}

Three arguments:

Angle[{(1/2), ((Sqrt[3])/2)}, {(1/2), 0}, {1, 0}]

Polar Ordering

Description

Given a list of points S = {(Subscript[p, 1]), …, (Subscript[p, n])} with numeric coordinates, a pole o, and a direction d, every point other than o has a well defined angle in the system of polar coordinates determined by o and d.

If two points p, p' have the same polar angle we say that p precedes p' if the norm of p is lower than that of p'.

NPolarOrdering[list, opts] gives the positions in list at which each element of NPolarSort[list, opts] appears.  opts is a sequence of options.

NPolarSort[{(Subscript[p, 1]), …, (Subscript[p, n])},  opts] sorts a list of numeric points by increasing polar angle using options opts.

The following options can be given to NPolarOrdering and NPolarSort

Origin->{x,y} orders w.r. to the pole {x, y} (numeric point defaulting to  {0, 0}). 
No point of list can coincide with the pole.

Direction->{x',y'} orders using the polar axis leading from the pole o to o + {x', y'}. {x', y'} defaults to  {1, 0} (horizontal axis).

Test

Let's consider a list s = {(Subscript[p, 1]), …, (Subscript[p, 10])} of ten random points (Fig.  - left) and the result of PolarSort[s] (right)

s = (SeedRandom[101]; 
        Table[{Random[Real, {-1, 1}], 
            Random[Real, {-1, 1}]}, {10}]);

Show[
      GraphicsArray[{Graphics[{Line[{{0, 0}, #}] & /@ #, 
                  MapIndexed[
                    Text[#2[[1]], #] &, #]}] &[s], Graphics[{Line[{{0, 0}, #}] & /@ #, 
                  MapIndexed[
                    Text[#2[[1]], #] &, #]}] &[NPolarSort[s]]}]];

Figure

The point with smallest polar angle is (Subscript[p, 2]), followed by (Subscript[p, 4]),…, and the point with largest polar angle is (Subscript[p, 7]). The sequence {2, 4, …, 7} is returned by

NPolarOrdering[s]

Changing the pole from {0, 0} to {-2, 2} one gets ( Fig. )

NPolarOrdering[s, Origin -> {-2, 2}]

With[{o = {-2, 2}}, 
      Show[GraphicsArray[{Graphics[{Line[{o, #}] & /@ #, 
                    MapIndexed[
                      Text[#2[[1]], #] &, #]} &[s]], 
            Graphics[{Line[{o, #}] & /@ #, 
                    MapIndexed[
                      Text[#2[[1]], #] &, #]} &[NPolarSort[s, Origin -> o]]]}]]];

Figure

In this example angles are measured from the negative vertical axis (Fig. ):

NPolarOrdering[s, Origin -> {-2, 2}, 
      Direction -> {0, -1}]

With[{o = {-2, 2}, d = {0, -1}}, 
      Show[GraphicsArray[{Graphics[{Line[{o, #}] & /@ #, 
                    MapIndexed[
                      Text[#2[[1]], #] &, #]}] &[s], Graphics[{Line[{o, #}] & /@ #, 
                    MapIndexed[
                      Text[#2[[1]], #] &, #]}] &[NPolarSort[s, Origin -> o, 
                Direction -> d]]}]]; ]

Figure

Lexicographic Ordering

Description

Two linearly independent vectors (directions) {(Subscript[d, 1]), (Subscript[d, 2])} totally order the plane since for every point p there exists a unique pair of scalars (Subscript[ξ, p]), (Subscript[η, p]) s.t.

\(TraditionalForm\`p = \(d\_1\) ξ\_p + \(\(d\_2\) \(η\_p\)\(\ \)\)\)

and one can define p ≺ q iff (Subscript[ξ, p]) < (Subscript[ξ, q]) or \(TraditionalForm\`ξ\_p = ξ\_q\ and\ \ η\_p < η\_q\) (generalized lexicographic ordering). The usual lexicographic ordering is a particular case with (Subscript[d, 1]) = {1, 0} and (Subscript[d, 2]) = {0, 1}.

NLexOrdering[{(Subscript[p, 1]), …, (Subscript[p, n])}] takes a set of numeric points and returns a list {(Subscript[i, 1]), …, (Subscript[i, n])} of indices s.t. if i precedes j then (Subscript[p, i]) lexicographically precedes (Subscript[p, j]).

NLexSort[{(Subscript[p, 1]), …, (Subscript[p, n])}] acts as NLexOrdering but it returns the sorted list of points (instead of indices).

NMinLexOrdering[{(Subscript[p, 1]), …, (Subscript[p, n])}] yields the position of a point at which the lexicographic minimum is attained.

NMinLexSort[{(Subscript[p, 1]), …, (Subscript[p, n])}] yields a point at which the lexicographic minimum is attained.

NMaxLexOrdering[{(Subscript[p, 1]), …, (Subscript[p, n])}] yields the position of a point at which the lexicographic maximum is attained.

NMaxLexSort[{(Subscript[p, 1]), …, (Subscript[p, n])}] yields a point at which the lexicographic maximum is attained.

Options: a direction or a pair of independent directions can be supplied to these functions,

Direction->{{x,y},{x',y'}} orders using {x, y} as (Subscript[d, 1]) and {x', y'} as (Subscript[d, 2]).

Direction->{x,y} orders using {x, y} as (Subscript[d, 1]) and the orthogonal direction {-y, x} as (Subscript[d, 2]).
The default is Direction->{{1,0},{0,1}} .

Test

Consider the vertices of a square and their lexicographic ordering

s = {{1, 1}, {-1, 1}, -{1, 1}, {1, -1}}; {NLexOrdering[s], 
      NMinLexOrdering[s], NMaxLexOrdering[s]}

The positions of the points inside s before and (bracketed) after the lexicographic ordering with default directions {{1, 0}, {0, 1}} are shown in Fig. .

\(With[{ii = NLexOrdering[s]}, 
      Show[Graphics[{Line[ClosePath[s]], 
            MapIndexed[
              Text[ToString[#1] <> "[" <> ToString[First@#2] <> "]", 
                  s[[#1]], \(-s[[#1]]\)] &, ii]}], 
        PlotRange -> All]]; \)

Figure

Comparing y-components first, then  x-components, if necessary (Fig. ):

{NLexOrdering[s, Direction -> {{0, 1}, {1, 0}}], 
      NMinLexOrdering[s, Direction -> {{0, 1}, {1, 0}}], 
      NMaxLexOrdering[s, Direction -> {{0, 1}, {1, 0}}]}

The same with sort:

{NLexSort[s, Direction -> {{0, 1}, {1, 0}}], 
      NMinLexSort[s, Direction -> {{0, 1}, {1, 0}}], 
      NMaxLexSort[s, Direction -> {{0, 1}, {1, 0}}]}

\(With[{ii = NLexOrdering[s, Direction -> {{0, 1}, {1, 0}}]}, 
      Show[Graphics[{Line[ClosePath[s]], 
            MapIndexed[
              Text[ToString[#1] <> "[" <> ToString[First@#2] <> "]", 
                  s[[#1]], \(-s[[#1]]\)] &, ii]}], 
        PlotRange -> All]]; \)

Figure

If the first direction is {1, 1} and the second {1, -1} one gets

NLexOrdering[s, Direction -> {{1, 1}, {1, -1}}]

Affine Transformations

Description

A plane affine transformations maps a point p as follows

p↦(({{(Subscript[t, 11]), (Subscript[t, 12])}, {(Subscript[t, 21]), (Subscript[t, 22])}}))p+(({{(Subscript[t, 31])}, {(Subscript[t, 32])}})).

In this section pts will denote a single point or a list of points.

Shifts

Shift[pts, q] shifts pts by the amount specified by a given point  q = {(Subscript[q, x]), (Subscript[q, y])};

Shift[pts, Offset[ q]] offsets pts horizontally and vertically by (Subscript[q, x]) and (Subscript[q, y])printer points, respectively.

Examples

Fig.  shows three black points and their shifts by {0.2,-0.5} (gray).

s = {{0, - .3}, {1, 0}, {0,  .3}}; 
    s1 = Shift[s, { .2, - .5}]

Show[
      Graphics[{AbsolutePointSize[5], Point /@ s, GrayLevel[ .75], 
          Point /@ s1}], AspectRatio -> Automatic, 
      Frame -> True];

Figure

The thin triangle in Fig.  is obtained from the thick one (list s) with a horizontal offset of  20 (printer's points)

Shift[s = {{0, - .3}, {1, 0}, {0,  .3}}, Offset@{20, 0}]

Show[
      Graphics[{{AbsoluteThickness[3], Line[ClosePath@s]}, 
          Line[Shift[ClosePath@s, Offset@{20, 0}]]}], 
      AspectRatio -> Automatic, PlotRange -> All];

Figure

Scaling and slanting

Scale[pts, {(Subscript[f, x]), (Subscript[f, y])}] scales pts horizontally by a factor (Subscript[f, x]) and vertically by (Subscript[f, y]);

ComplexScale[pts, {(Subscript[f, x]), (Subscript[f, y])}] implements {(Subscript[p, x]), 
          (Subscript[p, y])}  ↦  {(Subscript[p, x]) (Subscript[f, x]) - (Subscript[p, y]) (Subscript[f, y]), (Subscript[p, x]) 
              (Subscript[f, y]) - (Subscript[p, y]) (Subscript[f, x])} ;

Slant[pts, s] implements {(Subscript[p, x]), (Subscript[p, y])} ↦  {(Subscript[p, x]) + s (Subscript[p, y]), 
          (Subscript[p, y])}.

ComplexScale[p, f] gives the real and imaginary parts of the product of its arguments interpreted as complex numbers ((Subscript[p, x]) + i (Subscript[p, y])) ((Subscript[f, x]) + i (Subscript[f, y])). The same result is attained scaling both components by (Sqrt[((Subscript[f, x])^2) + ((Subscript[f, y])^2)]) and rotating by Angle[{(Subscript[f, x]), (Subscript[f, y])}].

Examples

Fig.  shows the effect of scaling both coordinates (same factor), only the x ,and  only the y coordinate:

s = {{0, - .3}, {1, 0}, {0,  .3}}; 
    Scale[s, 1.25 {1, 1}]
Scale[s, {1.25, 1}]
Scale[s, {1, 1.25}]

With[{cp = ClosePath@s}, 
      Show[GraphicsArray[{Graphics[{{Text["O", {0, 0}, {1, 0}], 
                  AbsolutePointSize[5], Point[{0, 0}], AbsoluteThickness[2], 
                  Line[cp]}, Line[Scale[cp, {1, 1} 1.25]]}, 
              AspectRatio -> Automatic], 
            Graphics[{{Text["O", {0, 0}, {1, 0}], AbsolutePointSize[5], 
                  Point[{0, 0}], AbsoluteThickness[2], Line[cp]}, 
                Line[Scale[cp, {1.25, 1}]]}, AspectRatio -> Automatic], 
            Graphics[{{Text["O", {0, 0}, {1, 0}], AbsolutePointSize[5], 
                  Point[{0, 0}], AbsoluteThickness[2], Line[s]}, 
                Line[Scale[cp, {1, 1.25}]]}, 
              AspectRatio -> Automatic]}]]];

Figure :scaling by {1.25,1.25}, by {1.25,1}, and  by {1,1.25}.

Some examples of complex scaling are shown in Fig. .

s = {{0, - .3}, {1, 
            0}, {0,  .3}};
ComplexScale[
      s, {1, 1}]
ComplexScale[s,UnitVector[{1, 1}]]
ComplexScale[s,{Sqrt[3]/2, -1/2}]

With[{cp = ClosePath[s]},Show[GraphicsArray[{Graphics[{{Text["O", {0, 0}, {1, 0}], 
                        AbsolutePointSize[5], Point[{0, 0}], 
                        AbsoluteThickness[2], Line[cp]}, 
                      Line[ComplexScale[cp, {1, 1}]]}, 
                    AspectRatio -> Automatic],Graphics[{{Text["O", {0, 0}, {1, 0}], 
                        AbsolutePointSize[5], Point[{0, 0}], 
                        AbsoluteThickness[2], Line[cp]}, 
                      Line[ComplexScale[cp, UnitVector@{1, 1}]]}, 
                    AspectRatio -> Automatic],Graphics[{{Text["O", {0, 0}, {1, 0}], 
                            AbsolutePointSize[5], Point[{0, 0}], 
                            AbsoluteThickness[2], Line[cp]},Line[ComplexScale[cp,{Sqrt[3]/2, -1/2}] ]},AspectRatio -> Automatic]}]]];

Figure :  complex scaling

Fig.  shows a slanted unit square

Slant[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}, ϵ]

With[{s = ClosePath@{{0, 0}, {1, 0}, {1, 1}, {0, 1}}}, 
        Show[Graphics[{Line[s], AbsoluteThickness [2], 
              Line[Slant[s,  .25]]}], Axes -> True, 
          AspectRatio -> Automatic]];

Figure

Rotations and reflections

Rotate[pts, θ, {(Subscript[c, x]), (Subscript[c, y])}] counterclockwise rotates by θ radians around {(Subscript[c, x]), (Subscript[c, y])}. The center of rotation defaults to {0, 0}.

Reflect[pts, {a,b}] reflects about the line determined by two distinct points a, b.

The coordinates of  p' (reflection of p about the line through points a and b) are given by

p'=a+(({{(Subscript[u, x]^2) - (Subscript[u, y]^2), 2  (Subscript[u, x]) (Subscript[u, y])}, {2  (Subscript[u, x]) (Subscript[u, y]), -(Subscript[u, x]^2) + (Subscript[u, y]^2)}}))(p - a)

where {(Subscript[u, x]), (Subscript[u, y])} = 
            UnitVector[{(Subscript[b, x]) - (Subscript[a, x]), (Subscript[b, y]) - (Subscript[a, y])}] .

Examples

Rotate a point p by 90° and check orthogonality

Rotate[p = {5, 3}, 90  °]
p . %

Rotate the unit square by -30° around the origin and around {1/2,1/2} (see Fig. )

u = {{0, 0}, {1, 0}, {1, 1}, {0, 1}}; 
    Rotate[u, -30 °]
Rotate[u, -30 °, {1/2, 1/2}]

With[{c = ClosePath@u}, 
      Show[GraphicsArray[{Graphics[{Line[c], AbsoluteThickness [2], 
                Line[Rotate[c, -30 °]], AbsolutePointSize[6], 
                Point[{0, 0}]}, AspectRatio -> Automatic], 
            Graphics[{Line[c], AbsoluteThickness [2], 
                Line[Rotate[c, -30 °, { .5,  .5}]], 
                AbsolutePointSize[4], Point[{ .5,  .5}]}, 
              AspectRatio -> Automatic]}]]];

Figure

A reflection about the axis determined by {0,1} and {1,2}

Reflect[{a, b}, {{0, 1}, {1, 2}}] // Expand

Fig.  shows an ∃ shape (thin) and its reflection about the line determined by points {1, 1.75} and {2, 1.5} (thick).

Module[{s = {{0, 0}, {2, 0}, {2, 2}, {1, 2}, {2, 2}, {2, 4}, {0, 4}}, 
        a = {1, 1.75}, b = {2, 1.5}, s1}, s1 = Reflect[s, {a, b}]; 
      Show[Graphics[{{RGBColor[1, 0, 0], Dashing[{0.02, 0.02, 0.01, 0.02}], 
              Line[{a - 2 (b - a), 
                  a + 2.2 (b - a)}]}, {AbsoluteThickness[2], 
              Line[s]}, {AbsoluteThickness[4], Line[s1]}, 
            Dashing[{0.02, 0.02}], MapThread[Line[{#1, #2}] &, {s, s1}]}], 
        AspectRatio -> Automatic, Frame -> True]];

Figure :  A reflected ∃ shape.

General transformation

Transform[pts, {{(Subscript[t, 11]), (Subscript[t, 12]), (Subscript[t, 13])},{(Subscript[t, 21]), (Subscript[t, 22]), (Subscript[t, 23])}}] implements  p↦(({{(Subscript[t, 11]), (Subscript[t, 12])}, {(Subscript[t, 21]), (Subscript[t, 22])}}))p+(({{(Subscript[t, 31])}, {(Subscript[t, 32])}}))

Examples

Apply the transformation T to the vertices u of the unit square (Fig. )

u = {{0, 0}, {1, 0}, {1, 1}, {0, 1}}; 
      T = {{2, 3, -1}, {4, 1, -2}}; 
u1 = 
        Transform[u, T]

Show[
      Graphics[{Line[ClosePath@u], AbsoluteThickness [2], 
          Line[ClosePath@u1]}], Axes -> True, 
      AspectRatio -> Automatic];

Figure :  The unit square transformed by  T={{2,3,-1}, {4,1,-2}}

Lines

Description

The line through two distinct points (Subscript[p, 1]), (Subscript[p, 2]) is the set of points p s.t.

(p - (Subscript[p, 1]))⋀((Subscript[p, 2]) - (Subscript[p, 1])) = 0,i.e.p⋀((Subscript[p, 2]) - (Subscript[p, 1])) = (Subscript[p, 1])⋀(Subscript[p, 2]).

With p = {x, y} one gets the cartesian equation

((Subscript[p, 2  y])  - (Subscript[p, 1  y])) 
          x - ((Subscript[p, 2  x]) - (Subscript[p, 1  x])) 
          y - ((Subscript[p, 1  x]) (Subscript[p, 2  y]) - (Subscript[p, 1  y]) 
              (Subscript[p, 2  x])) = 0.

The parametric equations are the components of p(t) = (Subscript[p, 1]) + t((Subscript[p, 2]) - (Subscript[p, 1])) (p(t) is the so called t-point of (Subscript[p, 1]) and (Subscript[p, 2])) i.e.

({{x(t) = (Subscript[p, 1  x]) + t((Subscript[p, 2  x]) - (Subscript[p, 1  x]))}, {y(t) = 
              (Subscript[p, 1  y]) + t((Subscript[p, 2  y]) - (Subscript[p, 1  y])).}})

CartesianLine[{(Subscript[p, 1]), (Subscript[p, 2])}]  yields the cartesian equation of the line through two  points (Subscript[p, 1]) ≠ (Subscript[p, 2]) as a pure function. 
CartesianLine[{(Subscript[p, 1]),(Subscript[p, 2])}][x,y] yields the traditional equation with x and y.

CartesianAxis[{(Subscript[p, 1]), (Subscript[p, 2])}]  yields the cartesian equation of the axis (orthogonal bisector) of a nondegenerate segment with endpoints (Subscript[p, 1]), (Subscript[p, 2]) as a pure function. 
CartesianAxis[{(Subscript[p, 1]),(Subscript[p, 2])}][x,y] yields the traditional equation with x and y.

ParametricLine[{(Subscript[p, 1]),(Subscript[p, 2])}] yields the parametric equation of the line through (Subscript[p, 1]) and (Subscript[p, 2]) as a pure function.
ParametricLine[{(Subscript[p, 1]),(Subscript[p, 2])}][t] yields the traditional equation with parameter t.

ParametricAxis[{(Subscript[p, 1]),(Subscript[p, 2])}] yields the parametric equation of the axis of a nondegenerate segment (Subscript[p, 1]), (Subscript[p, 2]) as a pure function.
ParametricAxis[{(Subscript[p, 1]),(Subscript[p, 2])}][t] yields  the traditional equation with parameter t.
When the parameter is zero one gets the midpoint of (Subscript[p, 1]) (Subscript[p, 2]) .

tPoint[t,{(Subscript[p, 1]),(Subscript[p, 2])}] yields the point (1 - t) (Subscript[p, 1]) + t (Subscript[p, 2]).

MidPoint[{(Subscript[p, 1]),(Subscript[p, 2])}] is the same as tPoint[1/2,{(Subscript[p, 1]),(Subscript[p, 2])}].

Test

The cartesian equation of a line defined by a pair of points as a pure function and as a function of x, y

CartesianLine[{{1, 2}, {0, 1}}]
CartesianLine[{{1, 2}, {0, 1}}][x, y]

The same for a pair of points with nonnumeric components

CartesianLine[{{px, py}, {qx, qy}}]
CartesianLine[{{px, py}, {qx, qy}}][x, y]

Parametric equations of a line

ParametricLine[{{1, 2}, {0, 1}}]
ParametricLine[{{1, 2}, {0, 1}}][t]
ParametricLine[{{px, py}, {qx, qy}}]
ParametricLine[{{px, py}, {qx, qy}}][t]

Segment axis:

CartesianAxis[{{1, 2}, {0, 1}}][x, y]
CartesianAxis[{{1, 2}, {1, 3}}][x, y]
ParametricAxis[{{1, 2}, {0, 1}}][t]
ParametricAxis[{{1, 2}, {1, 3}}][t]

tPoint[t,{p, q}] yields the same result as ParametricLine[{p, q}][t], and MidPoint[{p, q}] is equivalent to tPoint[1/2,{p, q}]

ParametricLine[{{px, py}, {qx, 
            qy}}][α]
tPoint[α, {{px, py}, {qx, qy}}]
Expand[% - %%] == {0, 0}

MidPoint[{{px, py}, {qx, qy}}] - tPoint[1/2, {{px, py}, {qx, qy}}] // 
      Simplify

Circular Arcs

The circumcircle of three points ((Subscript[p, 1]),({{(Subscript[p, 2]),, (Subscript[p, 3])}})) is returned as a pair {C, r} where C is the circumcenter and r the circumradius. A circle can be alternatively specified by two distinct points (Subscript[p, 1]), (Subscript[p, 2]) and an angle ϕ ≠ 0 formed by the directed line (Subscript[p, 1]) (Subscript[p, 2]) with the tangent to the circle at (Subscript[p, 1]).

CircumCircle[{(Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3])}] yields the center C and the radius r of the circumcircle of three points as {{(Subscript[C, x]), (Subscript[C, y])}, r}.

CircumCircle[(Subscript[p, 1]), ϕ, (Subscript[p, 2])] yields the center C and the radius r of the circle through (Subscript[p, 1]), (Subscript[p, 2])whose tangent at (Subscript[p, 1])forms an angle ϕ with the directed line from (Subscript[p, 1])to (Subscript[p, 2]).

CircumCircle[(Subscript[p, 1]), Absolute[ϕ], (Subscript[p, 2])] as above with the angle measured from to the horizontal axis.

Degenerate circles are returned as {{∞, ∞}, ∞}.

Test

Circumcircle of three noncollinear points and their distances from the computed circumcenter

{c, r} = 
      CircumCircle[s = {{0, 0}, {1, 0}, {1, 1}}]
VectorNorm[# - c] & /@ s

Three collinear points

CircumCircle[{{0, 1}, {2, 1}, {1, 1}}]

Circle through {0,0}, {0,1}, and with a 45° tangent measured from the directed line determined by (0, 0), (0, 1) or w.r.to the horizontal line (see Fig.  )

CircumCircle[{0, 0}, 45  °, {0, 1}]
CircumCircle[{0, 0}, Absolute[45  °], {0, 1}]

\(Show[
      GraphicsArray[{Graphics[{\(Circle[#[[1]], #[[2]]] &\)@
                CircumCircle[{0, 0}, 45  °, {0, 1}], 
              AbsolutePointSize[5], Point /@ {{0, 0}, {0, 1}}, 
              Text[\*""(Subscript[p, 1])"", {0, 0}, {1, 1}], 
              Text[\*""(Subscript[p, 2])"", {0, 1}, {1, \(-1\)}], 
              AbsoluteDashing[{2, 2}], 
              Line[{{0, 1}, {0, 0},  .5 {\(-1\), 1}}]}, 
            AspectRatio -> Automatic], 
          Graphics[{\(Circle[#[[1]], #[[2]]] &\)@
                CircumCircle[{0, 0}, Absolute[45  °], {0, 1}], 
              AbsolutePointSize[5], Point /@ {{0, 0}, {0, 1}}, 
              Text[\*""(Subscript[p, 1])"", {0, 0}, {1, \(-1\)}], 
              Text[\*""(Subscript[p, 2])"", {0, 1}, {1, 1}], 
              AbsoluteDashing[{2, 2}], 
              Line[{ .5 {01, 0}, {0, 0},  .5 {1, 1}}]}, 
            AspectRatio -> Automatic]}]]; \)

Figure :  CircumCircle[(Subscript[p, 1]), 45°, (Subscript[p, 2])] and CircumCircle[(Subscript[p, 1]), Absolute[ 45°], (Subscript[p, 2])]

Point Location

Description

We consider the location of a query point q w.r.to various geometric objects: line, halfline, segment, polygon, and circle.

The outcome of a point location w.r.to a line ℓ defined by distinct points (Subscript[p, 1]), (Subscript[p, 2]) and directed from (Subscript[p, 1])to (Subscript[p, 2]) is LeftPoint if q is to the left of ℓ, RightPoint if it is to the right of ℓ, and InPoint if q belongs to the relative interior of ℓ.

The outcome of a point location w.r.to halflines, segments, polygons, and circles is  InPoint if q belongs to the relative interior of the object, BorderPoint if q belongs to the boundary, and  OutPoint otherwise.

If obj denotes an object and q a query point then NPointLocation[q, obj] and NPointLocation[obj, q] give the same result.

In the sequel q,  (Subscript[p, 1]), (Subscript[p, 2]) will denote points with numeric components.

Lines:

NPointLocation[q,Line[ {(Subscript[p, 1]), (Subscript[p, 2])}]] takes a point q, a pair of distinct points (Subscript[p, 1]), (Subscript[p, 2]), and returns 
LeftPoint if q is to the left  of the line directed from (Subscript[p, 1])to (Subscript[p, 2]),
RightPoint if q is to the right of  (Subscript[p, 1]) (Subscript[p, 2]), 
InPoint if q is on the line (Subscript[p, 1]) (Subscript[p, 2]) .

Halflines and segments:

NPointLocation[q,obj[ {(Subscript[p, 1]), (Subscript[p, 2])}]] with obj equal to HalfLine or Segment, takes a point q, a pair of points (Subscript[p, 1]), (Subscript[p, 2]) defining a segment or a halfline emanating from the first point, and returns 
OutPoint if q is external to obj, 
InPoint if it belongs to the relative interior of obj, 
BorderPoint if q belongs to the relative boundary of obj.

Convex and simple polygons:

NPointLocation[q, Convex[𝒫]] takes a point q, a  convex polygon 𝒫,  and returns InPoint,  BorderPoint, or OutPoint if q is in the interior, on the boundary, or outside 𝒫. Both q and the vertices of 𝒫 must have numerical components.

NPointLocation[q, Simple[𝒫]] takes a point q and a  simple polygon 𝒫 and returns InPoint,  BorderPoint or OutPoint if q is in the interior, on the boundary, or outside 𝒫. Both q and the vertices of 𝒫 must have numeric components.

Circles:

PointLocation[q,Circle[{(Subscript[p, 0]), (Subscript[p, 1]), (Subscript[p, 2])}]] returns InPoint,  BorderPoint or OutPoint if the query point q lies in the interior, on the boundary, or in the exterior of circle C.

Test

Line

The following example illustrates the point location with respect to the line through {17,10} and {19,11} using a list qlist of query points (Fig. )

r = {{17, 10}, {19, 11}}; 
    qlist = {{10, 7}, {18, 8}, {21, 12}}; NPointLocation[#, Line[r]] & /@ 
      qlist

ParametricPlot[
      Evaluate[ParametricLine[{{1, 2}, {3, 3}}][t]], {t, 4, 11}, 
      Epilog -> {PointSize[0.025], Point /@ qlist, GrayLevel[0.5], 
          Point /@ r}];

Figure :  Point location w.r. to a line

HalfLine

Point location with respect to the halfline from  {17,10}through {19,11}: the query points are {10,7},{18, 8} (not on the halfline),{21,12} (in the relative interior) and {17,10} (on the relative boundary)

r = {{17, 10}, {19, 11}}; 
    qlist = {{10, 7}, {18, 8}, {21, 12}, {17, 10. }}; NPointLocation[#, 
          HalfLine[r]] & /@ qlist

Considering the halfline from {19,11} to {17,10} no point but {17,10}belongs to it:

r = Reverse@{{17, 10}, {19, 11}}; 
    qlist = {{10, 7}, {18, 8}, {21, 12}, {17, 10}}; NPointLocation[#, 
          HalfLine[r]] & /@ qlist

Segment

Point Location w.r.to the vertical segment seg defined by {1,2} and {1, 2.5} and a set of query points (qlist):  the first two points are external, the third coincides with an endpoint, and the fourth lies in the relative interior (Fig. )

seg = {{1, 2}, {1, 4}}; 
    qlist = {{0.5, 0.5}, {1, 0.5}, {1, 2}, {1, 2.5}}; NPointLocation[#, 
          Segment[seg]] & /@ qlist

Show[
      Graphics[{Line[seg], PointSize[0.025], 
          MapIndexed[{Point[#1], Text[First[#2], #1, {-1.5, 1.5}]} &, 
            qlist]}], PlotRange -> All];

Figure :  Point location w.r. to a vertical segment.

Non vertical segment seg (Fig. ):

seg = {{-1, -1}, {4, 4}}; 
    qlist = {{0.5, 0.5}, {1, 0.5}, {1, 2}, {1, 2.5}}; NPointLocation[#, 
          Segment[seg]] & /@ qlist

Show[
      Graphics[{Line[seg], PointSize[0.025], 
          MapIndexed[{Point[#1], Text[First[#2], #1, {-1.5, 1.5}]} &, 
            qlist]}], PlotRange -> All];

Figure :  Point location w.r. to a nonvertical segment.

Now the degenerate segments Segment[{{1,0.},{1.,0}}], {1.,0} and {} are considered

NPointLocation[#, Segment[{{1, 0. }, {1. , 0}}]] & /@ {{0, 0}, {1, 
          0}} 
NPointLocation[#, {1. , 0}] & /@ {{0, 0}, {1, 
          0}}
NPointLocation[#, {}] & /@ {{0, 0}, {1, 0}}

Convex polygon

A list of points s defining a convex polygon and a list qs of query points (Fig. )

s = {{0, 2}, {1, 1}, {1, 0}, {0, 0}}; 
    qs = {{1, 2}, { .5,  .8}, { .5, 1.5}}; 
NPointLocation[#, Convex[s]] & /@ qs

Show[
      Graphics[{Line[ClosePath@s], AbsolutePointSize[5], Point /@ qs, 
          MapIndexed[Text[First@#2, #, {0, 2}] &, qs]}]];

Figure

Simple polygon

A simple polygon s and a list qs of query points (Fig. )

s = {{0, 0}, {1, -1}, {2, 0}, {2.5, -1.1}, {3, 0}, {4, 0}, {4, 
          1}, {2.5, 2.5}, {0.5, 1}}; 
    qs = Table[{i, 0}, {i, -1/2, 4 + 1/2, 1/2}]; NPointLocation[#, 
          Simple[s]] & /@ qs

Show[
        Graphics[{Line[ClosePath[s]], PointSize[0.02], Point /@ qs, 
            MapIndexed[Text[First[#2], #1, {0, -1.5}] &, 
              qs]}]];

Figure

Circle

A circle and a sequence of query points qs (Fig. )

NPointLocation[#, 
          Circle[{{-1, 0}, {0, -1}, {1, 1}/(Sqrt[2])}]] & /@ (qs = 
          Table[{t, t/2 - 1}, {t, -1, 2, 1/3}])

Show[
        Graphics[{Circle @@ 
              CircumCircle[{{-1, 0}, {0, -1}, {1, 1}/(Sqrt[2])}], 
            PointSize[0.02], Point /@ qs, 
            MapIndexed[Text[First[#2], #1, {0, -1.5}] &, qs]}], 
        AspectRatio -> Automatic];

Figure

Intersections

Description

Intersection tests and intersection  computations are implemented for geometric objects defined by two pairs of points  (Subscript[p, 1]), (Subscript[p, 2]) and (Subscript[q, 1]), (Subscript[q, 2]). The result of an intersection test is

False if there is no intersection at all, otherwise it is Strong or Weak.

Weak intersections are intersections that can be destroyed by arbitrarily small suitably chosen perturbations;

Strong intersections can be destroyed only by means of sufficiently large perturbations.

For instance, the intersection between Segment[{{-1,0}, {1,0}}] and Segment[{{0,-1}, {0,1}}] is strong; that of Segment[{{-1,0}, {1,0}}] and Segment[{{0,-1}, {0,0}}] is weak.

NIntersection[obj1,obj2] computes the intersection of two objects defined by points with numeric coordinates.

NIntersection[obj1,obj2] is equivalent to NIntersection[obj2,obj1].

The following objects can be specified:

Line[{(Subscript[p, 1]),(Subscript[p, 2])}], HalfLine[{(Subscript[p, 1]),(Subscript[p, 2])}], Segment[{(Subscript[p, 1]),(Subscript[p, 2])}], Segment[{p, p}], p, {}, where p, (Subscript[p, 1]), (Subscript[p, 2]) are points with numeric components and (Subscript[p, 1]) ≠ (Subscript[p, 2]).
Segment[{p, p}], p, and {} denote segments degenerating to a point p or empty.

NIntersection[obj1,obj2] yields a pair {inttype, intset}:  inttype is the type of intersection (False, Weak, or Strong), and intset is the intersection set.

An empty intersection is returned as {False,{}},

A strong intersection at {x, y} is returned as {Strong,{x,y}}.

A weak intersection is returned as {Weak, intset}: here the intersection set intset can be a line, a halfline or a segment (degenerate or not). A weak intersection occurs when the relative interiors have an infinity of common points, or when the intersection contains a boundary point.

Test

Intersection of a line and…

The intersection of two lines may be empty, a point, or a line:

NIntersection[
          Line[{{0, 0}, {2, 1}}], #] & /@ {Line[{{0, 1}, {2, 2}}], 
        Line[{{0, 1}, {1, 0}}], Line[{{-2, -1}, {4, 2}}]}

Some examples of intersections of a line L and a list HLs of halflines (Fig ):

L = {{0, 0}, {2, 1}}; 
    HLs = {{{-1, 0}, {3, 8/3}}, {{0, 0}, {3, 3}}, {{0, -1}, {0, 
            3}}, {{2, 0}, {2, 3}}}; NIntersection[Line[L], 
          HalfLine[#]] & /@ HLs

Show[
      Graphics[{Line /@ 
            HLs, {AbsolutePointSize[5], Point[First[#]] & /@ 
              HLs}, {AbsoluteDashing[{3, 3}], 
            Line[ParametricLine[L][#] & /@ {- .5, 1.5}]}}], 
      Frame -> True];

Figure

A halfline contained in a line (weak intersection):

L = {{0, 0}, {2, 1}}; HL = ParametricLine[L] /@ {2, 3}; 
    NIntersection[Line[L], HalfLine[HL]]

Fig.  shows the intersections of a line r and a list sl of segments (only the last two  intersect  r)

r = {{1, 2}, {2, 4}}; 
sl = {{{3, 2}, {4, 2}}, {{3, 2}, {4, 8}}, {{3, 2}, {8, 
            17}}}; 
NIntersection[Line[r], Segment[#]] & /@ sl

Show[
      Graphics[{{RGBColor[1, 0, 1], 
            Line[{r[[1]], ParametricLine[r][8]}]}, Line /@ sl, 
          MapIndexed[Text[First[#2], Last[#1], {0, -1.2}] &, 
            sl]}]];

Figure

When a segment s contained in a line r  {Weak,s} is returned:

s = Segment[ParametricLine[r = {{1, 2}, {2, 4}}] /@ {0, 2}]; 
    NIntersection[Line[r], s]

A line r and some degenerate segments: {} and Segment[{{3/2,3},{3/2,3.0}}] on r

NIntersection[Line[r = {{1, 2}, {2, 4}}], #] & /@ {{}, 
        Segment[{{3/2, 3}, {3/2, 3.0}}]}

A segment degenerating to  a point {3/2,3.0} on r and one degenerating to a point outside r

NIntersection[Line[r], #] & /@ {{3/2, 3.0}, {5/2, 
          5. /2}}

Intersection of a halfline and…

A halfline h1 and a list of halflines h2L whose intersection with h1 if nonempty is a single point (Fig. )

h1 = {{-1, 1}, {1, 3/2}}; 
    h2L = {{0, 0}, {#, 1/2}} & /@ {-1, -1/2, 0, 
          1}; NIntersection[HalfLine[h1], HalfLine[#]] & /@ 
      h2L

Show[
      Graphics[{{RGBColor[0, 0, 1], 
            Line[ParametricLine[h1][#] & /@ {0, 2}], 
            AbsoluteDashing[{4, 3}], 
            Line[ParametricLine[h1][#] & /@ {2, 
                  4}]}, Line[{ParametricLine[#][
                    0], ParametricLine[#][8]}] & /@ h2L, 
          AbsolutePointSize[4], 
          Point /@ {{0, 0}, 
              h1[[1]]}}]];

Figure

Now the intersection of two halflines defined by four collinear points is considered. If the halflines have the same orientation their intersection is a halfline (Fig. ),

NIntersection[HalfLine[{{0, 0}, {1, 0}}], 
          HalfLine[#]] & /@ (t = {{{-2, 0}, {-1, 0}}, {{-1, 
                0}, {0, 0}}, {{0, 0}, {1, 0}}, {{1, 0}, {2, 
                0}}})

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], 
                Line[{{0, 0}, {5, 0}}], Dashing[{0.02, 0.03}], 
                Line[{{5, 0}, {6, 0}}]}, 
              Line[{Offset[{0, 
                      5}, #[[1]]], 
                  Offset[{0, 
                      5}, #[[2]]]}], {Dashing[{0.02, 0.03}], 
                Line[{Offset[{0, 
                        5}, #[[2]]], 
                    Offset[{0, 
                        5}, {1, 
                          0} + #[[2]]]}]}}, AspectRatio -> 1/8, Axes -> {True, False}]] &, t];

Figure

When the halflines have opposite orientations their intersection is either empty, or a point, or a nondegenerate segment (Fig. ).

NIntersection[HalfLine[{{0, 0}, {1, 0}}], 
          HalfLine[#]] & /@ (t = {{{-1, 0}, {-2, 0}}, {{0, 
                0}, {-1, 0}}, {{2, 0}, {1, 0}}})

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], 
                Line[{{0, 0}, {5, 0}}], Dashing[{0.02, 0.03}], 
                Line[{{5, 0}, {6, 0}}]}, Line[#]}, AspectRatio -> 1/8, 
            Axes -> {True, False}]] &, t];

Figure

The intersection of a halfline h and a nondegenerate segment s is either empty,  or a single point, or a segment. If the endpoints of s belong to different open halfspaces defined by h then the intersection can be empty, weak, or strong (Fig. )

NIntersection[HalfLine[{{0, 0}, {1, 0}}], 
          Segment[#]] & /@ (t = {{{-1, -1}, {-1, 
                1}}, {{0, -1}, {0, 1}}, {{1, -1}, {1, 
                1}}})

Scan[
        Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], 
                  Line[{{0, 0}, {5, 0}}], Dashing[{0.02, 0.03}], 
                  Line[{{5, 0}, {6, 0}}]}, 
                Line[Shift[#, Offset[{0, 5}]]], {Dashing[{0.02, 0.03}], 
                  Line[{Offset[{0, 
                          5}, #[[2]]], 
                      Offset[{0, 
                          5}, {1, 
                            0} + #[[2]]]}]}}, AspectRatio -> 1/8, Axes -> {True, False}]] &, 
        t];

Figure

Strong intersections cannot arise when s has an endpoint on the line determined by h (Fig.  and )

NIntersection[HalfLine[{{0, 0}, {1, 0}}], 
          Segment[#]] & /@ (t = {{{-1, 0}, {-1, 1}}, {{0, 0}, {0, 
                1}}, {{1, 0}, {1, 1}}})

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], 
                Line[{{0, 0}, {5, 0}}], Dashing[{0.02, 0.03}], 
                Line[{{5, 0}, {6, 0}}]}, Line[#]}, AspectRatio -> 1/8, 
            Axes -> {True, False}]] &, t];

Figure

NIntersection[HalfLine[{{0, 0}, {1, 0}}], 
          Segment[#]] & /@ (t = {{{-2, 0}, {-1, 0}}, {{-1, 
                0}, {0, 0}}, {{-1, 0}, {1, 0}}, {{1, 0}, {2, 
                0}}})

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], 
                Line[{{0, 0}, {5, 0}}], Dashing[{0.02, 0.03}], 
                Line[{{5, 0}, {6, 0}}]}, Line[#]}, AspectRatio -> 1/8, 
            Axes -> {True, False}]] &, t];

Figure

When the segment s is degenerate its intersection with h is empty or weak (Fig. )

NIntersection[
          HalfLine[{{0, 0}, {1, 0}}], #] & /@ (t = {{-2, 0}, {0, 
              0}, {1, 0}, {}})

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], 
                Line[{{0, 0}, {5, 0}}], Dashing[{0.02, 0.03}], 
                Line[{{5, 0}, {6, 0}}]}, AbsolutePointSize[5], Point[#]}, 
            AspectRatio -> 1/8, Axes -> {True, False}]] &, 
      Drop[t, -1]];

Figure

Intersection of a segment and…

The intersection of two nondegenerate segments can be empty, weak, or strong. Fig.  shows these possibilities when the segment s is kept fixed and the other  varies in a list sl .

s = {{1, 2}, {3, 4}}; 
    sl = {{{0, 1/2}, {1/2, 0}}, {{0, 4}, {2, 0}}, {{0, 4.2}, {2.1, 0}}, {{0, 
            6.8}, {6.8, 0}}, {{0, 7}, {7, 0}}}; NIntersection[Segment[s], 
          Segment[#]] & /@ sl

Show[
      Graphics[{{RGBColor[0, 0, 1], Line[s]}, Line /@ sl, 
          MapIndexed[Text[First[#2], First[#1], {1.5, 0}] &, sl]}], 
      PlotRange -> All];

Figure

Other intersections of two non degenerate segments are shown in Fig.  (the endpoints of the varying segment are on opposite sides with respect to the constant segment)

t = Table[{{i, -1}, {i, 1}}, {i, 0, 6}]; 
    s = {{1, 0}, {5, 0}}; NIntersection[Segment[s], Segment[#]] & /@ 
      t

Show[
      Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], Line[s]}, Line /@ t}], 
      AspectRatio -> 1/6];

Figure

In Fig.  the varying segment has an endpoint on the line determined by the segment that is kept fixed.

t = Table[{{i, 0}, {i, 1}}, {i, 0, 6}]; 
    s = {{1, 0}, {5, 0}}; NIntersection[Segment[s], Segment[#]] & /@ 
      t

Show[
      Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], Line[s]}, Line /@ t}], 
      AspectRatio -> 1/6];

Figure

In Fig.  the varying segment is contained in the line determined by the fixed segment.

t = Table[{{i, 0}, {i + 1, 0}}, {i, -1, 6}]; 
    s = {{1, 0}, {5, 0}}; NIntersection[Segment[s], Segment[#]] & /@ 
      t

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], Line[s]}, 
              Line[{Offset[{0, 
                      2}, #[[1]]], 
                  Offset[{0, 
                      2}, #[[2]]]}]}, 
            AspectRatio -> 1/8, Axes -> {True, False}]] &, 
      t];

Figure

The intersections of a nondegenerate segment with a varying degenerate one (the dot)

t = Table[{{i, 0}, {i, 0}}, {i, -1, 6}]; 
    s = {{1, 0}, {5, 0}}; NIntersection[Segment[s], Segment[#]] & /@ 
      t

Scan[
      Show[Graphics[{AbsoluteThickness[2], {GrayLevel[ .4], Line[s]}, 
              AbsolutePointSize[5], Point[First@#]}, AspectRatio -> 1/8, 
            Axes -> {True, False}]] &, t];

Figure

Some trivial intersections:

with the empty set

NIntersection[{}, #] & /@ {Segment[{{1, 0}, {5, 0}}], 
        Segment[{{1, 0}, {1, 0}}], {1, 0}, {}}

between degenerate segments

NIntersection[
          Segment[{{1, 2}, {1, 2}}], #] & /@ {Segment[{{1, 2}, {5, 0}}], 
        Segment[{{1, 2}, {1, 2}}], {1, 2}, {}}

as above with a point as first argument

NIntersection[Segment[{{1, 2}, {1, 2}}], #] & /@ {{1, 2}, {2, 
          1}}

intersection of points (degenerate segments):

NIntersection[#, {1, 2}] & /@ {{1, 2}, {2, 1}}

Polygons

Description

A path (i.e. a list of three or more points) is a polygon if the intersection of the segments determined by any pair of consecutive vertices coincides with the common endpoint. A polygon is simple if the intersection of any pair of nonconsecutive edges is empty.

In this section we consider lists of numeric points {(Subscript[p, 1]), (Subscript[p, 2]), …, (Subscript[p, n])} with the identifications (Subscript[p, 0]) = (Subscript[p, n]) and (Subscript[p, n + 1]) = (Subscript[p, 1]).

NPoligonQ[{(Subscript[p, 1]), (Subscript[p, 2]), …, (Subscript[p, n])}] yields True if the intersection of segments (Subscript[p, i - 1]) (Subscript[p, i]) and (Subscript[p, i]) (Subscript[p, i + 1]) is {(Subscript[p, i])} for every  1 ≤ i ≤ n. The (Subscript[p, i])'s must be numeric points.

NSimplePoligonQ[{(Subscript[p, 1]), (Subscript[p, 2]), …, (Subscript[p, n])}] yields True if NPoligonQ[{(Subscript[p, 1]), (Subscript[p, 2]), …, (Subscript[p, n])}] yields True and the intersection of any pair of nonconsecutive edges is empty (i.e. (Subscript[p, i]) (Subscript[p, i + 1])doesn't intersect  (Subscript[p, i + 2]) 
        (Subscript[p, i + 3]), …, (Subscript[p, n - 1]) (Subscript[p, n]),  (Subscript[p, n]) 
        (Subscript[p, 1]), …, (Subscript[p, i - 3]) (Subscript[p, i - 2]) for  1 ≤ i ≤ n).

The following list of points does not define a polygon since (Subscript[p, 1]) (Subscript[p, 2]) ⋂ (Subscript[p, 2]) 
            (Subscript[p, 3]) ≠ {(Subscript[p, 2])} (Fig. )

plist = {{0, 0}, {2, 0}, {1, 0}, {1/2, 1}}; NPolygonQ[plist]

Show[
      Graphics[{Line[ClosePath@plist], AbsolutePointSize[5], Point /@ plist, 
          MapIndexed[Text[First@#2, #1, {0, 2}] &, plist]}], 
      PlotRange -> All];

Figure

The following two lists define polygons which are not simple polygons (Fig. )

plist = {{0, 0}, {2, 0}, {1, -1}, {1/2, 1}}; 
    plist1 = {{0, 0}, {2, 0}, {3, 0}, {5/2, 1}, {2, 0}, {1/2, 
          1}}; {NPolygonQ[#], NSimplePolygonQ[#]} & /@ {plist, 
        plist1}

g = 
      Graphics[{Line[ClosePath@plist], AbsolutePointSize[5], Point /@ plist, 
          MapIndexed[Text[First@#2, #1, {0, 2}] &, plist]}, 
        PlotRange -> {All, {-1.5, 1.2}}]; 
    g1 = Graphics[{Line[ClosePath@plist1], AbsolutePointSize[5], 
          Point /@ plist1, 
          MapThread[
            Text[#2, #1, {0, 2}] &, {plist1, {1, "", 3, 4, "2,5", 
                6}}]}, PlotRange -> {All, {- .25, 1.2}}]; 
    Show[GraphicsArray[{g, g1}]];

Figure

A suitable small perturbation of vertex 5 makes the second polygon simple (Fig. )

plist1 = {{0, 0}, {2, 0}, {3, 0}, {5/2, 1}, {2, 1/3}, {1/2, 1}}; 
    NSimplePolygonQ[plist1]

g1 = 
      Graphics[{Line[ClosePath@plist1], AbsolutePointSize[5], 
          Point /@ plist1, MapIndexed[Text[First@#2, #1, {0, 2}] &, plist1]}, 
        PlotRange -> {All, {- .25, 1.2}}]; 
    Show[GraphicsArray[{g1}]];

Figure

Example I: Triangulation of a simple polygon

This example shows the implementation of a nonoptimal (O((n^2))) algorithm for the triangulation of a simple polygon. The algorithm is described in Ch. 1 of [].

The function Triangulate takes a list P of vertices of a simple polygon and returns a list of diagonals that triangulate the polygon. The diagonals are returned as pairs of positions of vertices in P.

Code

Triangulate calls TriangulateCCWPolygon which triangulates a simple polygon specified by a list of vertex labels (positions in P) inducing a counterclockwise (ccw) traversal of the boundary.

Recall that the vertex  \(TraditionalForm\`\(i\^*\)\) at the lexicographic minimum of P is convex and that a left turn occurs at \(TraditionalForm\`\(i\^*\)\) iff the list is ccw oriented.

Triangulate[P0_?NSimplePolygonQ] := 
        Block[{P = P0, 
            n = Length[
                P0]}, 
NMinLexOrdering[
                  P][[1]] // 
              Mod[# + {-1, 0, 1}, n, 1] & // 
            TriangulateCCWPolygon[If[
                    LeftTurnQ[P[[#]]], 
                    Identity, Reverse][Range@n]] &];

Triangles have no diagonal

TriangulateCCWPolygon[indices_] /; 
          Length[indices] === 3 :=  {};

Four or more vertices: find a diagonal and split the problem into two simpler subproblems.

A diagonal is found as follows:

TriangulateCCWPolygon[indices_] /; Length[indices] > 3 := 
        Module[{triple, inside, diagonal}, 
 (* 1 *)  
          triple = NMinLexOrdering[
                    P[[
                      indices]]][[1]] + {-1, 0, 1} // 
              indices[[Mod[#, Length@indices, 
                    1]]] &; 
 (* 
            2 *)  inside = 
            Select[Complement[indices, triple], 
              MatchQ[NPointLocation[
                    Convex[P[[triple]]], 
                    P[[#]]], 
                  InPoint | BorderPoint] &]; 
 (* 
            3 *) 
diagonal = 
            If[inside === {}, 
              triple[[{1, 
                  3}]], {triple[[2]], 
                Ordering[inside, -1, 
                    DoubleArea[
                        P[[{triple[[3]], 
                            triple[[1]], #}]]] &] // 
                  inside[[#[[1]]]] &}]; 
{Sequence @@ 
              TriangulateCCWPolygon[subpolygon[indices, diagonal]], diagonal, 
            Sequence @@ 
              TriangulateCCWPolygon[
                subpolygon[indices, Reverse@diagonal]]}];

subpolygon takes a list of vertex labels L, a diagonal diag, and returns the sublist determined by the ordered pair diag:

subpolygon = 
        Function[{L, 
            diag}, Take[#, Position[#, 
                    diag[[2]]][[1, 1]]] &[
            RotateLeft[
              L, Position[L, 
                    diag[[1]]][[1, 1]] - 1]]];

Test

Triangulation of a triangle: no diagonal

Triangulate[{{0, 1}, {1, 1}, {1, 0}}]

Triangulation of a polygon with 4 vertices: one diagonal

Triangulate[{{0, 1}, {1, 1}, {1, 0}, {0, 0}}]

A nonconvex polygon with 5 vertices: 2 diagonals

Triangulate[{{0, 1}, {1, 1}, {1/2, 1/3}, {1, 0}, {0, 0}}]

A polygon with n vertices has n - 3 diagonals

𝒫 = {{4, 10}, {5, 8}, {3, 7}, {5, 6}, {6, 8}, {8, 
            8}, {7, 7}, {7, 5}, {11, 5}, {11, 0}, {7, 0}, {7, 3}, {9.7, 
            2}, {8, 2}, {8, 1}, {10.5, 1}, {10, 4}, {6, 4}, {6, 0}, {4, 
            0}, {4, 2}, {3, 1}, {3, 4}, {2, 3}, {2, 6}, {0, 4}, {1, 8}, {3, 
            8}, {2, 9}};
diagonals = Triangulate[𝒫]
Length[𝒫] - Length[%] === 3

The triangulation:

Show[
        Graphics[{Line[
              ClosePath@𝒫], 
MapIndexed[
              Text[#2[[1]], #1] &, 𝒫], {AbsoluteDashing[{3, 4}], 
              Orange, Line[𝒫[[#]]] & /@ 
                diagonals}}], 
AspectRatio -> 
          Automatic];

Example II: Triangulation of a monotone polygon

This is a very classical topic in computational geometry (e.g. see Section 3.3 of [],  Section 2.1.2 of [], or Section 6.2.2.1 of []). Let's consider a y-monotone polygon, i.e. a polygon whose intersection with any horizontal line is a connected set (the empty set, a point, or a segment). The boundary of such a polygon is the union of two chains connecting the highest and the lowest verteices (left chain and right chain). In linear time one can find these chains and merge them so that the vertices are sorted from the highest to the lowest (the initial event queue).

The algorithm processes one vertex at a time starting from the highest one and generating diagonals that split the polygon into a triangle and a polygon with one less side, still to be triangulated. More precisely, the diagonals drawn from v  end at the vertices above it which are visible from v. A stack collects the vertices above the current vertex but not reached by a diagonal.

This implementation (see TriangulateMP) considers a generic polygon 𝒫 (a list of points), and a direction of monotonicity  (for instance y-monotonicity is expressed by the option Direction->{{0,-1},{1,0}}). The output is a list of diagonals triangulating 𝒫 (list of pairs of vertex positions).

Code

Two options can be specified to TriangulateMP(vertices, options) :

Verbose->True|False with the obvious meaning, and

Direction->{...} specifies a direction of monotonicity. The default is y-monotonicity. This choice amounts to sweeping the plane by means of a horizontal line descending from +∞

Options[TriangulateMP] = {Direction -> {{0, -1}, {1, 0}}, 
          Verbose -> False};

\(\(TriangulateMP[pts0_, opts___?OptionQ] := 
        Block[{pts = pts0, dir, n = Length@pts0, chainof, stack, 
            queue, δ, δstack = {}, MSG}, 
          dir = Rule[
              Direction, \(Direction /. {opts}\) /. 
                Options[
                  TriangulateMP]]; 
If[\(Verbose /. {opts}\) /. Options[TriangulateMP], MSG = Print, 
            MSG[___] := Null]; 
queue = 
            initqueue[step1[step0]]; MSG["initial queue: ", queue]; 
          stack = Table[
              0, {n}]; 
triangmpolygon];\)\)

step0

Get a CCW oriented list of vertex labels with the highest vertex first in the list. Uses the global variables  pts, n, dir

step0 := 
        With[{m = NMinLexOrdering[pts, 
                  dir][[1]]}, Mod[
                m + {-1, 0, 1}, n, 
                1] // If[
                    LeftTurnQ[pts[[#]]], 
                    Identity, Reverse][Range[n]] & // 
            RotateLeft[#, Position[#, m][[1, 
                    1]] - 1] &];

step1

This step uses the output of step0 (say inds) to generate the list {left chain, right chain}. Furthermore it defines the function chainof for each vertex other than the highest and the lowest.

step1 = 
        Function[inds, 
          With[{M = Position[
                    inds, NMaxLexOrdering[pts, 
                        dir][[1]]][[1, 
                  1]]}, {(Map[(chainof[#] = -1) &, Rest[#]]; #) &@
                Take[inds, 
                  M - 1], (Map[(chainof[#] = +1) &, 
                      Most[#]]; #) &@Reverse@
                  Drop[inds, M - 1]}]];

initqueue

Merge the left and right chains returned by step2 to get the initial queue of events

initqueue[{lch_, rch_}] := 
        Module[{il = 1, nl = Length@lch, ir = 1, nr = Length@rch, 
            stk = {}}, 
While[(il ≤ nl ∧ 
                ir ≤ nr), 
            If[(NMinLexOrdering[
                      pts[[{lch[[
                            il]], 
                          rch[[
                            ir]]}]], 
                      dir][[1]] === 
                  1), stk = {stk, 
                  lch[[il++]]}, 
              stk = {stk, 
                  rch[[ir++]]}]]; 
          While[il ≤ nl, 
            stk = {stk, 
                lch[[il++]]}]; 
While[ir ≤ nr, 
            stk = {stk, 
                rch[[ir++]]}]; 
          Flatten[stk]];

triangulate monotone polygon

Both the event queue and the vertex stack are implemented as linear lists (queue and stack;  the front of the queue and the top of the stack are pointed to by q and σ, respectively). The generated diagonals are efficiently collected as pairs δ[i, j] by δstack.

\(\(triangmpolygon := 
        Module[{σ = 2, q = 2, x, newδ}, 
 (*\ 
            store\ a\ new\ diagonal\ *) \ 
          newδ = \((δstack = {δstack, #}; 
                MSG["generated ", #])\) &; 
 (*\ 
            initialize\ the\ stack\ of\ vertices\ *) 
stack[[{1, 2}]] = 
            queue[[{1, 
                2}]]; 
 (*\ 
            sweep\ *) 
While[\((\(++q\))\) < 
              n, 
x = 
              chainof[queue[[q]]] =!= 
                chainof[
                  stack[[σ]]]; 
MSG["*** vertex ", 
              queue[[
                q]], " and top of stack on opposite chains?: ", x]; 

If[
              x, \  (*\ 
                draw\ diagonals\ *) 
\(newδ[δ[
                      queue[[q]], 
                      stack[[#]]]] &\) /@ 
                Range[2, σ]; 
              stack[[1]] = 
                stack[[σ]]; 
              stack[[2]] = 
                queue[[q]]; σ = 
                2; showstatus[q, σ], 
 (*\ 
                else\ *) 
While[σ ≥ 
                    2 ∧ \  (*\ \(visibility\)\(?\)\ *) 
                    chainof[queue[[q]]]*
                      DoubleArea[
                        pts[[{queue[[
                              q]], 
                            stack[[σ]], stack[[σ - 
                                1]]}]]] > 0, 
newδ[δ[
                    queue[[q]], 
                    stack[[\(--σ\)]]]]]; \(σ++\); 
              stack[[σ]] = 
                queue[[q]]; 
              showstatus[
                q, σ]];
]; 
MSG["*** vertex ", 
            queue[[q]], " (last vertex)"]; 
\(newδ[δ[
                  queue[[q]], 
                  stack[[σ]]]] &\) /@ Range[σ - 1, 
              2, \(-1\)]; 
MSG["*** END ***"]; 
          List @@@ Flatten@δstack];\)\)

showstatus [q_, σ_] :=  
        MSG["queue: ", Drop[queue, q], "\nstack: ", 
          Take[stack, σ]];

Test

ymp = {{0, 
          0}, {4, -2}, {-1, -7}, {0, -9}, {3, -10}, {2, -11}, {-2, -12}, {-4, -13}, {-1, -5}, {-3, -3}, {-3, -1}}; diags = TriangulateMP[ymp]

The triangulated polygon is as follows:

Show[
        Graphics[{Line[ClosePath@ymp], 
MapIndexed[
              Text[#2[[1]], #1] &, 
              ymp], {AbsoluteDashing[{3, 4}], 
              Orange, Line[
                    ymp[[#]]] & /@ 
                diags}}], 
AspectRatio -> 
          Automatic];

A call with the Verbose option gives

TriangulateMP[ymp, Verbose -> True]