diff --git a/Project.toml b/Project.toml
index fb0713aa3..527dd0464 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Anshul Singhvi <anshulsinghvi@gmail.com> and contributors"]
 version = "0.0.1-DEV"
 
 [deps]
+ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index fec67aab3..40f140f14 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -6,6 +6,7 @@ using GeoInterface
 using GeometryBasics
 import Proj
 using LinearAlgebra
+import ExactPredicates
 
 using GeoInterface.Extents: Extents
 
@@ -21,6 +22,10 @@ include("methods/signed_area.jl")
 include("methods/centroid.jl")
 include("methods/intersects.jl")
 include("methods/contains.jl")
+include("methods/crosses.jl")
+include("methods/disjoint.jl")
+include("methods/overlaps.jl")
+include("methods/within.jl")
 include("methods/polygonize.jl")
 include("methods/barycentric.jl")
 
diff --git a/src/methods/bools.jl b/src/methods/bools.jl
index 2591c1b1d..aac4f8075 100644
--- a/src/methods/bools.jl
+++ b/src/methods/bools.jl
@@ -13,11 +13,14 @@ export line_on_line, line_in_polygon, polygon_in_polygon
 
 Take a ring and return true or false whether or not the ring is clockwise or counter-clockwise.
 
-## Examples
+## Example
+
 ```jldoctest
 import GeoInterface as GI, GeometryOps as GO
-line = GI.LineString([(0, 0), (1, 1), (1, 0), (0, 0)])
-GO.isclockwise(line)
+
+ring = GI.LinearRing([(0, 0), (1, 1), (1, 0), (0, 0)])
+GO.isclockwise(ring)
+
 # output
 true
 ```
@@ -43,8 +46,10 @@ Take a polygon and return true or false as to whether it is concave or not.
 ## Examples
 ```jldoctest
 import GeoInterface as GI, GeometryOps as GO
+
 poly = GI.Polygon([[(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]])
 GO.isconcave(poly)
+
 # output
 false
 ```
@@ -53,8 +58,9 @@ function isconcave(poly)::Bool
     sign = false
 
     exterior = GI.getexterior(poly)
-    GI.npoint(exterior) <= 4 && return false
 
+    # FIXME handle not closed polygons
+    GI.npoint(exterior) <= 4 && return false
     n = GI.npoint(exterior) - 1
 
     for i in 1:n
@@ -82,33 +88,32 @@ function isconcave(poly)::Bool
     return false
 end
 
+equals(geo1, geo2) = _equals(trait(geo1), geo1, trait(geo2), geo2)
 
-function equals(geo1, geo2)
-    GI.geomtrait(geo1) !== GI.geomtrait(geo2) && return false
-
-    GI.geomtrait(geo1) isa PointTrait && return compare_points(geo1, geo2)
-    GI.geomtrait(geo1) isa LineStringTrait && return compare_lines(geo1, geo2)
-
-    error("Cant compare $(GI.trait(geo1)) and $(GI.trait(geo2)) yet")
-end
-
-function compare_points(p1, p2)
-    length(p1) !== length(p2) && return false
-
-    for i in eachindex(p1)
-        round(p1[i]; digits=10) !== round(p2[i]; digits=10) && return false
+_equals(::T, geo1, ::T, geo2) where T = error("Cant compare $T yet")
+function _equals(::T, p1, ::T, p2) where {T<:PointTrait}
+    GI.ncoord(p1) == GI.ncoord(p2) || return false
+    GI.x(p1) == GI.x(p2) || return false
+    GI.y(p1) == GI.y(p2) || return false
+    if GI.is3d(p1)
+        GI.z(p1) == GI.z(p2) || return false 
     end
-
     return true
 end
+function _equals(::T, l1, ::T, l2) where {T<:AbstractCurveTrait}
+    # Check line lengths match
+    GI.npoint(l1) == GI.npoint(l2) || return false
 
-function compare_lines(p1::Vector, p2::Vector)
-    # TODO: complete this
-    length(p1[1]) !== length(p2[1]) && return false
+    # Then check all points are the same
+    for (p1, p2) in zip(GI.getpoint(l1), GI.getpoint(l2))
+        equals(p1, p2) || return false
+    end
+    return true
 end
+_equals(t1, geo1, t2, geo2) = false
 
 # """
-#     parallel(line1::LineString, line2::LineString)::Bool
+#     isparallel(line1::LineString, line2::LineString)::Bool
 
 # Return `true` if each segment of `line1` is parallel to the correspondent segment of `line2`
 
@@ -148,17 +153,20 @@ end
 
 
 """
-    point_on_line(point::Point, line::LineString, ignoreEndVertices::Bool=false)::Bool
+    point_on_line(point::Point, line::LineString; ignore_end_vertices::Bool=false)::Bool
 
 Return true if a point is on a line. Accept a optional parameter to ignore the
 start and end vertices of the linestring.
 
 ## Examples
+
 ```jldoctest
 import GeoInterface as GI, GeometryOps as GO
-point = GI.Point(1, 1)
+
+point = (1, 1)
 line = GI.LineString([(0, 0), (3, 3), (4, 4)])
 GO.point_on_line(point, line)
+
 # output
 true
 ```
@@ -167,25 +175,25 @@ function point_on_line(point, line; ignore_end_vertices::Bool=false)::Bool
     line_points = tuple_points(line)
     n = length(line_points)
 
-    ignore = :none
+    exclude_boundary = :none
     for i in 1:n - 1
-        if ignore_end_vertices == true
+        if ignore_end_vertices
             if i === 1 
-                ignore = :start
+                exclude_boundary = :start
             elseif i === n - 2 
-                ignore = :end
+                exclude_boundary = :end
             elseif (i === 1 && i + 1 === n - 1) 
-                ignore = :both
+                exclude_boundary = :both
             end
         end
-        if point_on_segment(line_points[i], line_points[i + 1], point, ignore) 
+        if point_on_segment(point, (line_points[i], line_points[i + 1]); exclude_boundary) 
             return true
         end
     end
     return false
 end
 
-function point_on_segment(start, stop, point, exclude_boundary::Symbol=:none)::Bool
+function point_on_segment(point, (start, stop); exclude_boundary::Symbol=:none)::Bool
     x, y = GI.x(point), GI.y(point)
     x1, y1 = GI.x(start), GI.y(start)
     x2, y2 = GI.x(stop), GI.y(stop)
@@ -225,46 +233,72 @@ function point_on_segment(start, stop, point, exclude_boundary::Symbol=:none)::B
 end
 
 """
-    point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignoreBoundary::Bool=false)::Bool
+    point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignore_boundary::Bool=false)::Bool
 
 Take a Point and a Polygon and determine if the point
 resides inside the polygon. The polygon can be convex or concave. The function accounts for holes.
 
 ## Examples
+
 ```jldoctest
 import GeoInterface as GI, GeometryOps as GO
+
 point = (-77.0, 44.0)
-poly = GI.Polygon([[[-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]])
+poly = GI.Polygon([[(-81, 41), (-81, 47), (-72, 47), (-72, 41), (-81, 41)]])
 GO.point_in_polygon(point, poly)
+
 # output
 true
 ```
 """
-function point_in_polygon(p, polygon, ignore_boundary::Bool=false)::Bool
-    GI.trait(polygon) isa PolygonTrait || throw(ArgumentError("Not a polygon"))
-    
-    point_in_extent(p, GI.extent(polygon)) || return false
-    point_in_ring(p, GI.getexterior(polygon), ignore_boundary) || return false
-
-    for ring in GI.gethole(polygon)
-        point_in_ring(pt, ring, !ignore_boundary) && return false
+point_in_polygon(point, polygon; kw...)::Bool =
+    point_in_polygon(GI.trait(point), point, GI.trait(polygon), polygon; kw...)
+function point_in_polygon(
+    ::PointTrait, point, 
+    ::PolygonTrait, poly; 
+    ignore_boundary::Bool=false,
+    check_extent::Bool=false,
+)::Bool
+    # Cheaply check that the point is inside the polygon extent
+    if check_extent
+        point_in_extent(point, GI.extent(poly)) || return false
+    end
+
+    # Then check the point is inside the exterior ring
+    point_in_polygon(point, GI.getexterior(poly); ignore_boundary, check_extent=false) || return false
+
+    # Finally make sure the point is not in any of the holes,
+    # flipping the boundary condition
+    for ring in GI.gethole(poly)
+        point_in_polygon(point, ring; ignore_boundary=!ignore_boundary) && return false
     end
     return true
 end
+function point_in_polygon(
+    ::PointTrait, pt, 
+    ::Union{LineStringTrait,LinearRingTrait}, ring; 
+    ignore_boundary::Bool=false,
+    check_extent::Bool=false,
+)::Bool
+    # Cheaply check that the point is inside the ring extent
+    if check_extent
+        point_in_extent(point, GI.extent(ring)) || return false
+    end
 
-function point_in_ring(pt, ring, ignore_boundary::Bool=false)
-    GI.trait(ring) isa Union{LineStringTrait,LinearRingTrait} || throw(ArgumentError("Not a ring"))
+    # Then check the point is inside the ring
     inside = false
     n = GI.npoint(ring)
-    p1 = first(GI.getpoint(ring))
+    p_start = GI.getpoint(ring, 1)
     p_end = GI.getpoint(ring, n)
 
-    l = if GI.x(p1) == GI.x(p_end) && GI.y(p1) == GI.y(p_end) 
-        l = n -1
+    # Handle closed on non-closed rings
+    l = if GI.x(p_start) == GI.x(p_end) && GI.y(p_start) == GI.y(p_end) 
+        l = n - 1
     else
         n
     end
 
+    # Loop over all points in the ring
     for i in 1:l - 1
         j = i + 1
 
@@ -292,103 +326,58 @@ function point_in_ring(pt, ring, ignore_boundary::Bool=false)
 end
 
 function point_in_extent(p, extent::Extents.Extent)
-    extent.X[1] <= GI.x(p) && extent.Y[1] <= GI.y(p) &&
-        extent.X[2] >= GI.x(p) && extent.Y[2] >= GI.y(p)
-end
-
-function line_in_polygon(poly, line)
-    out = false
-
-    polybox = bbox(poly)
-    linebox = bbox(line)
-
-    !(bboxOverlap(polybox, linebox)) && return false
-
-    coords = line.coordinates
-
-    for i in 1:length(coords) - 1
-        mid = [(coords[i][1] + coords[i + 1][1]) / 2, (coords[i][2] + coords[i + 1][2]) / 2]
-        if point_in_polygon(Point(mid), poly, true)
-            out = true
-            break
-        end
-    end
-    return out
+    (x1, x2), (y1, y1) = extent.X, extent.Y
+    return x1 <= GI.x(p) && y1 <= GI.y(p) && x2 >= GI.x(p) && y2 >= GI.y(p)
 end
 
 line_on_line(line1, line2) = line_on_line(trait(line1), line1, trait(line2), line2)
 function line_on_line(t1::GI.AbstractCurveTrait, line1, t2::AbstractCurveTrait, line2)
     for p in GI.getpoint(line1)
+        # FIXME: all points being on the line doesn't
+        # actually mean the whole line is on the line...
         point_on_line(p, line2) || return false
     end
     return true
 end
 
 line_in_polygon(line, poly) = line_in_polygon(trait(line), line, trait(poly), poly)
-function line_in_polygon(::LineStringTrait, line, ::PolygonTrait, poly)
-    polybox = bbox(poly)
-    linebox = bbox(line)
-
-    !(bboxOverlap(polybox, linebox)) && return false
+function line_in_polygon(
+    ::AbstractCurveTrait, line, 
+    ::Union{AbstractPolygonTrait,LinearRingTrait}, poly
+)
+    Extents.intersects(GI.extent(poly), GI.extent(line)) || return false
 
-    coords = line.coordinates
     inside = false
-
-    for i in 1:length(coords) - 1
-        !(point_in_polygon(Point(coords[i]), poly)) && return false
-        !inside && (inside = point_in_polygon(Point(coords[i]), poly, true))
+    for i in 1:GI.npoint(line) - 1
+        p = GI.getpoint(line, i)
+        p2 = GI.getpoint(line, i + 1)
+        point_in_polygon(p, poly) || return false
+        if !inside 
+            inside = point_in_polygon(p, poly; ignore_boundary=true)
+        end
+        # FIXME This seems like a hack, we should check for intersections rather than midpoint??
         if !inside
-            mid = [(coords[i][1] + coords[i + 1][1]) / 2, (coords[i][2] + coords[i + 1][2]) / 2]
-            inside = point_in_polygon(Point(mid), poly, true)
+            mid = ((GI.x(p) + GI.x(p2)) / 2, (GI.y(p) + GI.y(p2)) / 2)
+            inside = point_in_polygon(mid, poly; ignore_boundary=true)
         end
     end
     return inside
 end
 
-# TODO - why were there two methods for this in Turf.jl?
-function polygon_in_polygon(ft1, ft2, reverse::Bool=false)
-    polybox1 = bbox(ft1)
-    polybox2 = bbox(ft2)
-    coords = []
-
-    if reverse
-        !(bbox_overlap(polybox2, polybox1)) && return false
-
-        for point in GI.getpoint(ft1)
-            !(point_in_polygon(point, ft2)) && return false
-        end
-    else
-        !(bbox_overlap(polybox1, polybox2)) && return false
-
-        for point in GI.getpoint(ft2)
-            !(point_in_polygon(point, ft1)) && return false
-        end
-    end
-
-    return true
-end
-function poly_in_poly(poly1, poly2)
-
-    for point in GI.getpoint(poly1)
-        (point_in_polygon(point, poly2)) && return true
-    end
-
-    for point in GI.getpoint(poly2)
-        (point_in_polygon(point, poly1)) && return true
-    end
-
-    inter = line_intersects(polygon_to_line(poly1), polygon_to_line(poly2))
-    inter != nothing && return true
+function polygon_in_polygon(poly1, poly2)
+     # edges1, edges2 = to_edges(poly1), to_edges(poly2)
+     # extent1, extent2 = to_extent(edges1), to_extent(edges2)
+     # Check the extents intersect
+     Extents.intersects(GI.extent(poly1), GI.extent(poly2)) || return false
 
-    return false
+     # Check all points in poly1 are in poly2
+     for point in GI.getpoint(poly1)
+         point_in_polygon(point, poly2) || return false
+     end
 
-end
-
-function bbox_overlap(box1::Vector{T}, box2::Vector{T}) where {T <: Real}
-    box1[1] > box2[1] && return false
-    box1[3] < box2[3] && return false
-    box1[2] > box2[2] && return false
-    box1[4] < box2[4] && return false
-    return true
-end
+     # Check the line of poly1 does not intersect the line of poly2
+     line_intersects(poly1, poly2) && return false
 
+     # poly1 must be in poly2
+     return true
+ end
diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index 855c46af0..cfd37098d 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -14,7 +14,7 @@ export centroid
 # However, it's totally fine to ignore this and not have this code path.  
 # We simply need to decide on this.
 
-function centroid(ls::LineString{2, T}) where T
+function centroid(ls::GB.LineString{2, T}) where T
     centroid = Point{2, T}(0)
     total_area = T(0)
     if length(ls) == 1
@@ -34,7 +34,7 @@ function centroid(ls::LineString{2, T}) where T
 end
 
 # a more optimized function, so we only calculate signed area once!
-function centroid_and_signed_area(ls::LineString{2, T}) where T
+function centroid_and_signed_area(ls::GB.LineString{2, T}) where T
     centroid = Point{2, T}(0)
     total_area = T(0)
     if length(ls) == 1
@@ -53,7 +53,7 @@ function centroid_and_signed_area(ls::LineString{2, T}) where T
     return (centroid ./ total_area, total_area)
 end
 
-function centroid(poly::GeometryBasics.Polygon{2, T}) where T
+function centroid(poly::GB.Polygon{2, T}) where T
     exterior_centroid, exterior_area = centroid_and_signed_area(poly.exterior)
 
     total_area = exterior_area
@@ -68,7 +68,7 @@ function centroid(poly::GeometryBasics.Polygon{2, T}) where T
 
 end
 
-function centroid(multipoly::MultiPolygon)
+function centroid(multipoly::GB.MultiPolygon)
     centroids = centroid.(multipoly.polygons)
     areas = signed_area.(multipoly.polygons)
     areas ./= sum(areas)
@@ -77,10 +77,10 @@ function centroid(multipoly::MultiPolygon)
 end
 
 
-function centroid(rect::Rect{N, T}) where {N, T}
+function centroid(rect::GB.Rect{N, T}) where {N, T}
     return Point{N, T}(rect.origin .- rect.widths ./ 2)
 end
 
-function centroid(sphere::HyperSphere{N, T}) where {N, T}
+function centroid(sphere::GB.HyperSphere{N, T}) where {N, T}
     return sphere.center
-end
\ No newline at end of file
+end
diff --git a/src/methods/contains.jl b/src/methods/contains.jl
index 090031d7d..a58e2a48b 100644
--- a/src/methods/contains.jl
+++ b/src/methods/contains.jl
@@ -2,80 +2,24 @@
 
 export contains
 
-# This currently works for point-in-linestring or point-in-polygon.
-
-# More GeometryBasics code
-
-_cross(p1, p2, p3) = (GI.x(p1) - GI.x(p3)) * (GI.y(p2) - GI.y(p3)) - (GI.x(p2) - GI.x(p3)) * (GI.y(p1) - GI.y(p3))
-
-"""
-    contains(pointlist, point)::Bool
-
-Returns `true` if `point` is contained in `pointlist` (geometrically, not as a set)
-,and  `false` otherwise.
 """
-contains(pointlist, point) = contains(GI.trait(pointlist), GI.trait(point), pointlist, point)
+    contains(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool
 
-# Implementation of a point-in-polygon algorithm
-# from Luxor.jl.  This is the Hormann-Agathos (2001) algorithm.
+Return true if the second geometry is completely contained by the first geometry.
+The interiors of both geometries must intersect and, the interior and boundary of the secondary (geometry b)
+must not intersect the exterior of the primary (geometry a).
+`contains` returns the exact opposite result of `within`.
 
-# For the source, see [the code from Luxor.jl](https://github.com/JuliaGraphics/Luxor.jl/blob/66d60fb51f6b1bb38690fe8dcc6c0084eeb80710/src/polygons.jl#L190-L229).
+## Examples
 
-function contains(::Union{GI.LineStringTrait, GI.LinearRingTrait}, ::GI.PointTrait, pointlist, point)
-    n = GI.npoint(pointlist)
-    c = false
-    q1 = GI.getpoint(pointlist, 1)
-    q2 = GI.getpoint(pointlist, 1)
-    @inbounds for (counter, current_point) in enumerate(Iterators.drop(GI.getpoint(pointlist), 1))
-        q1 = q2
-        ## if reached last point, set "next point" to first point.
-        ##
-        if counter == (n-1)
-            q2 = GI.getpoint(pointlist, 1)
-        else
-            q2 = current_point
-        end
-        if GI.x(q1) == GI.x(point) && GI.x(q1) == GI.y(point)
-            ## allowonedge || error("isinside(): VertexException a")
-            continue
-        end
-        if GI.y(q2) == GI.y(point)
-            if GI.x(q2) == GI.x(point)
-                ## allowonedge || error("isinside(): VertexException b")
-                continue
-            elseif (GI.y(q1) == GI.y(point)) && ((GI.x(q2) > GI.x(point)) == (GI.x(q1) < GI.x(point)))
-                ## allowonedge || error("isinside(): EdgeException")
-                continue
-            end
-        end
-        if (GI.y(q1) < GI.y(point)) != (GI.y(q2) < GI.y(point)) # crossing
-            if GI.x(q1) >= GI.x(point)
-                if GI.x(q2) > GI.x(point)
-                    c = !c
-                elseif ((_cross(q1, q2, point) > 0) == (GI.y(q2) > GI.y(q1)))
-                    c = !c
-                end
-            elseif GI.x(q2) > GI.x(point)
-                if ((_cross(q1, q2, point) > 0) == (GI.y(q2) > GI.y(q1)))
-                    c = !c
-                end
-            end
-        end
-    end
-    return c
+```jldoctest
+import GeometryOps as GO, GeoInterface as GI
+line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
+point = (1, 2)
 
-end
-
-function contains(poly::Polygon{2, T1}, point::Point{2, T2}) where {T1, T2}
-    c = contains(poly.exterior, point)
-    for interior in poly.interiors
-        if contains(interior, point)
-            return false
-        end
-    end
-    return c
-end
-
-# TODOs: implement contains for mesh, simplex, and 3d objects (eg rect, triangle, etc.)
-
-contains(mp::MultiPolygon{2, T1}, point::Point{2, T2}) where {T1, T2} = any((contains(poly, point) for poly in mp.polygons))
+GO.contains(line, point)
+# output
+true
+```
+"""
+contains(g1, g2)::Bool = within(g2, g1)
diff --git a/src/methods/crosses.jl b/src/methods/crosses.jl
index 113886852..7c215c857 100644
--- a/src/methods/crosses.jl
+++ b/src/methods/crosses.jl
@@ -1,51 +1,51 @@
 # # Crossing checks
 
 """
-     crosses(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool
+     crosses(geom1, geom2)::Bool
 
 Return `true` if the intersection results in a geometry whose dimension is one less than
 the maximum dimension of the two source geometries and the intersection set is interior to
 both source geometries.
 
-## Examples
-```jldoctest
-julia> line = LineString([[1, 1], [1, 2], [1, 3], [1, 4]])
-LineString(Array{Float64,1}[[1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]])
+TODO: broken
 
-julia> line2 = LineString([[-2, 2], [4, 2]])
-LineString(Array{Float64,1}[[-2.0, 2.0], [4.0, 2.0]])
+## Examples 
+```julia
+import GeoInterface as GI, GeometryOps as GO
 
-julia> crosses(line2, line)
+line1 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
+line2 = GI.LineString([(-2, 2), (4, 2)])
+
+GO.crosses(line1, line2)
+# output
 true
 ```
 """
 crosses(g1, g2)::Bool = crosses(trait(g1), g1, trait(g2), g2)::Bool
 crosses(t1::FeatureTrait, g1, t2, g2)::Bool = crosses(GI.geometry(g1), g2)
 crosses(t1, g1, t2::FeatureTrait, g2)::Bool = crosses(g1, geometry(g2))
-crosses(::MultiPointTrait, g1::LineStringTrait, , g2)::Bool = multipoint_cross_line(g1, g2)
-crosses(::MultiPointTrait, g1::PolygonTrait, , g2)::Bool = multipoint_cross_poly(g1, g2)
-crosses(::LineStringTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_cross_lines(g2, g1)
-crosses(::LineStringTrait, g1, ::PolygonTrait, g2)::Bool = line_cross_poly(g1, g2)
-crosses(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = line_cross_line(g1, g2)
-crosses(::PolygonTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_cross_poly(g2, g1)
-crosses(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = line_cross_poly(g2, g1)
-
-function multipoint_cross_line(geom1, geom2)
+crosses(::MultiPointTrait, g1, ::LineStringTrait, g2)::Bool = multipoint_crosses_line(g1, g2)
+crosses(::MultiPointTrait, g1, ::PolygonTrait, g2)::Bool = multipoint_crosses_poly(g1, g2)
+crosses(::LineStringTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_crosses_lines(g2, g1)
+crosses(::LineStringTrait, g1, ::PolygonTrait, g2)::Bool = line_crosses_poly(g1, g2)
+crosses(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = line_crosses_line(g1, g2)
+crosses(::PolygonTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_crosses_poly(g2, g1)
+crosses(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = line_crosses_poly(g2, g1)
+
+function multipoint_crosses_line(geom1, geom2)
     int_point = false
     ext_point = false
     i = 1
     np2 = GI.npoint(geom2)
 
-    while i < GI.npoint(geom1) && !intPoint && !extPoint
+    while i < GI.npoint(geom1) && !int_point && !ext_point
         for j in 1:GI.npoint(geom2) - 1
-            inc_vertices = (j === 1 || j === np2 - 2) ? :none : :both
-
-            if is_point_on_segment(GI.getpoint(geom2, j), GI.getpoint(geom2.coordinates, j + 1), GI.getpoint(geom1, i), inc_vertices)
+            exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both
+            if point_on_segment(GI.getpoint(geom1, i), (GI.getpoint(geom2, j), GI.getpoint(geom2, j + 1)); exclude_boundary)
                 int_point = true
             else
                 ext_point = true
             end
-
         end
         i += 1
     end
@@ -53,32 +53,30 @@ function multipoint_cross_line(geom1, geom2)
     return int_point && ext_point
 end
 
-function line_cross_line(line1, line2)
-    inter = intersection(line1, line2)
-
+function line_crosses_line(line1, line2)
     np2 = GI.npoint(line2)
-    if !isnothing(inter)
+    if line_intersects(line1, line2; meets=MEETS_CLOSED)
         for i in 1:GI.npoint(line1) - 1
             for j in 1:GI.npoint(line2) - 1
-                inc_vertices = (j === 1 || j === np2 - 2) ? :none : :both
+                exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both
                 pa = GI.getpoint(line1, i)
                 pb = GI.getpoint(line1, i + 1)
                 p = GI.getpoint(line2, j)
-                is_point_on_segment(pa, pb, p, inc_vertices) && return true
+                point_on_segment(p, (pa, pb); exclude_boundary) && return true
             end
         end
     end
     return false
 end
 
-function line_cross_poly(line, poly) = 
-
-    for line in flatten(AbstractCurveTrait, poly)
-        intersects(line)
+function line_crosses_poly(line, poly)
+    for l in flatten(AbstractCurveTrait, poly)
+        line_intersects(line, l) && return true
     end
+    return false
 end
 
-function multipoint_cross_poly(mp, poly)
+function multipoint_crosses_poly(mp, poly)
     int_point = false
     ext_point = false
 
@@ -88,7 +86,7 @@ function multipoint_cross_poly(mp, poly)
         else
             ext_point = true
         end
-        in_point && ext_point && return true
+        int_point && ext_point && return true
     end
     return false
 end
diff --git a/src/methods/disjoint.jl b/src/methods/disjoint.jl
index 5b7889f3d..b51e5ab66 100644
--- a/src/methods/disjoint.jl
+++ b/src/methods/disjoint.jl
@@ -6,25 +6,37 @@
 Return `true` if the intersection of the two geometries is an empty set.
 
 # Examples
+
 ```jldoctest
-julia> poly = Polygon([[[-1, 2], [3, 2], [3, 3], [-1, 3], [-1, 2]]])
-Polygon(Array{Array{Float64,1},1}[[[-1.0, 2.0], [3.0, 2.0], [3.0, 3.0], [-1.0, 3.0], [-1.0, 2.0]]])
+import GeometryOps as GO, GeoInterface as GI
 
-julia> point = Point([1, 1])
-Point([1.0, 1.0])
+poly = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]])
+point = (1, 1)
+GO.disjoint(poly, point)
 
-julia> disjoint(poly, point)
+# output
 true
 ```
 """
-disjoint(t1::FeatureTrait, g1, t2, g2)::Bool = disjoint(GI.geometry(g1), g2)
-disjoint(t1, g1, t2::FeatureTrait, g2)::Bool = disjoint(g1, geometry(g2))
-disjoint(t1::PointTrait, g1, t2::PointTrait, g2)::Bool = !point_equals_point(g1, g2)
-disjoint(t1::PointTrait, g1, t2::LineStringTrait, g2)::Bool = !point_on_line(g1, g2)
-disjoint(t1::PointTrait, g1, t2::PolygonTrait, g2)::Bool = !point_in_polygon(g1, g2)
-disjoint(t1::LineStringTrait, g1, t2::PointTrait, g2)::Bool = !point_on_line(g2, g1)
-disjoint(t1::LineStringTrait, g1, t2::LineStringTrait, g2)::Bool = !line_on_line(g1, g2)
-disjoint(t1::LineStringTrait, g1, t2::PolygonTrait, g2)::Bool = !line_in_polygon(g2, g1)
-disjoint(t1::PolygonTrait, g1, t2::PointTrait, g2)::Bool = !point_in_polygon(g2, g1)
-disjoint(t1::PolygonTrait, g1, t2::LineStringTrait, g2)::Bool = !line_in_polygon(g2, g1)
-disjoint(t1::PolygonTrait, g1, t2::PolygonTrait, g2)::Bool = !poly_in_poly(g2, g1)
+disjoint(g1, g2)::Bool = disjoint(trait(g1), g1, trait(g2), g2)
+disjoint(::FeatureTrait, g1, ::Any, g2)::Bool = disjoint(GI.geometry(g1), g2)
+disjoint(::Any, g1, t2::FeatureTrait, g2)::Bool = disjoint(g1, geometry(g2))
+disjoint(::PointTrait, g1, ::PointTrait, g2)::Bool = !point_equals_point(g1, g2)
+disjoint(::PointTrait, g1, ::LineStringTrait, g2)::Bool = !point_on_line(g1, g2)
+disjoint(::PointTrait, g1, ::PolygonTrait, g2)::Bool = !point_in_polygon(g1, g2)
+disjoint(::LineStringTrait, g1, ::PointTrait, g2)::Bool = !point_on_line(g2, g1)
+disjoint(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = !line_on_line(g1, g2)
+disjoint(::LineStringTrait, g1, ::PolygonTrait, g2)::Bool = !line_in_polygon(g2, g1)
+disjoint(::PolygonTrait, g1, ::PointTrait, g2)::Bool = !point_in_polygon(g2, g1)
+disjoint(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = !line_in_polygon(g2, g1)
+disjoint(::PolygonTrait, g1, ::PolygonTrait, g2)::Bool = polygon_disjoint(g2, g1)
+
+function polygon_disjoint(poly1, poly2)
+    for point in GI.getpoint(poly1)
+        point_in_polygon(point, poly2) && return false
+    end
+    for point in GI.getpoint(poly2)
+        point_in_polygon(point, poly1) && return false
+    end
+    return !line_intersects(poly1, poly2)
+end
diff --git a/src/methods/intersects.jl b/src/methods/intersects.jl
index 44add99e0..e0476bd81 100644
--- a/src/methods/intersects.jl
+++ b/src/methods/intersects.jl
@@ -7,69 +7,94 @@ export intersects, intersection
 # !!! note
 #     This does not compute intersections, only checks if they exist.
 
+const MEETS_OPEN = 1
+const MEETS_CLOSED = 0
 
 """
-    intersects(line_a, line_b)
+    line_intersects(line_a, line_b)
 
 Check if `line_a` intersects with `line_b`.
 
 These can be `LineTrait`, `LineStringTrait` or `LinearRingTrait`
+
+## Example
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
+line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
+GO.line_intersects(line1, line2)
+
+# output
+true
+```
 """
-intersects(a, b) = isnothing(intersection) # Probably faster ways to do this
+line_intersects(a, b; kw...) = line_intersects(trait(a), a, trait(b), b; kw...)
+# Skip to_edges for LineTrait
+function line_intersects(::GI.LineTrait, a, ::GI.LineTrait, b; meets=MEETS_OPEN)
+    a1 = _tuple_point(GI.getpoint(a, 1))
+    b1 = _tuple_point(GI.getpoint(b, 1))
+    a2 = _tuple_point(GI.getpoint(a, 2))
+    b2 = _tuple_point(GI.getpoint(b, 2))
+    return ExactPredicates.meet(a1, a2, b1, b2) == meets
+end
+function line_intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b; kw...)
+    edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
+    return line_intersects(edges_a, edges_b; kw...)
+end
+function line_intersects(edges_a::Vector{Edge}, edges_b::Vector{Edge}; meets=MEETS_OPEN)
+    # Extents.intersects(to_extent(edges_a), to_extent(edges_b)) || return false
+    for edge_a in edges_a
+        for edge_b in edges_b
+            ExactPredicates.meet(edge_a..., edge_b...) == meets && return true 
+        end
+    end
+    return false
+end
 
 """
-    intersection(line_a, line_b)
+    line_intersection(line_a, line_b)
 
 Find a point that intersects LineStrings with two coordinates each.
 
 Returns `nothing` if no point is found.
 
-# Examples
+## Example
 
 ```jldoctest
-import GeoInterface as GI
-import GeometryOps as GO
+import GeoInterface as GI, GeometryOps as GO
+
 line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
 line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
-GO.intersection(line1, line2)
+GO.line_intersection(line1, line2)
+
 # output
 (125.58375366067547, -14.83572303404496)
 ```
 """
-intersection(line_a, line_b) = intersection(trait(line_a), line_a, trait(line_b), line_b)
-function intersection(
-    ::Union{LineStringTrait,LinearRingTrait}, line_a, 
-    ::Union{LineStringTrait,LinearRingTrait}, line_b,
-)
-    result = Tuple{Float64,Float64}[] # TODO handle 3d, and other Real ?
-    a1 = GI.getpoint(line_a, 1)
-    b1 = GI.getpoint(line_b, 1)
-
-    # TODO we can check all of these against the extent 
-    # of line_b and continue the loop if theyre outside
-    for i in 1:GI.npoint(line_a) - 1
-        for j in 1:GI.npoint(line_b) - 1
-            a2 = GI.getpoint(line_a, i + 1)
-            b2 = GI.getpoint(line_b, j + 1)
-            inter = _intersection((a1, a2), (b1, b2))
-            isnothing(inter) || push!(result, inter)
-            a1 = a2
-            b1 = b2
+line_intersection(line_a, line_b) = line_intersection(trait(line_a), line_a, trait(line_b), line_b)
+function line_intersection(::GI.AbstractTrait, a, ::GI.AbstractTrait, b)
+    Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing
+    result = Tuple{Float64,Float64}[]
+    edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
+    for edge_a in edges_a
+        for edge_b in edges_b
+            x = _line_intersection(edge_a, edge_b)
+            isnothing(x) || push!(result, x)
         end
     end
-    return unique!(result)
+    return result
 end
-
-function intersection(::LineTrait, line_a, ::LineTrait, line_b)
+function line_intersection(::GI.LineTrait, line_a, ::GI.LineTrait, line_b)
     a1 = GI.getpoint(line_a, 1)
     b1 = GI.getpoint(line_b, 1)
     a2 = GI.getpoint(line_a, 2)
     b2 = GI.getpoint(line_b, 2)
 
-    return _intersection((a1, a2), (b1, b2))
+    return _line_intersection((a1, a2), (b1, b2))
 end
-
-function _intersection((p11, p12)::Tuple, (p21, p22)::Tuple)
+function _line_intersection((p11, p12)::Tuple, (p21, p22)::Tuple)
     # Get points from lines
     x1, y1 = GI.x(p11), GI.y(p11) 
     x2, y2 = GI.x(p12), GI.y(p12)
diff --git a/src/methods/overlaps.jl b/src/methods/overlaps.jl
index 00b7d471f..b846e43de 100644
--- a/src/methods/overlaps.jl
+++ b/src/methods/overlaps.jl
@@ -15,19 +15,20 @@ Multipoint/Multipoint, MultiLineString/MultiLineString and MultiPolygon/MultiPol
 
 ## Examples
 ```jldoctest
-julia> poly1 = Polygon([[[0,0],[0,5],[5,5],[5,0],[0,0]]])
-Polygon(Array{Array{Float64,1},1}[[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]])
+import GeometryOps as GO, GeoInterface as GI
+poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
+poly2 = GI.Polygon([[(1,1), (1,6), (6,6), (6,1), (1,1)]])
 
-julia> poly2 = Polygon([[[1,1],[1,6],[6,6],[6,1],[1,1]]])
-Polygon(Array{Array{Float64,1},1}[[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]])
-
-julia> overlap(poly1, poly2)
+GO.overlaps(poly1, poly2)
+# output
 true
 ```
 """
 overlaps(g1, g2)::Bool = overlaps(trait(g1), g1, trait(g2), g2)::Bool
 overlaps(t1::FeatureTrait, g1, t2, g2)::Bool = overlaps(GI.geometry(g1), g2)
 overlaps(t1, g1, t2::FeatureTrait, g2)::Bool = overlaps(g1, geometry(g2))
+overlaps(t1::FeatureTrait, g1, t2::FeatureTrait, g2)::Bool = overlaps(geometry(g1), geometry(g2))
+overlaps(::PolygonTrait, mp, ::MultiPolygonTrait, p)::Bool = overlaps(p, mp)
 function overlaps(::MultiPointTrait, g1, ::MultiPointTrait, g2)::Bool 
     for p1 in GI.getpoint(g1)
         for p2 in GI.getpoint(g2)
@@ -36,19 +37,15 @@ function overlaps(::MultiPointTrait, g1, ::MultiPointTrait, g2)::Bool
     end
 end
 function overlaps(::PolygonTrait, g1, ::PolygonTrait, g2)::Bool
-    line1 = polygon_to_line(geom1)
-    line2 = polygon_to_line(geom2)
-
-    intersection(line1, line2)
+    return line_intersects(g1, g2)
 end
-overlaps(::PolygonTrait, mp, ::MultiPointTrait, p)::Bool = overlaps(p, mp)
-function overlaps(t1::MultiPolygonTrait, mp, t2::Polygon, p1)::Bool
+function overlaps(t1::MultiPolygonTrait, mp, t2::PolygonTrait, p1)::Bool
     for p2 in GI.getgeom(mp)
-        overlaps(p1, p2)
+        overlaps(p1, thp2) && return true
     end
 end
 function overlaps(::MultiPolygonTrait, g1, ::MultiPolygonTrait, g2)::Bool
     for p1 in GI.getgeom(g1)
-        overlaps(PolygonTrait(), mp, PolygonTrait(), p1)
+        overlaps(PolygonTrait(), mp, PolygonTrait(), p1) && return true
     end
 end
diff --git a/src/methods/signed_area.jl b/src/methods/signed_area.jl
index 7e202701f..ecd13ae61 100644
--- a/src/methods/signed_area.jl
+++ b/src/methods/signed_area.jl
@@ -45,7 +45,7 @@ signed_area(x) = signed_area(GI.trait(x), x)
 # TODOS here:
 # 1. This could conceivably be multithreaded.  How to indicate that it should be so?
 # 2. What to do for corner cases (nan point, etc)?
-function signed_area(::Union{LineStringTrait, LinearRingTrait}, geom)
+function signed_area(::Union{GI.LineStringTrait,GI.LinearRingTrait}, geom)
     # Basically, we integrate the area under the line string, which gives us
     # the signed area.
     point₁ = GI.getpoint(geom, 1)
@@ -63,7 +63,7 @@ function signed_area(::Union{LineStringTrait, LinearRingTrait}, geom)
 end
 
 # This subtracts the 
-function signed_area(::PolygonTrait, geom)
+function signed_area(::GI.PolygonTrait, geom)
     s_area = signed_area(GI.getexterior(geom))
     area = abs(s_area)
     for hole in GI.gethole(geom)
@@ -72,7 +72,7 @@ function signed_area(::PolygonTrait, geom)
     return area * sign(s_area)
 end
 
-signed_area(::MultiPolygonTrait, geom) = sum((signed_area(poly) for poly in GI.getpolygon(geom)))
+signed_area(::GI.MultiPolygonTrait, geom) = sum((signed_area(poly) for poly in GI.getpolygon(geom)))
 
 # This should _theoretically_ work for anything, but I haven't actually tested yet!
 
@@ -104,4 +104,4 @@ signed_area(::MultiPolygonTrait, geom) = sum((signed_area(poly) for poly in GI.g
 # # WARNING: this may not do what you expect, since it's
 # # sensitive to winding order.  Use GeoInterface.area instead.
 # signed_area(mp::MultiPolygon) = sum(signed_area.(mp.polygons))
-# ```
\ No newline at end of file
+# ```
diff --git a/src/methods/within.jl b/src/methods/within.jl
index d4ffe1ebb..c930ce62f 100644
--- a/src/methods/within.jl
+++ b/src/methods/within.jl
@@ -13,21 +13,21 @@ must not intersect the exterior of the secondary (geometry b).
 
 ## Examples
 ```jldoctest setup=:(using GeometryOps, GeometryBasics)
-julia> line = LineString([[1, 1], [1, 2], [1, 3], [1, 4]])
-LineString(Array{Float64,1}[[1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]])
+import GeometryOps as GO, GeoInterface as GI
 
-julia> point = Point([1, 2])
-Point([1.0, 2.0])
+line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
+point = (1, 2)
+GO.within(point, line)
 
-julia> within(point, line)
+# output
 true
 ```
 """
 within(g1, g2)::Bool = within(trait(g1), g1, trait(g2), g2)::Bool
-within(t1::FeatureTrait, g1, t2, g2)::Bool = within(GI.geometry(g1), g2)
-within(t1, g1, t2::FeatureTrait, g2)::Bool = within(g1, geometry(g2))
-within(t1::PointTrait, g1::LineStringTrait, t2, g2)::Bool = point_on_line(ft1, ft2, true)
-within(t1::PointTrait, g1, t2::PolygonTrait, g2)::Bool = point_in_polygon(ft1, ft2, true)
-within(t1::LineStringTrait, g1, t2::PolygonTrait, g2)::Bool = line_in_polygon(ft1, ft2)
-within(t1::LineStringTrait, g1, t2::LineStringTrait, g2)::Bool = line_on_line(ft1, ft2)
-within(t1::PolygonTrait, g1, t2::PolygonTrait, g2)::Bool = polygon_in_polygon(ft1, ft2, true)
+within(::GI.FeatureTrait, g1, ::Any, g2)::Bool = within(GI.geometry(g1), g2)
+within(::Any, g1, t2::GI.FeatureTrait, g2)::Bool = within(g1, geometry(g2))
+within(::GI.PointTrait, g1, ::GI.LineStringTrait, g2)::Bool = point_on_line(g1, g2; ignore_end_vertices=true)
+within(::GI.PointTrait, g1, ::GI.PolygonTrait, g2)::Bool = point_in_polygon(g1, g2; ignore_boundary=true)
+within(::GI.LineStringTrait, g1, ::GI.PolygonTrait, g2)::Bool = line_in_polygon(g1, g2)
+within(::GI.LineStringTrait, g1, ::GI.LineStringTrait, g2)::Bool = line_on_line(g1, g2)
+within(::GI.PolygonTrait, g1, ::GI.PolygonTrait, g2)::Bool = polygon_in_polygon(g1, g2)
diff --git a/src/utils.jl b/src/utils.jl
index b59940a78..b1fbef8ba 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -5,3 +5,130 @@ _is3d(::GI.AbstractGeometryTrait, geom) = GI.is3d(geom)
 _is3d(::GI.FeatureTrait, feature) = _is3d(GI.geometry(feature))
 _is3d(::GI.FeatureCollectionTrait, fc) = _is3d(GI.getfeature(fc, 1))
 _is3d(::Nothing, geom) = _is3d(first(geom)) # Otherwise step into an itererable
+
+_npoint(x) = _npoint(trait(x), x)
+_npoint(::Nothing, xs::AbstractArray) = sum(_npoint, xs)
+_npoint(::GI.FeatureCollectionTrait, fc) = sum(_npoint, GI.getfeature(fc))
+_npoint(::GI.FeatureTrait, f) = _npoint(GI.geometry(f))
+_npoint(::GI.AbstractGeometryTrait, x) = GI.npoint(trait(x), x)
+
+_nedge(x) = _nedge(trait(x), x)
+_nedge(::Nothing, xs::AbstractArray) = sum(_nedge, xs)
+_nedge(::GI.FeatureCollectionTrait, fc) = sum(_nedge, GI.getfeature(fc))
+_nedge(::GI.FeatureTrait, f) = _nedge(GI.geometry(f))
+function _nedge(::GI.AbstractGeometryTrait, x)
+    n = 0
+    for g in GI.getgeom(x)
+        n += _nedge(g)
+    end
+    return n
+end
+_nedge(::GI.AbstractCurveTrait, x) = GI.npoint(x) - 1
+_nedge(::GI.PointTrait, x) = error("Cant get edges from points")
+
+
+"""
+    polygon_to_line(poly::Polygon)
+
+Converts a Polygon to LineString or MultiLineString
+
+# Examples
+
+```jldoctest
+import GeometryOps as GO, GeoInterface as GI
+
+poly = GI.Polygon([[(-2.275543, 53.464547), (-2.275543, 53.489271), (-2.215118, 53.489271), (-2.215118, 53.464547), (-2.275543, 53.464547)]])
+GO.polygon_to_line(poly)
+# output
+GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(-2.275543, 53.464547), (-2.275543, 53.489271), (-2.215118, 53.489271), (-2.215118, 53.464547), (-2.275543, 53.464547)], nothing, nothing)
+```
+"""
+function polygon_to_line(poly)
+    @assert GI.trait(poly) isa PolygonTrait
+    GI.ngeom(poly) > 1 && return GI.MultiLineString(collect(GI.getgeom(poly)))
+    return GI.LineString(collect(GI.getgeom(GI.getgeom(poly, 1))))
+end
+
+
+const TuplePoint = Tuple{Float64,Float64}
+const Edge = Tuple{TuplePoint,TuplePoint}
+
+"""
+    to_edges()
+
+Convert any geometry or collection of geometries into a flat 
+vector of `Tuple{Tuple{Float64,Float64},{Float64,Float64}}` edges.
+"""
+function to_edges(x)
+    edges = Vector{Edge}(undef, _nedge(x))
+    _to_edges!(edges, x, 1)
+    return edges
+end
+
+_to_edges!(edges::Vector, x, n) = _to_edges!(edges, trait(x), x, n)
+function _to_edges!(edges::Vector, ::GI.FeatureCollectionTrait, fc, n)
+    for f in GI.getfeature(fc)
+        n = _to_edges!(edges, f, n)
+    end
+end
+_to_edges!(edges::Vector, ::GI.FeatureTrait, f, n) = _to_edges!(edges, GI.geometry(f), n)
+function _to_edges!(edges::Vector, ::GI.AbstractGeometryTrait, fc, n)
+    for f in GI.getgeom(fc)
+        n = _to_edges!(edges, f, n)
+    end
+end
+function _to_edges!(edges::Vector, ::GI.AbstractCurveTrait, geom, n)
+    p1 = GI.getpoint(geom, 1) 
+    p1x, p1y = GI.x(p1), GI.y(p1)
+    for i in 2:GI.npoint(geom)
+        p2 = GI.getpoint(geom, i)
+        p2x, p2y = GI.x(p2), GI.y(p2)
+        edges[n] = (p1x, p1y), (p2x, p2y)
+        p1x, p1y = p2x, p2y
+        n += 1
+    end
+    return n
+end
+
+_tuple_point(p) = GI.x(p), GI.y(p)
+
+function to_extent(edges::Vector{Edge})
+    x, y = extrema(first, edges)
+    Extents.Extent(X=x, Y=y)
+end
+
+function to_extent(edges::Vector{Edge})
+    x, y = extrema(first, edges)
+    Extents.Extent(X=x, Y=y)
+end
+
+function to_points(xs)
+    points = Vector{TuplePoint}(undef, _npoint(x))
+    _to_points!(points, x, 1)
+    return points
+end
+
+_to_points!(points::Vector, x, n) = _to_points!(points, trait(x), x, n)
+function _to_points!(points::Vector, ::FeatureCollectionTrait, fc, n)
+    for f in GI.getfeature(fc)
+        n = _to_points!(points, f, n)
+    end
+end
+_to_points!(points::Vector, ::FeatureTrait, f, n) = _to_points!(points, GI.geometry(f), n)
+function _to_points!(points::Vector, ::AbstractGeometryTrait, fc, n)
+    for f in GI.getgeom(fc)
+        n = _to_points!(points, f, n)
+    end
+end
+function _to_points!(points::Vector, ::Union{AbstractCurveTrait,MultiPointTrait}, geom, n)
+    p1 = GI.getpoint(geom, 1) 
+    p1x, p1y = GI.x(p1), GI.y(p1)
+    for i in 2:GI.npoint(geom)
+        p2 = GI.getpoint(geom, i)
+        p2x, p2y = GI.x(p2), GI.y(p2)
+        points[n] = (p1x, p1y), (p2x, p2y)
+        p1 = p2
+        n += 1
+    end
+    return n
+end
diff --git a/test/methods/bools.jl b/test/methods/bools.jl
new file mode 100644
index 000000000..b7650cb87
--- /dev/null
+++ b/test/methods/bools.jl
@@ -0,0 +1,151 @@
+using Test, GeometryOps
+import GeoInterface as GI
+import GeometryOps as GO
+
+@testset "booleans" begin
+    line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]])
+    line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]])
+    line3 = GI.LineString([
+        (-111.544189453125, 24.186847428521244),
+        (-110.687255859375, 24.966140159912975),
+        (-110.4510498046875, 24.467150664739002),
+        (-109.9951171875, 25.180087808990645)
+    ])
+	line4 = GI.LineString([
+        (-111.4617919921875, 24.05148034322011),
+        (-110.8795166015625, 24.681961205014595),
+        (-110.841064453125, 24.14174098050432),
+        (-109.97863769531249, 24.617057340809524)
+    ])
+
+    # @test isparallel(line1, line2) == true
+	# @test isparallel(line3, line4) == false
+
+	poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]])
+	poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]])
+
+	@test isconcave(poly1) == true
+	@test isconcave(poly2) == false
+
+	l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]])
+	l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]])
+
+	@test isclockwise(l1) == true
+	@test isclockwise(l2) == false
+
+	l3 = GI.LineString([[0, 0], [3, 3], [4, 4]])
+	p1 = GI.Point([1,1])
+
+	l4 = GI.LineString([[0, 0], [3, 3]])
+	p2 = GI.Point([0, 0])
+
+	p3 = GI.Point([20, 20])
+	l5 = GI.LineString([[0, 0], [3, 3], [38.32, 5.96]])
+
+	@test GO.point_on_line(p2, l4; ignore_end_vertices=true) == false
+	@test GO.point_on_line(p3, l5; ignore_end_vertices=true) == false
+	@test GO.point_on_line(p1, l3) == true
+
+	pt = (-77, 44)
+	poly = GI.Polygon([[[-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]])
+
+	@test point_in_polygon(pt, poly) == true
+
+	poly3 = GI.Polygon([[(1, 1), (1, 10), (10, 10), (10, 1), (1, 1)]])
+	poly4 = GI.Polygon([[(1, 1), (2, 2), (3, 2), (1, 1)]])
+	line5 = GI.LineString([(1.0, 1.0), (2.0, 3.0), (2.0, 3.5)])
+
+	line6 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)])
+	poly5 = GI.Polygon([[(1.0, 1.0), (1.0, 20.0), (1.0, 3.0), (1.0, 4.0), (1.0, 1.0)]])
+	line7 = GI.LineString([(1.0, 2.0), (1.0, 3.0), (1.0, 3.5)])
+
+	@test GO.contains(poly3, poly4) == true
+	@test GO.contains(poly3, line5) == true
+	@test GO.contains(line6, (1, 2)) == true
+	@test GO.contains(poly3, poly5) == false
+	@test GO.contains(poly3 , line7) == false
+
+	@test GO.within(poly4, poly3) == true
+	@test GO.within(line5, poly3) == true
+	@test GO.within(poly5, poly3) == false
+	@test GO.within((1, 2), line6) == true
+	@test GO.within(line7, poly3) == false
+
+	poly6 = GI.Polygon([[(-11, -12), (-13, -12), (-13, -13), (-11, -13), (-11, -12)]])
+	poly7 = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]])
+	poly8 = GI.Polygon([[(-1, 2), (-13, -12), (-13, -13), (-11, -13), (-1, 2)]])
+
+	@test GO.disjoint(poly7, poly6) == true
+	@test GO.disjoint(poly7, (1, 1)) == true
+	@test GO.disjoint(poly7, GI.LineString([(0, 0), (12, 2), (12, 3), (12, 4)])) == true
+	@test GO.disjoint(poly8, poly7) == false
+
+	line8 = GI.LineString([(124.584961, -12.768946), (126.738281, -17.224758)])
+	line9 = GI.LineString([(123.354492, -15.961329), (127.22168, -14.008696)])
+
+    @test all(GO.line_intersection(line8, line9)[1] .≈ (125.583754, -14.835723))
+
+	line10 = GI.LineString([
+        (142.03125, -11.695273),
+        (138.691406, -16.804541),
+        (136.40625, -14.604847),
+        (135.966797, -12.039321),
+        (131.308594, -11.436955),
+        (128.232422, -15.36895),
+        (125.947266, -13.581921),
+        (121.816406, -18.729502),
+        (117.421875, -20.632784),
+        (113.378906, -23.402765),
+        (114.169922, -26.667096),
+    ])
+	line11 = GI.LineString([
+        (117.861328, -15.029686),
+        (122.124023, -24.886436),
+        (132.583008, -22.309426),
+        (132.890625, -7.754537),
+    ])
+
+	points = GO.line_intersection(line10, line11)
+    @test all(points[1] .≈ (119.832884, -19.58857))
+    @test all(points[2] .≈ (132.808697, -11.6309378))
+
+	@test GO.crosses(GI.LineString([(-2, 2), (4, 2)]), line6) == true
+	@test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), poly7) == true
+	@test GO.crosses(GI.MultiPoint([(1, 2), (12, 12)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == true
+	@test GO.crosses(GI.MultiPoint([(1, 0), (12, 12)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == false
+	@test GO.crosses(GI.LineString([(-2, 2), (-4, 2)]), poly7) == false
+
+    pl1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]])
+	pl2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]])
+
+	@test GO.overlaps(pl1, pl2) == true
+	@test_throws MethodError GO.overlaps(pl1, (1, 1))
+	@test_throws MethodError GO.overlaps((1, 1), pl2)
+
+	pl3 = pl4 = GI.Polygon([[
+        (-53.57208251953125, 28.287451910503744),
+        (-53.33038330078125, 28.29228897739706),
+        (-53.34136962890625, 28.430052892335723),
+        (-53.57208251953125, 28.287451910503744),
+    ]])
+	@test GO.overlaps(pl3, pl4) == false
+
+	mp1 = GI.MultiPoint([
+        (-36.05712890625, 26.480407161007275),
+        (-35.7220458984375, 27.137368359795584),
+        (-35.13427734375, 26.83387451505858),
+        (-35.4638671875, 27.254629577800063),
+        (-35.5462646484375, 26.86328062676624),
+        (-35.3924560546875, 26.504988828743404)
+    ])
+	mp2 = GI.MultiPoint([
+        (-35.4638671875, 27.254629577800063),
+        (-35.5462646484375, 26.86328062676624),
+        (-35.3924560546875, 26.504988828743404),
+        (-35.2001953125, 26.12091815959972),
+        (-34.9969482421875, 26.455820238459893)
+    ])
+
+	@test GO.overlaps(mp1, mp2) == true
+	@test GO.overlaps(mp1, mp2) == GO.overlaps(mp2, mp1)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index dd6f4d964..a270b6f31 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,6 +10,7 @@ const AG = ArchGDAL
 
 @testset "GeometryOps.jl" begin
     @testset "Primitives" begin include("primitives.jl") end
+    @testset "Bools" begin include("methods/bools.jl") end
     @testset "Signed Area" begin include("methods/signed_area.jl") end
     @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
     @testset "Reproject" begin include("transformations/reproject.jl") end