diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -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. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/docs/GEOS-DIFFERENCES.md b/docs/GEOS-DIFFERENCES.md new file mode 100644 --- /dev/null +++ b/docs/GEOS-DIFFERENCES.md @@ -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. diff --git a/geometry-simple.cabal b/geometry-simple.cabal new file mode 100644 --- /dev/null +++ b/geometry-simple.cabal @@ -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, diff --git a/src/Data/Geometry.hs b/src/Data/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry.hs @@ -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 diff --git a/src/Data/Geometry/Internal.hs b/src/Data/Geometry/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Internal.hs @@ -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 diff --git a/src/Data/Geometry/SimpleFeatures.hs b/src/Data/Geometry/SimpleFeatures.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/SimpleFeatures.hs @@ -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) diff --git a/src/Data/Geometry/Topology/Buffer.hs b/src/Data/Geometry/Topology/Buffer.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Buffer.hs @@ -0,0 +1,296 @@ +{- | Round planar buffers with polygonal approximations of circular arcs. + +The arc and corner rules follow GEOS 3.13.1. See +. +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 +. +-} +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 +. +-} +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 +. +-} +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] diff --git a/src/Data/Geometry/Topology/Measures.hs b/src/Data/Geometry/Topology/Measures.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Measures.hs @@ -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) diff --git a/src/Data/Geometry/Topology/Overlay.hs b/src/Data/Geometry/Topology/Overlay.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Overlay.hs @@ -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 +. +-} +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 diff --git a/src/Data/Geometry/Topology/Planar.hs b/src/Data/Geometry/Topology/Planar.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Planar.hs @@ -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 diff --git a/src/Data/Geometry/Topology/Relations.hs b/src/Data/Geometry/Topology/Relations.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Relations.hs @@ -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)) diff --git a/src/Data/Geometry/Topology/Snapping.hs b/src/Data/Geometry/Topology/Snapping.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Snapping.hs @@ -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)) diff --git a/src/Data/Geometry/Topology/Unary.hs b/src/Data/Geometry/Topology/Unary.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/Topology/Unary.hs @@ -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 _ = [] diff --git a/src/Data/Geometry/WKB.hs b/src/Data/Geometry/WKB.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/WKB.hs @@ -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)) diff --git a/src/Data/Geometry/WKT.hs b/src/Data/Geometry/WKT.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geometry/WKT.hs @@ -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 _ = ", "