BPCG

Basic Computational Geometry in the Plane

Adriano Pascoletti

Implementation

Definition of the package context

BeginPackage["BPCG`"]

Declaration of exported functions

wProd::usage = "wProd[v1,v2] computes the wedge product of two \
points in the plane. \nwProd[{x1,y1},{x2,y2}] yields x1 y2 - x2 y1.";
dArea::usage = "dArea[{v1,v2,...}] yields the signed double area of \
a plane polygon with three or more vertices.";
LexMin::usage = "LexMin[{v1,v2,...}] yields the lexicographic \
minimum of a nonempty set of points in the plane. \nSee also: \
ArgLexMin.";
ArgLexMin::usage = "ArgLexMin[{v1,v2,...}] yields the position of \
the lexicographic minimum of {v1,v2,...} inside {v1,v2,...}. \nSee also: \
LexMin.";
CCWOrientedQ::usage = "CCWOrientedQ[{v1,v2,...}] yields True if a \
list of vertices defining a simple polygon is counterclockwise \
oriented.";
CartesianLine::usage = "CartesianLine[{v1,v2}] yields the cartesian \
equation of the line through v1 and v2 as a pure function. The traditional \
equation in x and y is given by CartesianLine[{v1,v2}][x,y]. \nSee also: \
ParametricLine.";
ParametricLine::usage = "ParametricLine[{v1,v2}] yields the \
parametric equations of the line through v1 and v2 as a pure function. The \
traditional equations in t are given by ParametricLine[{v1,v2}][t]. \See \
also: CartesianLine.";
PLDLine::usage = "PLDLine[{v1,v2},q] considers the directed line r \
defined by the ordered pair of points {v1,v2}, a query point q, and yields \
+1/-1/0 if q is to the left of/to the right of/on r. \nSee also: PLSegment, \
PLCPolygon, PLSPolygon.";
PLCPolygon::usage = "PLCPolygon[{v1,v2,...},q] yields +1/-1/0 if \
the query point q is internal to/external to/on the boundary of a convex \
polygon with vertices {v1,v2,...}. \nSee also: PLDLine,PLSegment, PLSPolygon.\
";
PLSegment::usage = "PLSegment[{v1,v2},q] yields +1/-1 if a query \
point is internal to/external to the segment with endpoints v1,v2, or 0 if q \
coincides with one of the endpoints. \nSee also: PLDLine, PLCPolygon, \
PLSPolygon.";
PLSPolygon::usage = "PLSPolygon[{v1,v2,...},q] yields 1/0/-1 if the \
query point q is internal to/on the boundary of/external to a simple polygon \
with vertices v1,v2,... \nSee also: PLDLine, PLSegment, PLCPolygon.";
IntLL::usage = "IntLL[{v1,v2},{v3,v4}] yields True if the line \
through v3,v4 intersects the line through v1 and v2. \nSee also: wIntSS, \
sIntSS, IntSS, IntLS.";
IntLS::usage = "IntLS[{v1,v2},{v3,v4}] yields True if the segment \
{v3,v4} intersects the line through v1 and v2. \nSee also: wIntSS, sIntSS, \
IntSS, IntLL.";
sIntSS::usage = "sIntSS[{v1,v2},{v3,v4}]  yields True if the \
segments defined by two pairs of points intersect at some point other than \
v1,v2,v3,v4 (strong intersection). \nSee also: wIntSS, IntSS, IntLS, \
IntLL.";
wIntSS::usage = "IntSS[{v1,v2},{v3,v4}] yields True if the \
intersection of two segments defined by two pairs of points contains an \
endpoint (weak intersection). \nSee also: sIntSS, IntSS, IntLS, \
IntLL.";
IntSS::usage = "IntSS[{v1,v2},{v3,v4}] yields True if the \
intersection of two segments defined by two pairs of points is't empty. \nSee \
also: wIntSS, sIntSS, IntLS, IntLL.";
SimplePolygonQ::usage = "SimplePolygonQ[{v1,v2,...}] yields True if \
the polygon defined by a list of points is simple.";
PolarSort::usage = "PolarSort[{v1,v2,...}] yields the permutation \
which sorts a list of points by increasing  polar angle (polar angles in the \
interval [0,2π)). \nSee also: PolarSort1.";
PolarSort1::usage = "PolarSort1[{v1,v2,...}] yields the permutation \
which sorts a list of points in the plane by increasing  polar angle (polar \
angles in the interval [-π,π)). \nSee also: PolarSort.";
CComponents::usage = "CComponents[V,A] yields the list of the \
connected components of a graph with vertex set V and adjacency lists A. A \
connected component consisting of vertices i1, i2,... is represented by a \
list {i1,i2,...}.";

Begin the implementation part

Begin["`Private`"]

Definition of exported functions

Wedge product and signed double area

wProd[v1 : {x1_, y1_}, v2 : {x2_, y2_}] := Det[{v1, v2}];
dArea[ps : {({_, _}) .. }] /; Length[ps] ≥ 3 :=
        Plus @@ (wProd @@@ (Append[ps, First[ps]] //
                  Partition[#, 2, 1] &));

Lexicographic minimum

LexMin[l : {({_, _}) .. }] := First[Sort[l]];
ArgLexMin[
          l : {({_, _}) .. }] := Position[l,
            LexMin[l]][[1,
          1]];

Orientation of the boundary (of a simple polygon)

CCWOrientedQ[poly : {({_, _}) .. }] /;
          Length[poly] ≥ 3 :=
        With[{i = ArgLexMin[poly], n = Length[poly]},
            dArea[poly[[{Mod[i - 1, n, 1], i,
                  Mod[i + 1, n, 1]}]]]] > 0;

Equations of a straight line

CartesianLine[{v1 : {x1_, y1_},
            v2 : {x2_, y2_}}] := (wProd[v2 - v1, {#1, #2} - v1] ==
              0) &;
ParametricLine[{v1 : {x1_, y1_}, v2 : {x2_, y2_}}] :=
        v1 + # (v2 - v1) &;

Point location

point - directed line

PLDLine[dl : {v1 : {x1_, y1_}, v2 : {x2_, y2_}}, q : {x_, y_}] /;
          v1 - v2 ≠ {0, 0} :=
        Sign[dArea[{v1, v2, q}]];

point - segment

PLSegment[seg : {v1 : {x1_, y1_}, v2 : {x2_, y2_}}, q : {x_, y_}] /;
        q == v1 ∨ q == v2 := 0;
    PLSegment[seg : {v1 : {x1_, y1_}, v2 : {x2_, y2_}}, q : {x_, y_}] :=
      If[x1 - x2 == 0,
        If[x - x1 == 0 ∧ Min[y1, y2] < y < Max[y1, y2], 1, -1],
        If[Min[x1, x2] ≤ x ≤ Max[x1, x2] ∧
            dArea[{v1, v2, q}] == 0, 1, -1]];

point - convex polygon

PLCPolygon[poly : {({_, _}) .. }, q : {x_, y_}] /;
          Length[poly] ≥ 3 :=
        If[CCWOrientedQ[poly], PLCP1[poly, q],
          PLCP1[Reverse[poly], q]];
PLCP1[poly : {({_, _}) .. }, q : {x_, y_}] :=
        Module[{edges = (Append[poly, First[poly]] //
                  Partition[#, 2, 1] &), signs},
          signs = Sign[dArea[#]] & @@@ Map[{Append[#, q]} &, edges];
          Which[MemberQ[signs, -1], -1, FreeQ[signs, 0], 1, True,
            0]];

point - simple polygon

PLSPolygon[poly : {({_, _}) .. }, q : {x_, y_}] /;
          Length[poly] ≥ 3 :=
        Module[{edges = (Append[poly, First[poly]] //
                  Partition[#, 2, 1] &)},
          If[Apply[
              Or, (PLSegment[#, q] ≥ 0) & /@ edges], 0,
            PLSPoygon1[edges, q]]];
PLSPoygon1[sp_, q : {qx_, qy_}] :=
        Module[{t},
          t = DeleteCases[sp, {{x1_, y_}, {x2_, y_}}]; 
t =
            DeleteCases[
              t, {{x1_, y1_}, {x2_, y2_}} /;
                Min[y1, y2] ≥ qy ∨
                  Max[y1, y2] < qy]; 
t =
            Map[If[#[[1,
                      2]] > #[[2,
                      2]], #, Reverse[#]] &,
              t]; 
If[
            OddQ[Count[t, e_ /; Sign[dArea[Prepend[e, {qx, qy}]]] == 1]],
            1, -1]];

Intersections

Line - line

IntLL[s1 : {v1 : {x1_, y1_}, v2 : {x2_, y2_}},
            s2 : {v3 : {x3_, y3_}, v4 : {x4_, y4_}}] /;
          v1 - v2 ≠ {0, 0} ∧
            v3 - v4 ≠ {0, 0} := (dArea[{v1, v2, v3}] -
                dArea[{v1, v2, v4}] ≠
              0) ∨ (dArea[{v1, v2, v3}] ==
              dArea[{v1, v2, v4}] == 0);

Line - segment

IntLS[r : {v1 : {x1_, y1_}, v2 : {x2_, y2_}},
            seg : {v3 : {x3_, y3_}, v4 : {x4_, y4_}}] /;
          v1 - v2 ≠ {0, 0} :=
        dArea[{v1, v2, v3}] dArea[{v1, v2, v4}] ≤ 0;

Segment - segment (strong intersection)

sIntSS[s1 : {v1 : {x1_, y1_}, v2 : {x2_, y2_}},
            s2 : {v3 : {x3_, y3_}, v4 : {x4_, y4_}}] /;
          v1 - v2 ≠ {0, 0} ∧ v3 - v4 ≠ {0, 0} :=
        dArea[{v1, v2, v3}] dArea[{v1, v2, v4}] < 0 ∧
          dArea[{v3, v4, v1}] dArea[{v3, v4, v2}] < 0;

Segment - segment (weak intersection)

wIntSS[s1 : {v1 : {x1_, y1_}, v2 : {x2_, y2_}},
            s2 : {v3 : {x3_, y3_}, v4 : {x4_, y4_}}] /;
          v1 - v2 ≠ {0, 0} ∧ v3 - v4 ≠ {0, 0} :=
        PLSegment[s1, v3] ≥ 0 ∨
          PLSegment[s1, v4] ≥ 0 ∨
          PLSegment[s2, v1] ≥ 0 ∨
          PLSegment[s2, v2] ≥ 0;

Segment - segment (strong or weak intersection)

IntSS[s1 : {v1 : {x1_, y1_}, v2 : {x2_, y2_}},
            s2 : {v3 : {x3_, y3_}, v4 : {x4_, y4_}}] /;
          v1 - v2 ≠ {0, 0} ∧
            v3 - v4 ≠ {0,
                0} := (dArea[{v1, v2, v3}] dArea[{v1, v2, v4}] < 0 ∧
              dArea[{v3, v4, v1}] dArea[{v3, v4, v2}] <
                0) ∨ (PLSegment[s1, v3] ≥ 0 ∨
              PLSegment[s1, v4] ≥ 0 ∨
              PLSegment[s2, v1] ≥ 0 ∨
              PLSegment[s2, v2] ≥ 0);

Simplicity test

SimplePolygonQ[poly : {({_, _}) .. }] /;
          Length[poly] ≥
            3 := (Sort[poly] == Union[poly]) ∧
          edgesOK[poly];
edgesOK[poly_] :=         Module[{n = Length[poly],
            edges = Partition[Append[poly, First[poly]], 2, 1], nonadjacent},
          nonadjacent =
            DeleteCases[Range[n],
                j_ /; MemberQ[{Mod[# - 1, n, 1], #, Mod[# + 1, n, 1]},
                    j]] &; ¬ Or @@
              Table[Or @@ (IntSS[
                          edges[[
                            i]], #] & /@
                      edges[[nonadjacent[
                          i]]]), {i, n}]];

Polar sorting

For polar angles in [0,2π)

PolarSort[s : {({_, _}) ... }] /; FreeQ[s, {0, 0}] :=
      With[{poslabels =
            Select[Range[Length[s]],
              s[[#, 2]] > 0 ∨
                  s[[#, 2]] ==
                      0  ∧
                    s[[#, 1]] > 0 &],
          neglabels =
            Select[Range[Length[s]],
              s[[#, 2]] < 0 ∨
                  s[[#, 2]] ==
                      0  ∧
                    s[[#, 1]] <
                      0 &]}, 
Join[
          Sort[poslabels,
            wProd[s[[#1]],
                  s[[#2]]] ≥ 0 &],
          Sort[neglabels,
            wProd[s[[#1]],
                  s[[#2]]] ≥ 0 &]]]

For polar angles in [-π,π)

PolarSort1[s : {({_, _}) ... }] /; FreeQ[s, {0, 0}] :=
      With[{poslabels =
            Select[Range[Length[s]],
              s[[#, 2]] > 0 ∨
                  s[[#, 2]] ==
                      0  ∧
                    s[[#, 1]] > 0 &],
          neglabels =
            Select[Range[Length[s]],
              s[[#, 2]] < 0 ∨
                  s[[#, 2]] ==
                      0  ∧
                    s[[#, 1]] <
                      0 &]}, 
Join[
          Sort[neglabels,
            wProd[s[[#1]],
                  s[[#2]]] ≥ 0 &],
          Sort[poslabels,
            wProd[s[[#1]],
                  s[[#2]]] ≥ 0 &]]]

Connected components

CComponents[V_, A_] :=
        Block[{ mark = Table[False, {Length[V]}],  cccounter = 0, s},
          Scan[If[mark[[#]] ===
                  False, ++cccounter; s[cccounter] = {}; DFS[A, #]] &,
            Range[Length[V]]];
          Table[Flatten[s[i]], {i,
              cccounter}]];

SetAttributes[DFS, HoldFirst];
DFS[a_, i_] := (  (*1*)           mark[[i]] = True;  (*2*)           s[cccounter] = {s[cccounter], i};  (*3*)           Scan[If[¬ mark[[#]],
                DFS[a, #]] &,
            a[[i]]]);

Epilogue

End[]

EndPackage[]