BPCGtest

Test package BPCG.m

Adriano Pascoletti

Dipartimento di Matematica e Informatica
Via delle Scienze, 206
I-33100 Udine (Italy)
pascolet@dimi.uniud.it
http://www.dimi.uniud.it/~pascolet

The package BPCG.m provides a set of very basic operations needed to solve more complex problems of computational geometry in the plane. Occasionally some Mathematica based solutions to some problems of computational geometry have appeared in MathSource, and a nice set of important functions is available in the Computational Geometry directory (contained in the Add-ons/Standard Packages/Discrete Mathematics/ directory), but it seems that no reasonable Mathematica package defining the basic operations is available. My belief is that Mathematica’ s possibility of dealing with symbolic, numeric, high precision numerics, and graphics offers an ideal environment for understanding, exploring, and prototyping computational geometry algorithms.  Here is what I tried to offer to my students: they have used a preliminary draft of this package to successfully display, solve, and animate many algorithms.

The package has been devoloped using Mathematica v. 4.

The functions can be grouped as follows:
Intersections: IntLL, IntLS, IntSS, wIntSS, sIntSS; Point location: PLDLine, PLSegment, PLCPolygon, PLSPolygon; Straight lines: CartesianLine, ParametricLine; Polygons: CCWOrientedQ, SimplePolygonQ; Sorting by polar angle: PolarSort, PolarSort1; Connected components: CComponents; Signed areas: wProd, dArea; Lexicographic order: LexMin, ArgLexMin.

References

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

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

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

Loading the package

If the file BPCG.m doesn't reside inside a directory known to Mathematica (a good place is the Mathematica 4.0 Files/AddOns/Applications/), setup your path to the directory containing BPCG.m (Input menu: Get File Path...) and load it:

AppendTo[$Path,DirectoryName[<your path to BPCG . m >]];

<< BPCG`

The therein defined functions are

?BPCG`*

Testing the package

Points, segments, polygons

The basic entities involved in the description of a geometric object in the plane are those of point, segment and polygon. Points are pairs of coordinates ({x, y}), segments are pairs of points ({{(Subscript[x, i]), (Subscript[y, i])}, {(Subscript[x, j]), (Subscript[y, j])}}). A set of points {(Subscript[v, 1]), …, (Subscript[v, n])} defines a polygon with vertices {(Subscript[v, 1]), …, (Subscript[v, n])} and  edges \(TraditionalForm\`e\_1 = {v\_1, v\_2},       e\_2 = {v\_2, v\_3}, …, e\_n = {v\_n, v\_1}\). A simple polygon is a polygon s.t. the intersection of two consecutive edges is the common vertex ((Subscript[e, i]) ⋂
          (Subscript[e, i + 1]) = {(Subscript[v, i + 1])}), and the intersection of any pair of non consecutive edges is empty ((Subscript[e, i]) ⋂ (Subscript[e, j]) = ∅ for all j ≠ i - 1, i + 1).

The function SimplePolygonQ[poly] takes a list poly consisting of three or more points, interprets it as defining a polygon and returns True if the polygon is simple.

Examples

Here poly defines a simple polygon:

poly = {{0, 0}, {1, -1}, {2, 0}, {2.5`, -1.1`}, {3, 0}, {4,
          0}, {4, 1}, {2.5`, 2.5`}, {0.5`, 1}};
    SimplePolygonQ[poly]

Show[Graphics[{Thickness[0.01], Line[Append[poly, First[poly]]]}],
        Axes -> True];

Inserting a new vertex  {4.5, 0} after vertex {4, 1} destroys the simplicity:

poly = {{0, 0}, {1, -1}, {2, 0}, {2.5`, -1.1`}, {3, 0}, {4,
          0}, {4, 1}, {4.5, 0}, {2.5`, 2.5`}, {0.5`, 1}};
    SimplePolygonQ[poly]

Show[Graphics[{Thickness[0.01], Line[Append[poly, First[poly]]]}],
        Axes -> True];

The following isn’t simple because vertex {0, 2} belongs to segment {{0, 0}, {0, 3}}

poly = {{0, 0}, {1, 1}, {0, 2}, {1, 3}, {0, 3}};
    SimplePolygonQ[poly]

Show[
        Graphics[{Thickness[0.010], Line[Append[poly, First[poly]]],
            PointSize[0.025], Point /@ poly}],
        Axes -> True];

A famous theorem by C. Jordan’s states that a simple polygon divides the plane into two regions, one bounded (interior) and one unbounded (exterior). Traversing the boundary counterclockwise the interior is always to the left and the exterior to the right. CCWOrientedQ[poly] returns True if poly determines a counterclockwise traversal of the boundary.

Example

poly = {{0, 0}, {0, 1}, {1, 0}}; CCWOrientedQ[poly]

CCWOrientedQ[poly // Reverse]

Wedge product and signed area

The wedge product of two points in the plane (Subscript[v, 1]) = {(Subscript[x, 1]), (Subscript[y, 1])},  (Subscript[v, 2]) = {(Subscript[x, 2]), (Subscript[y, 2])}, denoted (Subscript[v, 1])⋀(Subscript[v, 2]), is det[({{(Subscript[x, 1]), (Subscript[y, 1])}, {(Subscript[x, 2]), (Subscript[y, 2])}})]=(Subscript[x, 1]) (Subscript[y, 2]) - (Subscript[x, 2]) (Subscript[y, 1])and it is twice the area of the triangle O (Subscript[v, 1]) (Subscript[v, 2]). These three points are  collinear iff (Subscript[v, 1])⋀(Subscript[v, 2]) = 0. If (Subscript[v, 1])⋀(Subscript[v, 2]) ≠ 0 then (Subscript[v, 1])⋀(Subscript[v, 2])is positive iff the triple  O (Subscript[v, 1]) (Subscript[v, 2]) determines a counterclockwise path. The signed double area of the triangle O (Subscript[v, 1]) (Subscript[v, 2])is given by (Subscript[v, 1])⋀(Subscript[v, 2]). This formula can be extended to an arbitrary simple polygon 𝒫 with vertices {(Subscript[v, 1]), …, (Subscript[v, n])} because (Subscript[v, 1])⋀(Subscript[v, 2]) +
        (Subscript[v, 2])⋀(Subscript[v, 3]) + ⋯ + (Subscript[v, n])⋀(Subscript[v, 1]) is twice the area of 𝒫 and it is positive iff the list {(Subscript[v, 1]), …, (Subscript[v, n])} defines a counterclockwise traversal of the boundary of 𝒫. The wedge product is given by wProd[{x1,y1},{x2,y2}], and the signed double area of a simple polygon by dArea[{{x1,y1},...,{xn,yn}}].

Examples
The triangle defined by O and the two points {1, 3}, {-1, 4}has positive double area 7

dArea[{{0, 0}, {1, 3}, {-1, 4}}]

Show[
        Graphics[{{Thickness[0.010],
              Line[{{0, 0}, {1, 3}, {-1, 4}, {0, 0}}]}, PointSize[0.025],
            Point /@ {{0, 0}, {1, 3}, {-1, 4}},
            Line[{{-2, 0}, {2, 0}}], Line[{{0, -1}, {0, 5}}],
            Line[{{-1, 4}, {-1, 0}}], Line[{{1, 3}, {1, 0}}]}],
        Axes -> True];

The signed double area of a generic triangle is

dArea[{{x1, y1}, {x2, y2}, {x3, y3}}]

Point Location

Various queries are possible when we are given a pair of distinct points (Subscript[v, 1]), (Subscript[v, 2])in the plane and a query point q (in the plane, of course) [2].

If r denotes the line through (Subscript[v, 1]) and (Subscript[v, 2]) directed from (Subscript[v, 1]) to (Subscript[v, 2]) then there are three possibilities for q: it is to the left of r, on r, to the right of r. The function PLDLine[{v1,v2},q] returns +1, 0 and -1 respectively; PLSegment[{v1,v2},q]  returns +1 when q belongs to the (relative) interior of the segment s defined by (Subscript[v, 1])and (Subscript[v, 2]), 0 if q coincides with (Subscript[v, 1])or (Subscript[v, 2]), -1 if q is in the exterior of s.

Examples

The following example illustrates the point location problem with respect to the straight line through {1, 2} and {3, 3}; the query points are {10, 7}, {21,1 2},{18, 8}

r = {{1, 2}, {3, 3}};
    qlist = {{10, 7}, {21, 12}, {18, 8}}; PLDLine[r, #] & /@
      qlist

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

Now the location of a set of query points (qlist) w.r.to the vertical segment defined by {1, 2} and {1, 2.5} is considered:

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

This answer tells us that the first two points are external, the third coincides with an endpoint and the fourth lies in the relative interior.

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

A convex polygon with a counterclockwise boundary is the intersection of the left half-spaces defined by the (directed) lines through the edges. The function PLCPolygon[poly, q] uses this idea and returns +1 / -1 / 0 if a query point q is internal , on the boundary, or external to a convex polygon poly, respectively. The position of a query point w.r. to a simple polygon 𝒫 is handled by PLSPolygon[𝒫, q] and returns +1 / 0 / -1 if q is internal / on the boundary / external to 𝒫 [3].

Examples

This list  cp specifies a clockwise oriented boundary of a convex polygon and qs a sequence of queries:

cp = {{0, 2}, {1, 1}, {1, 0}, {0, 0}};
    qs = {{2, 2}, { .5,  .8}, { .5, 1.5}}; PLCPolygon[cp, #] & /@
      qs

Show[
        Graphics[{Thickness[0.008], Line[Append[#, First[#]] &[cp]],
            Point /@ cp, PointSize[0.025], Point /@ qs,
            MapIndexed[Text[First[#2], #1, {1, 1.5}] &, qs]}]];

The following example takes a simple polygon sp and a sequence of query points qs

sp = {{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}]; PLSPolygon[sp, #]
& /@
      qs

Points 1 and 11 are external, 3,4 and 5 internal, and 2, 6, 8, 9, 10 on the boundary.

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

Intersections

The intersection test between a line r determined by two points (Subscript[v, 1]), (Subscript[v, 2]) and the segment s defined by two points (Subscript[v, 3]), (Subscript[v, 4]) is handled by IntLS[{v1,v2},{v3,v4}]: the result is True if r ⋂ s ≠ ∅. If the segment s' defined by (Subscript[v, 1]), (Subscript[v, 2]) is considered, instead two kinds of intersection are tested. Strong intersections: s' strongly intersects s if their relative interiors have nonempty intersection (sIntSS[s', s] ⟶ True); weak intersection occurs if s' ⋂ s contains an endpoint (wIntSS[s', s] ⟶ True). The function IntSS[s', s] yields True if strong or weak intersection occurs, i.e. if s' ⋂ s ≠ ∅.
The function IntLL[{v1,v2},{v3,v4}] yields True if the line through (Subscript[v, 1]), (Subscript[v, 2]) intersects the line through (Subscript[v, 3]), (Subscript[v, 4]).

Examples

Line-Segment: straight line r through {1, 2} and {2, 4} and a sequence of segments (slist). Only the last two segments intersect r

r = {{1, 2}, {2, 4}};
slist = {{{3, 2}, {4, 2}}, {{3, 2}, {4, 8}}, {{3, 2}, {8,
              17}}};
IntLS[r, #] & /@ slist

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

Segment - segment: the first segment s1 is kept fixed, the second varies in a list (slist)

s1 = {{1, 2}, {3, 4}};
    slist = {{{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}}};
    Print["weak: ", wIntSS[s1, #] & /@
        slist, "\nstrong: ", sIntSS[s1, #] & /@
        slist, "\nweak or strong: ", IntSS[s1, #] & /@
        slist];

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

Line - line: consider the lines determined by  s1 and a generic element of slist. All these pairs of lines are non parallel:

IntLL[s1, #] & /@ slist

If s2 is parallel to s1 the corresponfing lines are either distinct

s2 = {0, 1/2} + # & /@ s1; IntLL[s1, s2]

or coincident

s2 =  .1 (s1[[2]] -
                  s1[[1]]) + # & /@
        s1; IntLL[s1, s2]

Lexicographic minimum

The lexicographic minimum of a set of points s is returned by LexMin[s] and the point where this minimum is attained is computed by ArgLexMin[s]. Of course {(Subscript[x, 1]), (Subscript[y, 1])} precedes {(Subscript[x, 2]), (Subscript[y, 2])} if (Subscript[x, 1]) < (Subscript[x, 2]) or (Subscript[x, 1]) = (Subscript[x, 2])and (Subscript[y, 1]) < (Subscript[y, 2]).

Example

The lexicographic minimum of this set of points is {0, 0} and it is attained at the first point of the list, so  ArgLexMin returns 1

p = {{0, 0}, {1, -1}, {2, 0}, {2.5, -1.1}, {3, 0}, {4, 0}, {4,
          1}, {2.5, 2.5}, {0.5, 1}}; LexMin[p]
ArgLexMin[p]

Straight lines

Let (Subscript[v, 1]), (Subscript[v, 2])denote two distinct points of the plane.  The ‘pure form’ cartesian equation of the line through them is given by CartesianLine[v1,v2] and the more familiar equation with x and y by CartesianLine[v1,v2][x,y]. Similarly the ‘pure’ parametric equation is given by ParametricLine[v1,v2], and its form with a parameter t by ParametricLine[v1,v2][t].

Examples

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

This is another example (a vertical line)

CartesianLine[{{1, 2}, {1, 3}}]
CartesianLine[{{1, 2}, {1, 3}}][x, y]
ParametricLine[{{1, 2}, {1, 3}}]
ParametricLine[{{1, 2}, {1, 3}}][t]

Sort by polar angle

The polar angle of a point v ≠ 0 is a real number in the interval [ 0, 2  π) or in [ -π, +π). In the framework of computational geometry both criteria are used, and this package provides two functions for sorting a set s  of nonzero points of the plane by increasing polar angle (PolarSort[s] and PolarSort1[s] consider the first and second interval, respectively). Both functions return the permutation that yields the sorting. The sorted version of s is given by s[[PolarSort[s]]] or s[[PolarSort1[s]]]

Examples

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

This is the desired permutation of {1, 2,...}: point 2 has minimum polar angle, point 7 maximum

With[{g1 =
            Graphics[{Line[{{0, 0}, #}] & /@ s,
                MapIndexed[Text[First[#2], #1] &, s]}],
          g2 = Graphics[{Line[{{0, 0}, #}] & /@ s,
                MapIndexed[
                  Text[ToString[First[#2]] <> "'",
                      s[[#1]]] &,
                  PolarSort[s]]}]},
        Show[GraphicsArray[{g1, g2}]]];

When [ -π,  +π ) is considered the minimum polar angle is that of point 6, the maximum that of 8

PolarSort1[s]

With[{g1 =
            Graphics[{Line[{{0, 0}, #}] & /@ s,
                MapIndexed[Text[First[#2], #1] &, s]}],
          g2 = Graphics[{Line[{{0, 0}, #}] & /@ s,
                MapIndexed[
                  Text[ToString[First[#2]] <> "'",
                      s[[#1]]] &,
                  PolarSort1[s]]}]},
        Show[GraphicsArray[{g1, g2}]]];

Connected components

An undirected graph is connected if for every pair of vertices there exists a path joining them. A connected component C of an undirected graph G is a subset of vertices s.t. (1) for every pair of vertices in C there exists a path in C  joining them, (2) C is maximal, i.e. no connected subset of G properly contains C. [1]
CComponents[V,A] computes the connected components of a graph described by a vertex set V and a set of adjacency lists A, where the i-th adjacency list A[[ i ]] collects the labels of vertices j s.t. there exists an edge with endpoints i and j. CComponents[V,A] yields a list of sublists, where each sublist collects the (labels of) vertices belonging to the same connected component.

Examples

A graph on 16 vertices and 5 connected components described by its adjacency lists

G = {V = {{1, 0}, {2, 0}, {3, 0}, {1, 1}, {2, 1}, {3, 1}, {1, 2}, {2,
                2}, {3, 2}, {5, 0}, {7, 0}, {5, 2}, {7, 2}, {2.5, 1.5}, {4,
                1}, {6, 1}},
          A = {{2, 4}, {1, 5, 3}, {2, 6}, {1, 5, 7}, {2, 6, 4, 8}, {5, 3,
                9}, {4, 8}, {7, 5, 9}, {8, 6}, {11, 12}, {10, 13}, {10,
                13}, {12, 11}, {}, {}, {}}};
CComponents[V, A]

With[{f =
            Function[{i}, Line[{V[[
                        i]],
                      V[[#]]}] & /@
                Select[A[[
                    i]], (# > i) &]]},
        Show[Graphics[{f /@ Range[Length[V]], PointSize[0.02], Point /@ V,
              MapIndexed[Text[First[#2], #1, {1.5, -1. }] &, V]}],
          PlotRange -> All]];

The empty graph has no component:

CComponents[{}, {}]

Using the package

The functions supplied by BPCG can used to perform more complex tasks. As an example we develop a non optimal  O((n^2)) algorithm for triangulating a simple polygon (Triangulate) (see [2]). It takes a list of points defining the boundary of a simple polygon and returns a list of diagonals of a triangulation. The diagonal connecting the i-th to the j-th vertex is returned as  {i, j}.

Triangulation of a simple polygon

Triangulate[p] prepares what the core of the algorithm (function Triang2) needs: a counterclockwise list L of vertex labels, and a global list of vertices (polygon)

Triangulate[p : {({_, _}) .. }] /; Length[p] ≥ 3 :=
        Block[{polygon = p,
            L = If[CCWOrientedQ[p], #, Reverse[#]] &[Range[Length[p]]]},
          Triang2[L]];

Triang2 works as follows:
1 - It finds the vertex attaining the lexicographic minimum ( i ) its predecessor ( (Subscript[i, 1])) and successor ( (Subscript[i, 2]));
2 - finds the set S of vertices other than (Subscript[i, 1]), i, (Subscript[i, 2]) in the interior of the triangle ((Subscript[i, 2]), i, (Subscript[i, 1])) (note the counterclockwise boundary) or on the segment (Subscript[i, 1]) (Subscript[i, 2])(see getS).
3 - If S is empty then (Subscript[i, 1]) (Subscript[i, 2]) is a diagonal, else find a vertex j of S maximizing the signed area of the triangle ((Subscript[i, 2]), (Subscript[i, 1]), j) (see maximizer) and the diagonal to be returned is i j . 
4 - Use the diagonal just found to split the boundary into two counterclockwise oriented boundaries (see splitpolygon) and recursively triangulate the subpolygons.

Triang2[L_] /; Length[L] == 3 := {};   (*      no diagonal in a triangle *) 
Triang2[L_] :=
      Module[{n = Length[L], i, i1, i2, S, diag, subpolygons},  (*          step 1 *) 
i =
          ArgLexMin[
            polygon[[
              L]]]; 
i1 =
          Mod[i - 1, n, 1];
        i2 = Mod[i + 1, n, 1]; 
 (*          step 2 *)  
S =
          getS[L, {i1, i, i2}]; 
 (* step 3 *)         diag = If[S == {},
            L[[{i1,
                i2}]], {L[[
                i]],
              maximizer[S,
                L[[{i2,
                    i1}]]]}]; 
 (*          step 4 *) 
subpolygons =
          splitpolygon[L, diag]; 
 (*          recursion *) 
Join[{diag},
          Triang2[First[subpolygons]],
          Triang2[Last[subpolygons]]]];

getS[L,{i1,i,i2}]
1 - Compute the list Lc of vertex labels other than i, (Subscript[i, 1]), (Subscript[i, 2])
2 - Select vertices of Lc in the relative interior of segment (Subscript[i, 1]) (Subscript[i, 2]) or in the interior of the triangle (Subscript[i, 1]), i, (Subscript[i, 2])

getS[L_, indices : {i1_, i_, i2_}] :=
        With[{ (* step 1 *)             Lc = Delete[L, {{i1}, {i}, {i2}}] }, 
 (*            step 2 *) 
Select[Lc,
            PLSegment[
                    polygon[[
                      L[[{i1,
                          i2}]]]],
                    polygon[[#]]] >
                  0 ∨
                PLCPolygon[
                    polygon[[
                      L[[
                        indices]]]],
                    polygon[[#]]] >
                  0 &]];

maximizer[S, {i2, i1}]
Let S(≠∅) denote the set of labels returned by getS:
1 - Compute the max area M of triangles (Subscript[i, 2]) (Subscript[i, 1]) j (for all j∈S)
2 - return a label \(TraditionalForm\`\(i\^*\)\) s.t. \(TraditionalForm\`area(i\_2, i\_1, \(i\^*\)) = M\)

maximizer[S_, {i2_, i1_}] :=
        With[{M =
              Max[dArea[
                      polygon[[{i2,
                          i1, #}]]] & /@ S]},
          Select[S, (dArea[
                      polygon[[{i2,
                          i1, #}]]] == M) &] //
            First];

splitpolygon[L, {d1, d2}]
Let L denote a counterclockwise list of vertex labels of a polygon , d1, d2 the labels of the endpoints of a diagonal. In order to split the polygon with vertices in L into two (counterclockwise) subpolygons
1 - Find the positions i1, i2 of d1, d2 inside L;
2 - The first subpolygon is given by the labels from i1 to i2, the second by those from i2 to i1

splitpolygon[L_, {d1_, d2_}] :=
      With[{ (* step 1 *)           i1 = Position[L, d1][[1,
              1]],
          i2 = Position[L, d2][[1,
              1]], n = Length[L]}, 
 (*          step 2 *) 
{NestWhileList[Mod[# + 1, n, 1] &,
              i1, (L[[#]] ≠
                    d2) &] //
            L[[#]] &,
          NestWhileList[Mod[# + 1, n, 1] &,
              i2, (L[[#]] ≠
                    d1) &] //
            L[[#]] &}]

Example: triangulate the following simple polygon

p = {{4, 10}, {5, 8}, { .85, 7}, {5, 6}, {4, 2}, {3, 1}, {3, 4}, {2,
            3}, {2, 6}, {0, 4}, {1, 8}, {3, 8}, {2, 9}};

Show[
        Graphics[{Line[Append[p, First[p]]],
            MapIndexed[Text[First[#2], #1] &, p]}]];

The list of diagonals is returned by

diagonals  = Triangulate[p]

Show[
        Graphics[{Line[Append[p, First[p]]], {Dashing[{0.01, 0.01}],
              RGBColor[1, 0,
                0], Line[p[[#]]] & /@
                diagonals},
            MapIndexed[Text[First[#2], #1] &, p]}]];