geometry-simple (empty) → 0.1.0.0

raw patch · 17 files changed

+4162/−0 lines, 17 filesdep +basedep +binarydep +bytestring

Dependencies added: base, binary, bytestring, containers, deepseq, text, transformers, vector

Files

• CHANGELOG.md +21/−0
• LICENSE +21/−0
• README.md +164/−0
• docs/GEOS-DIFFERENCES.md +177/−0
• geometry-simple.cabal +70/−0
• src/Data/Geometry.hs +39/−0
• src/Data/Geometry/Internal.hs +472/−0
• src/Data/Geometry/SimpleFeatures.hs +564/−0
• src/Data/Geometry/Topology/Buffer.hs +296/−0
• src/Data/Geometry/Topology/Measures.hs +145/−0
• src/Data/Geometry/Topology/Overlay.hs +311/−0
• src/Data/Geometry/Topology/Planar.hs +421/−0
• src/Data/Geometry/Topology/Relations.hs +260/−0
• src/Data/Geometry/Topology/Snapping.hs +80/−0
• src/Data/Geometry/Topology/Unary.hs +307/−0
• src/Data/Geometry/WKB.hs +340/−0
• src/Data/Geometry/WKT.hs +474/−0

@@ -0,0 +1,21 @@+# 0.1.0.0++Initial release.++- Seven Simple Features geometry families with XY, XYZ, XYM, and XYZM layouts.+  Points, coordinate sequences, and polygon rings retain their own layouts.+  Coordinates and multipoints use unboxed storage.+- ISO WKB and WKT codecs with construction checks and exact finite-coordinate+  round trips.+- Coordinate and point accessors, measurements, envelopes, and convex hulls.+  Indices start at zero. Point observers preserve stored ordinates and layouts.+- Pure Haskell topology checks, spatial predicates, DE-9IM relations, distance,+  overlays, and round buffers. Overlays and buffers return+  `Either TopologyException Geometry`. They use bounded snapping when output+  rounding changes topology; exhausted retries return `Left PrecisionFailure`.+  Unrepresentable buffer offsets return `Left CoordinateOverflow`.+  The retry schedule follows GEOS 3.13.1, without its self-union and+  precision-grid attempts, so results and precision failures can differ. See the [overlay precision policy](https://github.com/Tritlo/geometry-simple/blob/main/docs/GEOS-DIFFERENCES.md#overlay-precision).+- Constructed planar results use XY. Polygon exteriors run counterclockwise+  and holes clockwise. Hull vertices use a deterministic XY order.+- Measured-location queries with linear interpolation along segments.

@@ -0,0 +1,21 @@+MIT License++Copyright (c) 2026 Matthias Pall Gissurarson++Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all+copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE+SOFTWARE.

@@ -0,0 +1,164 @@+# geometry-simple++A pure Haskell library for the seven core Simple Features geometry families:+points, line strings, polygons, their three multi-geometry forms, and geometry+collections. It provides ISO WKB and WKT codecs, planar measurements, spatial+predicates, overlays, buffers, and measured-location queries.++Coordinates use unboxed vectors. The library has no database or native-library+dependency. DuckDB interchange through WKB or WKT requires a+[common coordinate layout](https://duckdb.org/docs/current/sql/data_types/geometry)+across all members.++## Example++```haskell+import Data.Geometry+import Data.Geometry.WKB (decodeWKB, encodeWKB)+import Data.Geometry.WKT (decodeWKT, encodeWKT)+import qualified Data.Geometry.SimpleFeatures as SF+import qualified Data.Vector.Unboxed as U++line = LineString (CoordinatesXY (U.fromList [XY 0 0, XY 2 2]))+point = PointGeometry (PointXY (XY 1 1))++intersectsLine = SF.intersects line point+firstX = SF.startPoint line >>= SF.pointX+binaryRoundTrip = encodeWKB line >>= decodeWKB+textRoundTrip = encodeWKT line >>= decodeWKT+```++Add `geometry-simple` and `vector` to your component's `build-depends`:++```cabal+build-depends:+  geometry-simple >=0.1 && <0.2,+  vector >=0.13 && <0.14,+```++WKB uses `ByteString` from `bytestring`. WKT uses `Text` from `text`.++## Geometry values++Each `Point` and `Coordinates` value has an XY, XYZ, XYM, or XYZM layout.+Empty points and coordinate sequences retain their layout, such as+`EmptyPoint DimXYZ`. Empty multi-geometries and empty collections have no+layout of their own. Writers give them their containing collection's layout,+or XY when no containing layout is available.+`PolygonRings` stores an exterior ring and a boxed vector of interior rings.+Each ring has its own layout. Collection members can have different layouts.++Coordinate sequences and multipoints use unboxed storage. Multilines,+multipolygons, and geometry collections use boxed vectors for their members.+Vector slices share their source buffers; use `U.force` to copy a slice.++`withPoint` and `withCoordinates` apply an operation without matching every+coordinate constructor. The `pointX`, `pointY`, `pointZ`, and `pointM` accessors+return `Nothing` for an empty point or an absent ordinate. Stored NaN ordinates+remain present. The `x`, `y`, `z`, and `m` functions operate on coordinate values.++The constructors do not validate geometry. Codecs check line lengths, ring+closure, and empty shells. Use `SF.isValid` to check polygon topology.+The derived `Eq` compares storage; `SF.equals` compares XY point sets.+As with `Double`, positive and negative zero compare equal under `Eq`.++CRS and SRID metadata belongs to the caller. `Data.Geometry.Internal` exposes+implementation details and can change without following the package versioning+policy.++## Operations++| Group | Functions in `Data.Geometry.SimpleFeatures` |+| --- | --- |+| Properties | `geometryType`, `dimension`, `coordinateDimension`, `spatialDimension`, `is3D`, `isMeasured`, `isEmpty` |+| Ordinates | `x`, `y`, `z`, `m`, `pointX`, `pointY`, `pointZ`, `pointM` |+| Members | `numGeometries`, `geometryN` |+| Line points | `numPoints`, `pointN`, `startPoint`, `endPoint`, `isClosed` |+| Rings | `exteriorRing`, `numInteriorRings`, `interiorRingN` |+| Measurements | `envelope`, `area`, `geometryLength`, `curveLength`, `perimeter`, `centroid`, `convexHull` |+| Topology | `boundary`, `isSimple`, `isRing`, `isValid`, `pointOnSurface` |+| Relations | `relate`, `relatePattern`, `equals`, `disjoint`, `intersects`, `touches`, `crosses`, `within`, `contains`, `overlaps`, `covers`, `coveredBy` |+| Construction | `distance`, `intersection`, `union`, `difference`, `symmetricDifference`, `buffer`, `bufferWithSegments` |+| Measures | `locateAlong`, `locateBetween` |++Indices start at zero. Member counts include empty members. Accessors return+`Nothing` for an index out of range or an unsupported geometry family.+A geometry that is not a collection is its own sole member.++Planar calculations require finite X and Y. Polygon measurements and binary+spatial operations require valid topology. Lengths use coordinate units; areas+use square units. Longitude and latitude inputs therefore give planar lengths+in degrees. These are not geodesic calculations.++Constructed planar results use XY coordinates. Polygon exteriors run+counterclockwise and holes clockwise. Accessors preserve stored Z/M; measured+queries use M and interpolate Z when present. These rules follow the map+geometry model in OGC Simple Feature Access 1.2.1, section 6.1.2.5.++`intersects` and `disjoint` reject separated envelopes and stop at the first+contact. Point containment uses direct point-location tests. Full relation+matrices, overlays, and validity checks use exact spatial indexes. Point-location+queries reuse ring indexes and aggregated winding counts. These indexes can+still take quadratic time when many segment bounds overlap, even for disjoint+shapes such as slanted interleaved combs. Many intersections also increase the+cost of exact rational arithmetic. Measure the shapes and sizes used by your+application.+For large indexed workloads, use a native library such as+[`geos`](https://hackage.haskell.org/package/geos).++Overlays and buffers return `Either SF.TopologyException Geometry` and validate+rounded output. If rounding changes topology, they retry with bounded snapping;+thin regions can collapse. Exhausted precision retries return+`Left SF.PrecisionFailure`. See the+[precision policy](https://github.com/Tritlo/geometry-simple/blob/main/docs/GEOS-DIFFERENCES.md#overlay-precision).++```haskell+clippedArea :: Geometry -> Geometry -> Either SF.TopologyException Double+clippedArea a b = SF.area <$> SF.intersection a b+```++See [Simple Features and GEOS](https://github.com/Tritlo/geometry-simple/blob/main/docs/GEOS-DIFFERENCES.md) for numerical limits,+empty-value rules, format conversions, and deliberate differences from GEOS.+The package implements the seven-family core; it does not claim full OGC SFA+conformance or implement SQL, CRS metadata, Triangle, TIN, or PolyhedralSurface.++## Codecs++`encodeWKB` writes little-endian ISO WKB. `decodeWKB` accepts either byte order,+including mixed byte orders in collections. The WKT decoder accepts explicit+or inferred coordinate layouts, lowercase keywords, and scientific notation.+All four codec functions return `Either String` and reject invalid construction.+The decoders also reject trailing input.++Finite ordinates retain their exact bits through WKB and through text produced+by `encodeWKT`, including negative zero and subnormals. Writers can promote+layouts where a format requires one layout, such as WKT multi-geometries and+WKB polygon rings. Missing Z/M values become NaN. A structural round trip is+therefore not guaranteed for every mixed-layout value.++WKT collections use a parent dimension tag when their members share one output+layout. Mixed-layout collections omit it and retain child tags. That form is+an extension accepted by this library and GEOS; DuckDB rejects it.++## Development++Tests, benchmarks, the Shapely driver, and standard-derived fixtures live in+`dev/`. They are excluded from the published library archive. A repository+checkout can run them with:++```sh+cabal build all+cabal test all --test-show-details=direct+```++CI compares results with Shapely/GEOS and runs independent properties and OGC+examples. It covers GHC 9.6.7 through 9.14.1 on Linux and GHC 9.14.1 on macOS.+See [CONTRIBUTING.md](https://github.com/Tritlo/geometry-simple/blob/main/CONTRIBUTING.md)+for the pinned Nix environment, benchmarks, and release checks.++## License++The published library is MIT licensed. See+[LICENSE](https://github.com/Tritlo/geometry-simple/blob/main/LICENSE).+Standard-derived development fixtures retain their separate OGC notices in+`dev/` and are not part of the release package.

@@ -0,0 +1,177 @@+# Simple Features and GEOS++The target is the planar Simple Features core for seven geometry families.+GEOS is an independent comparison reference. It does not define every Haskell+API or output convention in this package.++## Deliberate differences++| Behavior | geometry-simple | GEOS 3.13.1 |+| --- | --- | --- |+| Constructed planar coordinates | Hulls, centroids, representative points, overlays, and buffers use XY, including empty results. | Some operations retain, interpolate, or discard Z/M according to operation-specific rules. |+| Polygon construction | Exterior rings run counterclockwise; holes run clockwise. | Hulls and many constructed polygon exteriors run clockwise. |+| Hull ordering | The first vertex and line endpoints follow lexicographic XY order. | Polygon starts use Y then X; two-point hulls can retain input order. |+| Overlay line components | Join consecutive edges through vertices with two neighbors. | Can retain separate lines at source vertices and intersection nodes. |+| Representative points | Use the first available polygon, interior line vertex, endpoint, or point, in that order. Empty components are skipped. | Selection can depend on centroid distance, interval width, and empty members. |+| Point observers | Preserve the stored layout and every ordinate, including NaN Z/M. | Extracted points can lose dimensions whose ordinate is NaN. |+| Mixed WKT collections | Each child carries its own dimension tag; the parent has none. | A parent tag can conflict with a child and make the writer's output unreadable. |+| WKT numbers | Shortest scientific notation that decodes to the same `Double`. | Decimal formatting differs. |+| Numerical reductions | Compensated centroid sums retain small contributions during cancellation. | Evaluation order and final floating-point digits can differ. |++OGC 06-103r4 section 6.1.2.5 distinguishes observers, which retain stored Z/M,+from map operations, which ignore them when computing new geometry.+Section 6.1.11.1 describes exterior rings as counterclockwise when viewed from+the top, with opposite winding for holes. We use that convention in the XY+plane. Lexicographic start vertices are our deterministic convention; the+standard does not prescribe them.++The SQL test annex also permits either polygon winding when checking its+expected answers. A different winding alone does not establish a topology bug.++Section 6.1.2.4 defines overlays by their point sets. It does not require a+particular grouping into line components. Sections 6.1.10.2 and 6.1.13.2 require+`PointOnSurface` to return a point on the surface. They do not prescribe which+point to select. The Haskell method also supports points, curves, and mixed+collections, preferring a nonempty component of the highest dimension.++Sources: [OGC common architecture](https://docs.ogc.org/is/06-103r4/06-103r4.pdf),+[OGC SQL test examples](https://docs.ogc.org/is/06-104r4/06-104r4.pdf).++## Geometry contracts++- Point, line, and polygon empties retain their family's topological dimension.+  A collection with no atomic members has dimension -1.+- Coordinate dimensions combine stored member metadata. XYZ and XYM members+  together have coordinate dimension 3, while both Z and M flags are present.+- Empty multi-geometries and empty collections have no stored coordinate layout.+  Writers inherit the containing collection's output tag, or use XY at the top+  level. An empty container does not force a collection to use mixed layouts.+- Polygon area subtracts holes regardless of input winding. Area and perimeter+  close rings supplied directly through Haskell constructors.+- `geometryLength` includes lines and polygon boundaries. `curveLength` counts+  only lines; `perimeter` counts only polygon rings.+- A centroid weights polygons by area, then falls back to segment lengths, then+  points. Each collapsed line or ring contributes its first coordinate once.+  Lower-dimensional components do not affect a higher-dimensional centroid.+- Envelopes use the bounding rectangle specified in section 6.1.2.2. Empty+  input gives an empty point; one XY location gives a point. Flat bounds retain+  repeated rectangle corners and can form a topologically invalid polygon.+- Line boundaries use the mod-2 endpoint rule. `boundary` returns `Nothing` for+  a geometry collection. Stored polygon rings retain their coordinate layouts.+- `relate` returns a nine-character DE-9IM matrix. `relatePattern` accepts+  `T`, `F`, `*`, `0`, `1`, and `2`; invalid patterns return `False`.+  Distance to an empty geometry is NaN.+- Buffers have round caps and joins, with eight segments per quadrant by+  default. Negative distances erode polygons. For points and lines, they give+  an empty polygon.+  Zero distance extracts polygonal regions. Invalid input can lose regions,+  such as one lobe of a self-crossing bowtie. Circular arcs are approximations.+- Overlay line results join through vertices with exactly two neighbors.+  They stop at endpoints and branches. The point set is preserved, but line+  component counts, component order, and polygon hole order can differ from GEOS.+- `pointOnSurface` returns an XY point on a nonempty component. Polygons use+  a horizontal interior interval. Lines prefer stored interior vertices over+  endpoints. A collection of empty components gives an empty XY point.+- Measured queries select points and curve portions using M. Polygon queries+  select their boundary positions, as permitted by the implementation-defined+  surface rule. Empty input gives `Nothing`; no match gives an empty point.++## Numerical limits++Planar operations require finite XY coordinates. Measurements use `Double`.+Polygon cross products use coordinates relative to a ring vertex. Centroid+sums keep polygon positions separate from local moments and compensate for+rounding during addition. Products can still overflow or lose precision.++Centroids can overflow for polygon spans above about 1e100 units or line spans+above about 1e150 units. Underflow can occur below about 1e-100 and 1e-150,+respectively. Compensation cannot recover rounding already lost in products.++Hull orientation tests use a bounded `Double` calculation with an exact+`Rational` fallback. Topology uses exact rational intersections and rounds+constructed output coordinates to `Double`. Exact bounding-box indexes prune+segment pairs, validity checks, point locations, and distance candidates.+Representative points use exact scanline crossings. Their rounded X ordinate+must stay within the selected interval. A shell vertex is used when no interval+contains a representable point. Distance retains exact scaling until the final+conversion, including subnormal projected distances.++Buffers use scaled norms when intermediate products overflow or underflow.+An offset outside the finite `Double` range returns `Left CoordinateOverflow`.+These rules can differ from GEOS at extreme coordinate scales.++Repeated winding queries use aggregated crossing counts. Arrangements in which+many segment bounds overlap can still take quadratic time, even when the shapes+do not intersect, such as slanted interleaved combs. Predicates can reject incompatible dimensions+or bounds before constructing a full relation matrix.++### Overlay precision++Overlays and buffers first compute an exact arrangement and round its result+to `Double`. Buffer offsets approximate circular arcs with floating-point coordinates.+Rounding can move all vertices of a small ring onto one line. The operation+removes such a ring and does not retry because of it.+A valid result keeps its coordinates. If rounding makes it invalid, the operation+retries with vertex and segment snapping. Separate groups of overlapping input+bounds use separate tolerances. A distant component therefore does not set the+precision of a local operation. Positive buffers expand these bounds by their+distance before grouping, because offsets from disjoint inputs can overlap.++The first tolerance is the group's largest absolute ordinate divided by 10^12,+with a floor of the smallest positive `Double`. Five attempts increase this+tolerance tenfold, each starting from the original inputs. Nearby vertices and+intersection nodes share coordinates. Edges also snap to nearby vertices.+Narrow regions can collapse. These functions return+`Either TopologyException Geometry`. If all precision attempts fail, they return+`Left PrecisionFailure`. Offsets beyond the finite `Double` range return+`Left CoordinateOverflow`. `Left OpenBoundary` and `Left UncontainedHole` indicate+a library defect. `Data.Geometry.SimpleFeatures` exports these+constructors. Applications can handle failures without catching exceptions.+The combined output is checked again after separate groups are processed.++The retry count and tolerance schedule follow+[GEOS OverlayNGRobust](https://github.com/libgeos/geos/blob/3.13.1/include/geos/operation/overlayng/OverlayNGRobust.h).+This implementation retains exact noding and deterministic XY representatives.+Buffers use this same retry schedule on their source coordinates; this differs+from GEOS buffer-specific precision reduction.+It does not reproduce GEOS's additional self-union and precision-grid attempts.+It can therefore report a precision failure before GEOS exhausts its retry+strategies. The two implementations do not have the same failure behavior.+Valid results can retain tiny and thin regions that GEOS discards. At subnormal scales,+GEOS's floating-point validity checks can also disagree with exact orientation.+The precision tests include an independent rational check for such a region.++## Format rules++Untagged WKT infers XY, XYZ, or XYZM from two, three, or four ordinates. XYM+requires M. Untagged collection children infer their layouts independently.+In multi-geometries, empties before the first coordinate remain XY; later+empties use the inferred layout. Explicit parent tags require matching child+layouts. Writers tag collections when all members share one output layout,+including nested collections. Empty containers inherit the containing layout;+standalone empty containers use XY.++Mixed-layout collection WKT omits the parent tag and retains each child's tag.+This extends the OGC grammar in section 7. GEOS accepts this form, but DuckDB+requires one layout across all members for both WKT and WKB. Use a common+layout for DuckDB interchange.++The WKT decoder accepts attached tags such as `POINTZ`, both multipoint+syntaxes, signed numbers, fractions, and exponents. Whitespace between ordinates+is required: space, tab, CR, or LF. The decoder rounds each decimal ordinate+to the nearest `Double`. Underflow produces signed zero; overflow produces infinity.++Both codecs accept nonfinite ordinates for storage. NaN in both WKB point XY+ordinates denotes an empty point. Nonfinite values do not have the finite-bit+round-trip guarantee. Use `EmptyPoint` for an empty point in Haskell.++Lines must have zero or at least two coordinates. Rings must have zero or at+least three coordinates and close in XY. An empty exterior cannot contain a+nonempty hole. Writers normalize polygons with only empty rings to an empty+polygon, so empty-hole counts can change after writing. These construction+checks do not establish valid topology.++WKB rejects unknown type tags, incorrect child families, and counts exceeding+the remaining input before allocation. Counts must fit in 32 bits. The codecs+have no nesting limit, and WKT has no count limit. EWKB, EWKT, and embedded SRIDs+are outside the package's scope.

@@ -0,0 +1,70 @@+cabal-version: 3.4+name: geometry-simple+version: 0.1.0.0+license: MIT+license-file: LICENSE+author: Matthias Pall Gissurarson+maintainer: mpg@mpg.is+category: Data, Geometry+build-type: Simple+tested-with:+  ghc ==9.10.3+  ghc ==9.12.4+  ghc ==9.14.1+  ghc ==9.6.7+  ghc ==9.8.4++homepage: https://github.com/Tritlo/geometry-simple+bug-reports: https://github.com/Tritlo/geometry-simple/issues+synopsis: Simple Features geometries, codecs, and pure planar operations+description:+  Geometry values for the seven OGC Simple Features families: points,+  linestrings, polygons, their multi-geometry forms, and geometry collections.+  Coordinates are XY, XYZ, XYM, or XYZM, and coordinate sequences are stored+  in unboxed vectors.++  "Data.Geometry.WKB" and "Data.Geometry.WKT" read and write ISO WKB and WKT.+  The codecs check input structure, coordinate dimensions, line lengths,+  and ring closure. "Data.Geometry.SimpleFeatures" has accessors, planar+  measurements, spatial predicates, validity checks, overlays, round buffers,+  and measured-location queries.++  The package needs no native library or database.++extra-doc-files:+  CHANGELOG.md+  docs/GEOS-DIFFERENCES.md+  README.md++source-repository head+  type: git+  location: https://github.com/Tritlo/geometry-simple.git++library+  exposed-modules:+    Data.Geometry+    Data.Geometry.Internal+    Data.Geometry.SimpleFeatures+    Data.Geometry.WKB+    Data.Geometry.WKT++  hs-source-dirs: src+  other-modules:+    Data.Geometry.Topology.Buffer+    Data.Geometry.Topology.Measures+    Data.Geometry.Topology.Overlay+    Data.Geometry.Topology.Planar+    Data.Geometry.Topology.Relations+    Data.Geometry.Topology.Snapping+    Data.Geometry.Topology.Unary++  default-language: Haskell2010+  build-depends:+    base >=4.18 && <5,+    binary >=0.8.4 && <0.9,+    bytestring >=0.11.2 && <0.13,+    containers >=0.6 && <0.9,+    deepseq >=1.4.8 && <1.6,+    text >=2.0 && <2.2,+    transformers >=0.6 && <0.7,+    vector >=0.13 && <0.14,

@@ -0,0 +1,39 @@+{- | Simple Features geometry values with XY, XYZ, XYM, or XYZM coordinates.+Each point and coordinate sequence has its own layout. Collection members+and polygon rings can use different layouts.++@+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++square :: Geometry+square = Polygon (PolygonRings+    (CoordinatesXY (U.fromList [XY 0 0, XY 1 0, XY 1 1, XY 0 1, XY 0 0]))+    V.empty)+@++A geometry does not store a coordinate reference system. Keep the CRS or SRID+next to the value.++"Data.Geometry.WKB" and "Data.Geometry.WKT" convert geometries to and from+ISO WKB and WKT. "Data.Geometry.SimpleFeatures" provides accessors,+measurements, spatial predicates, and planar geometry operations.+"Data.Geometry.Internal" exposes the 'Coordinate' methods. That module can+change between releases without following the package versioning policy.+-}+module Data.Geometry (+    XY (..),+    XYZ (..),+    XYM (..),+    XYZM (..),+    Dimensions (..),+    Coordinate,+    Point (..),+    Coordinates (..),+    PolygonRings (..),+    Geometry (..),+    withPoint,+    withCoordinates,+) where++import Data.Geometry.Internal

@@ -0,0 +1,472 @@+{-# LANGUAGE DerivingVia #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# OPTIONS_HADDOCK not-home #-}++{- | Geometry types, coordinate conversion, and shared codec validation.++This module is internal. It does not follow the PVP, and any release can+change it. Import "Data.Geometry" and the codec modules for a stable API.+-}+module Data.Geometry.Internal where++import Control.DeepSeq (NFData (..), rwhnf)+import Control.Monad (unless, when)+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Unboxed as U+import Data.Word (Word8)++-- | A coordinate in the XY plane.+data XY+    = -- | @XY x y@.+      XY !Double !Double+    deriving (Eq, Show, Read)++-- | A coordinate with X, Y, and elevation Z.+data XYZ+    = -- | @XYZ x y z@.+      XYZ !Double !Double !Double+    deriving (Eq, Show, Read)++-- | A coordinate with X, Y, and a measure M.+data XYM+    = -- | @XYM x y m@.+      XYM !Double !Double !Double+    deriving (Eq, Show, Read)++-- | A coordinate with X, Y, elevation Z, and a measure M.+data XYZM+    = -- | @XYZM x y z m@.+      XYZM !Double !Double !Double !Double+    deriving (Eq, Show, Read)++-- | The ordinates stored by a point or coordinate sequence.+data Dimensions+    = -- | X and Y.+      DimXY+    | -- | X, Y, and Z.+      DimXYZ+    | -- | X, Y, and M.+      DimXYM+    | -- | X, Y, Z, and M.+      DimXYZM+    deriving (Eq, Ord, Show, Read, Enum, Bounded)++-- | The dimension of a point set, ordered from empty collections to surfaces.+data TopologicalDimension+    = -- | A collection with no atomic members.+      NoDimension+    | -- | A point or multipoint, including an empty value.+      PointDimension+    | -- | A line or multiline, including an empty value.+      CurveDimension+    | -- | A polygon or multipolygon, including an empty value.+      SurfaceDimension+    deriving (Eq, Ord, Show, Read)++{- | The coordinate types t'XY', t'XYZ', t'XYM', and t'XYZM'. Other instances+are not supported.+-}+class (Eq c, Show c, Read c, NFData c, U.Unbox c) => Coordinate c where+    -- | The ordinates available in this coordinate type.+    coordinateDimensions :: proxy c -> Dimensions++    -- | The X, Y, Z, and M ordinates. An ordinate that the type does not have is zero.+    coordinateComponents :: c -> (Double, Double, Double, Double)++    -- | Construct a coordinate from X, Y, Z, and M. Ignore ordinates absent from the type.+    coordinateFromComponents :: (Double, Double, Double, Double) -> c++instance Coordinate XY where+    coordinateDimensions _ = DimXY+    coordinateComponents (XY x y) = (x, y, 0, 0)+    coordinateFromComponents (x, y, _, _) = XY x y++instance Coordinate XYZ where+    coordinateDimensions _ = DimXYZ+    coordinateComponents (XYZ x y z) = (x, y, z, 0)+    coordinateFromComponents (x, y, z, _) = XYZ x y z++instance Coordinate XYM where+    coordinateDimensions _ = DimXYM+    coordinateComponents (XYM x y m) = (x, y, 0, m)+    coordinateFromComponents (x, y, _, m) = XYM x y m++instance Coordinate XYZM where+    coordinateDimensions _ = DimXYZM+    coordinateComponents (XYZM x y z m) = (x, y, z, m)+    coordinateFromComponents (x, y, z, m) = XYZM x y z m++-- | A point with its own coordinate layout. Empty points retain their layout.+data Point+    = -- | No coordinate, with an explicit layout for serialization.+      EmptyPoint !Dimensions+    | -- | A point in the XY plane.+      PointXY !XY+    | -- | A point with elevation Z.+      PointXYZ !XYZ+    | -- | A point with measure M.+      PointXYM !XYM+    | -- | A point with elevation Z and measure M.+      PointXYZM !XYZM+    deriving (Eq, Show, Read)++-- | An unboxed coordinate sequence. Empty sequences retain their layout.+data Coordinates+    = -- | An XY sequence, including an empty XY sequence.+      CoordinatesXY !(U.Vector XY)+    | -- | An XYZ sequence, including an empty XYZ sequence.+      CoordinatesXYZ !(U.Vector XYZ)+    | -- | An XYM sequence, including an empty XYM sequence.+      CoordinatesXYM !(U.Vector XYM)+    | -- | An XYZM sequence, including an empty XYZM sequence.+      CoordinatesXYZM !(U.Vector XYZM)+    deriving (Eq, Show, Read)++-- | An exterior ring and its holes. Each ring has its own coordinate layout.+data PolygonRings+    = -- | @PolygonRings shell holes@. Ring orientation does not affect area.+      PolygonRings+        -- | Exterior ring. An empty polygon has an empty exterior.+        !Coordinates+        -- | Interior rings, in stored order.+        !(V.Vector Coordinates)+    deriving (Eq, Show, Read)++{- | The seven Simple Features geometry families. Each point, line, and ring+retains its coordinate layout. Collection members can have different layouts.+The constructors do not check minimum lengths, ring closure, or topology.+-}+data Geometry+    = -- | One point, which may be empty.+      PointGeometry !Point+    | -- | A sequence joined by straight segments.+      LineString !Coordinates+    | -- | An exterior ring and any holes.+      Polygon !PolygonRings+    | -- | Points with independent layouts and empty values.+      MultiPoint !(U.Vector Point)+    | -- | Line strings with independent coordinate layouts.+      MultiLineString !(V.Vector Coordinates)+    | -- | Polygons with independent ring layouts.+      MultiPolygon !(V.Vector PolygonRings)+    | -- | Geometries of any family, including nested collections.+      GeometryCollection !(V.Vector Geometry)+    deriving (Eq, Show, Read)++-- Coordinates have strict fields, so weak head normal form is normal form.+instance NFData XY where rnf = rwhnf+instance NFData XYZ where rnf = rwhnf+instance NFData XYM where rnf = rwhnf+instance NFData XYZM where rnf = rwhnf++instance NFData Dimensions where rnf = rwhnf++-- Points and unboxed sequences have strict fields with no lazy parts.+instance NFData Point where rnf = rwhnf+instance NFData Coordinates where rnf = rwhnf++instance NFData PolygonRings where+    rnf (PolygonRings shell holes) = rnf shell `seq` rnf holes++-- Boxed vectors of rings and members can contain unevaluated elements.+instance NFData Geometry where+    rnf geometry = case geometry of+        PointGeometry point -> rnf point+        LineString points -> rnf points+        Polygon rings -> rnf rings+        MultiPoint points -> rnf points+        MultiLineString lineStrings -> rnf lineStrings+        MultiPolygon polygons -> rnf polygons+        GeometryCollection children -> rnf children++-- | Apply an operation to the typed buffer of a coordinate sequence.+withCoordinates :: (forall c. (Coordinate c) => U.Vector c -> a) -> Coordinates -> a+withCoordinates f coordinates = case coordinates of+    CoordinatesXY values -> f values+    CoordinatesXYZ values -> f values+    CoordinatesXYM values -> f values+    CoordinatesXYZM values -> f values++-- | Apply an operation to a nonempty point's coordinate.+withPoint :: (forall c. (Coordinate c) => c -> a) -> Point -> Maybe a+withPoint f point = case point of+    EmptyPoint _ -> Nothing+    PointXY value -> Just (f value)+    PointXYZ value -> Just (f value)+    PointXYM value -> Just (f value)+    PointXYZM value -> Just (f value)++-- | The layout stored by a coordinate sequence.+dimensionsOf :: Coordinates -> Dimensions+dimensionsOf (CoordinatesXY _) = DimXY+dimensionsOf (CoordinatesXYZ _) = DimXYZ+dimensionsOf (CoordinatesXYM _) = DimXYM+dimensionsOf (CoordinatesXYZM _) = DimXYZM++-- | The layout stored by a point, including an empty point.+pointDimensions :: Point -> Dimensions+pointDimensions (EmptyPoint dimensions) = dimensions+pointDimensions (PointXY _) = DimXY+pointDimensions (PointXYZ _) = DimXYZ+pointDimensions (PointXYM _) = DimXYM+pointDimensions (PointXYZM _) = DimXYZM++-- | Construct a nonempty point. Ignore tuple fields absent from the chosen layout.+{-# INLINE pointFromComponents #-}+pointFromComponents :: Dimensions -> (Double, Double, Double, Double) -> Point+pointFromComponents dimensions values = case dimensions of+    DimXY -> PointXY (coordinateFromComponents values)+    DimXYZ -> PointXYZ (coordinateFromComponents values)+    DimXYM -> PointXYM (coordinateFromComponents values)+    DimXYZM -> PointXYZM (coordinateFromComponents values)++-- | Combine Z and M flags independently.+unionDimensions :: Dimensions -> Dimensions -> Dimensions+unionDimensions DimXY b = b+unionDimensions a DimXY = a+unionDimensions a b+    | a == b = a+    | otherwise = DimXYZM++-- | The number of ordinates stored by a layout.+dimensionCount :: Dimensions -> Int+dimensionCount DimXY = 2+dimensionCount DimXYZM = 4+dimensionCount _ = 3++-- | The union of the layouts stored by all polygon rings.+polygonDimensions :: PolygonRings -> Dimensions+polygonDimensions (PolygonRings shell holes) = V.foldl' (\acc ring -> unionDimensions acc (dimensionsOf ring)) (dimensionsOf shell) holes++-- | The union of all stored Z and M flags. Empty collections report XY.+geometryDimensions :: Geometry -> Dimensions+geometryDimensions geometry = case geometry of+    PointGeometry point -> pointDimensions point+    LineString points -> dimensionsOf points+    Polygon rings -> polygonDimensions rings+    MultiPoint points -> U.foldl' (\acc point -> unionDimensions acc (pointDimensions point)) DimXY points+    MultiLineString lineStrings -> V.foldl' (\acc points -> unionDimensions acc (dimensionsOf points)) DimXY lineStrings+    MultiPolygon polygons -> V.foldl' (\acc rings -> unionDimensions acc (polygonDimensions rings)) DimXY polygons+    GeometryCollection children -> V.foldl' (\acc child -> unionDimensions acc (geometryDimensions child)) DimXY children++{- | The topological dimension of the geometry family.+Empty values keep their family's dimension. Collections use the greatest+member dimension, or 'NoDimension' when they have no atomic members.+-}+topologicalDimension :: Geometry -> TopologicalDimension+topologicalDimension geometry = case geometry of+    PointGeometry _ -> PointDimension+    MultiPoint _ -> PointDimension+    LineString _ -> CurveDimension+    MultiLineString _ -> CurveDimension+    Polygon _ -> SurfaceDimension+    MultiPolygon _ -> SurfaceDimension+    GeometryCollection children -> V.foldl' (\n child -> max n (topologicalDimension child)) NoDimension children++-- | The greatest coordinate count among members. XYZ and XYM together give 3.+geometryCoordinateDimension :: Geometry -> Int+geometryCoordinateDimension geometry = case geometry of+    PointGeometry point -> dimensionCount (pointDimensions point)+    LineString points -> dimensionCount (dimensionsOf points)+    Polygon (PolygonRings shell holes) -> V.foldl' (\acc ring -> max acc (dimensionCount (dimensionsOf ring))) (dimensionCount (dimensionsOf shell)) holes+    MultiPoint points -> U.foldl' (\acc point -> max acc (dimensionCount (pointDimensions point))) 2 points+    MultiLineString lineStrings -> V.foldl' (\acc points -> max acc (dimensionCount (dimensionsOf points))) 2 lineStrings+    MultiPolygon polygons -> V.foldl' (\acc rings -> max acc (geometryCoordinateDimension (Polygon rings))) 2 polygons+    GeometryCollection children -> V.foldl' (\acc child -> max acc (geometryCoordinateDimension child)) 2 children++-- | Construct an empty sequence with the given layout.+emptyCoordinates :: Dimensions -> Coordinates+emptyCoordinates DimXY = CoordinatesXY U.empty+emptyCoordinates DimXYZ = CoordinatesXYZ U.empty+emptyCoordinates DimXYM = CoordinatesXYM U.empty+emptyCoordinates DimXYZM = CoordinatesXYZM U.empty++-- | Whether a sequence contains no coordinates.+coordinatesEmpty :: Coordinates -> Bool+coordinatesEmpty = withCoordinates U.null++-- | Test whether an ordinate is neither NaN nor infinity.+finite :: Double -> Bool+{-# INLINE finite #-}+finite value = not (isNaN value || isInfinite value)++{- | Check line lengths, ring closure, and polygon emptiness.+Apply the supplied count check to each sequence and collection before checking+its contents. WKB uses it to enforce its 32-bit count limit.+-}+validateGeometry :: (Int -> Either String ()) -> Geometry -> Either String ()+validateGeometry checkLength geometry = case geometry of+    PointGeometry _ -> pure ()+    LineString points -> checkLineLength points >> validateLine points+    Polygon rings@(PolygonRings shell holes) -> do+        checkLength (1 + V.length holes)+        checkLineLength shell+        V.mapM_ checkLineLength holes+        validatePolygon rings+    MultiPoint points -> checkLength (U.length points)+    MultiLineString lineStrings -> do+        checkLength (V.length lineStrings)+        V.mapM_ (validateGeometry checkLength . LineString) lineStrings+    MultiPolygon polygons -> do+        checkLength (V.length polygons)+        V.mapM_ (validateGeometry checkLength . Polygon) polygons+    GeometryCollection children -> do+        checkLength (V.length children)+        V.mapM_ (validateGeometry checkLength) children+  where+    checkLineLength = checkLength . withCoordinates U.length++-- | A line has zero coordinates or at least two coordinates.+validateLine :: Coordinates -> Either String ()+validateLine points = when (withCoordinates U.length points == 1) (Left "Geometry line must have zero or at least two coordinates")++-- | Rings have zero or at least three coordinates and close in X and Y.+validatePolygon :: PolygonRings -> Either String ()+validatePolygon (PolygonRings shell holes) = do+    validateRing shell+    V.mapM_ validateRing holes+    when (coordinatesEmpty shell && not (V.all coordinatesEmpty holes)) $+        Left "Geometry polygon has an empty shell and nonempty holes"+  where+    validateRing = withCoordinates $ \points -> case U.length points of+        0 -> pure ()+        count | count < 3 -> Left "Geometry ring must have zero or at least three coordinates"+        _ ->+            let (x, y, _, _) = coordinateComponents (U.head points)+                (x', y', _, _) = coordinateComponents (U.last points)+             in unless (x == x' && y == y') (Left "Geometry ring is not closed")++-- | Test whether a geometry has no stored coordinates.+geometryEmpty :: Geometry -> Bool+geometryEmpty geometry = case geometry of+    PointGeometry point -> emptyPoint point+    LineString points -> coordinatesEmpty points+    Polygon (PolygonRings shell holes) -> coordinatesEmpty shell && V.all coordinatesEmpty holes+    MultiPoint points -> U.all emptyPoint points+    MultiLineString lineStrings -> V.all coordinatesEmpty lineStrings+    MultiPolygon polygons -> V.all (geometryEmpty . Polygon) polygons+    GeometryCollection children -> V.all geometryEmpty children+  where+    emptyPoint (EmptyPoint _) = True+    emptyPoint _ = False++-- Unboxed vectors store coordinates as tuples, so each ordinate has its own buffer.++instance U.IsoUnbox XY (Double, Double) where+    toURepr (XY x y) = (x, y)+    {-# INLINE toURepr #-}+    fromURepr (x, y) = XY x y+    {-# INLINE fromURepr #-}++-- | Mutable unboxed storage for XY values.+newtype instance U.MVector s XY+    = -- | Wrap the unboxed buffers for X and Y.+      MVXY (U.MVector s (Double, Double))++-- | Immutable unboxed storage for XY values.+newtype instance U.Vector XY+    = -- | Wrap the unboxed buffers for X and Y.+      VXY (U.Vector (Double, Double))++deriving via (U.As XY (Double, Double)) instance M.MVector U.MVector XY+deriving via (U.As XY (Double, Double)) instance G.Vector U.Vector XY+instance U.Unbox XY++instance U.IsoUnbox XYZ (Double, Double, Double) where+    toURepr (XYZ x y z) = (x, y, z)+    {-# INLINE toURepr #-}+    fromURepr (x, y, z) = XYZ x y z+    {-# INLINE fromURepr #-}++-- | Mutable unboxed storage for XYZ values.+newtype instance U.MVector s XYZ+    = -- | Wrap the unboxed buffers for X, Y, and Z.+      MVXYZ (U.MVector s (Double, Double, Double))++-- | Immutable unboxed storage for XYZ values.+newtype instance U.Vector XYZ+    = -- | Wrap the unboxed buffers for X, Y, and Z.+      VXYZ (U.Vector (Double, Double, Double))++deriving via (U.As XYZ (Double, Double, Double)) instance M.MVector U.MVector XYZ+deriving via (U.As XYZ (Double, Double, Double)) instance G.Vector U.Vector XYZ+instance U.Unbox XYZ++instance U.IsoUnbox XYM (Double, Double, Double) where+    toURepr (XYM x y m) = (x, y, m)+    {-# INLINE toURepr #-}+    fromURepr (x, y, m) = XYM x y m+    {-# INLINE fromURepr #-}++-- | Mutable unboxed storage for XYM values.+newtype instance U.MVector s XYM+    = -- | Wrap the unboxed buffers for X, Y, and M.+      MVXYM (U.MVector s (Double, Double, Double))++-- | Immutable unboxed storage for XYM values.+newtype instance U.Vector XYM+    = -- | Wrap the unboxed buffers for X, Y, and M.+      VXYM (U.Vector (Double, Double, Double))++deriving via (U.As XYM (Double, Double, Double)) instance M.MVector U.MVector XYM+deriving via (U.As XYM (Double, Double, Double)) instance G.Vector U.Vector XYM+instance U.Unbox XYM++instance U.IsoUnbox XYZM (Double, Double, Double, Double) where+    toURepr (XYZM x y z m) = (x, y, z, m)+    {-# INLINE toURepr #-}+    fromURepr (x, y, z, m) = XYZM x y z m+    {-# INLINE fromURepr #-}++-- | Mutable unboxed storage for XYZM values.+newtype instance U.MVector s XYZM+    = -- | Wrap the unboxed buffers for X, Y, Z, and M.+      MVXYZM (U.MVector s (Double, Double, Double, Double))++-- | Immutable unboxed storage for XYZM values.+newtype instance U.Vector XYZM+    = -- | Wrap the unboxed buffers for X, Y, Z, and M.+      VXYZM (U.Vector (Double, Double, Double, Double))++deriving via (U.As XYZM (Double, Double, Double, Double)) instance M.MVector U.MVector XYZM+deriving via (U.As XYZM (Double, Double, Double, Double)) instance G.Vector U.Vector XYZM+instance U.Unbox XYZM++-- One tag buffer stores the layout and presence. Unused ordinates are zero.+instance U.IsoUnbox Point (Word8, Double, Double, Double, Double) where+    toURepr (EmptyPoint dimensions) = (fromIntegral (fromEnum dimensions), 0, 0, 0, 0)+    toURepr (PointXY (XY x y)) = (4, x, y, 0, 0)+    toURepr (PointXYZ (XYZ x y z)) = (5, x, y, z, 0)+    toURepr (PointXYM (XYM x y m)) = (6, x, y, 0, m)+    toURepr (PointXYZM (XYZM x y z m)) = (7, x, y, z, m)+    {-# INLINE toURepr #-}+    fromURepr (tag, x, y, z, m) = case tag of+        0 -> EmptyPoint DimXY+        1 -> EmptyPoint DimXYZ+        2 -> EmptyPoint DimXYM+        3 -> EmptyPoint DimXYZM+        4 -> PointXY (XY x y)+        5 -> PointXYZ (XYZ x y z)+        6 -> PointXYM (XYM x y m)+        _ -> PointXYZM (XYZM x y z m)+    {-# INLINE fromURepr #-}++-- | Mutable unboxed storage for Point values.+newtype instance U.MVector s Point+    = -- | Wrap the unboxed buffers for point tags and X, Y, Z, and M.+      MVPoint (U.MVector s (Word8, Double, Double, Double, Double))++-- | Immutable unboxed storage for Point values.+newtype instance U.Vector Point+    = -- | Wrap the unboxed buffers for point tags and X, Y, Z, and M.+      VPoint (U.Vector (Word8, Double, Double, Double, Double))++deriving via (U.As Point (Word8, Double, Double, Double, Double)) instance M.MVector U.MVector Point+deriving via (U.As Point (Word8, Double, Double, Double, Double)) instance G.Vector U.Vector Point+instance U.Unbox Point

@@ -0,0 +1,564 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}++{- | Simple Features accessors, planar topology, and measurements.++Most names follow OGC Simple Feature Access. Indices start at zero, as in+GEOS. The module exports short names such as 'x' and 'area', so import it+qualified:++> import qualified Data.Geometry.SimpleFeatures as SF++Measurements use only X and Y. Z and M stay available through the accessors.+Constructed planar results use t'XY' coordinates. Observers retain stored layouts.++Planar operations require finite X and Y. Measurements use 'Double'+arithmetic and can overflow or underflow. Polygon measurements and binary+spatial operations assume valid topology; use 'isValid' to check it.+Area and perimeter close open rings. Topology uses exact rational segment+intersections and rounds constructed coordinates to 'Double'.++Overlays and buffers return @Either TopologyException Geometry@. If rounding+makes a result invalid, they retry with bounded snapping, and thin regions can+collapse. 'Left' 'PrecisionFailure' means that every retry failed; see+'intersection' for the tolerance schedule. Buffers approximate circular arcs+with straight segments.++The module covers the seven Simple Features geometry families. The surface+types and reference systems of OGC SFA are outside its scope.+-}+module Data.Geometry.SimpleFeatures (+    -- * Geometry properties+    geometryType,+    dimension,+    coordinateDimension,+    spatialDimension,+    is3D,+    isMeasured,+    isEmpty,++    -- * Coordinate ordinates+    x,+    y,+    z,+    m,+    pointX,+    pointY,+    pointZ,+    pointM,++    -- * Members, points, and rings+    numGeometries,+    geometryN,+    numPoints,+    pointN,+    startPoint,+    endPoint,+    isClosed,+    exteriorRing,+    numInteriorRings,+    interiorRingN,++    -- * Measurements and bounds+    envelope,+    area,+    geometryLength,+    curveLength,+    perimeter,+    centroid,+    convexHull,++    -- * Topology+    boundary,+    isSimple,+    isRing,+    isValid,+    pointOnSurface,++    -- * Spatial relations+    relate,+    relatePattern,+    equals,+    disjoint,+    intersects,+    touches,+    crosses,+    within,+    contains,+    overlaps,+    covers,+    coveredBy,++    -- * Distance and geometry construction+    TopologyException (..),+    distance,+    intersection,+    union,+    difference,+    symmetricDifference,+    buffer,+    bufferWithSegments,++    -- * Measured locations+    locateAlong,+    locateBetween,+) where++import Control.Applicative ((<|>))+import Control.Monad (join)+import Data.Geometry.Internal+import Data.Geometry.Topology.Buffer (buffer, bufferWithSegments)+import Data.Geometry.Topology.Measures (locateAlong, locateBetween)+import Data.Geometry.Topology.Overlay (TopologyException (..), difference, intersection, symmetricDifference, union)+import Data.Geometry.Topology.Relations (contains, coveredBy, covers, crosses, disjoint, distance, equals, intersects, overlaps, relate, relatePattern, touches, within)+import Data.Geometry.Topology.Unary (boundary, isRing, isSimple, isValid, pointOnSurface)+import qualified Data.List as List+import Data.Maybe (fromMaybe)+import Data.Proxy (Proxy (..))+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++-- | The uppercase WKT family name, such as @POLYGON@, without a dimension tag.+geometryType :: Geometry -> String+geometryType geometry = case geometry of+    PointGeometry _ -> "POINT"+    LineString _ -> "LINESTRING"+    Polygon _ -> "POLYGON"+    MultiPoint _ -> "MULTIPOINT"+    MultiLineString _ -> "MULTILINESTRING"+    MultiPolygon _ -> "MULTIPOLYGON"+    GeometryCollection _ -> "GEOMETRYCOLLECTION"++{- | The topological dimension: 0 for points, 1 for lines, and 2 for polygons.+Empty values keep their family's dimension. A collection has the largest+dimension of its atomic members, or -1 when it has none.+-}+dimension :: Geometry -> Int+dimension geometry = case topologicalDimension geometry of+    NoDimension -> -1+    PointDimension -> 0+    CurveDimension -> 1+    SurfaceDimension -> 2++{- | The largest number of ordinates in any stored point or sequence: 2, 3,+or 4. XYZ and XYM members together give 3, even though both+'is3D' and 'isMeasured' are true. Atomic empty geometries retain their layout.+Collections without atomic members report 2.+-}+coordinateDimension :: Geometry -> Int+coordinateDimension = geometryCoordinateDimension++-- | The number of spatial ordinates: 3 when 'is3D' is true, or 2 otherwise.+spatialDimension :: Geometry -> Int+spatialDimension geometry = if is3D geometry then 3 else 2++-- | Whether the geometry's coordinate layout has Z. See 'coordinateDimension'.+is3D :: Geometry -> Bool+is3D geometry = geometryDimensions geometry `elem` [DimXYZ, DimXYZM]++-- | Whether the geometry's coordinate layout has M. See 'coordinateDimension'.+isMeasured :: Geometry -> Bool+isMeasured geometry = geometryDimensions geometry `elem` [DimXYM, DimXYZM]++-- | Whether the geometry has no coordinates. A collection of empty members is empty.+isEmpty :: Geometry -> Bool+isEmpty = geometryEmpty++-- | The X ordinate.+x :: (Coordinate c) => c -> Double+x coordinate = let (value, _, _, _) = coordinateComponents coordinate in value++-- | The Y ordinate.+y :: (Coordinate c) => c -> Double+y coordinate = let (_, value, _, _) = coordinateComponents coordinate in value++-- | The Z ordinate, or 'Nothing' when the coordinate type has no Z.+z :: forall c. (Coordinate c) => c -> Maybe Double+z coordinate = case coordinateDimensions (Proxy :: Proxy c) of+    DimXYZ -> Just value+    DimXYZM -> Just value+    _ -> Nothing+  where+    (_, _, value, _) = coordinateComponents coordinate++-- | The M ordinate, or 'Nothing' when the coordinate type has no M.+m :: forall c. (Coordinate c) => c -> Maybe Double+m coordinate = case coordinateDimensions (Proxy :: Proxy c) of+    DimXYM -> Just value+    DimXYZM -> Just value+    _ -> Nothing+  where+    (_, _, _, value) = coordinateComponents coordinate++-- | The X ordinate of a point, or 'Nothing' for an empty point.+pointX :: Point -> Maybe Double+pointX = withPoint x++-- | The Y ordinate of a point, or 'Nothing' for an empty point.+pointY :: Point -> Maybe Double+pointY = withPoint y++-- | The Z ordinate, or 'Nothing' for an empty point or a layout without Z.+pointZ :: Point -> Maybe Double+pointZ = join . withPoint z++-- | The M ordinate, or 'Nothing' for an empty point or a layout without M.+pointM :: Point -> Maybe Double+pointM = join . withPoint m++{- | The number of direct members of a multi-geometry or collection, including+empty members. Other geometries count as one member, also when empty.+-}+numGeometries :: Geometry -> Int+numGeometries geometry = case geometry of+    MultiPoint points -> U.length points+    MultiLineString lineStrings -> V.length lineStrings+    MultiPolygon polygons -> V.length polygons+    GeometryCollection children -> V.length children+    _ -> 1++{- | The direct member at a zero-based index, or 'Nothing' when the index is out+of range. A geometry that is not a collection is its own first member.+-}+geometryN :: Int -> Geometry -> Maybe Geometry+geometryN index geometry = case geometry of+    MultiPoint points -> PointGeometry <$> points U.!? index+    MultiLineString lineStrings -> LineString <$> lineStrings V.!? index+    MultiPolygon polygons -> Polygon <$> polygons V.!? index+    GeometryCollection children -> children V.!? index+    _ -> if index == 0 then Just geometry else Nothing++-- | The number of coordinates in a 'LineString'. Other families give 'Nothing'.+numPoints :: Geometry -> Maybe Int+numPoints (LineString points) = Just (withCoordinates U.length points)+numPoints _ = Nothing++{- | The 'LineString' point at a zero-based index. Return 'Nothing' for an+index out of range or another geometry family. Preserve the stored layout+and all ordinates, including NaN Z or M values.+-}+pointN :: Int -> Geometry -> Maybe Point+pointN index (LineString points) = withCoordinates (\values -> coordinatePoint <$> values U.!? index) points+pointN _ _ = Nothing++-- | Copy a coordinate into a point with the same layout and ordinates.+coordinatePoint :: forall c. (Coordinate c) => c -> Point+coordinatePoint coordinate = pointFromComponents (coordinateDimensions (Proxy :: Proxy c)) (coordinateComponents coordinate)++-- | The first point of a 'LineString', or 'Nothing' for an empty line or another family.+startPoint :: Geometry -> Maybe Point+startPoint = pointN 0++-- | The last point of a 'LineString', or 'Nothing' for an empty line or another family.+endPoint :: Geometry -> Maybe Point+endPoint geometry@(LineString points) = pointN (withCoordinates U.length points - 1) geometry+endPoint _ = Nothing++{- | Whether a nonempty 'LineString' starts and ends at the same XY position.+A 'MultiLineString' is closed when it has lines and all of them are closed.+Other families give 'False'. The test does not check whether a line is simple.+-}+isClosed :: Geometry -> Bool+isClosed geometry = case geometry of+    LineString points -> withCoordinates closed points+    MultiLineString lineStrings -> not (V.null lineStrings) && V.all (withCoordinates closed) lineStrings+    _ -> False+  where+    closed points = not (U.null points) && xy (U.head points) == xy (U.last points)++-- | The exterior ring of a 'Polygon', including its empty layout. Other families give 'Nothing'.+exteriorRing :: Geometry -> Maybe Coordinates+exteriorRing (Polygon (PolygonRings shell _)) = Just shell+exteriorRing _ = Nothing++-- | The number of holes in a 'Polygon'. Other families give 'Nothing'.+numInteriorRings :: Geometry -> Maybe Int+numInteriorRings (Polygon (PolygonRings _ holes)) = Just (V.length holes)+numInteriorRings _ = Nothing++{- | The 'Polygon' hole at a zero-based index. Index 0 is the first hole.+Return 'Nothing' for an index out of range or another geometry family.+-}+interiorRingN :: Int -> Geometry -> Maybe Coordinates+interiorRingN index (Polygon (PolygonRings _ holes)) = holes V.!? index+interiorRingN _ _ = Nothing++{- | The smallest XY bounding rectangle, as a counterclockwise 'Polygon'.+Empty input gives an empty point. A single XY location gives a point.+Horizontal and vertical bounds give a polygon with repeated corners.+-}+envelope :: Geometry -> Geometry+envelope geometry+    | minX > maxX = PointGeometry (EmptyPoint DimXY)+    | minX == maxX && minY == maxY = PointGeometry (PointXY (XY minX minY))+    | otherwise = Polygon (PolygonRings (CoordinatesXY (U.fromList [XY minX minY, XY maxX minY, XY maxX maxY, XY minX maxY, XY minX minY])) V.empty)+  where+    -- An inverted infinite box remains inverted when there are no coordinates.+    (minX, minY, maxX, maxY) = foldCoordinates extend (infinity, infinity, -infinity, -infinity) geometry+    infinity = 1 / 0+    extend (!left, !bottom, !right, !top) coordinate =+        (min left (x coordinate), min bottom (y coordinate), max right (x coordinate), max top (y coordinate))++{- | The total polygon area in square coordinate units. The first ring of each+polygon is the exterior, and the other rings are holes. Ring orientation does+not matter. Other families add zero. The cross products of each ring use its+first vertex as the origin, which limits cancellation far from zero.+-}+area :: Geometry -> Double+area geometry = let (weight, _, _) = surfaceMoments (0, 0) geometry in weight / 2++-- | The total XY length of lines and polygon boundaries, including holes.+geometryLength :: Geometry -> Double+geometryLength geometry = curveLength geometry + perimeter geometry++{- | The total length of all lines, including lines in collections, in coordinate+units. Polygon boundaries and points add zero.+-}+curveLength :: Geometry -> Double+curveLength geometry = case geometry of+    LineString points -> withCoordinates (pathLength False) points+    MultiLineString lineStrings -> V.foldl' (\total points -> total + withCoordinates (pathLength False) points) 0 lineStrings+    GeometryCollection children -> V.foldl' (\total child -> total + curveLength child) 0 children+    _ -> 0++-- | The total length of all polygon rings, including holes. Lines and points add zero.+perimeter :: Geometry -> Double+perimeter geometry = case geometry of+    Polygon (PolygonRings shell holes) -> V.foldl' (\total points -> total + withCoordinates (pathLength True) points) (withCoordinates (pathLength True) shell) holes+    MultiPolygon polygons -> V.foldl' (\total rings -> total + perimeter (Polygon rings)) 0 polygons+    GeometryCollection children -> V.foldl' (\total child -> total + perimeter child) 0 children+    _ -> 0++{- | The XY centroid, or an empty XY point for empty input. Polygons are+weighted by area. If the total area is zero, segments are weighted by length. If all+segments have zero length, the result is the mean of the points, and each+line or ring counts as one point at its first coordinate, as in GEOS.+Lower-dimensional parts do not affect a higher-dimensional centroid.+The centroid can be outside the geometry, for example in a hole.++Compensated sums retain small contributions when large moments cancel.+Polygon positions and local moments stay separate until summation. Moments+can overflow when a polygon spans more than about 1e100 units or a line more+than about 1e150 units. They can underflow below about 1e-100 units for polygons+or 1e-150 units for lines. Products still round to 'Double'. Strong cancellation+between products can reduce accuracy even when all intermediate values are finite.+-}+centroid :: Geometry -> Point+centroid geometry = fromMaybe (EmptyPoint DimXY) (weightedMean SurfaceDimension <|> weightedMean CurveDimension <|> weightedMean PointDimension)+  where+    weightedMean dimensionToMeasure+        | weight == 0 = Nothing+        | finite mx && finite my = Just (PointXY (XY (mx / weight) (my / weight)))+        -- Normalize before summation only when raw moments overflow.+        | otherwise = Just (PointXY (XY (if finite mx then mx / weight else normalizedX) (if finite my then my / weight else normalizedY)))+      where+        (weight, mx, my) = moments dimensionToMeasure 1+        (_, normalizedX, normalizedY) = moments dimensionToMeasure weight+    moments dimensionToMeasure divisor =+        let CentroidMoments w wc mx mxc my myc = foldCentroidMoments dimensionToMeasure divisor (CentroidMoments 0 0 0 0 0 0) geometry+         in (compensatedValue (w, wc), compensatedValue (mx, mxc), compensatedValue (my, myc))++{- | The XY convex hull, using Andrew's monotone chain algorithm.+Return an empty collection, a point, a line, or a counterclockwise polygon.+Line endpoints and the first polygon vertex use lexicographic XY order.+Ignore Z and M. The orientation tests are exact.+-}+convexHull :: Geometry -> Geometry+convexHull geometry = hullGeometry hull+  where+    coordinates = foldCoordinates (\rest coordinate -> xy coordinate : rest) [] geometry+    points = [point | point : _ <- List.group (List.sort coordinates)]+    hull = case points of+        [] -> []+        [_] -> points+        [_, _] -> points+        _ -> init (chain points) ++ init (chain (reverse points))+    chain = reverse . List.foldl' push []+    push (b : a : rest) point+        | orientation a b point /= GT = push (a : rest) point+    push rest point = point : rest++-- | Construct an XY hull from its extreme vertices.+hullGeometry :: [(Double, Double)] -> Geometry+hullGeometry [] = GeometryCollection V.empty+hullGeometry points@(first : _) = case points of+    [(a, b)] -> PointGeometry (PointXY (XY a b))+    [_, _] -> LineString (coordinates points)+    _ -> Polygon (PolygonRings (coordinates (points ++ [first])) V.empty)+  where+    coordinates = CoordinatesXY . U.fromList . map (uncurry XY)++-- | Extract the planar coordinate pair.+xy :: (Coordinate c) => c -> (Double, Double)+xy coordinate = (x coordinate, y coordinate)++-- | Fold coordinates in stored order. Skip explicit empty points.+foldCoordinates :: (forall c. (Coordinate c) => a -> c -> a) -> a -> Geometry -> a+foldCoordinates step initial geometry = case geometry of+    PointGeometry point -> fromMaybe initial (withPoint (step initial) point)+    LineString points -> sequenceFold initial points+    Polygon rings -> polygonFold initial rings+    MultiPoint points -> U.foldl' (\total point -> fromMaybe total (withPoint (step total) point)) initial points+    MultiLineString lineStrings -> V.foldl' sequenceFold initial lineStrings+    MultiPolygon polygons -> V.foldl' polygonFold initial polygons+    GeometryCollection children -> V.foldl' (foldCoordinates step) initial children+  where+    sequenceFold total = withCoordinates (U.foldl' step total)+    polygonFold total (PolygonRings shell holes) = V.foldl' sequenceFold (sequenceFold total shell) holes++-- | Fold adjacent pairs and optionally close the path.+foldSegments :: (U.Unbox c) => Bool -> (a -> c -> c -> a) -> a -> U.Vector c -> a+foldSegments close step initial points+    | U.null points = initial+    | otherwise =+        let result = U.foldl' (\total (a, b) -> step total a b) initial (U.zip points (U.tail points))+         in if close then step result (U.last points) (U.head points) else result++{- | The XY distance. Scale before the square root so that the squares cannot+overflow or underflow.+-}+segmentLength :: (Coordinate c) => c -> c -> Double+segmentLength a b+    | large == 0 = 0+    | otherwise = large * sqrt (1 + ratio * ratio)+  where+    deltaX = abs (x b - x a)+    deltaY = abs (y b - y a)+    large = max deltaX deltaY+    ratio = min deltaX deltaY / large++-- | Sum XY segment lengths. Optionally include the final closing segment.+pathLength :: (Coordinate c) => Bool -> U.Vector c -> Double+pathLength close = foldSegments close (\total a b -> total + segmentLength a b) 0++-- | A weight and the weighted X and Y offsets from an origin.+type Moments = (Double, Double, Double)++-- | Add contributions to one centroid.+addMoments :: Moments -> Moments -> Moments+addMoments (weight, mx, my) (otherWeight, otherX, otherY) =+    let !total = weight + otherWeight+        !totalX = mx + otherX+        !totalY = my + otherY+     in (total, totalX, totalY)++{- | Calculate twice the ring area and its moments about an origin, whatever+the ring winding. The cross products use the first ring vertex as a local+origin, which limits cancellation far from zero.+-}+ringMoments :: (Coordinate c) => (Double, Double) -> U.Vector c -> Moments+ringMoments (originX, originY) ring+    | U.null ring = (0, 0, 0)+    | otherwise =+        let (weight, mx, my) = foldSegments True step (0, 0, 0) ring+            size = abs weight+            direction = signum weight+         in (size, size * (baseX - originX) + direction * mx / 3, size * (baseY - originY) + direction * my / 3)+  where+    (baseX, baseY) = xy (U.head ring)+    step (!weight, !mx, !my) a b =+        let ax = x a - baseX+            ay = y a - baseY+            bx = x b - baseX+            by = y b - baseY+            cross = ax * by - bx * ay+         in (weight + cross, mx + (ax + bx) * cross, my + (ay + by) * cross)++-- | Add exterior ring moments and subtract hole moments, regardless of winding.+surfaceMoments :: (Double, Double) -> Geometry -> Moments+surfaceMoments origin geometry = case geometry of+    Polygon (PolygonRings shell holes) -> V.foldl' subtractRing (withCoordinates (ringMoments origin) shell) holes+    MultiPolygon polygons -> V.foldl' (\total rings -> addMoments total (surfaceMoments origin (Polygon rings))) (0, 0, 0) polygons+    GeometryCollection children -> V.foldl' (\total child -> addMoments total (surfaceMoments origin child)) (0, 0, 0) children+    _ -> (0, 0, 0)+  where+    subtractRing total ring =+        let (weight, mx, my) = withCoordinates (ringMoments origin) ring+         in addMoments total (-weight, -mx, -my)++-- | A sum and the rounding error retained by Neumaier summation.+type Compensated = (Double, Double)++-- | Strict sums and corrections for the weight, X moment, and Y moment.+data CentroidMoments = CentroidMoments !Double !Double !Double !Double !Double !Double++-- | Add one value without discarding small terms when larger terms cancel.+addCompensated :: Compensated -> Double -> Compensated+{-# INLINE addCompensated #-}+addCompensated (total, correction) value =+    let !next = total + value+        !errorTerm = if abs total >= abs value then (total - next) + value else (value - next) + total+        !nextCorrection = correction + errorTerm+     in (next, nextCorrection)++-- | Include the retained rounding error in the result.+compensatedValue :: Compensated -> Double+compensatedValue (total, correction) = total + correction++{- | Fold weights and separate moment contributions for one centroid dimension.+Keep polygon bases separate from local moments, and keep segment endpoints+separate. Rounding their absolute centroids first would discard small offsets.+The divisor scales moments on the overflow retry. Weights stay unscaled.+-}+foldCentroidMoments :: TopologicalDimension -> Double -> CentroidMoments -> Geometry -> CentroidMoments+{-# INLINE foldCentroidMoments #-}+foldCentroidMoments dimensionToMeasure divisor = go+  where+    go initial geometry = case geometry of+        PointGeometry point | dimensionToMeasure == PointDimension -> fromMaybe initial (withPoint (single initial) point)+        LineString points -> case dimensionToMeasure of+            CurveDimension -> withCoordinates (foldSegments False segment initial) points+            PointDimension -> first initial points+            _ -> initial+        Polygon rings@(PolygonRings shell holes) -> case dimensionToMeasure of+            SurfaceDimension -> surface initial rings+            CurveDimension -> V.foldl' path (path initial shell) holes+            _ -> V.foldl' first (first initial shell) holes+        MultiPoint points | dimensionToMeasure == PointDimension -> U.foldl' (\total point -> fromMaybe total (withPoint (single total) point)) initial points+        MultiLineString lineStrings -> V.foldl' (\total points -> go total (LineString points)) initial lineStrings+        MultiPolygon polygons -> V.foldl' (\total rings -> go total (Polygon rings)) initial polygons+        GeometryCollection children -> V.foldl' go initial children+        _ -> initial+    step (CentroidMoments w wc mx mxc my myc) a b c =+        let (w', wc') = addCompensated (w, wc) a+            (mx', mxc') = addCompensated (mx, mxc) b+            (my', myc') = addCompensated (my, myc) c+         in CentroidMoments w' wc' mx' mxc' my' myc'+    single total coordinate = step total 1 (x coordinate / divisor) (y coordinate / divisor)+    -- A zero-length line or ring counts as one point in the point fallback.+    first total = withCoordinates (\points -> maybe total (single total) (points U.!? 0))+    path total = withCoordinates (foldSegments True segment total)+    segment total a b+        | weight == 0 = total+        | otherwise = step (step total weight (factor * x a) (factor * y a)) 0 (factor * x b) (factor * y b)+      where+        weight = segmentLength a b+        factor = (weight / divisor) / 2+    surface total rings@(PolygonRings shell _) = case withCoordinates (\points -> xy <$> points U.!? 0) shell of+        Nothing -> total+        Just (baseX, baseY) ->+            -- Subtract holes locally before multiplying by the polygon position.+            let (weight, mx, my) = surfaceMoments (baseX, baseY) (Polygon rings)+                factor = weight / divisor+             in if weight == 0 then total else step (step total weight (factor * baseX) (factor * baseY)) 0 (mx / divisor) (my / divisor)++{- | The turn of three XY points: 'GT' for counterclockwise, 'LT' for clockwise,+and 'EQ' for collinear. Use the Double determinant when its error bound+decides the sign (Shewchuk 1997). Otherwise use exact Rational arithmetic.+-}+orientation :: (Double, Double) -> (Double, Double) -> (Double, Double) -> Ordering+orientation (ax, ay) (bx, by) (cx, cy)+    | magnitude >= 9.332636185032189e-302 && abs determinant > 4.440892098500626e-16 * magnitude = compare determinant 0+    | otherwise = compare exact 0+  where+    left = (bx - ax) * (cy - ay)+    right = (by - ay) * (cx - ax)+    determinant = left - right+    -- The bound 2^-51 * magnitude exceeds the rounding error of the+    -- determinant. The limit 2^-1000 also covers products that underflow.+    -- Overflow gives NaN or infinity, which fail both tests.+    magnitude = abs left + abs right+    exact =+        (toRational bx - toRational ax) * (toRational cy - toRational ay)+            - (toRational by - toRational ay) * (toRational cx - toRational ax)

@@ -0,0 +1,296 @@+{- | Round planar buffers with polygonal approximations of circular arcs.++The arc and corner rules follow GEOS 3.13.1. See+<https://github.com/libgeos/geos/blob/3.13.1/src/operation/buffer/OffsetSegmentGenerator.cpp OffsetSegmentGenerator>.+The relative distance thresholds below match that implementation.+-}+module Data.Geometry.Topology.Buffer (buffer, bufferWithSegments) where++import Data.Geometry.Internal+import Data.Geometry.Topology.Overlay+import Data.Geometry.Topology.Planar+import Data.List (group, groupBy)+import Data.Maybe (catMaybes)+import qualified Data.Vector as V++{- | Buffer by a distance in coordinate units, with eight segments per quadrant.+Positive distances expand geometry. Negative distances erode polygons and give+empty polygons for points and lines. All results use XY coordinates.+A zero distance extracts polygonal regions. Invalid input can lose regions,+such as one lobe of a self-crossing bowtie. It is not a general validity repair.+The distance and XY coordinates must be finite. Rounded results use the same+validation and bounded snapping as 'intersection'. Return 'Left' when construction+fails; a 'Right' result has valid topology. Return 'Left' 'CoordinateOverflow'+when an offset coordinate cannot fit in a finite Double.+-}+buffer :: Double -> Geometry -> Either TopologyException Geometry+buffer = bufferWithSegments 8++{- | Construct a round buffer with the requested number of segments per quadrant.+Values below one use one segment. The distance and XY coordinates must be finite.+The distance and coordinate-layout rules are the same as for 'buffer'.+-}+bufferWithSegments :: Int -> Double -> Geometry -> Either TopologyException Geometry+bufferWithSegments quadrants radius geometry =+    robustOperation SurfaceDimension (max 0 (toRational radius)) (\source _ -> bufferPlanar (max 1 quadrants) radius source) (planar geometry) (Planar [] [] [])++-- | Build buffer boundaries from source coordinates at one precision attempt.+bufferPlanar :: Int -> Double -> Planar -> Either TopologyException Geometry+bufferPlanar count radius source = do+    lines' <- if radius > 0 then concat <$> traverse (lineBuffer count width) (planarLines source) else pure []+    rings <- traverse checkedCurve (lines' ++ [circle count width (toDouble p) | radius > 0, p <- planarPoints source])+    offsetRings <- if radius == 0 then pure [] else concat <$> traverse (offsetPolygon count radius) (planarPolygons source)+    let surfaces = Planar [] [] (planarPolygons source)+        bands = Planar [] [] [[map toExact ring] | ring <- rings]+        edges = nodeSegments ((if radius /= 0 then concatMap ringSegments offsetRings else segments surfaces) ++ segments bands) []+        queryEdges = segmentQuery edges+        depths = map prepareDepth (planarPolygons source)+        offsetWinding = prepareWinding (concatMap ringSegments (offsetRings ++ map (map toExact) rings))+        selected p+            | radius == 0 = any (\depth -> depth p > 0) depths+            | otherwise = offsetWinding p < 0+        boundary =+            [ if leftInside then (a, b) else (b, a)+            | edge@(a, b) <- edges+            , let (left, right) = sidePoints queryEdges edge+            , let leftInside = selected left+            , leftInside /= selected right+            ]+    polygons <- polygonize boundary+    pure (assemble SurfaceDimension polygons [] [])+  where+    width = abs radius++-- | Reject nonfinite offset coordinates before conversion to exact positions.+checkedCurve :: [FloatingPosition] -> Either TopologyException [FloatingPosition]+checkedCurve points+    | all (\(x, y) -> finite x && finite y) points = Right points+    | otherwise = Left CoordinateOverflow++-- | Prepare ring winding and orientation for zero-buffer queries.+prepareDepth :: [[Position]] -> Position -> Int+prepareDepth [] = const 0+prepareDepth (shell : holes) = \point -> shellDepth point - sum (map ($ point) holeDepths)+  where+    shellDepth = depth shell+    holeDepths = map depth holes+    depth ring =+        let sign = if ringArea ring > 0 then 1 else -1+            winding = prepareWinding (ringSegments ring)+         in \point -> sign * winding point++{- | Offset shells clockwise and holes counterclockwise, with interior on the right.+GEOS discards inverted erosion curves whose samples all lie within 0.99 of+the requested radius from the source boundary. See+<https://github.com/libgeos/geos/blob/3.13.1/src/operation/buffer/BufferCurveSetBuilder.cpp#L418 BufferCurveSetBuilder>.+-}+offsetPolygon :: Int -> Double -> [[Position]] -> Either TopologyException [[Position]]+offsetPolygon _ _ [] = Right []+offsetPolygon count distance (shell : holes)+    | distance < 0 && eroded shell = Right []+    | otherwise = do+        shellCurve <- offset (distance < 0) shell+        if distance < 0 && inverted shell shellCurve+            then pure []+            else (shellCurve :) . catMaybes <$> traverse holeCurve holes+  where+    radius = abs distance+    holeCurve hole+        | distance >= 0 && eroded hole = Right Nothing+        | otherwise = do+            curve <- offset (distance > 0) hole+            pure (if distance >= 0 && inverted hole curve then Nothing else Just curve)+    offset ccw ring = do+        let flipped = (ringArea ring > 0) /= ccw+            simplified = simplify (if flipped then LT else GT) (radius / 100) (map toDouble ring)+            input = if flipped then reverse simplified else simplified+        curve <- map toExact <$> (offsetRing count radius flipped input >>= checkedCurve)+        pure (if distance < 0 then reverse curve else curve)+    eroded [] = True+    eroded ring+        | length (unique ring) == 3 = toRational radius * sum [vectorMagnitude (subtractPosition b a) | (a, b) <- ringSegments ring] > abs (ringArea ring)+        | otherwise =+            let xs = map fst ring+                ys = map snd ring+             in 2 * toRational radius > min (maximum xs - minimum xs) (maximum ys - minimum ys)+    -- GEOS tests only rings with 4 to 8 coordinates after it removes repeated+    -- points, and skips curves with more than four times as many coordinates.+    -- Without these limits, a radius below the coordinate precision makes every+    -- ring look inverted.+    inverted ring curve =+        coordinateCount > 3 && coordinateCount < 9 && length curve <= 4 * coordinateCount && not (any farEnough (curve ++ map midpoint (ringSegments curve)))+      where+        coordinateCount = length (group ring)+        query = segmentQuery (ringSegments ring)+        reach = toRational (0.99 * radius)+        farEnough point@(x, y) = all ((> reach ^ (2 :: Int)) . squaredLength . segmentOffset point) (query ((x - reach, y - reach), (x + reach, y + reach)))++-- | Build a connected left offset curve before resolving its self-intersections.+offsetRing :: Int -> Double -> Bool -> [FloatingPosition] -> Either TopologyException [FloatingPosition]+offsetRing count radius flipped original = case points of+    [] -> Right []+    [p] -> Right (circle count radius p)+    _ -> close . separate radius . concat <$> traverse (offsetCorner count radius flipped) (zip3 (last points : init points) points (drop 1 points ++ take 1 points))+  where+    uniquePoints = [p | p : _ <- group original]+    points = case uniquePoints of+        first : _ | first == last uniquePoints -> init uniquePoints+        _ -> uniquePoints++{- | Join two left offset segments, rounding convex turns and trimming concave turns.+The 1e-3 relative separation and 80:1 inside-corner connector match GEOS.+-}+offsetCorner :: Int -> Double -> Bool -> (FloatingPosition, FloatingPosition, FloatingPosition) -> Either TopologyException [FloatingPosition]+offsetCorner count radius flipped (a, b, c) = do+    _ <- checkedCurve [end, start]+    case turn of+        LT | separation < radius * 1e-3 -> pure [end]+        LT -> checkedCurve (arc count radius b before after True)+        EQ -> checkedCurve (if (a < b) /= (b < c) then arc count radius b before after (not flipped) else [])+        GT -> do+            _ <- checkedCurve [plus a before, plus c after]+            checkedCurve $ case segmentIntersection (toExact (plus a before), toExact end) (toExact start, toExact (plus c after)) of+                p : _ -> [toDouble p]+                [] -> if separation < radius * 1e-3 then [end] else [end, toward end, toward start, start]+  where+    turn = orientation (toExact a) (toExact b) (toExact c)+    before = normal radius a b+    after = normal radius b c+    end = plus b before+    start = plus b after+    separation = pointDistance end start+    factor = if count >= 8 then 80 else 1+    toward (x, y) = let (u, v) = b in (weighted x u, weighted y v)+    weighted x y =+        let result = (factor * x + y) / (factor + 1)+         in if finite result then result else fromRational ((toRational factor * toRational x + toRational y) / toRational (factor + 1))++-- | Construct one connected outline for an open line, or two offsets for a ring.+lineBuffer :: Int -> Double -> [Position] -> Either TopologyException [[FloatingPosition]]+lineBuffer count radius original = case points of+    [] -> Right []+    [point] -> Right [circle count radius point]+    first : _ : _+        | first == last points -> map (map toDouble) <$> offsetPolygon count radius [map toExact points, map toExact points]+        | otherwise -> do+            left <- leftSide leftPath+            right <- leftSide rightPath+            outline <- checkedCurve (left ++ cap leftPath ++ right ++ cap rightPath)+            pure [close (separate radius outline)]+      where+        leftPath = simplify GT (radius / 100) points+        rightPath = reverse (simplify LT (radius / 100) points)+        cap path =+            let final = last path+                penultimate = path !! (length path - 2)+                endNormal = normal radius penultimate final+             in arc count radius final endNormal (opposite endNormal) True+  where+    points = [point | point : _ <- group (map toDouble original)]+    leftSide path@(_ : _ : _) = do+        corners <- concat <$> traverse (offsetCorner count radius False) (zip3 path (drop 1 path) (drop 2 path))+        pure (corners ++ [plus (last path) (normal radius (path !! (length path - 2)) (last path))])+    leftSide _ = Right []++-- | Match GEOS's minimum vertex separation of 1e-4 times the buffer radius.+separate :: Double -> [FloatingPosition] -> [FloatingPosition]+separate radius points = [first | first : _ <- groupBy (\a b -> pointDistance a b < radius * 1e-4) points]++-- | The Euclidean distance used by native buffer vertex thresholds.+pointDistance :: FloatingPosition -> FloatingPosition -> Double+pointDistance a@(x, y) b@(u, v)+    | finite squared && not (isDenormalized squared) && (squared > 0 || a == b) = sqrt squared+    | otherwise = vectorLength (subtractPosition (toExact a) (toExact b))+  where+    squared = (x - u) * (x - u) + (y - v) * (y - v)++{- | Remove shallow concave vertices in the same directed passes as GEOS.+GEOS uses a tolerance of radius / 100. See+<https://github.com/libgeos/geos/blob/3.13.1/src/operation/buffer/BufferInputLineSimplifier.cpp BufferInputLineSimplifier>.+-}+simplify :: Ordering -> Double -> [FloatingPosition] -> [FloatingPosition]+simplify direction tolerance points = map (input V.!) (repeatPass [0 .. V.length input - 1])+  where+    input = V.fromList points+    exact index = toExact (input V.! index)+    threshold = toRational tolerance ^ (2 :: Int)+    shallow middle first final = squaredLength (segmentOffset (exact middle) (exact first, exact final)) < threshold+    deletable a b c = orientation (exact a) (exact b) (exact c) == direction && shallow b a c && all (shallow b a) [a, a + max 1 ((c - a) `div` 10) .. c - 1]+    pass (first : rest) = first : scan rest+    pass [] = []+    scan (a : b : c : rest)+        | deletable a b c = a : scan (c : rest)+        | otherwise = a : scan (b : c : rest)+    scan rest = rest+    repeatPass indexes = let next = pass indexes in if next == indexes then indexes else repeatPass next++-- | A floating-point position used to approximate circular arcs.+type FloatingPosition = (Double, Double)++-- | Round an exact position once for metric calculations.+toDouble :: Position -> FloatingPosition+toDouble (x, y) = (fromRational x, fromRational y)++-- | Preserve the exact binary coordinates of an approximate arc.+toExact :: FloatingPosition -> Position+toExact (x, y) = (toRational x, toRational y)++-- | Add a displacement to a position.+plus :: FloatingPosition -> FloatingPosition -> FloatingPosition+plus (x, y) (u, v) = (x + u, y + v)++-- | Negate a displacement.+opposite :: FloatingPosition -> FloatingPosition+opposite (x, y) = (-x, -y)++{- | The left perpendicular displacement at the requested distance.+Keep the usual evaluation order. Use exact scaling when a product overflows+or enters the subnormal range and loses precision.+-}+normal :: Double -> FloatingPosition -> FloatingPosition -> FloatingPosition+normal radius a@(x, y) b@(u, v)+    | not (isDenormalized distance || isDenormalized scaledX || isDenormalized scaledY) && finite nx && finite ny && (nx /= 0 || dy == 0 || radius == 0) && (ny /= 0 || dx == 0 || radius == 0) = (nx, ny)+    | otherwise = toDouble (scalePosition (toRational radius / vectorMagnitude direction) (-ey, ex))+  where+    dx = u - x+    dy = v - y+    distance = pointDistance a b+    scaledX = radius * dx+    scaledY = radius * dy+    nx = -(scaledY / distance)+    ny = scaledX / distance+    direction@(ex, ey) = subtractPosition (toExact b) (toExact a)++-- | Close a nonempty polygon ring.+close :: [a] -> [a]+close [] = []+close points@(first : _) = points ++ [first]++-- | Approximate a clockwise circle with a fixed quadrant resolution.+circle :: Int -> Double -> FloatingPosition -> [FloatingPosition]+circle count radius center = close [plus center (radial radius angle) | i <- [0 .. 4 * count - 1], let angle = negate (fromIntegral i * (2 * pi / fromIntegral (4 * count)))]++{- | Match GEOS snapping of trigonometric values at the coordinate axes.+The 5e-16 threshold comes from+<https://github.com/libgeos/geos/blob/3.13.1/include/geos/algorithm/Angle.h#L232 Angle.sinCosSnap>.+-}+radial :: Double -> Double -> FloatingPosition+radial radius angle = (radius * snap (cos angle), radius * snap (sin angle))+  where+    snap value = if abs value < 5e-16 then 0 else value++-- | Approximate the directed arc between two radial displacement vectors.+arc :: Int -> Double -> FloatingPosition -> FloatingPosition -> FloatingPosition -> Bool -> [FloatingPosition]+arc count radius center first last' clockwise = plus center first : interior ++ [plus center last']+  where+    direction endpoint = let (x, y) = plus center endpoint; (cx, cy) = center in atan2 (y - cy) (x - cx)+    initial = direction first+    end = direction last'+    start+        | clockwise && initial <= end = initial + 2 * pi+        | not clockwise && initial >= end = initial - 2 * pi+        | otherwise = initial+    total = abs (start - end)+    steps = max 1 (floor (total / ((pi / 2) / fromIntegral count) + 0.5) :: Int)+    increment = (if clockwise then -1 else 1) * (total / fromIntegral steps)+    interior = [plus center (radial radius angle) | i <- [1 .. steps - 1], let angle = start + fromIntegral i * increment]

@@ -0,0 +1,145 @@+-- | Select points and curve portions by their M coordinate.+module Data.Geometry.Topology.Measures (locateAlong, locateBetween) where++import Data.Geometry.Internal+import qualified Data.List as List+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++{- | Select the points and curve portions whose measure equals the argument.+Constant-measure curve portions remain curves. See 'locateBetween' for the+empty-input, layout, and polygon rules.+-}+locateAlong :: Double -> Geometry -> Maybe Geometry+locateAlong value = locateBetween value value++{- | Select the inclusive interval between two M values. Interpolate X, Y,+and Z linearly within each segment. Retain the source coordinate layouts.+Consecutive selected portions stay connected within each input curve.++Return 'Nothing' for an empty input. Return an empty Point for no match or a+reversed interval. An input without M returns an empty XY Point.+Point-only results are MultiPoints. Curve-only results are MultiLineStrings.+A mixed result contains a MultiLineString and a MultiPoint.++For polygons, apply the operation to each boundary ring. This is the+implementation-defined surface rule permitted by OGC SFA 1.2.1, 6.1.2.6.5.+Do not interpolate measures across a polygon interior or between members.+Segments with nonfinite M values do not interpolate. Their matching endpoints+can still contribute points. A NaN interval bound matches nothing.+-}++{- HLINT ignore locateBetween "Use >" -}+locateBetween :: Double -> Double -> Geometry -> Maybe Geometry+locateBetween lower upper geometry+    | geometryEmpty geometry = Nothing+    | not (measured layout) = Just (PointGeometry (EmptyPoint DimXY))+    -- A NaN bound fails this comparison, so it selects nothing.+    | not (lower <= upper) = Just empty+    | otherwise = Just result+  where+    layout = geometryDimensions geometry+    empty = PointGeometry (EmptyPoint layout)+    (points, curves) = select lower upper geometry+    result = case (points, curves) of+        ([], []) -> empty+        (_, []) -> MultiPoint (U.fromList points)+        ([], _) -> MultiLineString (V.fromList curves)+        _ -> GeometryCollection (V.fromList [MultiLineString (V.fromList curves), MultiPoint (U.fromList points)])++-- | Whether a coordinate layout includes M.+measured :: Dimensions -> Bool+measured layout = layout == DimXYM || layout == DimXYZM++-- | Select each atomic member without joining distinct curves.+select :: Double -> Double -> Geometry -> ([Point], [Coordinates])+select lower upper geometry = case geometry of+    PointGeometry point -> ([point | selectedPoint point], [])+    LineString coordinates -> selectCoordinates lower upper coordinates+    Polygon (PolygonRings shell holes) -> selectCoordinates lower upper shell <> foldMap (selectCoordinates lower upper) holes+    MultiPoint points -> (filter selectedPoint (U.toList points), [])+    MultiLineString curves -> foldMap (selectCoordinates lower upper) curves+    MultiPolygon polygons -> foldMap (select lower upper . Polygon) polygons+    GeometryCollection children -> foldMap (select lower upper) children+  where+    selectedPoint point = measured (pointDimensions point) && maybe False (inRange lower upper . measure) (withPoint coordinateComponents point)++-- | Read the M ordinate from the common coordinate representation.+measure :: (Double, Double, Double, Double) -> Double+measure (_, _, _, value) = value++-- | Test the closed measure interval. NaN values do not match.+inRange :: Double -> Double -> Double -> Bool+inRange lower upper value = lower <= value && value <= upper++-- | Clip a measured sequence and retain its concrete coordinate layout.+selectCoordinates :: Double -> Double -> Coordinates -> ([Point], [Coordinates])+selectCoordinates lower upper coordinates = case coordinates of+    CoordinatesXYM values -> build PointXYM CoordinatesXYM values+    CoordinatesXYZM values -> build PointXYZM CoordinatesXYZM values+    _ -> ([], [])+  where+    build wrapPoint wrapCurve values =+        let runs = selectedRuns lower upper (U.toList values)+         in ( [wrapPoint value | [value] <- runs]+            , [wrapCurve (U.fromList run) | run@(_ : _ : _) <- runs]+            )++-- | Join consecutive selected pieces. An excluded segment separates runs.+selectedRuns :: (Coordinate c) => Double -> Double -> [c] -> [[c]]+selectedRuns _ _ [] = []+selectedRuns lower upper [value] = [[value] | inRange lower upper (measure (coordinateComponents value))]+selectedRuns lower upper values@(start : _) = joinEnds (reverse (finish (List.foldl' append ([], []) pieces)))+  where+    pieces = concat [clipSegment lower upper a b | (a, b) <- zip values (drop 1 values)]+    finish ([], completed) = completed+    finish (current, completed) = reverse current : completed+    append state [] = ([], finish state)+    append ([], completed) piece = (reverse piece, completed)+    append state@(current@(lastValue : _), completed) piece@(firstValue : rest)+        | sameCoordinate lastValue firstValue = (reverse rest ++ current, completed)+        | otherwise = (reverse piece, finish state)+    joinEnds runs@(initial@(first : _) : rest)+        | sameCoordinate start (last values) = case reverse rest of+            final : middle | sameCoordinate first (last final) -> (final ++ drop 1 initial) : reverse middle+            _ -> runs+    joinEnds runs = runs++-- | Clip one segment in measure space and retain its original direction.+clipSegment :: (Coordinate c) => Double -> Double -> c -> c -> [[c]]+clipSegment lower upper a b+    | not (finite start && finite end) = [[a] | inRange lower upper start] ++ [[]] ++ [[b] | inRange lower upper end]+    | start == end && not (inRange lower upper start) = [[]]+    | start == end = [if sameCoordinate a b then [a] else [a, b]]+    | selectedStart > selectedEnd = [[]]+    | selectedStart == selectedEnd = [[at selectedStart]]+    | start < end = [[at selectedStart, at selectedEnd]]+    | otherwise = [[at selectedEnd, at selectedStart]]+  where+    start = measure (coordinateComponents a)+    end = measure (coordinateComponents b)+    selectedStart = max lower (min start end)+    selectedEnd = min upper (max start end)+    at value+        | value == start = a+        | value == end = b+        | otherwise = interpolate ((toRational value - toRational start) / (toRational end - toRational start)) value a b++-- | Interpolate coordinates exactly before their final conversion to Double.+interpolate :: (Coordinate c) => Rational -> Double -> c -> c -> c+interpolate fraction value a b = coordinateFromComponents (along x u, along y v, along z w, value)+  where+    (x, y, z, _) = coordinateComponents a+    (u, v, w, _) = coordinateComponents b+    along first second+        | first == second = first+        | finite first && finite second = fromRational ((1 - fraction) * toRational first + fraction * toRational second)+        | otherwise = (1 - fromRational fraction) * first + fromRational fraction * second++-- | Compare shared endpoints while retaining unknown Z values.+sameCoordinate :: (Coordinate c) => c -> c -> Bool+sameCoordinate a b = equal x u && equal y v && equal z w && equal m n+  where+    (x, y, z, m) = coordinateComponents a+    (u, v, w, n) = coordinateComponents b+    equal first second = first == second || (isNaN first && isNaN second)

@@ -0,0 +1,311 @@+{- | Planar set operations over exact segment arrangements.++Precision retries use five snapping tolerances. The first is the largest+absolute ordinate in a group divided by 10^12. Use at least the smallest+positive Double. Each later attempt multiplies this tolerance by ten. All attempts+start from the original inputs. The schedule follows+<https://github.com/libgeos/geos/blob/3.13.1/include/geos/operation/overlayng/OverlayNGRobust.h GEOS OverlayNGRobust>.+-}+module Data.Geometry.Topology.Overlay where++import Control.DeepSeq (NFData (..), rwhnf)+import Control.Exception (Exception)+import Data.Geometry.Internal+import Data.Geometry.Topology.Planar+import Data.Geometry.Topology.Snapping (snapPlanars)+import Data.Geometry.Topology.Unary (isValid)+import qualified Data.Graph as Graph+import Data.List (maximumBy, minimumBy)+import qualified Data.List as List+import qualified Data.Map.Strict as Map+import Data.Maybe (fromMaybe)+import Data.Ord (comparing)+import qualified Data.Set as Set+import Data.Tree (flatten)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++{- | A topology construction failure returned by overlays and buffers.+It is an 'Exception', so callers can rethrow it, for example with+@either throwIO pure@.+-}+data TopologyException+    = -- | Rounded output remains invalid after all precision retries.+      PrecisionFailure+    | -- | A constructed coordinate exceeds the finite Double range.+      CoordinateOverflow+    | -- | A selected boundary does not close. This is a library defect; please report it.+      OpenBoundary+    | -- | A hole has no containing exterior ring. This is a library defect; please report it.+      UncontainedHole+    deriving (Eq, Show, Read)++instance Exception TopologyException++instance NFData TopologyException where+    rnf = rwhnf++{- | The points common to both geometries. Coordinates must have finite XY values.+Line results join consecutive edges through vertices with exactly two neighbors.+Their component count and order can differ from other implementations.+If rounding to Double makes the result invalid, retry with snapping tolerances+of 10^-12 to 10^-8 times the largest absolute ordinate, as GEOS does. Thin+regions can collapse. Return 'Left' 'PrecisionFailure' if every attempt remains+invalid.+-}+intersection :: Geometry -> Geometry -> Either TopologyException Geometry+intersection a b = overlay (&&) (min (topologicalDimension a) (topologicalDimension b)) a b++-- | The points in either geometry. The precision and failure rules match 'intersection'.+union :: Geometry -> Geometry -> Either TopologyException Geometry+union a b = overlay (||) (max (topologicalDimension a) (topologicalDimension b)) a b++{- | The closure of the points in the first geometry but not the second.+The precision and failure rules match 'intersection'.+-}+difference :: Geometry -> Geometry -> Either TopologyException Geometry+difference a = overlay (\x y -> x && not y) (topologicalDimension a) a++{- | The closure of the points in exactly one geometry.+The precision and failure rules match 'intersection'.+-}+symmetricDifference :: Geometry -> Geometry -> Either TopologyException Geometry+symmetricDifference a b = overlay (/=) (max (topologicalDimension a) (topologicalDimension b)) a b++{- | Try exact overlay first. If output rounding breaks topology, retry each+connected group with five snapping tolerances from magnitude / 10^12 to+magnitude / 10^8. Each attempt starts from the original inputs. These bounds+follow GEOS OverlayNGRobust; separate groups keep their own precision scale.+-}+overlay :: (Bool -> Bool -> Bool) -> TopologicalDimension -> Geometry -> Geometry -> Either TopologyException Geometry+overlay select emptyDimension first second+    | Just unchanged <- trivialOverlay select emptyDimension a b = Right unchanged+    | otherwise = robustOperation emptyDimension 0 (overlayPlanar select emptyDimension) a b+  where+    a = planar first+    b = planar second++{- | Validate construction before and after precision retries. Retry connected+input groups separately so distant components do not set the local tolerance.+Expand group bounds by the positive buffer distance when constructing buffers.+Reuse the first result when all inputs belong to one group.+-}+robustOperation :: TopologicalDimension -> Rational -> (Planar -> Planar -> Either TopologyException Geometry) -> Planar -> Planar -> Either TopologyException Geometry+robustOperation emptyDimension padding operation a b = case operation a b >>= validateResult of+    Right result -> Right result+    Left CoordinateOverflow -> Left CoordinateOverflow+    Left _ -> do+        parts <- case overlayGroups padding a b of+            [_] -> (: []) <$> firstValid (snappedAttempts a b)+            groups -> traverse (\(x, y) -> firstValid (operation x y : snappedAttempts x y)) groups+        let result = combinePlanar (map planar parts)+        validateResult (assemble emptyDimension (planarPolygons result) (planarLines result) (planarPoints result))+  where+    -- The Either semigroup keeps the first valid result, or else the last failure.+    firstValid = foldr1 (<>) . map (>>= validateResult)+    snappedAttempts x y =+        [ let (snappedA, snappedB) = snapPlanars (tolerance * 10 ^ attemptIndex) x y+           in operation snappedA snappedB+        | attemptIndex <- [0 :: Int .. 4]+        ]+      where+        magnitude = maximum (0 : [max (abs u) (abs v) | (u, v) <- allPositions x ++ allPositions y])+        minimumSpacing = toRational (encodeFloat 1 (fst (floatRange (0 :: Double)) - floatDigits (0 :: Double)) :: Double)+        tolerance = max minimumSpacing (magnitude / 10 ^ (12 :: Int))++-- | Reject invalid rounded output without raising an exception from pure code.+validateResult :: Geometry -> Either TopologyException Geometry+validateResult geometry = if isValid geometry then Right geometry else Left PrecisionFailure++-- | Group atomic inputs with overlapping bounds before a precision retry.+overlayGroups :: Rational -> Planar -> Planar -> [(Planar, Planar)]+overlayGroups padding first second = map group (Graph.components graph)+  where+    parts shape = [Planar [p] [] [] | p <- planarPoints shape] ++ [Planar [] [line] [] | line <- planarLines shape] ++ [Planar [] [] [polygon] | polygon <- planarPolygons shape]+    inputs = [(side, part) | (side, shape) <- [(False, first), (True, second)], part <- parts shape, not (null (allPositions part))]+    pairs = overlappingPairs [(expand bounds, i) | (i, (_, part)) <- zip [0 :: Int ..] inputs, Just bounds <- [pointBounds (allPositions part)]]+    expand ((x, y), (u, v)) = ((x - padding, y - padding), (u + padding, v + padding))+    adjacent = Map.fromListWith (++) [(a, [b]) | (i, j) <- pairs, (a, b) <- [(i, j), (j, i)]]+    (graph, entry, _) = Graph.graphFromEdges [(part, i, Map.findWithDefault [] i adjacent) | (i, part) <- zip [0 :: Int ..] inputs]+    group tree =+        let members = [part | vertex <- flatten tree, let (part, _, _) = entry vertex]+         in (combinePlanar [part | (False, part) <- members], combinePlanar [part | (True, part) <- members])++-- | Handle empty or disjoint single components without new intersection coordinates.+trivialOverlay :: (Bool -> Bool -> Bool) -> TopologicalDimension -> Planar -> Planar -> Maybe Geometry+trivialOverlay select emptyDimension a b+    | earlyEmpty = Just (emptyGeometry emptyDimension)+    | separate && singleComponent a && singleComponent b =+        Just (assemble emptyDimension (concatMap planarPolygons retained) (concatMap planarLines retained) (concatMap planarPoints retained))+    | otherwise = Nothing+  where+    separate = disjointBounds (allPositions a) (allPositions b)+    retained = [shape | (keep, shape) <- [(select True False, a), (select False True, b)], keep]+    singleComponent shape = length (componentPoints shape) <= 1+    firstEmpty = null (allPositions a)+    secondEmpty = null (allPositions b)+    earlyEmpty = (firstEmpty && not (select False True)) || (secondEmpty && not (select True False)) || (firstEmpty && secondEmpty) || (separate && not (select True False) && not (select False True))++-- | Evaluate a Boolean set operation on exact faces, edges, and vertices.+overlayPlanar :: (Bool -> Bool -> Bool) -> TopologicalDimension -> Planar -> Planar -> Either TopologyException Geometry+overlayPlanar select emptyDimension a b = do+    polygons <- polygonize boundaryEdges+    let surfaces = Planar [] [] polygons+        (locateSurfaces, _) = prepareLocations surfaces+        lineEdges = [edge | edge <- edges, selected (midpoint edge), locateSurfaces (midpoint edge) == Exterior]+        paths = linePaths lineEdges+        curves = Planar [] paths polygons+        (locateCurves, _) = prepareLocations curves+        nodes = unique (vertices a ++ vertices b ++ concatMap (\(u, v) -> [u, v]) edges)+        points = [p | p <- nodes, selected p, locateCurves p == Exterior]+    pure (assemble emptyDimension polygons paths points)+  where+    edges = nodeSegments (segments a ++ segments b) (vertices a ++ vertices b)+    queryEdges = segmentQuery edges+    (locateA, insideA) = prepareLocations a+    (locateB, insideB) = prepareLocations b+    selected p = select (locateA p /= Exterior) (locateB p /= Exterior)+    selectedFace p = select (insideA p) (insideB p)+    boundaryEdges =+        [ if leftInside then (u, v) else (v, u)+        | edge@(u, v) <- edges+        , let (left, right) = sidePoints queryEdges edge+        , let leftInside = selectedFace left+        , leftInside /= selectedFace right+        ]++-- | Collect directed cycles with their selected region on the left.+boundaryRings :: [Segment] -> Either TopologyException [[Position]]+boundaryRings edges = concatMap splitRing <$> collect (Set.fromList edges)+  where+    outgoingEdges = Map.fromListWith Set.union [(a, Set.singleton (a, b)) | (a, b) <- edges]+    collect remaining = case Set.minView remaining of+        Nothing -> Right []+        Just (edge@(start, _), _) -> do+            (ring, rest) <- walk start edge [] remaining+            (ring :) <$> collect rest+    walk start edge@(a, b) accumulated remaining+        | b == start = Right (reverse (b : a : accumulated), rest)+        | otherwise = do+            next <- case if null preceding then outgoing else preceding of+                [] -> Left OpenBoundary+                candidates -> Right (maximumBy (\(_, x) (_, y) -> compareDirection (subtractPosition x b) (subtractPosition y b)) candidates)+            walk start next (a : accumulated) rest+      where+        rest = Set.delete edge remaining+        outgoing = Set.toList (Set.intersection (Map.findWithDefault Set.empty b outgoingEdges) rest)+        reverseDirection = subtractPosition a b+        preceding = filter (\(_, c) -> compareDirection (subtractPosition c b) reverseDirection == LT) outgoing++-- | Separate rings that touch at one vertex without joining their interiors.+splitRing :: [Position] -> [[Position]]+splitRing = walk Set.empty []+  where+    walk _ _ [] = []+    walk seen path (point : rest)+        | Set.notMember point seen = walk (Set.insert point seen) (point : path) rest+        | otherwise = case break (== point) path of+            (_, []) -> walk (Set.insert point seen) (point : path) rest+            (before, _ : after) -> (point : reverse before ++ [point]) : walk (seen `Set.difference` Set.fromList before) (point : after) rest++-- | Twice the signed area of a closed ring.+ringArea :: [Position] -> Rational+ringArea = sum . map (uncurry cross) . ringSegments++-- | Find a simple ring's orientation at its lexicographically smallest vertex.+ringOrientation :: [Position] -> Ordering+ringOrientation ring = case ringSegments ring of+    [] -> EQ+    edges ->+        let (before, at) = minimumBy (comparing snd) edges+         in case [after | (start, after) <- edges, start == at] of+                after : _ -> orientation before at after+                [] -> EQ++-- | Group each clockwise hole with its innermost containing shell.+polygonize :: [Segment] -> Either TopologyException [[[Position]]]+polygonize edges = do+    rings <- map (\ring -> (ringOrientation ring, ring)) <$> boundaryRings edges+    let shells = [(ring, prepareRing ring) | (GT, ring) <- rings]+        holes = [ring | (LT, ring) <- rings]+    assignments <- traverse (\hole -> do shell <- containingShell shells hole; pure (shell, [hole])) holes+    let groupedHoles = Map.fromListWith (++) assignments+    pure [shell : Map.findWithDefault [] shell groupedHoles | (shell, _) <- shells]+  where+    containingShell shells hole = case ringSegments hole of+        [] -> Left UncontainedHole+        edge : _ -> case filter (\(_, locate) -> locate (midpoint edge) == Interior) shells of+            [] -> Left UncontainedHole+            first : rest -> Right (fst (List.foldl' innermost first rest))+    innermost current@(_, locate) candidate@(ring, _) = case ringSegments ring of+        edge : _ | locate (midpoint edge) == Interior -> candidate+        _ -> current++{- | Join edges through degree-two vertices. Stop at original endpoints and+junctions, even after some of their edges have been removed.+-}+linePaths :: [Segment] -> [[Position]]+linePaths edges = collect (neighbors, starts)+  where+    neighbors = Map.fromListWith Set.union [(a, Set.singleton b) | (p, q) <- edges, (a, b) <- [(p, q), (q, p)]]+    starts = Map.keysSet (Map.filter ((/= 2) . Set.size) neighbors)+    collect remaining@(graph, ends) = case Map.lookupMin graph of+        Nothing -> []+        Just (first, _) ->+            let start = fromMaybe first (Set.lookupMin ends)+                (path, rest) = walk start start [] remaining+             in path : collect rest+    walk start current accumulated remaining@(graph, _) = case Set.lookupMin adjacent of+        Nothing -> (reverse (current : accumulated), remaining)+        _ | not (null accumulated) && (current == start || Set.member current starts) -> (reverse (current : accumulated), remaining)+        Just next -> walk start next (current : accumulated) (removeNeighbor current next (removeNeighbor next current remaining))+      where+        adjacent = Map.findWithDefault Set.empty current graph+    removeNeighbor neighbor point (graph, ends) =+        let adjacent = Set.delete neighbor (Map.findWithDefault Set.empty point graph)+         in if Set.null adjacent+                then (Map.delete point graph, Set.delete point ends)+                else (Map.insert point adjacent graph, ends)++-- | Construct the smallest XY family that holds the selected components.+assemble :: TopologicalDimension -> [[[Position]]] -> [[Position]] -> [Position] -> Geometry+assemble emptyDimension polygons lines' points = case parts of+    [] -> emptyGeometry emptyDimension+    [part] -> part+    _ -> GeometryCollection (V.fromList (concatMap atomic parts))+  where+    parts = pointParts ++ lineParts ++ polygonParts+    roundedLines = [[point | point : _ <- List.group (map rounded line)] | line <- lines']+    pointParts = case [PointXY (XY x y) | (x, y) <- unique (map rounded points ++ [p | [p] <- roundedLines])] of+        [] -> []+        [point] -> [PointGeometry point]+        values -> [MultiPoint (U.fromList values)]+    lineParts = case [CoordinatesXY (U.fromList [XY x y | (x, y) <- line]) | line@(_ : _ : _) <- roundedLines] of+        [] -> []+        [line] -> [LineString line]+        values -> [MultiLineString (V.fromList values)]+    polygonParts = case [PolygonRings (coordinates (oriented GT shell)) (V.fromList (map (coordinates . oriented LT) (filter spansArea holes))) | shell : holes <- map (map roundedRing) polygons, spansArea shell] of+        [] -> []+        [rings] -> [Polygon rings]+        values -> [MultiPolygon (V.fromList values)]+    rounded (x, y) = (fromRational x :: Double, fromRational y :: Double)+    -- Rounding can merge neighbouring vertices. Remove the repeats and drop rings+    -- whose vertices all lie on one line, so one collapsed sliver cannot+    -- invalidate the whole result. Validation rejects other degenerate rings.+    roundedRing ring = [(toRational x, toRational y) | (x, y) : _ <- List.group (map rounded ring)]+    coordinates = CoordinatesXY . U.fromList . map (\(x, y) -> XY (fromRational x) (fromRational y))+    oriented direction ring = if ringOrientation ring == direction then ring else reverse ring++    atomic geometry = case geometry of+        MultiPoint values -> map PointGeometry (U.toList values)+        MultiLineString values -> map LineString (V.toList values)+        MultiPolygon values -> map Polygon (V.toList values)+        _ -> [geometry]++-- | Construct an empty XY geometry of the requested topological dimension.+emptyGeometry :: TopologicalDimension -> Geometry+emptyGeometry dimension = case dimension of+    PointDimension -> PointGeometry (EmptyPoint DimXY)+    CurveDimension -> LineString (CoordinatesXY U.empty)+    SurfaceDimension -> Polygon (PolygonRings (CoordinatesXY U.empty) V.empty)+    NoDimension -> GeometryCollection V.empty

@@ -0,0 +1,421 @@+{-# LANGUAGE BangPatterns #-}++-- | Exact planar primitives shared by the topology operations.+module Data.Geometry.Topology.Planar where++import Data.Geometry.Internal+import Data.List (sortBy, sortOn)+import qualified Data.List as List+import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++-- | An exact position in the XY plane.+type Position = (Rational, Rational)++-- | A closed straight segment.+type Segment = (Position, Position)++-- | A balanced tree of segments with exact bounding boxes on its branches.+data SegmentIndex+    = -- | One segment, which can have coincident endpoints.+      SegmentLeaf !Segment+    | -- | The enclosing box and two nonempty subtrees.+      SegmentBranch !Segment !SegmentIndex !SegmentIndex++-- | Enclose the indexed segments. Leaf endpoints also specify their bounds.+indexBounds :: SegmentIndex -> Segment+indexBounds (SegmentLeaf edge) = edge+indexBounds (SegmentBranch bounds _ _) = bounds++-- | Divide segments at the median, alternating X and Y at each level.+indexSegments :: [Segment] -> Maybe SegmentIndex+indexSegments = build fst snd+  where+    build _ _ [] = Nothing+    build _ _ [edge] = Just (SegmentLeaf edge)+    build coordinate other edges = do+        let (first, second) = splitAt (length edges `div` 2) (sortOn (\(a, b) -> coordinate a + coordinate b) edges)+        left <- build other coordinate first+        right <- build other coordinate second+        let ((ax, ay), (bx, by)) = indexBounds left+            ((cx, cy), (dx, dy)) = indexBounds right+            bounds = ((minimum [ax, bx, cx, dx], minimum [ay, by, cy, dy]), (maximum [ax, bx, cx, dx], maximum [ay, by, cy, dy]))+        pure (SegmentBranch bounds left right)++-- | Build one index and select segments whose closed bounds meet each query box.+segmentQuery :: [Segment] -> Segment -> [Segment]+segmentQuery edges = maybe (const []) querySegments (indexSegments edges)++-- | Select leaves whose closed bounds intersect the supplied box.+querySegments :: SegmentIndex -> Segment -> [Segment]+querySegments tree bounds+    | not (overlapsBounds bounds (indexBounds tree)) = []+    | otherwise = case tree of+        SegmentLeaf edge -> [edge]+        SegmentBranch _ left right -> querySegments left bounds ++ querySegments right bounds++-- | Attach values to indexed bounds, retaining values with equal boxes.+boundsQuery :: [(Segment, a)] -> Segment -> [a]+boundsQuery entries = concatMap (table Map.!) . query+  where+    table = Map.fromListWith (++) [(bounds, [value]) | (bounds, value) <- entries]+    query = segmentQuery (Map.keys table)++-- | Enumerate unordered pairs with intersecting bounds, including equal boxes.+overlappingPairs :: [(Segment, a)] -> [(a, a)]+overlappingPairs entries = [(a, b) | (bounds, (i, a)) <- indexed, (j, b) <- query bounds, i < j]+  where+    indexed = zipWith (\i (bounds, value) -> (bounds, (i, value))) [0 :: Int ..] entries+    query = boundsQuery indexed++-- | Prepare exact winding queries over directed edges, retaining duplicate edges.+prepareWinding :: [Segment] -> Position -> Int+prepareWinding edges = maybe (const 0) windingIndex (indexSegments edges)++{- | Count crossings without visiting branches wholly to the right of the point.+Each branch stores cumulative changes at endpoint Y values. An upward edge+adds one between its endpoints; a downward edge subtracts one. Shared vertices+cancel when branches combine. Branches that contain the query X still use exact+orientation tests at their leaves.+-}+windingIndex :: SegmentIndex -> Position -> Int+windingIndex = snd . build+  where+    build (SegmentLeaf (a@(_, ay), b@(_, by))) =+        ( Map.filter (/= 0) (Map.fromListWith (+) [(ay, 1), (by, -1)])+        , \point@(_, y) ->+            if ay <= y && by > y && orientation a b point == GT+                then 1+                else if by <= y && ay > y && orientation a b point == LT then -1 else 0+        )+    build (SegmentBranch ((ax, ay), (bx, by)) left right) = (changes, classify)+      where+        (leftChanges, leftWinding) = build left+        (rightChanges, rightWinding) = build right+        changes = Map.filter (/= 0) (Map.unionWith (+) leftChanges rightChanges)+        cumulative = snd (Map.mapAccum (\total delta -> let next = total + delta in (next, next)) 0 changes)+        classify point@(x, y)+            | x >= bx || y < ay || y >= by = 0+            | x < ax = maybe 0 snd (Map.lookupLE y cumulative)+            | otherwise = leftWinding point + rightWinding point++-- | A point's location relative to a geometry.+data Location = Exterior | Boundary | Interior deriving (Eq, Ord, Show, Read)++-- | The atomic components of a geometry, projected into the XY plane.+data Planar = Planar+    { planarPoints :: [Position]+    , planarLines :: [[Position]]+    , planarPolygons :: [[[Position]]]+    }+    deriving (Eq, Show, Read)++-- | Convert finite coordinates without rounding their binary values.+positions :: Coordinates -> [Position]+positions = withCoordinates (map position . U.toList)++-- | Project one finite coordinate into the plane.+position :: (Coordinate c) => c -> Position+position coordinate = let (x, y, _, _) = coordinateComponents coordinate in (toRational x, toRational y)++-- | Flatten collections while retaining polygon and line boundaries.+planar :: Geometry -> Planar+planar geometry =+    let Planar points lines' polygons = collect (Planar [] [] []) geometry+     in Planar (reverse points) (reverse lines') (reverse polygons)+  where+    collect rest@(Planar points lines' polygons) shape = case shape of+        PointGeometry point -> maybe rest (\p -> Planar (p : points) lines' polygons) (withPoint position point)+        LineString line -> Planar points (positions line : lines') polygons+        Polygon (PolygonRings shell holes) -> Planar points lines' (map positions (shell : V.toList holes) : polygons)+        MultiPoint values -> U.foldl' (\acc point -> collect acc (PointGeometry point)) rest values+        MultiLineString values -> V.foldl' (\acc line -> collect acc (LineString line)) rest values+        MultiPolygon values -> V.foldl' (\acc rings -> collect acc (Polygon rings)) rest values+        GeometryCollection children -> V.foldl' collect rest children++-- | Concatenate atomic components without changing their coordinates.+combinePlanar :: [Planar] -> Planar+combinePlanar parts = Planar (concatMap planarPoints parts) (concatMap planarLines parts) (concatMap planarPolygons parts)++{- | Find the dimension of a valid point set, ignoring empty components.+Collapsed lines and collinear rings contribute only their stored point set.+-}+planarDimension :: Planar -> TopologicalDimension+planarDimension shape+    | any (any spansArea) (planarPolygons shape) = SurfaceDimension+    | not (null (segments shape)) = CurveDimension+    | not (null (allPositions shape)) = PointDimension+    | otherwise = NoDimension++-- | Test whether the positions do not all lie on one line.+spansArea :: [Position] -> Bool+spansArea (first : rest) = case dropWhile (== first) rest of+    second : remaining -> any ((/= EQ) . orientation first second) remaining+    [] -> False+spansArea [] = False++-- | Remove duplicates and return values in ascending order.+unique :: (Ord a) => [a] -> [a]+unique = Set.toAscList . Set.fromList++-- | List adjacent pairs, excluding segments with zero length.+lineSegments :: [Position] -> [Segment]+lineSegments points = [(a, b) | (a, b) <- zip points (drop 1 points), a /= b]++-- | List the segments of a closed ring.+ringSegments :: [Position] -> [Segment]+ringSegments [] = []+ringSegments points = lineSegments (points ++ take 1 points)++-- | List every line segment and polygon boundary segment.+segments :: Planar -> [Segment]+segments shape = concatMap lineSegments (planarLines shape) ++ concatMap ringSegments (concat (planarPolygons shape))++-- | List distinct input coordinates, including collapsed lines.+vertices :: Planar -> [Position]+vertices = unique . allPositions++-- | Traverse coordinates without sorting or removing duplicates.+allPositions :: Planar -> [Position]+allPositions shape = planarPoints shape ++ concat (planarLines shape) ++ concat (concat (planarPolygons shape))++-- | Enclose a nonempty point set in an exact axis-aligned rectangle.+pointBounds :: [Position] -> Maybe Segment+pointBounds [] = Nothing+pointBounds (first : rest) = Just (List.foldl' extend (first, first) rest)+  where+    extend ((!ax, !ay), (!bx, !by)) (x, y) = ((min ax x, min ay y), (max bx x, max by y))++-- | Test strict separation of the input envelopes. Empty inputs are disjoint.+disjointBounds :: [Position] -> [Position] -> Bool+disjointBounds a b = case (pointBounds a, pointBounds b) of+    (Just first, Just second) -> not (overlapsBounds first second)+    _ -> True++-- | Test whether the first nonempty envelope contains the second one.+enclosesBounds :: [Position] -> [Position] -> Bool+enclosesBounds a b = case (pointBounds a, pointBounds b) of+    (Just ((ax, ay), (bx, by)), Just ((cx, cy), (dx, dy))) -> ax <= cx && ay <= cy && bx >= dx && by >= dy+    _ -> False++-- | Choose one point from each connected component before testing its edges.+componentPoints :: Planar -> [Position]+componentPoints shape =+    planarPoints shape+        ++ [point | point : _ <- planarLines shape]+        ++ [point | (point : _) : _ <- planarPolygons shape]++-- | Test segment contacts using both coordinate bounds and stop at the first hit.+segmentsIntersect :: [Segment] -> [Segment] -> Bool+segmentsIntersect first second = any (\edge -> not (all (null . segmentIntersection edge) (query edge))) first+  where+    query = segmentQuery second++-- | The displacement from the second position to the first.+subtractPosition :: Position -> Position -> Position+subtractPosition (x, y) (u, v) = (x - u, y - v)++-- | Translate a position by an XY displacement.+addPosition :: Position -> Position -> Position+addPosition (x, y) (u, v) = (x + u, y + v)++-- | Scale both XY components by an exact factor.+scalePosition :: Rational -> Position -> Position+scalePosition t (x, y) = (t * x, t * y)++-- | The signed area of the parallelogram spanned by two vectors.+cross :: Position -> Position -> Rational+cross (x, y) (u, v) = x * v - y * u++-- | The turn from the first segment to the second: GT is counterclockwise.+orientation :: Position -> Position -> Position -> Ordering+orientation a b c = compare (cross (subtractPosition b a) (subtractPosition c a)) 0++-- | Test membership of a closed segment with exact arithmetic.+pointOnSegment :: Position -> Segment -> Bool+pointOnSegment point@(x, y) (a@(ax, ay), b@(bx, by)) =+    x >= min ax bx && x <= max ax bx && y >= min ay by && y <= max ay by && orientation a b point == EQ++-- | Return no point, one intersection, or the two ends of an overlap.+segmentIntersection :: Segment -> Segment -> [Position]+segmentIntersection first@(a, b) second@(c, d)+    | not (overlapsBounds first second) = []+    | denominator == 0 = case unique [p | p <- [a, b, c, d], pointOnSegment p first, pointOnSegment p second] of+        [] -> []+        [p] -> [p]+        firstPoint : rest -> [firstPoint, last rest]+    | t >= 0 && t <= 1 && u >= 0 && u <= 1 = [addPosition a (scalePosition t ab)]+    | otherwise = []+  where+    ab = subtractPosition b a+    cd = subtractPosition d c+    ac = subtractPosition c a+    denominator = cross ab cd+    t = cross ac cd / denominator+    u = cross ac ab / denominator++-- | Test whether two segment bounding boxes meet.+overlapsBounds :: Segment -> Segment -> Bool+overlapsBounds ((ax, ay), (bx, by)) ((cx, cy), (dx, dy)) =+    max ax bx >= min cx dx && max cx dx >= min ax bx && max ay by >= min cy dy && max cy dy >= min ay by++-- | Split segments at every intersection and supplied point. Return unique edges.+nodeSegments :: [Segment] -> [Position] -> [Segment]+nodeSegments input points = unique (concatMap split edges)+  where+    edges = unique [if a < b then (a, b) else (b, a) | (a, b) <- input, a /= b]+    query = segmentQuery edges+    -- Segment intersections already account for every endpoint.+    endpoints = Set.fromList (concatMap (\(a, b) -> [a, b]) edges)+    extraPoints = Set.toList (Set.fromList points `Set.difference` endpoints)+    queryPoints = segmentQuery [(p, p) | p <- extraPoints]+    split edge@(a, b) = lineSegments (unique (a : b : [p | (p, _) <- queryPoints edge, pointOnSegment p edge] ++ concatMap (segmentIntersection edge) (query edge)))++-- | The exact midpoint of a segment.+midpoint :: Segment -> Position+midpoint (a, b) = scalePosition (1 / 2) (addPosition a b)++-- | Build one edge index for repeated exact location queries in a ring.+prepareRing :: [Position] -> Position -> Location+prepareRing ring = case indexSegments (ringSegments ring) of+    Nothing -> const Exterior+    Just tree ->+        let winding = windingIndex tree+         in \point ->+                if any (pointOnSegment point) (querySegments tree (point, point))+                    then Boundary+                    else if odd (winding point) then Interior else Exterior++-- | Index the shell and holes once. Skip holes whose bounds exclude the point.+preparePolygon :: [[Position]] -> Position -> Location+preparePolygon [] = const Exterior+preparePolygon [shell] = prepareRing shell+preparePolygon (shell : holes) = classify+  where+    locateShell = prepareRing shell+    queryHoles = boundsQuery [(bounds, prepareRing hole) | hole <- holes, Just bounds <- [pointBounds hole]]+    classify point = case locateShell point of+        Exterior -> Exterior+        shellLocation+            | Interior `elem` holesHere -> Exterior+            | Boundary `elem` holesHere -> Boundary+            | otherwise -> shellLocation+          where+            holesHere = map ($ point) (queryHoles (point, point))++-- | Index polygon components and retain the boundary rules for their union.+prepareSurface :: [[[Position]]] -> Position -> Location+prepareSurface [] = const Exterior+prepareSurface [polygon] = preparePolygon polygon+prepareSurface polygons = classify+  where+    queryPolygons = boundsQuery [(bounds, preparePolygon polygon) | polygon <- polygons, Just bounds <- [pointBounds (concat polygon)]]+    queryEdges = segmentQuery (concatMap ringSegments (concat polygons))+    locations point = map ($ point) (queryPolygons (point, point))+    classify point+        | Interior `elem` here = Interior+        | otherwise = case filter (== Boundary) here of+            _ : _ : _ | all (elem Interior . locations) (sectorSamples queryEdges point) -> Interior+            _ : _ -> Boundary+            [] -> Exterior+      where+        here = locations point++-- | Prepare point locations and surface membership for repeated arrangement queries.+prepareLocations :: Planar -> (Position -> Location, Position -> Bool)+prepareLocations shape = (classify, (== Interior) . locateSurface)+  where+    locateSurface = prepareSurface (planarPolygons shape)+    queryLines = segmentQuery (concatMap lineSegments lines')+    lines' = planarLines shape+    linePoints = Set.fromList (concat lines')+    points = Set.fromList (planarPoints shape)+    endpoints = Map.fromListWith (+) [(p, 1 :: Int) | line@(first : _) <- lines', p <- [first, last line]]+    classify point = case locateSurface point of+        Exterior+            | Set.member point linePoints || any (pointOnSegment point) (queryLines (point, point)) ->+                if odd (Map.findWithDefault 0 point endpoints) then Boundary else Interior+            | Set.member point points -> Interior+            | otherwise -> Exterior+        result -> result++-- | Select a nearby point without crossing any segment after the start point.+nearPoint :: (Segment -> [Segment]) -> Position -> Position -> Position+nearPoint query origin direction = addPosition origin (scalePosition step direction)+  where+    -- Only crossings before t=2 can reduce the initial step of one.+    reach = (origin, addPosition origin (scalePosition 2 direction))+    step = minimum (1 : [t / 2 | edge <- query reach, t <- rayParameters edge, t > 0])+    rayParameters (a, b)+        | determinant /= 0 = [t | u >= 0 && u <= 1]+        | cross offset direction /= 0 = []+        | otherwise = map parameter [a, b]+      where+        edge = subtractPosition b a+        offset = subtractPosition a origin+        determinant = cross direction edge+        t = cross offset edge / determinant+        u = cross offset direction / determinant+    parameter point =+        let (dx, dy) = direction+            (x, y) = subtractPosition point origin+         in if dx /= 0 then x / dx else y / dy++-- | Sample the faces immediately to the left and right of a noded edge.+sidePoints :: (Segment -> [Segment]) -> Segment -> (Position, Position)+sidePoints query segment@(a, b) = (nearPoint query middle normal, nearPoint query middle (scalePosition (-1) normal))+  where+    middle = midpoint segment+    (dx, dy) = subtractPosition b a+    normal = (-dy, dx)++-- | Sort nonzero directions counterclockwise from the positive X axis.+compareDirection :: Position -> Position -> Ordering+compareDirection a@(x, y) b@(u, v) = case compare (half x y) (half u v) of+    EQ -> compare 0 (cross a b)+    result -> result+  where+    half p q = not (q > 0 || (q == 0 && p >= 0))++-- | Sample boundary sectors using an existing edge index.+sectorSamples :: (Segment -> [Segment]) -> Position -> [Position]+sectorSamples query point = [nearPoint query point (addPosition a b) | (a, b) <- zip directions (drop 1 directions ++ take 1 directions)]+  where+    directions = sortBy compareDirection (unique ([(1, 0), (0, 1), (-1, 0), (0, -1)] ++ incident))+    incident = [normalize direction | edge@(a, b) <- query (point, point), pointOnSegment point edge, q <- [a, b], q /= point, let direction = subtractPosition q point]+    normalize (x, y) = let size = abs x + abs y in (x / size, y / size)++-- | Return the vector from the closest point on a segment to a point.+segmentOffset :: Position -> Segment -> Position+segmentOffset point (a, b)+    | a == b = offset+    | otherwise = subtractPosition offset (scalePosition fraction direction)+  where+    offset = subtractPosition point a+    direction = subtractPosition b a+    fraction = max 0 (min 1 (dot offset direction / squaredLength direction))++-- | The exact scalar product of two vectors.+dot :: Position -> Position -> Rational+dot (x, y) (u, v) = x * u + y * v++-- | The exact squared length of a vector.+squaredLength :: Position -> Rational+squaredLength vector = dot vector vector++{- | Scale before the square root, then round the complete length to Double.+Keep the scale exact so subnormal projections do not round to zero too early.+-}+vectorLength :: Position -> Double+vectorLength = fromRational . vectorMagnitude++-- | Approximate a norm with exact scaling and a Double square root.+vectorMagnitude :: Position -> Rational+vectorMagnitude (x, y)+    | scale == 0 = 0+    | otherwise = scale * toRational (sqrt (1 + fromRational (ratio * ratio)) :: Double)+  where+    scale = max (abs x) (abs y)+    ratio = min (abs x) (abs y) / scale

@@ -0,0 +1,260 @@+-- | Planar spatial relations and distance. Coordinates must have finite XY values.+module Data.Geometry.Topology.Relations (+    relate,+    relatePattern,+    equals,+    disjoint,+    intersects,+    touches,+    crosses,+    within,+    contains,+    overlaps,+    covers,+    coveredBy,+    distance,+) where++import Data.Geometry.Internal (Geometry (..), TopologicalDimension (..), withPoint)+import Data.Geometry.Topology.Planar+import qualified Data.List as List+import qualified Data.Set as Set+import qualified Data.Vector.Unboxed as U++{- | Return the nine-character DE-9IM intersection matrix in row-major order.+Rows and columns denote interior, boundary, and exterior. Each character is+@F@ for an empty intersection, or @0@, @1@, or @2@ for its dimension.+Z and M ordinates do not affect the result. Line boundaries use the mod-2 rule.+-}+relate :: Geometry -> Geometry -> String+relate first second = map symbol (U.toList matrix)+  where+    a = planar first+    b = planar second+    inputVertices = vertices a ++ vertices b+    edges = nodeSegments (segments a ++ segments b) inputVertices+    queryEdges = segmentQuery edges+    nodes = unique (inputVertices ++ concatMap (\(p, q) -> [p, q]) edges)+    (locateA, insideA) = prepareLocations a+    (locateB, insideB) = prepareLocations b+    samples =+        [(locateA p, locateB p, 0) | p <- nodes]+            ++ [(locateA p, locateB p, 1) | edge <- edges, let p = midpoint edge]+            ++ [ (areaLocation insideA p, areaLocation insideB p, 2)+               | edge <- edges+               , let (left, right) = sidePoints queryEdges edge+               , p <- [left, right]+               ]+    areaLocation containsPoint p = if containsPoint p then Interior else Exterior+    matrix :: U.Vector Int+    matrix = U.accum max (U.fromList (replicate 8 (-1) ++ [2])) [(locationIndex row * 3 + locationIndex column, dimension) | (row, column, dimension) <- samples]+    symbol dimension = case dimension of+        0 -> '0'+        1 -> '1'+        2 -> '2'+        _ -> 'F'++-- | Convert the location order to the DE-9IM row and column order.+locationIndex :: Location -> Int+locationIndex location = case location of+    Interior -> 0+    Boundary -> 1+    Exterior -> 2++{- | Test a DE-9IM pattern. The pattern must contain nine characters from+@TF*012@. Invalid patterns return 'False'.+-}+relatePattern :: String -> Geometry -> Geometry -> Bool+relatePattern patternText first second = matches patternText (relate first second)++-- | Validate a pattern and match it against a complete intersection matrix.+matches :: String -> String -> Bool+matches patternText matrix = length patternText == 9 && all (`elem` "TF*012") patternText && and (zipWith match patternText matrix)+  where+    match '*' _ = True+    match 'T' actual = actual /= 'F'+    match expected actual = expected == actual++{- | Test whether both geometries contain the same XY point set.+All empty geometries are spatially equal, regardless of their family.+-}+equals :: Geometry -> Geometry -> Bool+equals first second =+    planarDimension a == planarDimension b+        && pointBounds (allPositions a) == pointBounds (allPositions b)+        && (matrix == "FFFFFFFF2" || matches "T*F**FFF*" matrix)+  where+    a = planar first+    b = planar second+    matrix = relate first second++-- | Test whether the geometries have no point in common.+disjoint :: Geometry -> Geometry -> Bool+disjoint first second = not (intersects first second)++{- | Test whether the geometries have at least one point in common.+Reject disjoint envelopes, then test component points and segment contacts.+Stop when an intersection is found without constructing a relation matrix.+-}+intersects :: Geometry -> Geometry -> Bool+intersects first second = planarIntersects (planar first) (planar second)++-- | Reuse projected components when testing contact before a distance search.+planarIntersects :: Planar -> Planar -> Bool+planarIntersects a b =+    not (disjointBounds (allPositions a) (allPositions b))+        && ( any ((/= Exterior) . locateB) (componentPoints a)+                || any ((/= Exterior) . locateA) (componentPoints b)+                || segmentsIntersect (segments a) (segments b)+           )+  where+    (locateA, _) = prepareLocations a+    (locateB, _) = prepareLocations b++-- | Test whether the geometries meet but their interiors do not intersect.+touches :: Geometry -> Geometry -> Bool+touches first second+    | not (planarIntersects a b) = False+    | any insideA (componentPoints b) || any insideB (componentPoints a) = False+    | otherwise = any (`matches` matrix) ["FT*******", "F**T*****", "F***T****"]+  where+    a = planar first+    b = planar second+    matrix = relate first second+    (_, insideA) = prepareLocations a+    (_, insideB) = prepareLocations b++{- | Whether the geometries cross in their interiors.+For different dimensions, the interiors must intersect and the geometry with+the lower dimension must extend outside the other. Two lines cross when+their interiors meet at points. Other equal-dimension pairs return 'False'.+-}+crosses :: Geometry -> Geometry -> Bool+crosses first second+    | dimensionA == dimensionB && dimensionA /= CurveDimension = False+    | disjointBounds (allPositions a) (allPositions b) = False+    | dimensionA < dimensionB = matches "T*T******" matrix+    | dimensionA > dimensionB = matches "T*****T**" matrix+    | otherwise = matches "0********" matrix+  where+    matrix = relate first second+    a = planar first+    b = planar second+    dimensionA = planarDimension a+    dimensionB = planarDimension b++{- | Test whether the first geometry lies in the second geometry and their+interiors intersect. A geometry on only the second boundary is not within it.+-}+within :: Geometry -> Geometry -> Bool+within first second = contains second first++{- | Whether the second geometry lies in the first and their interiors intersect.+Contact confined to the boundary does not count. Use 'covers' to include it.+-}+contains :: Geometry -> Geometry -> Bool+contains first (PointGeometry point) = maybe False ((== Interior) . fst (prepareLocations (planar first))) (withPoint position point)+contains first second =+    planarDimension a >= planarDimension b+        && enclosesBounds (allPositions a) (allPositions b)+        && relatePattern "T*****FF*" first second+  where+    a = planar first+    b = planar second++{- | Test whether geometries of the same dimension share an interior part+of that dimension and each has a part outside the other.+-}+overlaps :: Geometry -> Geometry -> Bool+overlaps first second+    | dimensionA /= dimensionB || dimensionA == NoDimension = False+    | disjointBounds (allPositions a) (allPositions b) = False+    | dimensionA == CurveDimension = matches "1*T***T**" matrix+    | otherwise = matches "T*T***T**" matrix+  where+    matrix = relate first second+    a = planar first+    b = planar second+    dimensionA = planarDimension a+    dimensionB = planarDimension b++{- | Test whether the second geometry has no point outside the first.+Return 'False' if either geometry is empty. Boundary points are included.+-}+covers :: Geometry -> Geometry -> Bool+covers first (PointGeometry point) = maybe False ((/= Exterior) . fst (prepareLocations (planar first))) (withPoint position point)+covers first second =+    planarDimension a >= planarDimension b+        && enclosesBounds (allPositions a) (allPositions b)+        && matches "******FF*" matrix+        && not (matches "FF*FF****" matrix)+  where+    a = planar first+    b = planar second+    matrix = relate first second++-- | Whether the first geometry is covered by the second, including its boundary.+coveredBy :: Geometry -> Geometry -> Bool+coveredBy first second = covers second first++{- | Return the minimum Euclidean distance in the XY plane.+Empty geometries return NaN. Intersecting geometries return zero.+Exact projections avoid overflow and cancellation before the final square root.+-}+distance :: Geometry -> Geometry -> Double+distance first second = case (verticesA, verticesB) of+    (p : _, q : _)+        | planarIntersects a b -> 0+        | otherwise ->+            let offset = subtractPosition p q+                initial = (offset, squaredLength offset)+                closest = case (treeA, treeB) of+                    (Just firstTree, Just secondTree) -> nearestPair firstTree secondTree initial+                    _ -> initial+             in vectorLength (fst closest)+    _ -> 0 / 0+  where+    a = planar first+    b = planar second+    verticesA = vertices a+    verticesB = vertices b+    edgesA = segments a+    edgesB = segments b+    treeA = distanceIndex verticesA edgesA+    treeB = distanceIndex verticesB edgesB++-- | Index segments and retain isolated vertices as zero-length segments.+distanceIndex :: [Position] -> [Segment] -> Maybe SegmentIndex+distanceIndex points edges = indexSegments (edges ++ [(p, p) | p <- points, Set.notMember p endpoints])+  where+    endpoints = Set.fromList (concatMap (\(a, b) -> [a, b]) edges)++-- | Search pairs of subtrees, pruning whole groups by their exact distance bounds.+nearestPair :: SegmentIndex -> SegmentIndex -> (Position, Rational) -> (Position, Rational)+nearestPair first second = visit (bound first second) first second+  where+    bound a b = boxesDistanceSquared (indexBounds a) (indexBounds b)+    visit lowerBound a b best@(_, bestSquared)+        | lowerBound >= bestSquared = best+        | otherwise = case (a, b) of+            (SegmentLeaf firstEdge@(p, q), SegmentLeaf secondEdge@(r, s)) ->+                List.foldl' closer best [segmentOffset p secondEdge, segmentOffset q secondEdge, segmentOffset r firstEdge, segmentOffset s firstEdge]+            (SegmentBranch _ left right, SegmentLeaf _) -> descend left b right b best+            (SegmentLeaf _, SegmentBranch _ left right) -> descend a left a right best+            (SegmentBranch _ leftA rightA, SegmentBranch _ leftB rightB)+                | extent a >= extent b -> descend leftA b rightA b best+                | otherwise -> descend a leftB a rightB best+    descend a b c d best+        | firstBound <= secondBound = visit secondBound c d (visit firstBound a b best)+        | otherwise = visit firstBound a b (visit secondBound c d best)+      where+        firstBound = bound a b+        secondBound = bound c d+    extent tree = let ((x, y), (u, v)) = indexBounds tree in max (abs (u - x)) (abs (v - y))+    closer best@(_, current) candidate = let squared = squaredLength candidate in if squared < current then (candidate, squared) else best++-- | The exact squared distance between two closed axis-aligned bounding boxes.+boxesDistanceSquared :: Segment -> Segment -> Rational+boxesDistanceSquared ((ax, ay), (bx, by)) ((cx, cy), (dx, dy)) = squaredLength (gap ax bx cx dx, gap ay by cy dy)+  where+    gap a b c d = max 0 (max (min a b - max c d) (min c d - max a b))

@@ -0,0 +1,80 @@+{- | Tolerance-based noding for precision retries.++Nearby vertices share one stored position. Segments receive shared intersection+nodes and nearby vertices before exact overlay. The algorithm follows the+vertex and near-segment snapping strategy used by GEOS SnappingNoder.+-}+module Data.Geometry.Topology.Snapping (snapPlanars) where++import Data.Geometry.Internal (TopologicalDimension (..))+import Data.Geometry.Topology.Planar+import Data.List (minimumBy, sortOn)+import qualified Data.List as List+import qualified Data.Map.Strict as Map+import Data.Ord (comparing)++{- | Snap two inputs together within a positive distance tolerance.+Process stored vertices before generated intersections. Select representatives+in XY order so operand order does not determine the coordinate choice.+-}+snapPlanars :: Rational -> Planar -> Planar -> (Planar, Planar)+snapPlanars tolerance first second = (finish a, finish b)+  where+    stored = representatives tolerance (unique (allPositions first ++ allPositions second))+    a = mapPositions (stored Map.!) first+    b = mapPositions (stored Map.!) second+    edges = segments a ++ segments b+    queryEdges = segmentQuery edges+    noded = nodeSegments edges []+    nodes = representatives tolerance (unique (allPositions a ++ allPositions b) ++ unique (concatMap (\(p, q) -> [p, q]) noded))+    candidates = segmentQuery [(p, p) | p <- unique (Map.elems nodes)]+    squaredTolerance = tolerance * tolerance+    finish shape = Planar (planarPoints shape ++ collapsed) curves polygons+      where+        paths = map path (planarLines shape)+        collapsed = [p | [p] <- paths]+        curves = [line | line@(_ : _ : _) <- paths]+        polygons = [shell : filter hasArea holes | rings <- planarPolygons shape, shell : holes <- [map path rings], hasArea shell]+    hasArea ring = planarDimension (Planar [] [] [[ring]]) == SurfaceDimension+    path [] = []+    path values@(start : _) = start : concatMap insert (lineSegments values)+    insert edge@(p@(px, py), q@(qx, qy)) = sortOn parameter (filter between (unique (crossings ++ nearby))) ++ [q]+      where+        direction = subtractPosition q p+        parameter point = dot (subtractPosition point p) direction+        between point = parameter point > 0 && parameter point < squaredLength direction+        crossings = [nodes Map.! point | other <- queryEdges edge, point <- segmentIntersection edge other]+        expanded = ((min px qx - tolerance, min py qy - tolerance), (max px qx + tolerance, max py qy + tolerance))+        nearby =+            [ point+            | (point, _) <- candidates expanded+            , squaredLength (subtractPosition point p) >= squaredTolerance+            , squaredLength (subtractPosition point q) >= squaredTolerance+            , squaredLength (segmentOffset point edge) < squaredTolerance+            ]++-- | Transform every position without changing the component structure.+mapPositions :: (Position -> Position) -> Planar -> Planar+mapPositions f (Planar points lines' polygons) = Planar (map f points) (map (map f) lines') (map (map (map f)) polygons)++{- | Select an existing nearby representative or store a rounded position.+The spatial buckets contain only representatives. A chain of nearby inputs+therefore cannot move an endpoint farther than the tolerance.+-}+representatives :: Rational -> [Position] -> Map.Map Position Position+representatives tolerance = snd . List.foldl' add (Map.empty, Map.empty)+  where+    squaredTolerance = tolerance * tolerance+    bucket (x, y) = (floor (x / tolerance) :: Integer, floor (y / tolerance) :: Integer)+    add state@(buckets, assigned) point+        | Map.member point assigned = state+        | otherwise = case nearby of+            [] ->+                let chosen = rounded point+                 in (Map.insertWith (++) (bucket chosen) [chosen] buckets, Map.insert point chosen assigned)+            _ -> (buckets, Map.insert point (minimumBy (comparing (\p -> (distance p, p))) nearby) assigned)+      where+        (i, j) = bucket point+        distance p = squaredLength (subtractPosition p point)+        nearby = [p | x <- [i - 1 .. i + 1], y <- [j - 1 .. j + 1], p <- Map.findWithDefault [] (x, y) buckets, distance p <= squaredTolerance]+    rounded (x, y) = (toRational (fromRational x :: Double), toRational (fromRational y :: Double))

@@ -0,0 +1,307 @@+-- | Boundary, topology validation, and representative interior points.+module Data.Geometry.Topology.Unary (+    boundary,+    isSimple,+    isRing,+    isValid,+    pointOnSurface,+) where++import Control.Applicative ((<|>))+import Data.Geometry.Internal+import Data.Geometry.Topology.Planar+import qualified Data.Graph as Graph+import Data.List (group, sort)+import qualified Data.List as List+import qualified Data.Map.Strict as Map+import Data.Maybe (fromMaybe, listToMaybe, mapMaybe)+import Data.Ord (comparing)+import qualified Data.Set as Set+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++{- | The topological boundary. Geometry collections have no boundary operation+in GEOS and return 'Nothing'. Line endpoints use the mod-2 boundary rule.+-}+boundary :: Geometry -> Maybe Geometry+boundary geometry = case geometry of+    PointGeometry _ -> Just (GeometryCollection V.empty)+    MultiPoint _ -> Just (GeometryCollection V.empty)+    LineString line -> Just (MultiPoint (U.fromList (lineBoundary line)))+    MultiLineString lines' -> Just (MultiPoint (U.fromList endpoints))+      where+        counts = Map.fromListWith (\(_, n) (old, m) -> (old, n + m)) [(pointXY point, (point, 1 :: Int)) | point <- foldMap lineEnds lines']+        selected = [point | (point, count) <- Map.elems counts, odd count]+        withZ = not (all (isNaN . elevationOrNaN) selected)+        endpoints = map (boundaryPoint withZ) selected+    Polygon rings@(PolygonRings shell holes)+        | geometryEmpty (Polygon rings) -> Just (MultiLineString V.empty)+        | V.null holes -> Just (LineString shell)+        | otherwise -> Just (MultiLineString (V.cons shell holes))+    MultiPolygon polygons -> Just (MultiLineString (V.fromList (foldMap polygonBoundary polygons)))+    GeometryCollection _ -> Nothing++-- | Return the rings of a nonempty polygon.+polygonBoundary :: PolygonRings -> [Coordinates]+polygonBoundary rings@(PolygonRings shell holes)+    | geometryEmpty (Polygon rings) = []+    | otherwise = shell : V.toList holes++-- | Return the two endpoints, also for a closed line.+lineEnds :: Coordinates -> [Point]+lineEnds line = withCoordinates endpoints line+  where+    endpoints points = case U.uncons points of+        Nothing -> []+        Just (first, _) -> map (pointFromComponents (dimensionsOf line) . coordinateComponents) [first, U.last points]++-- | Return the distinct endpoints of an open line.+lineBoundary :: Coordinates -> [Point]+lineBoundary line = case lineEnds line of+    [a, b] | pointXY a /= pointXY b -> [a, b]+    _ -> []++-- | Copy coordinates without changing their layout.+coordinatePoints :: Coordinates -> [Point]+coordinatePoints coordinates = withCoordinates (map (pointFromComponents (dimensionsOf coordinates) . coordinateComponents) . U.toList) coordinates++-- | Read the XY ordinates of a nonempty point.+pointXY :: Point -> (Double, Double)+pointXY point = case withPoint coordinateComponents point of+    Just (x, y, _, _) -> (x, y)+    Nothing -> (0 / 0, 0 / 0)++-- | Read elevation, with NaN for absent Z.+elevationOrNaN :: Point -> Double+elevationOrNaN (PointXYZ (XYZ _ _ z)) = z+elevationOrNaN (PointXYZM (XYZM _ _ z _)) = z+elevationOrNaN _ = 0 / 0++-- | Multi-line boundaries discard M and share their output coordinate layout.+boundaryPoint :: Bool -> Point -> Point+boundaryPoint withZ point =+    let (x, y) = pointXY point+     in if withZ then PointXYZ (XYZ x y (elevationOrNaN point)) else PointXY (XY x y)++{- | Whether curves have no self-intersections except their closing endpoint,+and multi-points have no duplicate XY positions. Polygon rings are checked+separately. Geometry collection members are checked separately, as in GEOS.+The planar tests require finite XY coordinates.+-}+isSimple :: Geometry -> Bool+isSimple geometry = case geometry of+    PointGeometry _ -> True+    MultiPoint points -> let xy = mapMaybe (withPoint position) (U.toList points) in length xy == Set.size (Set.fromList xy)+    LineString line -> simpleLines [positions line]+    MultiLineString lines' -> simpleLines (map positions (V.toList lines'))+    Polygon (PolygonRings shell holes) -> all (simpleLines . (: []) . positions) (shell : V.toList holes)+    MultiPolygon polygons -> V.all (isSimple . Polygon) polygons+    GeometryCollection children -> V.all isSimple children++-- | Whether a line string is both closed and simple. Other families return false.+isRing :: Geometry -> Bool+isRing geometry@(LineString line) = case lineEnds line of+    [a, b] -> pointXY a == pointXY b && isSimple geometry+    _ -> False+isRing _ = False++-- | Compare segment intersections after removing adjacent repeated positions.+simpleLines :: [[Position]] -> Bool+simpleLines lines' = all allowed (overlappingPairs [(edge, value) | value@(_, _, _, _, edge) <- indexed])+  where+    indexed = [(lineIndex, edgeIndex, length points - 2, first == last points, edge) | (lineIndex, input) <- zip [0 :: Int ..] lines', let points = [p | p : _ <- group input], first : _ <- [points], (edgeIndex, edge) <- zip [0 :: Int ..] (lineSegments points)]+    allowed ((lineA, indexA, lastA, closedA, edgeA), (lineB, indexB, lastB, closedB, edgeB)) = case segmentIntersection edgeA edgeB of+        [] -> True+        [p]+            | p `notElem` ends edgeA || p `notElem` ends edgeB -> False+            | lineA == lineB && abs (indexA - indexB) == 1 -> True+            | otherwise -> endpoint p indexA lastA edgeA && endpoint p indexB lastB edgeB && (lineA == lineB || not (closedA || closedB))+        _ -> False+    endpoint p i lastIndex (a, b) = (i == 0 && p == a) || (i == lastIndex && p == b)+    ends (a, b) = [a, b]++{- | Check finite XY positions and the Simple Features topology rules.+Polygon holes must lie inside the shell. Ring contacts must leave the+polygon interior connected. Multi-polygon interiors must be disjoint.+-}+isValid :: Geometry -> Bool+isValid geometry = finiteGeometry geometry && valid geometry+  where+    valid shape = case shape of+        PointGeometry _ -> True+        MultiPoint _ -> True+        LineString line -> validLine (positions line)+        MultiLineString lines' -> V.all (validLine . positions) lines'+        Polygon rings -> validPolygon (polygonPositions rings)+        MultiPolygon polygons -> all validPolygon rings && all disjointPolygons (overlappingPairs [(bounds, polygon) | polygon <- rings, Just bounds <- [pointBounds (concat polygon)]])+          where+            rings = map polygonPositions (V.toList polygons)+        GeometryCollection children -> V.all valid children++-- | Check XY only. Elevations and measures do not affect topology.+finiteGeometry :: Geometry -> Bool+finiteGeometry geometry = case geometry of+    PointGeometry point -> maybe True finiteCoordinate (withPoint coordinateComponents point)+    MultiPoint points -> U.all (finiteGeometry . PointGeometry) points+    LineString line -> finiteLine line+    MultiLineString lines' -> V.all finiteLine lines'+    Polygon (PolygonRings shell holes) -> finiteLine shell && V.all finiteLine holes+    MultiPolygon polygons -> V.all (finiteGeometry . Polygon) polygons+    GeometryCollection children -> V.all finiteGeometry children+  where+    finiteCoordinate (x, y, _, _) = finite x && finite y+    finiteLine = withCoordinates (U.all (finiteCoordinate . coordinateComponents))++-- | Read the shell followed by holes in the XY plane.+polygonPositions :: PolygonRings -> [[Position]]+polygonPositions (PolygonRings shell holes) = map positions (shell : V.toList holes)++-- | Nonempty valid lines have at least two distinct XY positions.+validLine :: [Position] -> Bool+validLine [] = True+validLine (first : rest) = any (/= first) rest++-- | Nonempty valid rings are closed, simple, and have three distinct vertices.+validRing :: [Position] -> Bool+validRing [] = True+validRing points@(first : _) = first == last points && Set.size (Set.fromList points) >= 3 && simpleLines [points]++-- | A ring with shared edge and point-location indexes for validity checks.+data IndexedRing = IndexedRing+    { ringPoints :: [Position]+    -- ^ The original ring positions.+    , ringSize :: Int+    -- ^ The number of stored positions.+    , queryRing :: Segment -> [Segment]+    -- ^ Select ring edges whose bounds meet a query box.+    , locateRing :: Position -> Location+    -- ^ Classify a position against the ring.+    }++-- | Prepare a ring once for all comparisons with the other rings.+indexRing :: [Position] -> IndexedRing+indexRing points = IndexedRing points (length points) (segmentQuery (ringSegments points)) (prepareRing points)++-- | Validate rings, containment, hole separation, and contact cycles.+validPolygon :: [[Position]] -> Bool+validPolygon [] = True+validPolygon (shell : holes)+    | null shell = all null holes+    | otherwise = all (validRing . ringPoints) rings && all noOverlap pairs && all insideShell nonemptyHoles && all separateHoles holePairs && acyclic contacts+  where+    shellIndex = indexRing shell+    nonemptyHoles = map indexRing (filter (not . null) holes)+    rings = shellIndex : nonemptyHoles+    candidates = overlappingPairs [(bounds, (i, ring)) | (i, ring) <- zip [0 :: Int ..] rings, Just bounds <- [pointBounds (ringPoints ring)]]+    pairs = [(i, a, j, b, intersections a b) | ((i, a), (j, b)) <- candidates]+    intersections a b+        | ringSize a <= ringSize b = ringIntersections (ringPoints a) (queryRing b)+        | otherwise = ringIntersections (ringPoints b) (queryRing a)+    noOverlap (_, _, _, _, crossings) = all ((<= 1) . length) crossings+    holePairs = [(a, b) | (i, a, j, b, _) <- pairs, i > 0, j > 0]+    insideShell hole = all ((/= Exterior) . locateRing shellIndex) (ringSamplesAgainst (ringPoints hole) (queryRing shellIndex))+    separateHoles (a, b) = all ((/= Interior) . locateRing b) (ringSamplesAgainst (ringPoints a) (queryRing b)) && all ((/= Interior) . locateRing a) (ringSamplesAgainst (ringPoints b) (queryRing a))+    contacts = Set.toList (Set.fromList [(Left k, Right p) | (i, _, j, _, crossings) <- pairs, p <- concat crossings, k <- [i, j]])++-- | Intersections with an existing edge index. Two points denote a shared edge.+ringIntersections :: [Position] -> (Segment -> [Segment]) -> [[Position]]+ringIntersections ring query = [segmentIntersection x y | x <- ringSegments ring, y <- query x]++-- | Whether a ring shares a segment of positive length with indexed boundaries.+sharedEdge :: [Position] -> (Segment -> [Segment]) -> Bool+sharedEdge ring query = any ((> 1) . length) (ringIntersections ring query)++-- | Split one ring against indexed boundaries before sampling its open edges.+ringSamplesAgainst :: [Position] -> (Segment -> [Segment]) -> [Position]+ringSamplesAgainst ring query = concatMap sample (ringSegments ring)+  where+    sample edge@(a, b) = map midpoint (lineSegments (unique (a : b : concatMap (segmentIntersection edge) (query edge))))++{- | A cycle through distinct contact positions disconnects a polygon interior.+An undirected graph is acyclic when its edge count is its vertex count minus+its connected-component count. Each tree in the DFS forest is one component.+-}+acyclic :: [(Either Int Position, Either Int Position)] -> Bool+acyclic contacts = length contacts == Map.size adjacent - length (Graph.dff graph)+  where+    adjacent = Map.fromListWith (++) [(a, [b]) | (first, second) <- contacts, (a, b) <- [(first, second), (second, first)]]+    (graph, _, _) = Graph.graphFromEdges [((), point, neighbors) | (point, neighbors) <- Map.toList adjacent]++-- | Multi-polygons can touch at isolated points but cannot share interior area.+disjointPolygons :: ([[Position]], [[Position]]) -> Bool+disjointPolygons (a, b) =+    not (any (`sharedEdge` queryB) a)+        && all ((/= Interior) . preparePolygon b) (samples a queryB)+        && all ((/= Interior) . preparePolygon a) (samples b queryA)+  where+    queryA = segmentQuery (concatMap ringSegments a)+    queryB = segmentQuery (concatMap ringSegments b)+    samples rings query = concatMap (`ringSamplesAgainst` query) rings++{- | Choose a point on a nonempty component, preferring polygons, then lines,+then points. For polygons, use an interior horizontal interval. For lines,+prefer an interior stored vertex, then an endpoint. Empty components are skipped.+Results use XY coordinates. Empty input gives an empty XY point.+The selected point can differ from other implementations.+-}+pointOnSurface :: Geometry -> Point+pointOnSurface geometry = fromMaybe (EmptyPoint DimXY) (surfacePoint <|> storedPoint)+  where+    surfacePoint = listToMaybe (mapMaybe polygonInterior (polygonMembers geometry))+    storedPoint = boundaryPoint False <$> listToMaybe (interiors ++ concatMap endpoints lines' ++ mapMaybe pointCoordinate (pointMembers geometry))+    lines' = map coordinatePoints (lineMembers geometry)+    interiors = concatMap (drop 1 . takeInterior) lines'+    takeInterior [] = []+    takeInterior points = init points+    endpoints [] = []+    endpoints points@(first : _) = [first, last points]+    pointCoordinate (EmptyPoint _) = Nothing+    pointCoordinate point = Just point++-- | List atomic point members in input order.+pointMembers :: Geometry -> [Point]+pointMembers (PointGeometry point) = [point]+pointMembers (MultiPoint points) = U.toList points+pointMembers (GeometryCollection children) = foldMap pointMembers children+pointMembers _ = []++-- | List atomic line members in input order.+lineMembers :: Geometry -> [Coordinates]+lineMembers (LineString line) = [line]+lineMembers (MultiLineString lines') = V.toList lines'+lineMembers (GeometryCollection children) = foldMap lineMembers children+lineMembers _ = []++-- | List atomic polygon members in input order.+polygonMembers :: Geometry -> [PolygonRings]+polygonMembers (Polygon rings) = [rings]+polygonMembers (MultiPolygon polygons) = V.toList polygons+polygonMembers (GeometryCollection children) = foldMap polygonMembers children+polygonMembers _ = []++{- | Select the widest horizontal interval with a representable point.+Compute crossings exactly and check that rounding stays within the interval.+Use a shell vertex if no interval contains a representable coordinate.+-}+polygonInterior :: PolygonRings -> Maybe Point+polygonInterior (PolygonRings shell holes) = case map pointXY (coordinatePoints shell) of+    [] -> Nothing+    shellPoints@((x, y) : _) -> Just (snd (List.maximumBy (comparing fst) ((0, PointXY (XY x y)) : intervals)))+      where+        rings = shellPoints : map (map pointXY . coordinatePoints) (V.toList holes)+        ys = map snd (concat rings)+        lo = minimum (map snd shellPoints)+        hi = maximum (map snd shellPoints)+        center = mean lo hi+        below = maximum (lo : filter (<= center) ys)+        above = minimum (hi : filter (> center) ys)+        scanY = mean below above+        crossings = sort [crossing a b | ring <- rings, (a@(_, ay), b@(_, by)) <- zip ring (drop 1 ring), ay /= by, min ay by <= scanY, max ay by >= scanY, not (ay == scanY && by < scanY), not (by == scanY && ay < scanY)]+        mean a b = fromRational ((toRational a + toRational b) / 2)+        crossing (ax, ay) (bx, by) = toRational ax + (toRational scanY - toRational ay) * (toRational bx - toRational ax) / (toRational by - toRational ay)+        intervals = [(b - a, PointXY (XY midpointX scanY)) | (a, b) <- adjacentPairs crossings, a < b, let midpointX = fromRational ((a + b) / 2), toRational midpointX >= a, toRational midpointX <= b]++-- | Pair sorted crossings into interior intervals.+adjacentPairs :: [a] -> [(a, a)]+adjacentPairs (a : b : rest) = (a, b) : adjacentPairs rest+adjacentPairs _ = []

@@ -0,0 +1,340 @@+{-# LANGUAGE ScopedTypeVariables #-}++{- | Checked ISO WKB decoding and encoding.++Collection members retain their own coordinate layouts. Polygon rings share+one header; encoding pads them to their combined layout. Children can use+different byte orders and layouts. The codecs check+line lengths, ring closure, counts, and type codes. They accept non-finite+ordinates and do not validate polygon topology. EWKB and SRIDs are not supported.+-}+module Data.Geometry.WKB (decodeWKB, encodeWKB) where++import Control.Monad (unless, when)+import Control.Monad.ST (runST)+import Data.Binary.Get (Get, bytesRead, getByteString, getDoublebe, getDoublele, getWord32be, getWord32le, getWord8, lookAhead, runGetOrFail, skip)+import Data.Bits (shiftL, (.|.))+import Data.ByteString (ByteString)+import qualified Data.ByteString as BS+import Data.ByteString.Builder (Builder)+import qualified Data.ByteString.Builder as Builder+import Data.ByteString.Builder.Prim ((>$<), (>*<))+import qualified Data.ByteString.Builder.Prim as Prim+import qualified Data.ByteString.Lazy as BL+import Data.Geometry.Internal+import Data.Int (Int64)+import Data.List (stripPrefix)+import Data.Maybe (fromMaybe)+import Data.Proxy (Proxy (..))+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM+import Data.Word (Word32, Word64)+import GHC.Float (castDoubleToWord64, castWord64ToDouble)++{- | Decode one complete ISO WKB geometry. Each child uses its own header.+A point with NaN in both X and Y becomes an empty point with the same layout.+Return 'Left' for malformed WKB, invalid construction, or trailing bytes.+-}+decodeWKB :: ByteString -> Either String Geometry+decodeWKB bytes = runDecoder (getGeometry (fromIntegral (BS.length bytes))) bytes++{- | Encode little-endian ISO WKB. Child headers retain their layouts.+Polygon rings use their combined layout, with NaN for absent Z or M ordinates.+Finite ordinates retain their exact bits, including negative zero.+Return 'Left' for invalid construction or a count that exceeds 32 bits.+-}+encodeWKB :: Geometry -> Either String ByteString+encodeWKB geometry = do+    validateGeometry checkedLength geometry+    pure (BL.toStrict (Builder.toLazyByteString (putGeometry geometry)))++-- | Run a decoder and require complete input consumption.+runDecoder :: Get a -> ByteString -> Either String a+runDecoder parser bytes = case runGetOrFail parser (BL.fromStrict bytes) of+    Left (_, _, message) -> Left $ case stripPrefix "Geometry WKB " message of+        Just _ -> message+        Nothing -> "Geometry WKB " ++ fromMaybe message (stripPrefix "Geometry " message)+    Right (remaining, _, value)+        | BL.null remaining -> Right value+        | otherwise -> Left "Geometry WKB has trailing bytes"++-- | Read one geometry's byte order, dimensions, and family.+getHeader :: Get (Bool, Dimensions, Word32)+getHeader = do+    marker <- getWord8+    little <- case marker of+        0 -> pure False+        1 -> pure True+        _ -> fail "Geometry WKB has an invalid byte order"+    tag <- getWord little+    let (dimensionTag, family) = tag `quotRem` 1000+    unless (dimensionTag <= 3 && family >= 1 && family <= 7) $+        fail "Geometry WKB has an unsupported type"+    pure (little, toEnum (fromIntegral dimensionTag), family)++-- | Read a word in the current geometry's byte order.+getWord :: Bool -> Get Word32+getWord little = if little then getWord32le else getWord32be++-- | Check a count against remaining input before allocating a vector.+getCount :: Int64 -> Bool -> Int64 -> Get Int+getCount total little minimumBytes = do+    count <- getWord little+    consumed <- bytesRead+    when (fromIntegral count > (total - consumed) `div` minimumBytes) $+        fail "Geometry WKB count exceeds the remaining bytes"+    pure (fromIntegral count)++-- | Read a complete geometry, including each collection member's own header.+getGeometry :: Int64 -> Get Geometry+getGeometry total = do+    (little, dimensions, family) <- getHeader+    case family of+        1 -> PointGeometry <$> getPoint little dimensions+        2 -> LineString <$> getCoordinates total little dimensions+        3 -> Polygon <$> getPolygon total little dimensions+        4 -> MultiPoint <$> getMultiPoints total little+        5 -> do+            count <- getCount total little 9+            MultiLineString <$> V.replicateM count (getChild 2 (getCoordinates total))+        6 -> do+            count <- getCount total little 9+            MultiPolygon <$> V.replicateM count (getChild 3 (getPolygon total))+        _ -> do+            count <- getCount total little 9+            GeometryCollection <$> V.replicateM count (do child <- getGeometry total; pure $! child)++-- | Check a multi-geometry child header before reading its typed body.+getChild :: Word32 -> (Bool -> Dimensions -> Get a) -> Get a+getChild expected body = do+    (little, dimensions, family) <- getHeader+    unless (family == expected) (fail "Geometry WKB multi child has the wrong family")+    child <- body little dimensions+    pure $! child++-- | Read a polygon's ring count and validate its shell and holes.+getPolygon :: Int64 -> Bool -> Dimensions -> Get PolygonRings+getPolygon total little dimensions = do+    count <- getCount total little 4+    rings <-+        if count == 0+            then pure (PolygonRings (emptyCoordinates dimensions) V.empty)+            else PolygonRings <$> getCoordinates total little dimensions <*> V.replicateM (count - 1) (getCoordinates total little dimensions)+    either fail pure (validatePolygon rings)+    pure rings++-- | Read point ordinates before applying WKB's XY-NaN empty convention.+getPoint :: Bool -> Dimensions -> Get Point+getPoint little dimensions = do+    let number = if little then getDoublele else getDoublebe+    x <- number+    y <- number+    (z, m) <- case dimensions of+        DimXY -> pure (0, 0)+        DimXYZ -> do z <- number; pure (z, 0)+        DimXYM -> do m <- number; pure (0, m)+        DimXYZM -> (,) <$> number <*> number+    pure (if isNaN x && isNaN y then EmptyPoint dimensions else pointFromComponents dimensions (x, y, z, m))++-- | Read one checked coordinate block into its typed unboxed buffer.+getCoordinates :: Int64 -> Bool -> Dimensions -> Get Coordinates+getCoordinates total little dimensions = do+    let stride = 8 * dimensionCount dimensions+    count <- getCount total little (fromIntegral stride)+    when (count == 1) (fail "Geometry line must have zero or at least two coordinates")+    if count == 0+        then pure (emptyCoordinates dimensions)+        else do+            bytes <- getByteString (count * stride)+            let values :: (Coordinate c) => U.Vector c+                values = U.generate count (\i -> coordinateAt little bytes (i * stride))+            pure $ case dimensions of+                DimXY -> CoordinatesXY values+                DimXYZ -> CoordinatesXYZ values+                DimXYM -> CoordinatesXYM values+                DimXYZM -> CoordinatesXYZM values++-- | Read point children without a separate Get action for each child.+getMultiPoints :: Int64 -> Bool -> Get (U.Vector Point)+getMultiPoints total little = do+    count <- getCount total little 21+    if count == 0+        then pure U.empty+        else do+            consumed <- bytesRead+            bytes <- lookAhead (getByteString (fromIntegral (total - consumed)))+            (points, size) <- either fail pure (multiPointsAt count bytes)+            skip size+            pure points++-- | Read checked variable-width point records into one unboxed vector.+multiPointsAt :: Int -> ByteString -> Either String (U.Vector Point, Int)+multiPointsAt count bytes = runST $ do+    target <- UM.new count+    let go index offset+            | index == count = do+                points <- U.unsafeFreeze target+                pure (Right (points, offset))+            | BS.length bytes - offset < 5 = pure (Left "Geometry WKB point header exceeds the remaining bytes")+            | otherwise = case BS.index bytes offset of+                0 -> child index offset False+                1 -> child index offset True+                _ -> pure (Left "Geometry WKB has an invalid byte order")+        child index offset little+            | dimensionTag > 3 || family < 1 || family > 7 = pure (Left "Geometry WKB has an unsupported type")+            | family /= 1 = pure (Left "Geometry WKB multi child has the wrong family")+            | stride > BS.length bytes - offset = pure (Left "Geometry WKB point exceeds the remaining bytes")+            | otherwise = UM.write target index value >> go (index + 1) (offset + stride)+          where+            (dimensionTag, family) = word32At little bytes (offset + 1) `quotRem` 1000+            dimensions = toEnum (fromIntegral dimensionTag)+            stride = 5 + 8 * dimensionCount dimensions+            value =+                let (x, y, z, m) = componentsAt dimensions little bytes (offset + 5)+                 in if isNaN x && isNaN y then EmptyPoint dimensions else pointFromComponents dimensions (x, y, z, m)+    go 0 0++-- | Read the ordinates present in a checked coordinate record.+{-# INLINE componentsAt #-}+componentsAt :: Dimensions -> Bool -> ByteString -> Int -> (Double, Double, Double, Double)+componentsAt dimensions little bytes offset =+    let number i = castWord64ToDouble (word64At little bytes (offset + i))+        x = number 0+        y = number 8+     in case dimensions of+            DimXY -> (x, y, 0, 0)+            DimXYZ -> (x, y, number 16, 0)+            DimXYM -> (x, y, 0, number 16)+            DimXYZM -> (x, y, number 16, number 24)++-- | Read a typed coordinate from a block with a checked byte length.+coordinateAt :: forall c. (Coordinate c) => Bool -> ByteString -> Int -> c+coordinateAt little bytes offset = coordinateFromComponents (componentsAt (coordinateDimensions (Proxy :: Proxy c)) little bytes offset)++-- | Read four checked bytes in either byte order without alignment assumptions.+word32At :: Bool -> ByteString -> Int -> Word32+word32At little bytes offset =+    let byte i = fromIntegral (BS.index bytes (offset + i))+     in if little+            then byte 0 .|. shiftL (byte 1) 8 .|. shiftL (byte 2) 16 .|. shiftL (byte 3) 24+            else shiftL (byte 0) 24 .|. shiftL (byte 1) 16 .|. shiftL (byte 2) 8 .|. byte 3++-- | Read eight checked bytes and preserve the exact IEEE-754 representation.+word64At :: Bool -> ByteString -> Int -> Word64+word64At little bytes offset =+    let byte i = fromIntegral (BS.index bytes (offset + i))+     in if little+            then+                byte 0+                    .|. shiftL (byte 1) 8+                    .|. shiftL (byte 2) 16+                    .|. shiftL (byte 3) 24+                    .|. shiftL (byte 4) 32+                    .|. shiftL (byte 5) 40+                    .|. shiftL (byte 6) 48+                    .|. shiftL (byte 7) 56+            else+                shiftL (byte 0) 56+                    .|. shiftL (byte 1) 48+                    .|. shiftL (byte 2) 40+                    .|. shiftL (byte 3) 32+                    .|. shiftL (byte 4) 24+                    .|. shiftL (byte 5) 16+                    .|. shiftL (byte 6) 8+                    .|. byte 7++-- | Check that a vector length fits the WKB unsigned 32-bit count.+checkedLength :: Int -> Either String ()+checkedLength count = when (toInteger count > toInteger (maxBound :: Word32)) (Left "Geometry count exceeds Word32")++-- | Write a geometry with its own aggregate header and each child's own layout.+putGeometry :: Geometry -> Builder+putGeometry geometry = snd (geometryBuilder geometry) DimXY++-- | Compute layouts once. Empty containers inherit their containing header.+geometryBuilder :: Geometry -> (Maybe Dimensions, Dimensions -> Builder)+geometryBuilder geometry = case geometry of+    PointGeometry point -> tagged (Just (pointDimensions point)) (const (putPoint point))+    LineString points -> tagged (Just (dimensionsOf points)) (const (withCoordinates (putLine (dimensionsOf points)) points))+    Polygon rings@(PolygonRings shell holes) ->+        let dimensions = polygonDimensions rings+         in tagged (Just dimensions)+                $ const+                $ if coordinatesEmpty shell+                    then putLength 0+                    else putLength (1 + V.length holes) <> withCoordinates (putLine dimensions) shell <> V.foldMap (withCoordinates (putLine dimensions)) holes+    MultiPoint points -> tagged (if U.null points then Nothing else Just (geometryDimensions geometry)) (const (putLength (U.length points) <> U.foldMap (putGeometry . PointGeometry) points))+    MultiLineString lines' -> tagged (if V.null lines' then Nothing else Just (geometryDimensions geometry)) (const (putLength (V.length lines') <> V.foldMap (putGeometry . LineString) lines'))+    MultiPolygon polygons -> tagged (if V.null polygons then Nothing else Just (geometryDimensions geometry)) (const (putLength (V.length polygons) <> V.foldMap (putGeometry . Polygon) polygons))+    GeometryCollection children ->+        let parts = V.map geometryBuilder children+            dimensions = V.foldl' (\acc (layout, _) -> combine acc layout) Nothing parts+         in tagged dimensions (\layout -> putLength (V.length children) <> V.foldMap (\(_, render) -> render layout) parts)+  where+    combine Nothing second = second+    combine first Nothing = first+    combine (Just first) (Just second) = Just (unionDimensions first second)+    tagged declared body =+        ( declared+        , \inherited ->+            let dimensions = fromMaybe inherited declared+             in Builder.word8 1 <> Builder.word32LE (geometryFamily geometry + 1000 * fromIntegral (fromEnum dimensions)) <> body dimensions+        )++-- | The ISO WKB family tag for a geometry.+geometryFamily :: Geometry -> Word32+geometryFamily (PointGeometry _) = 1+geometryFamily (LineString _) = 2+geometryFamily (Polygon _) = 3+geometryFamily (MultiPoint _) = 4+geometryFamily (MultiLineString _) = 5+geometryFamily (MultiPolygon _) = 6+geometryFamily (GeometryCollection _) = 7++-- | Write a vector length that validation has checked.+putLength :: Int -> Builder+putLength = Builder.word32LE . fromIntegral++-- | Empty points use canonical quiet NaNs in every declared ordinate.+putPoint :: Point -> Builder+putPoint point = fromMaybe empty (withPoint (putCoordinate dimensions) point)+  where+    dimensions = pointDimensions point+    empty = mconcat (replicate (dimensionCount dimensions) (Builder.word64LE 0x7ff8000000000000))++-- | Write exact ordinate bits, padding absent Z or M ordinates with NaN.+putCoordinate :: forall c. (Coordinate c) => Dimensions -> c -> Builder+putCoordinate target coordinate =+    let (x, y, z, m) = coordinateComponents coordinate+        source = coordinateDimensions (Proxy :: Proxy c)+        number = Builder.word64LE . castDoubleToWord64+        nan = castWord64ToDouble 0x7ff8000000000000+        zValue = if source == DimXYZ || source == DimXYZM then z else nan+        mValue = if source == DimXYM || source == DimXYZM then m else nan+        extra = case target of+            DimXY -> mempty+            DimXYZ -> number zValue+            DimXYM -> number mValue+            DimXYZM -> number zValue <> number mValue+     in number x <> number y <> extra++-- | Write a line or ring. Homogeneous buffers use the fixed-width builder.+putLine :: forall c. (Coordinate c) => Dimensions -> U.Vector c -> Builder+putLine dimensions points =+    putLength (U.length points)+        <> if dimensions == coordinateDimensions (Proxy :: Proxy c)+            then Prim.primUnfoldrFixed coordinatePrim next 0+            else U.foldMap (putCoordinate dimensions) points+  where+    next i+        | i >= U.length points = Nothing+        | otherwise = Just (points U.! i, i + 1)++-- | Write a complete coordinate with one buffer-size check.+coordinatePrim :: forall c. (Coordinate c) => Prim.FixedPrim c+coordinatePrim = case coordinateDimensions (Proxy :: Proxy c) of+    DimXY -> (\c -> let (x, y, _, _) = coordinateComponents c in (x, y)) >$< (Prim.doubleLE >*< Prim.doubleLE)+    DimXYZ -> (\c -> let (x, y, z, _) = coordinateComponents c in (x, (y, z))) >$< (Prim.doubleLE >*< Prim.doubleLE >*< Prim.doubleLE)+    DimXYM -> (\c -> let (x, y, _, m) = coordinateComponents c in (x, (y, m))) >$< (Prim.doubleLE >*< Prim.doubleLE >*< Prim.doubleLE)+    DimXYZM -> (\c -> let (x, y, z, m) = coordinateComponents c in ((x, y), (z, m))) >$< ((Prim.doubleLE >*< Prim.doubleLE) >*< (Prim.doubleLE >*< Prim.doubleLE))

@@ -0,0 +1,474 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-}++{- | WKT codecs for the seven Simple Features families.++The decoder retains each collection member's layout. For untagged input,+it infers XY, XYZ, or XYZM from the number of ordinates. XYM requires an M tag.+Within a multi-geometry, the first coordinate sets the layout for the remaining+coordinates. Empty members before that coordinate retain XY.+The codecs check line lengths and ring closure, but not polygon topology.+NaN and infinity are accepted. EWKT and SRIDs are not supported.+-}+module Data.Geometry.WKT (decodeWKT, encodeWKT) where++import Control.Monad (unless, when)+import Control.Monad.ST (runST)+import Control.Monad.Trans.Class (lift)+import Control.Monad.Trans.State.Strict (StateT (..), get, modify', put)+import Data.ByteString.Builder (Builder)+import qualified Data.ByteString.Builder as Builder+import qualified Data.ByteString.Builder.RealFloat as RealFloat+import qualified Data.ByteString.Lazy as BL+import Data.Char (isAsciiLower, isAsciiUpper, isDigit)+import Data.Geometry.Internal+import Data.Maybe (fromMaybe)+import Data.Proxy (Proxy (..))+import Data.Ratio ((%))+import Data.Text (Text)+import qualified Data.Text as Text+import qualified Data.Text.Encoding as TextEncoding+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM+import GHC.Float (castWord64ToDouble)++-- | The remaining text and a controlled parse error.+type Parser = StateT Text (Either String)++-- | The geometry family selected by a WKT keyword.+data Family+    = -- | A point or an empty point.+      PointFamily+    | -- | One coordinate sequence.+      LineFamily+    | -- | An exterior ring and holes.+      PolygonFamily+    | -- | A collection of points.+      MultiPointFamily+    | -- | A collection of line strings.+      MultiLineFamily+    | -- | A collection of polygons.+      MultiPolygonFamily+    | -- | A collection of arbitrary geometries.+      CollectionFamily++{- | Decode one complete geometry and preserve each member's coordinate layout.+Return 'Left' for malformed WKT or trailing input. See the module documentation+for layout inference and the construction checks.+-}+decodeWKT :: Text -> Either String Geometry+decodeWKT input = do+    ((geometry, _), remaining) <- runStateT (geometryParser <* spaces) input+    if Text.null remaining then Right geometry else Left "Geometry WKT has trailing input"++{- | Write dimension tags and shortest scientific decimal ordinates.+Multi-geometries and polygon rings pad absent Z or M ordinates with NaN.+Geometry collections use a parent tag when all children share one output layout.+Mixed collections omit that tag and retain each child's layout and ordinates.+Mixed-layout collections extend the standard WKT grammar.+Return 'Left' for invalid line lengths, ring closure, or polygon emptiness.+-}+encodeWKT :: Geometry -> Either String Text+encodeWKT geometry = do+    validateGeometry (const (Right ())) geometry+    pure (TextEncoding.decodeUtf8 (BL.toStrict (Builder.toLazyByteString (snd (geometryWKT geometry) DimXY))))++-- | Stop parsing with a geometry-specific error.+failure :: String -> Parser a+failure message = lift (Left ("Geometry WKT " ++ message))++-- | Consume whitespace before a structural token.+spaces :: Parser ()+spaces = modify' (Text.dropWhile whitespace)++-- | Accept the ASCII whitespace that WKT readers use: space, tab, LF, and CR.+whitespace :: Char -> Bool+whitespace c = c == ' ' || c == '\t' || c == '\n' || c == '\r'++-- | Keywords and named numbers use ASCII letters only.+letter :: Char -> Bool+letter c = isAsciiLower c || isAsciiUpper c++-- | Require a punctuation character, with optional leading whitespace.+symbol :: Char -> Parser ()+symbol expected = do+    spaces+    input <- get+    case Text.uncons input of+        Just (actual, rest) | actual == expected -> put rest+        _ -> failure ("expected " ++ show expected)++-- | Read an ASCII keyword without consuming following whitespace.+word :: Parser Text+word = do+    spaces+    input <- get+    let (name, rest) = Text.span letter input+    when (Text.null name) (failure "expected a keyword")+    put rest+    pure (Text.toUpper name)++-- | Consume EMPTY when it is the next complete keyword.+emptyKeyword :: Parser Bool+emptyKeyword = do+    spaces+    input <- get+    let (name, rest) = Text.span letter input+    if Text.toUpper name == "EMPTY"+        then put rest >> pure True+        else pure False++-- | Read the geometry family and an attached or separate dimension suffix.+header :: Parser (Family, Maybe Dimensions)+header = do+    name <- word+    let families = [("POINT", PointFamily), ("LINESTRING", LineFamily), ("POLYGON", PolygonFamily), ("MULTIPOINT", MultiPointFamily), ("MULTILINESTRING", MultiLineFamily), ("MULTIPOLYGON", MultiPolygonFamily), ("GEOMETRYCOLLECTION", CollectionFamily)]+        suffixes = [("ZM", DimXYZM), ("Z", DimXYZ), ("M", DimXYM)]+        attached = [(family, dimensions) | (suffix, dimensions) <- suffixes, Just base <- [Text.stripSuffix suffix name], Just family <- [lookup base families]]+    case lookup name families of+        Just family -> do+            spaces+            input <- get+            let (tag, rest) = Text.span letter input+            case lookup (Text.toUpper tag) suffixes of+                Just dimensions -> put rest >> pure (family, Just dimensions)+                Nothing -> pure (family, Nothing)+        Nothing -> case attached of+            [(family, dimensions)] -> pure (family, Just dimensions)+            _ -> failure "has an unsupported geometry type"++-- | Read a geometry and retain the parsed layout for an explicit parent tag.+geometryParser :: Parser (Geometry, LayoutSummary)+geometryParser = do+    (family, declared) <- header+    case family of+        CollectionFamily -> do+            (members, ()) <- boxedSequence () (stateless geometryParser)+            case declared of+                Just dimensions -> unless (V.all ((== Uniform dimensions) . snd) members) (failure "has mixed coordinate dimensions")+                Nothing -> pure ()+            let layout = if V.null members then Uniform (fromMaybe DimXY declared) else V.foldl' (\acc (_, child) -> combineLayout acc child) Inherited members+            pure (GeometryCollection (V.map fst members), layout)+        PointFamily -> known PointGeometry (point True declared)+        LineFamily -> known LineString (line declared)+        PolygonFamily -> known Polygon (polygon declared)+        MultiPointFamily -> known MultiPoint (multiPoint declared)+        MultiLineFamily -> known MultiLineString (boxedSequence declared line)+        MultiPolygonFamily -> known MultiPolygon (boxedSequence declared polygon)+  where+    known wrap parser = do+        (value, dimensions) <- parser+        pure (wrap value, Uniform (fromMaybe DimXY dimensions))+    line current = do+        (values, dimensions) <- coordinates current+        lift (validateLine values)+        pure (values, dimensions)++-- | Inspect one coordinate without consuming its text.+lookAheadParser :: Parser a -> Parser a+lookAheadParser parser = StateT $ \input -> do+    (value, _) <- runStateT parser input+    pure (value, input)++-- | Infer an untagged coordinate's two, three, or four ordinates.+inferDimensions :: Parser Dimensions+inferDimensions = do+    _ <- number+    _ <- nextNumber+    third <- more+    if not third+        then pure DimXY+        else do+            _ <- nextNumber+            fourth <- more+            if fourth then nextNumber >> pure DimXYZM else pure DimXYZ+  where+    more = do+        remaining <- get+        pure $ case Text.uncons (Text.dropWhile whitespace remaining) of+            Nothing -> False+            Just (c, _) -> c /= ',' && c /= ')'++-- | Read a point body. Bare MULTIPOINT coordinates cannot contain EMPTY.+point :: Bool -> Maybe Dimensions -> Parser (Point, Maybe Dimensions)+point parenthesized current = do+    empty <- if parenthesized then emptyKeyword else pure False+    if empty+        then pure (EmptyPoint (fromMaybe DimXY current), current)+        else do+            when parenthesized (symbol '(')+            spaces+            dimensions <- maybe (lookAheadParser inferDimensions) pure current+            x <- number+            y <- nextNumber+            (z, m) <- case dimensions of+                DimXY -> pure (0, 0)+                DimXYZ -> do z <- nextNumber; pure (z, 0)+                DimXYM -> do m <- nextNumber; pure (0, m)+                DimXYZM -> (,) <$> nextNumber <*> nextNumber+            when parenthesized (symbol ')')+            pure (pointFromComponents dimensions (x, y, z, m), Just dimensions)++-- | Read one sequence into the unboxed buffer for its inferred layout.+coordinates :: Maybe Dimensions -> Parser (Coordinates, Maybe Dimensions)+coordinates current = do+    empty <- emptyKeyword+    if empty+        then pure (emptyCoordinates (fromMaybe DimXY current), current)+        else do+            dimensions <- maybe (lookAheadParser (symbol '(' *> spaces *> inferDimensions)) pure current+            values <- case dimensions of+                DimXY -> CoordinatesXY . fst <$> unboxedSequence () (stateless coordinate)+                DimXYZ -> CoordinatesXYZ . fst <$> unboxedSequence () (stateless coordinate)+                DimXYM -> CoordinatesXYM . fst <$> unboxedSequence () (stateless coordinate)+                DimXYZM -> CoordinatesXYZM . fst <$> unboxedSequence () (stateless coordinate)+            pure (values, Just dimensions)++-- | Preserve empty rings and their layouts when reading a polygon.+polygon :: Maybe Dimensions -> Parser (PolygonRings, Maybe Dimensions)+polygon current = do+    (rings, dimensions) <- boxedSequence current coordinates+    let values =+            if V.null rings+                then PolygonRings (emptyCoordinates (fromMaybe DimXY current)) V.empty+                else PolygonRings (V.head rings) (V.tail rings)+    lift (validatePolygon values)+    pure (values, dimensions)++-- | Use one MULTIPOINT spelling throughout its body.+multiPoint :: Maybe Dimensions -> Parser (U.Vector Point, Maybe Dimensions)+multiPoint current = do+    spaces+    input <- get+    -- Inspect only the first token. Uppercasing the rest of the input would be quadratic.+    let first = Text.dropWhile whitespace (Text.drop 1 input)+        parenthesized = Text.isPrefixOf "(" first || Text.toUpper (Text.takeWhile letter first) == "EMPTY"+    unboxedSequence current (point parenthesized)++{- | Read EMPTY or a parenthesized sequence into a boxed vector. Carry the+inferred layout from each element to the next. A list keeps deep nesting+linear: boxed mutable buffers for every open level would be rescanned by+each garbage collection.+-}+boxedSequence :: s -> (s -> Parser (a, s)) -> Parser (V.Vector a, s)+boxedSequence initialState element = do+    empty <- emptyKeyword+    if empty+        then pure (V.empty, initialState)+        else symbol '(' >> go 1 [] initialState+  where+    go !count values current = do+        (value, next) <- element current+        finished <- delimiter+        let values' = value `seq` value : values+        if finished then pure (V.fromListN count (reverse values'), next) else go (count + 1) values' next++{- | Read EMPTY or a parenthesized sequence into an unboxed vector. Write the+elements into a growable buffer, which avoids an intermediate list for long+coordinate sequences.+-}+unboxedSequence :: (U.Unbox a) => s -> (s -> Parser (a, s)) -> Parser (U.Vector a, s)+unboxedSequence initialState element = do+    empty <- emptyKeyword+    if empty+        then pure (U.empty, initialState)+        else do+            symbol '('+            StateT $ \input -> runST $ do+                initial <- UM.new 16+                let go !count buffer current remaining = case runStateT (element current) remaining of+                        Left message -> pure (Left message)+                        Right ((value, next), afterElement) -> case runStateT delimiter afterElement of+                            Left message -> pure (Left message)+                            Right (finished, rest) -> do+                                target <- if count == UM.length buffer then UM.grow buffer (UM.length buffer) else pure buffer+                                UM.write target count value+                                if finished+                                    then do+                                        result <- U.freeze (UM.slice 0 (count + 1) target)+                                        pure (Right ((result, next), rest))+                                    else go (count + 1) target next rest+                go 0 initial initialState input++-- | Read elements that do not share layout state.+stateless :: Parser a -> () -> Parser (a, ())+stateless element () = do+    value <- element+    pure (value, ())++-- | Consume a comma or a closing parenthesis.+delimiter :: Parser Bool+delimiter = do+    spaces+    input <- get+    case Text.uncons input of+        Just (',', rest) -> put rest >> pure False+        Just (')', rest) -> put rest >> pure True+        _ -> failure "expected ',' or ')'"++-- | Parse exactly the ordinates required by the coordinate type.+coordinate :: forall c. (Coordinate c) => Parser c+coordinate = do+    spaces+    x <- number+    y <- nextNumber+    (z, m) <- case coordinateDimensions (Proxy :: Proxy c) of+        DimXY -> pure (0, 0)+        DimXYZ -> do z <- nextNumber; pure (z, 0)+        DimXYM -> do m <- nextNumber; pure (0, m)+        DimXYZM -> (,) <$> nextNumber <*> nextNumber+    pure (coordinateFromComponents (x, y, z, m))++-- | Require whitespace between ordinates so adjacent numbers cannot be split.+nextNumber :: Parser Double+nextNumber = do+    input <- get+    unless (maybe False (whitespace . fst) (Text.uncons input)) $+        failure "requires whitespace between coordinates"+    spaces+    number++{- | Read named IEEE values or an exactly rounded decimal ordinate.+Bound extreme exponents so integer powers stay proportional to input length.+Use 'fromRational' for rounding. @Data.Text.Read.double@ and+@Data.Text.Read.rational@ can underflow intermediate powers, including the+power in @5e-324@.+WKT also permits @.5@ and @1.@, which those readers do not consume fully.+-}+number :: Parser Double+number = do+    input <- get+    let (negative, unsigned) = case Text.uncons input of+            Just ('-', rest) -> (True, rest)+            Just ('+', rest) -> (False, rest)+            _ -> (False, input)+        (keyword, afterKeyword) = Text.span letter unsigned+        special = lookup (Text.toUpper keyword) [("NAN", castWord64ToDouble 0x7ff8000000000000), ("INF", 1 / 0), ("INFINITY", 1 / 0)]+    case special of+        Just value -> put afterKeyword >> pure (if negative then negate value else value)+        Nothing -> decimalNumber negative unsigned++-- | Round a decimal once. Clamp only exponents whose values must be zero or infinity.+decimalNumber :: Bool -> Text -> Parser Double+decimalNumber negative unsigned = do+    let (whole, afterWhole) = Text.span isDigit unsigned+        (fraction, afterFraction) = case Text.uncons afterWhole of+            Just ('.', rest) -> Text.span isDigit rest+            _ -> (Text.empty, afterWhole)+        -- Allow all mantissa digits to compensate for the exponent.+        -- The extra 400 exceeds Double's decimal range (-324 to 308).+        exponentLimit = toInteger (Text.length whole) + toInteger (Text.length fraction) + 400+    when (Text.null whole && Text.null fraction) (failure "expected a decimal number")+    (power, rest) <- case Text.uncons afterFraction of+        Just (marker, afterMarker) | marker == 'e' || marker == 'E' -> do+            let (negativeExponent, afterSign) = case Text.uncons afterMarker of+                    Just ('-', tailText) -> (True, tailText)+                    Just ('+', tailText) -> (False, tailText)+                    _ -> (False, afterMarker)+                (digits, afterDigits) = Text.span isDigit afterSign+            when (Text.null digits) (failure "expected exponent digits")+            let magnitude = Text.foldl' (\n c -> min (exponentLimit + 1) (10 * n + toInteger (fromEnum c - fromEnum '0'))) 0 digits+            pure (if negativeExponent then negate magnitude else magnitude, afterDigits)+        _ -> pure (0, afterFraction)+    if (Text.all (== '0') whole && Text.all (== '0') fraction) || power < negate exponentLimit+        then put rest >> pure (if negative then -0.0 else 0.0)+        else do+            let coefficient = digitsValue (whole <> fraction)+                adjustedPower = power - toInteger (Text.length fraction)+                magnitude+                    | power > exponentLimit = 1 / 0+                    | adjustedPower >= 0 = fromInteger (coefficient * 10 ^ adjustedPower)+                    | otherwise = fromRational (coefficient % (10 ^ negate adjustedPower))+                value = if negative then negate magnitude else magnitude+            put rest >> pure value++{- | Read decimal digits. Split long input in halves, because a digit-by-digit+loop over a large Integer takes quadratic time.+-}+digitsValue :: Text -> Integer+digitsValue digits+    | size <= 64 = Text.foldl' (\value c -> 10 * value + toInteger (fromEnum c - fromEnum '0')) 0 digits+    | otherwise = digitsValue high * 10 ^ (size - half) + digitsValue low+  where+    size = Text.length digits+    half = size `div` 2+    (high, low) = Text.splitAt half digits++-- | Whether a WKT container has one layout, no stored layout, or mixed layouts.+data LayoutSummary = Uniform Dimensions | Inherited | Mixed+    deriving (Eq)++-- | Combine child layouts. Containers without a stored layout are neutral.+combineLayout :: LayoutSummary -> LayoutSummary -> LayoutSummary+combineLayout Inherited second = second+combineLayout first Inherited = first+combineLayout (Uniform first) (Uniform second) | first == second = Uniform first+combineLayout _ _ = Mixed++-- | Compute layouts once and give empty containers their parent's output tag.+geometryWKT :: Geometry -> (LayoutSummary, Dimensions -> Builder)+geometryWKT geometry =+    ( layout+    , \inherited ->+        let dimensions = case layout of Uniform value -> value; Inherited -> inherited; Mixed -> DimXY+            suffix = case layout of+                Mixed -> ""+                _ -> case dimensions of DimXYZ -> " Z"; DimXYM -> " M"; DimXYZM -> " ZM"; DimXY -> ""+         in name <> suffix <> " " <> body dimensions+    )+  where+    sourceDimensions = geometryDimensions geometry+    (name, layout, body) = case geometry of+        PointGeometry value -> ("POINT", Uniform sourceDimensions, (`pointWKT` value))+        LineString points -> ("LINESTRING", Uniform sourceDimensions, (`coordinatesWKT` points))+        Polygon rings -> ("POLYGON", Uniform sourceDimensions, (`polygonWKT` rings))+        MultiPoint points -> ("MULTIPOINT", if U.null points then Inherited else Uniform sourceDimensions, \d -> sequenceWKT (pointWKT d) points)+        MultiLineString lineStrings -> ("MULTILINESTRING", if V.null lineStrings then Inherited else Uniform sourceDimensions, \d -> sequenceWKT (coordinatesWKT d) lineStrings)+        MultiPolygon polygons -> ("MULTIPOLYGON", if V.null polygons then Inherited else Uniform sourceDimensions, \d -> sequenceWKT (polygonWKT d) polygons)+        GeometryCollection members ->+            let children = V.map geometryWKT members+                common = V.foldl' (\acc (childLayout, _) -> combineLayout acc childLayout) Inherited children+             in ("GEOMETRYCOLLECTION", common, \d -> sequenceWKT (\(_, render) -> render d) children)++-- | Empty points have no ordinates in WKT. Nonempty points use the writer's layout.+pointWKT :: Dimensions -> Point -> Builder+pointWKT dimensions = fromMaybe "EMPTY" . withPoint (\value -> "(" <> coordinateWKT dimensions value <> ")")++-- | Empty polygons discard their empty holes only during writing.+polygonWKT :: Dimensions -> PolygonRings -> Builder+polygonWKT dimensions (PolygonRings shell holes)+    | coordinatesEmpty shell = "EMPTY"+    | otherwise = "(" <> coordinatesWKT dimensions shell <> V.foldMap (\ring -> ", " <> coordinatesWKT dimensions ring) holes <> ")"++-- | Render a sequence in the dimensions selected by its containing geometry.+coordinatesWKT :: Dimensions -> Coordinates -> Builder+coordinatesWKT dimensions = withCoordinates (sequenceWKT (coordinateWKT dimensions))++-- | Preserve finite bits with scientific notation and pad absent extra ordinates.+coordinateWKT :: forall c. (Coordinate c) => Dimensions -> c -> Builder+coordinateWKT target coordinateValue =+    let (x, y, z, m) = coordinateComponents coordinateValue+        source = coordinateDimensions (Proxy :: Proxy c)+        ordinate = RealFloat.formatDouble RealFloat.scientific+        zValue = if source == DimXYZ || source == DimXYZM then z else 0 / 0+        mValue = if source == DimXYM || source == DimXYZM then m else 0 / 0+        extra = case target of+            DimXY -> mempty+            DimXYZ -> " " <> ordinate zValue+            DimXYM -> " " <> ordinate mValue+            DimXYZM -> " " <> ordinate zValue <> " " <> ordinate mValue+     in ordinate x <> " " <> ordinate y <> extra++-- | Render a vector as EMPTY or a parenthesized sequence.+sequenceWKT :: (G.Vector v a) => (a -> Builder) -> v a -> Builder+sequenceWKT render values+    | G.null values = "EMPTY"+    | otherwise = "(" <> G.ifoldr (\i value rest -> separator i <> render value <> rest) mempty values <> ")"++-- | Separate elements after the first element.+separator :: Int -> Builder+separator 0 = mempty+separator _ = ", "