diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 85b1594e5..eaefcae5f 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,8 +18,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
+          - '1.6'
           - '1'
-          - 'nightly'
+          # - 'nightly'
         os:
           - ubuntu-latest
         arch:
diff --git a/Project.toml b/Project.toml
index 887eb4b81..fb0713aa3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,9 @@ version = "0.0.1-DEV"
 [deps]
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 GeometryBasics = "0.4.7"
diff --git a/README.md b/README.md
new file mode 100644
index 000000000..a18bd67a0
--- /dev/null
+++ b/README.md
@@ -0,0 +1,27 @@
+## GeometryOps.jl
+
+
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://asinghvi17.github.io/GeometryOps.jl/stable/)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://asinghvi17.github.io/GeometryOps.jl/dev/)
+[![Build Status](https://github.com/asinghvi17/GeometryOps.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/asinghvi17/GeometryOps.jl/actions/workflows/CI.yml?query=branch%3Amain)
+
+<img src="docs/src/assets/logo.png" alt="GeometryOps logo" width="250">
+
+GeometryOps.jl is a package for geometric calculations on (primarily 2D) geometries.
+
+The driving idea behind this package is to unify all the disparate packages for geometric calculations in Julia, and make them GeoInterface.jl-compatible. We seem to be focusing primarily on 2/2.5D geometries for now.
+
+Most of the usecases are driven by GIS and similar Earth data workflows, so this might be a bit specialized towards that, but methods should always be general to any coordinate space.
+
+## Methods 
+
+- Signed area, centroid, distance, etc
+- Iteration into geometries (`apply`)
+- Line and polygon simplification
+- Generalized barycentric coordinates in polygons
+
+### Planned additions
+
+- OGC methods (crosses, contains, intersects, etc)
+- Polygon union, intersection and clipping
+- Arclength interpolation (absolute and relative)
diff --git a/docs/make.jl b/docs/make.jl
index fc76d61db..d1ad2ed52 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,9 +2,20 @@ using GeometryOps
 using Documenter
 using Literate
 using Makie, CairoMakie
-CairoMakie.activate!(px_per_unit = 2, type = "png") # TODO: make this svg
+CairoMakie.activate!(px_per_unit = 2, type = "png", inline = true) # TODO: make this svg
 
-DocMeta.setdocmeta!(GeometryOps, :DocTestSetup, :(using GeometryOps); recursive=true)
+DocMeta.setdocmeta!(GeometryOps, :DocTestSetup, :(using GeometryOps; using GeometryOps.GeometryBasics); recursive=true)
+
+# First, remove any codecov files that may have been generated by the CI run
+for (root, dirs, files) in walkdir(dirname(@__DIR__)) # walk through `GeometryOps/*`
+    # Iterate through all files in the current directory
+    for file in files
+        # If the file is a codecov file, remove it
+        if splitext(file)[2] == ".cov"
+            rm(joinpath(root, file))
+        end
+    end
+end
 
 source_path = joinpath(dirname(@__DIR__), "src")
 output_path = joinpath(@__DIR__, "src", "source")
@@ -12,13 +23,15 @@ mkpath(output_path)
 
 literate_pages = String[]
 
-for (root_path, dirs, files) in walkdir(source_path)
-    output_dir = joinpath(output_path, relpath(root_path, source_path))
-    for file in files
-        # convert the source file to Markdown
-        Literate.markdown(joinpath(root_path, file), output_dir; documenter = false)
-        # TODO: make this respect nesting somehow!
-        push!(literate_pages, joinpath("source", relpath(root_path, source_path), splitext(file)[1] * ".md"))
+withenv("JULIA_DEBUG" => "Literate") do # allow Literate debug output to escape to the terminal!
+    for (root_path, dirs, files) in walkdir(source_path)
+        output_dir = joinpath(output_path, relpath(root_path, source_path))
+        for file in files
+            # convert the source file to Markdown
+            Literate.markdown(joinpath(root_path, file), output_dir; documenter = false)
+            # TODO: make this respect nesting somehow!
+            push!(literate_pages, joinpath("source", relpath(root_path, source_path), splitext(file)[1] * ".md"))
+        end
     end
 end
 
@@ -36,7 +49,8 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "Source code" => literate_pages
-    ]
+    ],
+    strict=false,
 )
 
 deploydocs(;
diff --git a/docs/src/assets/logo.png b/docs/src/assets/logo.png
new file mode 100644
index 000000000..8a3120119
Binary files /dev/null and b/docs/src/assets/logo.png differ
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 4031c384d..983b6a99e 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -1,18 +1,24 @@
+# # GeometryOps.jl
+
 module GeometryOps
 
 using GeoInterface
 using GeometryBasics
 import Proj
+using LinearAlgebra
 
 const GI = GeoInterface
 const GB = GeometryBasics
 
 include("primitives.jl")
 include("utils.jl")
+
 include("methods/signed_distance.jl")
 include("methods/signed_area.jl")
 include("methods/centroid.jl")
 include("methods/contains.jl")
+include("methods/barycentric.jl")
+
 include("transformations/flip.jl")
 include("transformations/simplify.jl")
 include("transformations/reproject.jl")
diff --git a/src/methods/barycentric.jl b/src/methods/barycentric.jl
new file mode 100644
index 000000000..8a91a701c
--- /dev/null
+++ b/src/methods/barycentric.jl
@@ -0,0 +1,454 @@
+# # Barycentric coordinates
+
+export barycentric_coordinates, barycentric_coordinates!, barycentric_interpolate
+export MeanValue
+
+# Generalized barycentric coordinates are a generalization of barycentric coordinates, 
+# which are typically used in triangles, to arbitrary polygons. 
+
+# They provide a way to express a point within a polygon as a weighted average 
+# of the polygon's vertices.
+
+# In the case of a triangle, barycentric coordinates are a set of three numbers 
+# $(λ_1, λ_2, λ_3)$, each associated with a vertex of the triangle. Any point within 
+# the triangle can be expressed as a weighted average of the vertices, where the 
+# weights are the barycentric coordinates. The weights sum to 1, and each is non-negative.
+
+# For a polygon with $n$ vertices, generalized barycentric coordinates are a set of 
+# $n$ numbers $(λ_1, λ_2, ..., λ_n)$, each associated with a vertex of the polygon. 
+# Any point within the polygon can be expressed as a weighted average of the vertices, 
+# where the weights are the generalized barycentric coordinates. 
+
+# As with the triangle case, the weights sum to 1, and each is non-negative.
+
+
+# ## Example
+# This example was taken from [this page of CGAL's documentation](https://doc.cgal.org/latest/Barycentric_coordinates_2/index.html).
+
+# ```@example barycentric
+# using GeometryOps, Makie
+# using GeometryOps.GeometryBasics
+# # Define a polygon
+# polygon_points = Point3f[
+# (0.03, 0.05, 0.00), (0.07, 0.04, 0.02), (0.10, 0.04, 0.04),
+# (0.14, 0.04, 0.06), (0.17, 0.07, 0.08), (0.20, 0.09, 0.10),
+# (0.22, 0.11, 0.12), (0.25, 0.11, 0.14), (0.27, 0.10, 0.16),
+# (0.30, 0.07, 0.18), (0.31, 0.04, 0.20), (0.34, 0.03, 0.22),
+# (0.37, 0.02, 0.24), (0.40, 0.03, 0.26), (0.42, 0.04, 0.28),
+# (0.44, 0.07, 0.30), (0.45, 0.10, 0.32), (0.46, 0.13, 0.34),
+# (0.46, 0.19, 0.36), (0.47, 0.26, 0.38), (0.47, 0.31, 0.40),
+# (0.47, 0.35, 0.42), (0.45, 0.37, 0.44), (0.41, 0.38, 0.46),
+# (0.38, 0.37, 0.48), (0.35, 0.36, 0.50), (0.32, 0.35, 0.52),
+# (0.30, 0.37, 0.54), (0.28, 0.39, 0.56), (0.25, 0.40, 0.58),
+# (0.23, 0.39, 0.60), (0.21, 0.37, 0.62), (0.21, 0.34, 0.64),
+# (0.23, 0.32, 0.66), (0.24, 0.29, 0.68), (0.27, 0.24, 0.70),
+# (0.29, 0.21, 0.72), (0.29, 0.18, 0.74), (0.26, 0.16, 0.76),
+# (0.24, 0.17, 0.78), (0.23, 0.19, 0.80), (0.24, 0.22, 0.82),
+# (0.24, 0.25, 0.84), (0.21, 0.26, 0.86), (0.17, 0.26, 0.88),
+# (0.12, 0.24, 0.90), (0.07, 0.20, 0.92), (0.03, 0.15, 0.94),
+# (0.01, 0.10, 0.97), (0.02, 0.07, 1.00)]
+# # Plot it!
+# # First, we'll plot the polygon using Makie's rendering:
+# f, a1, p1 = poly(
+#     polygon_points; 
+#     color = last.(polygon_points), colormap = cgrad(:jet, 18; categorical = true), 
+#     axis = (; 
+#         aspect = DataAspect(), title = "Makie mesh based polygon rendering", subtitle = "CairoMakie"
+#     ), 
+#     figure = (; resolution = (800, 400),)
+# )
+#
+# Makie.update_state_before_display!(f) # We have to call this explicitly, to get the axis limits correct
+# # Now that we've plotted the first polygon,
+# # we can render it using barycentric coordinates.
+# a1_bbox = a1.finallimits[] # First we get the extent of the axis
+# ext = GeometryOps.GI.Extent(NamedTuple{(:X, :Y)}(zip(minimum(a1_bbox), maximum(a1_bbox))))
+#
+# a2, p2box = poly( # Now, we plot a cropping rectangle around the axis so we only show the polygon
+#     f[1, 2], 
+#     GeometryOps.GeometryBasics.Polygon( # This is a rectangle with an internal hole shaped like the polygon.
+#         Point2f[(ext.X[1], ext.Y[1]), (ext.X[2], ext.Y[1]), (ext.X[2], ext.Y[2]), (ext.X[1], ext.Y[2]), (ext.X[1], ext.Y[1])], 
+#         [reverse(Point2f.(polygon_points))]
+#     ); 
+#     color = :white, xautolimits = false, yautolimits = false,
+#     axis = (; 
+#         aspect = DataAspect(), title = "Barycentric coordinate based polygon rendering", subtitle = "GeometryOps",
+#         limits = (ext.X, ext.Y), 
+#     )
+# )
+# hidedecorations!(a1)
+# hidedecorations!(a2)
+# cb = Colorbar(f[2, :], p1.plots[1]; vertical = false, flipaxis = true)
+# # Finally, we perform barycentric interpolation on a grid,
+# xrange = LinRange(ext.X..., widths(a2.scene.px_area[])[1] * 4) # 2 rendered pixels per "physical" pixel
+# yrange = LinRange(ext.Y..., widths(a2.scene.px_area[])[2] * 4) # 2 rendered pixels per "physical" pixel
+# @time mean_values = barycentric_interpolate.(
+#     (MeanValue(),), # The barycentric coordinate algorithm (MeanValue is the only one for now)
+#     (Point2f.(polygon_points),), # The polygon points as `Point2f`
+#     (last.(polygon_points,),),   # The values per polygon point - can be anything which supports addition and division
+#     Point2f.(xrange, yrange')    # The points at which to interpolate
+# )
+# # and render!
+# hm = heatmap!(
+#     a2, xrange, yrange, mean_values;
+#     colormap = p1.colormap, # Use the same colormap as the original polygon plot
+#     colorrange = p1.plots[1].colorrange[], # Access the rendered mesh plot's colorrange directly
+#     transformation = (; translation = Vec3f(0,0,-1)), # This gets the heatmap to render "behind" the previously plotted polygon
+#     xautolimits = false, yautolimits = false
+# )
+# f
+# ```
+
+# ## Barycentric-coordinate API
+# In some cases, we actually want barycentric interpolation, and have no interest 
+# in the coordinates themselves.  
+# 
+# However, the coordinates can be useful for debugging, and when performing 3D rendering,
+# multiple barycentric values (depth, uv) are needed for depth buffering.
+#
+
+const _VecTypes = Union{Tuple{Vararg{T, N}}, GeometryBasics.StaticArraysCore.StaticArray{Tuple{N}, T, 1}} where {N, T}
+
+"""
+    abstract type AbstractBarycentricCoordinateMethod
+
+Abstract supertype for barycentric coordinate methods.  
+The subtypes may serve as dispatch types, or may cache 
+some information about the target polygon.  
+
+## API
+The following methods must be implemented for all subtypes:
+- `barycentric_coordinates!(λs::Vector{<: Real}, method::AbstractBarycentricCoordinateMethod, exterior::Vector{<: Point{2, T1}}, point::Point{2, T2})`
+- `barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, exterior::Vector{<: Point{2, T1}}, values::Vector{V}, point::Point{2, T2})::V`
+- `barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, exterior::Vector{<: Point{2, T1}}, interiors::Vector{<: Vector{<: Point{2, T1}}} values::Vector{V}, point::Point{2, T2})::V`
+The rest of the methods will be implemented in terms of these, and have efficient dispatches for broadcasting.
+"""
+abstract type AbstractBarycentricCoordinateMethod end
+
+
+Base.@propagate_inbounds function barycentric_coordinates!(λs::Vector{<: Real}, method::AbstractBarycentricCoordinateMethod, polypoints::AbstractVector{<: Point{N1, T1}}, point::Point{N2, T2}) where {N1, N2, T1 <: Real, T2 <: Real}
+    @boundscheck @assert length(λs) == length(polypoints)
+    @boundscheck @assert length(polypoints) >= 3
+
+    @error("Not implemented yet for method $(method).")
+end
+Base.@propagate_inbounds barycentric_coordinates!(λs::Vector{<: Real}, polypoints::AbstractVector{<: Point{N1, T1}}, point::Point{N2, T2}) where {N1, N2, T1 <: Real, T2 <: Real} = barycentric_coordinates!(λs, MeanValue(), polypoints, point)
+
+Base.@propagate_inbounds function barycentric_coordinates(method::AbstractBarycentricCoordinateMethod, polypoints::AbstractVector{<: Point{N1, T1}}, point::Point{N2, T2}) where {N1, N2, T1 <: Real, T2 <: Real}
+    λs = zeros(promote_type(T1, T2), length(polypoints))
+    barycentric_coordinates!(λs, method, polypoints, point)
+    return λs
+end
+Base.@propagate_inbounds barycentric_coordinates(polypoints::AbstractVector{<: Point{N1, T1}}, point::Point{N2, T2}) where {N1, N2, T1 <: Real, T2 <: Real} = barycentric_coordinates(MeanValue(), polypoints, point)
+
+Base.@propagate_inbounds function barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, polypoints::AbstractVector{<: Point{N, T1}}, values::AbstractVector{V}, point::Point{N, T2}) where {N, T1 <: Real, T2 <: Real, V}
+    @boundscheck @assert length(values) == length(polypoints)
+    @boundscheck @assert length(polypoints) >= 3
+    λs = barycentric_coordinates(method, polypoints, point)
+    return sum(λs .* values)
+end
+Base.@propagate_inbounds barycentric_interpolate(polypoints::AbstractVector{<: Point{N, T1}}, values::AbstractVector{V}, point::Point{N, T2}) where {N, T1 <: Real, T2 <: Real, V} = barycentric_interpolate(MeanValue(), polypoints, values, point)
+
+Base.@propagate_inbounds function barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, exterior::AbstractVector{<: Point{N, T1}}, interiors::AbstractVector{<: Point{N, T1}}, values::AbstractVector{V}, point::Point{N, T2}) where {N, T1 <: Real, T2 <: Real, V}
+    @boundscheck @assert length(values) == length(exterior) + isempty(interiors) ? 0 : sum(length.(interiors))
+    @boundscheck @assert length(exterior) >= 3
+    λs = barycentric_coordinates(method, exterior, interiors, point)
+    return sum(λs .* values)
+end
+Base.@propagate_inbounds barycentric_interpolate(exterior::AbstractVector{<: Point{N, T1}}, interiors::AbstractVector{<: Point{N, T1}}, values::AbstractVector{V}, point::Point{N, T2}) where {N, T1 <: Real, T2 <: Real, V} = barycentric_interpolate(MeanValue(), exterior, interiors, values, point)
+
+Base.@propagate_inbounds function barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, polygon::Polygon{2, T1}, values::AbstractVector{V}, point::Point{2, T2}) where {T1 <: Real, T2 <: Real, V}
+    exterior = decompose(Point{2, promote_type(T1, T2)}, polygon.exterior)
+    if isempty(polygon.interiors)
+        @boundscheck @assert length(values) == length(exterior)
+        return barycentric_interpolate(method, exterior, values, point)
+    else # the poly has interiors
+        interiors = reverse.(decompose.((Point{2, promote_type(T1, T2)},), polygon.interiors))
+        @boundscheck @assert length(values) == length(exterior) + sum(length.(interiors))
+        return barycentric_interpolate(method, exterior, interiors, values, point)
+    end
+end
+Base.@propagate_inbounds barycentric_interpolate(polygon::Polygon{2, T1}, values::AbstractVector{V}, point::Point{2, T2}) where {T1 <: Real, T2 <: Real, V} = barycentric_interpolate(MeanValue(), polygon, values, point)
+
+# 3D polygons are considered to have their vertices in the XY plane, 
+# and the Z coordinate must represent some value.  This is to say that
+# the Z coordinate is interpreted as an M coordinate.
+Base.@propagate_inbounds function barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, polygon::Polygon{3, T1}, point::Point{2, T2}) where {T1 <: Real, T2 <: Real}
+    exterior_point3s = decompose(Point{3, promote_type(T1, T2)}, polygon.exterior)
+    exterior_values = getindex.(exterior_point3s, 3)
+    exterior_points = Point2f.(exterior_point3s)
+    if isempty(polygon.interiors)
+        return barycentric_interpolate(method, exterior_points, exterior_values, point)
+    else # the poly has interiors
+        interior_point3s = decompose.((Point{3, promote_type(T1, T2)},), polygon.interiors)
+        interior_values = collect(Iterators.flatten((getindex.(point3s, 3) for point3s in interior_point3s)))
+        interior_points = map(point3s -> Point2f.(point3s), interior_point3s)
+        return barycentric_interpolate(method, exterior_points, interior_points, vcat(exterior_values, interior_values), point)
+    end
+end
+Base.@propagate_inbounds barycentric_interpolate(polygon::Polygon{3, T1}, point::Point{2, T2}) where {T1 <: Real, T2 <: Real} = barycentric_interpolate(MeanValue(), polygon, point)
+
+# This method is the one which supports GeoInterface.
+Base.@propagate_inbounds function barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, polygon, values::AbstractVector{V}, point) where V
+    @assert GeoInterface.trait(polygon) isa GeoInterface.PolygonTrait
+    @assert GeoInterface.trait(point) isa GeoInterface.PointTrait
+    passable_polygon = GeoInterface.convert(GeometryBasics, polygon)
+    @assert passable_polygon isa GeometryBasics.Polygon "The polygon was converted to a $(typeof(passable_polygon)), which is not a `GeometryBasics.Polygon`."
+    ## first_poly_point = GeoInterface.getpoint(GeoInterface.getexterior(polygon))
+    passable_point = GeoInterface.convert(GeometryBasics, point)
+    return barycentric_interpolate(method, passable_polygon, Point2(passable_point))
+end
+Base.@propagate_inbounds barycentric_interpolate(polygon, values::AbstractVector{V}, point) where V = barycentric_interpolate(MeanValue(), polygon, values, point)
+
+"""
+    weighted_mean(weight::Real, x1, x2)
+
+Returns the weighted mean of `x1` and `x2`, where `weight` is the weight of `x1`.
+
+Specifically, calculates `x1 * weight + x2 * (1 - weight)`.
+
+!!! note
+    The idea for this method is that you can override this for custom types, like Color types, in extension modules.
+"""
+function weighted_mean(weight::WT, x1, x2) where {WT <: Real}
+    return muladd(x1, weight, x2 * (oneunit(WT) - weight))
+end
+
+
+"""
+    MeanValue() <: AbstractBarycentricCoordinateMethod
+
+This method calculates barycentric coordinates using the mean value method.
+
+## References
+
+"""
+struct MeanValue <: AbstractBarycentricCoordinateMethod 
+end
+
+# Before we go to the actual implementation, there are some quick and simple utility functions
+# that we need to implement.  These are mainly for convenience and code brevity.
+
+"""
+    _det(s1::Point2{T1}, s2::Point2{T2}) where {T1 <: Real, T2 <: Real}
+
+Returns the determinant of the matrix formed by `hcat`'ing two points `s1` and `s2`.
+
+Specifically, this is: 
+```julia
+s1[1] * s2[2] - s1[2] * s2[1]
+```
+"""
+function _det(s1::_VecTypes{2, T1}, s2::_VecTypes{2, T2}) where {T1 <: Real, T2 <: Real}
+    return s1[1] * s2[2] - s1[2] * s2[1]
+end
+
+"""
+    t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)
+
+Returns the "T-value" as described in Hormann's presentation [^HormannPresentation] on how to calculate
+the mean-value coordinate.  
+
+Here, `sᵢ` is the vector from vertex `vᵢ` to the point, and `rᵢ` is the norm (length) of `sᵢ`.
+`s` must be `Point` and `r` must be real numbers.
+
+```math
+tᵢ = \\frac{\\mathrm{det}\\left(sᵢ, sᵢ₊₁\\right)}{rᵢ * rᵢ₊₁ + sᵢ ⋅ sᵢ₊₁}
+```
+
+[^HormannPresentation]: K. Hormann and N. Sukumar. Generalized Barycentric Coordinates in Computer Graphics and Computational Mechanics. Taylor & Fancis, CRC Press, 2017.
+```
+
+"""
+function t_value(sᵢ::_VecTypes{N, T1}, sᵢ₊₁::_VecTypes{N, T1}, rᵢ::T2, rᵢ₊₁::T2) where {N, T1 <: Real, T2 <: Real}
+    return _det(sᵢ, sᵢ₊₁) / muladd(rᵢ, rᵢ₊₁, dot(sᵢ, sᵢ₊₁))
+end
+
+
+function barycentric_coordinates!(λs::Vector{<: Real}, ::MeanValue, polypoints::AbstractVector{<: Point{2, T1}}, point::Point{2, T2}) where {T1 <: Real, T2 <: Real}
+    @boundscheck @assert length(λs) == length(polypoints)
+    @boundscheck @assert length(polypoints) >= 3
+    n_points = length(polypoints)
+    ## Initialize counters and register variables
+    ## Points - these are actually vectors from point to vertices
+    ##  polypoints[i-1], polypoints[i], polypoints[i+1]
+    sᵢ₋₁ = polypoints[end] - point
+    sᵢ   = polypoints[begin] - point
+    sᵢ₊₁ = polypoints[begin+1] - point
+    ## radius / Euclidean distance between points.
+    rᵢ₋₁ = norm(sᵢ₋₁) 
+    rᵢ   = norm(sᵢ  )
+    rᵢ₊₁ = norm(sᵢ₊₁)
+    ## Perform the first computation explicitly, so we can cut down on 
+    ## a mod in the loop.
+    λs[1] = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+    ## Loop through the rest of the vertices, compute, store in λs
+    for i in 2:n_points
+        ## Increment counters + set variables
+        sᵢ₋₁ = sᵢ
+        sᵢ   = sᵢ₊₁
+        sᵢ₊₁ = polypoints[mod1(i+1, n_points)] - point
+        rᵢ₋₁ = rᵢ
+        rᵢ   = rᵢ₊₁
+        rᵢ₊₁ = norm(sᵢ₊₁) # radius / Euclidean distance between points.
+        λs[i] = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ
+    end
+    ## Normalize λs to the 1-norm (sum=1)
+    λs ./= sum(λs) 
+    return λs
+end
+
+# ```julia
+# function barycentric_coordinates(::MeanValue, polypoints::NTuple{N, Point{2, T2}}, point::Point{2, T1},) where {N, T1, T2}
+#     ## Initialize counters and register variables
+#     ## Points - these are actually vectors from point to vertices
+#     ##  polypoints[i-1], polypoints[i], polypoints[i+1]
+#     sᵢ₋₁ = polypoints[end] - point
+#     sᵢ   = polypoints[begin] - point
+#     sᵢ₊₁ = polypoints[begin+1] - point
+#     ## radius / Euclidean distance between points.
+#     rᵢ₋₁ = norm(sᵢ₋₁) 
+#     rᵢ   = norm(sᵢ  )
+#     rᵢ₊₁ = norm(sᵢ₊₁)
+#     λ₁ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+#     λs = ntuple(N) do i
+#         if i == 1
+#             return λ₁
+#         end
+#         ## Increment counters + set variables
+#         sᵢ₋₁ = sᵢ
+#         sᵢ   = sᵢ₊₁
+#         sᵢ₊₁ = polypoints[mod1(i+1, N)] - point
+#         rᵢ₋₁ = rᵢ
+#         rᵢ   = rᵢ₊₁
+#         rᵢ₊₁ = norm(sᵢ₊₁) # radius / Euclidean distance between points.
+#         return (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ
+#     end
+#
+#     ∑λ = sum(λs)
+#
+#     return ntuple(N) do i
+#         λs[i] / ∑λ
+#     end
+# end
+# ```
+
+# This performs an inplace accumulation, using less memory and is faster.
+# That's particularly good if you are using a polygon with a large number of points...
+function barycentric_interpolate(::MeanValue, polypoints::AbstractVector{<: Point{2, T1}}, values::AbstractVector{V}, point::Point{2, T2}) where {T1 <: Real, T2 <: Real, V}
+    @boundscheck @assert length(values) == length(polypoints)
+    @boundscheck @assert length(polypoints) >= 3
+    
+    n_points = length(polypoints)
+    ## Initialize counters and register variables
+    ## Points - these are actually vectors from point to vertices
+    ##  polypoints[i-1], polypoints[i], polypoints[i+1]
+    sᵢ₋₁ = polypoints[end] - point
+    sᵢ   = polypoints[begin] - point
+    sᵢ₊₁ = polypoints[begin+1] - point
+    ## radius / Euclidean distance between points.
+    rᵢ₋₁ = norm(sᵢ₋₁) 
+    rᵢ   = norm(sᵢ  )
+    rᵢ₊₁ = norm(sᵢ₊₁)
+    ## Now, we set the interpolated value to the first point's value, multiplied
+    ## by the weight computed relative to the first point in the polygon. 
+    wᵢ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+    wₜₒₜ = wᵢ
+    interpolated_value = values[begin] * wᵢ
+    for i in 2:n_points
+        ## Increment counters + set variables
+        sᵢ₋₁ = sᵢ
+        sᵢ   = sᵢ₊₁
+        sᵢ₊₁ = polypoints[mod1(i+1, n_points)] - point
+        rᵢ₋₁ = rᵢ
+        rᵢ   = rᵢ₊₁
+        rᵢ₊₁ = norm(sᵢ₊₁) 
+        ## Now, we calculate the weight:
+        wᵢ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+        ## perform a weighted sum with the interpolated value:
+        interpolated_value += values[i] * wᵢ
+        ## and add the weight to the total weight accumulator.
+        wₜₒₜ += wᵢ
+    end
+    ## Return the normalized interpolated value.
+    return interpolated_value / wₜₒₜ
+end
+
+# When you have holes, then you have to be careful 
+# about the order you iterate around points.
+
+# Specifically, you have to iterate around each linear ring separately 
+# and ensure there are no degenerate/repeated points at the start and end!
+
+function barycentric_interpolate(::MeanValue, exterior::AbstractVector{<: Point{N, T1}}, interiors::AbstractVector{<: AbstractVector{<: Point{N, T1}}}, values::AbstractVector{V}, point::Point{N, T2}) where {N, T1 <: Real, T2 <: Real, V}
+    ## @boundscheck @assert length(values) == (length(exterior) + isempty(interiors) ? 0 : sum(length.(interiors)))
+    ## @boundscheck @assert length(exterior) >= 3
+
+    current_index = 1
+    l_exterior = length(exterior)
+
+    sᵢ₋₁ = exterior[end] - point
+    sᵢ   = exterior[begin] - point
+    sᵢ₊₁ = exterior[begin+1] - point
+    rᵢ₋₁ = norm(sᵢ₋₁) # radius / Euclidean distance between points.
+    rᵢ   = norm(sᵢ  ) # radius / Euclidean distance between points.
+    rᵢ₊₁ = norm(sᵢ₊₁) # radius / Euclidean distance between points.
+    # Now, we set the interpolated value to the first point's value, multiplied
+    # by the weight computed relative to the first point in the polygon. 
+    wᵢ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+    wₜₒₜ = wᵢ
+    interpolated_value = values[begin] * wᵢ
+
+    for i in 2:l_exterior
+        # Increment counters + set variables
+        sᵢ₋₁ = sᵢ
+        sᵢ   = sᵢ₊₁
+        sᵢ₊₁ = exterior[mod1(i+1, l_exterior)] - point
+        rᵢ₋₁ = rᵢ
+        rᵢ   = rᵢ₊₁
+        rᵢ₊₁ = norm(sᵢ₊₁) # radius / Euclidean distance between points.
+        wᵢ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+        # Updates - first the interpolated value,
+        interpolated_value += values[current_index] * wᵢ
+        # then the accumulators for total weight and current index.
+        wₜₒₜ += wᵢ
+        current_index += 1
+
+    end
+    for hole in interiors
+        l_hole = length(hole)
+        sᵢ₋₁ = hole[end] - point
+        sᵢ   = hole[begin] - point
+        sᵢ₊₁ = hole[begin+1] - point
+        rᵢ₋₁ = norm(sᵢ₋₁) # radius / Euclidean distance between points.
+        rᵢ   = norm(sᵢ  ) # radius / Euclidean distance between points.
+        rᵢ₊₁ = norm(sᵢ₊₁) # radius / Euclidean distance between points.
+        ## Now, we set the interpolated value to the first point's value, multiplied
+        ## by the weight computed relative to the first point in the polygon. 
+        wᵢ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ
+
+        interpolated_value += values[current_index] * wᵢ
+
+        wₜₒₜ += wᵢ
+        current_index += 1
+    
+        for i in 2:l_hole
+            ## Increment counters + set variables
+            sᵢ₋₁ = sᵢ
+            sᵢ   = sᵢ₊₁
+            sᵢ₊₁ = hole[mod1(i+1, l_hole)] - point
+            rᵢ₋₁ = rᵢ
+            rᵢ   = rᵢ₊₁
+            rᵢ₊₁ = norm(sᵢ₊₁) ## radius / Euclidean distance between points.
+            wᵢ = (t_value(sᵢ₋₁, sᵢ, rᵢ₋₁, rᵢ) + t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)) / rᵢ 
+            interpolated_value += values[current_index] * wᵢ
+            wₜₒₜ += wᵢ
+            current_index += 1
+        end
+    end
+    return interpolated_value / wₜₒₜ
+
+end
+
+struct Wachspress <: AbstractBarycentricCoordinateMethod
+end
diff --git a/src/methods/centroid.jl b/src/methods/centroid.jl
index df78f8360..855c46af0 100644
--- a/src/methods/centroid.jl
+++ b/src/methods/centroid.jl
@@ -1,3 +1,5 @@
+# # Centroid
+
 export centroid
 
 # These are all GeometryBasics.jl methods so far.
diff --git a/src/methods/contains.jl b/src/methods/contains.jl
index 2fe1f929b..5d978e89e 100644
--- a/src/methods/contains.jl
+++ b/src/methods/contains.jl
@@ -1,44 +1,52 @@
+# # Containment
+
+# This currently works for point-in-linestring or point-in-polygon.
+
 # More GeometryBasics code
 
-_cross(p1, p2, p3) = (p1[1] - p3[1]) * (p2[2] - p3[2]) - (p2[1] - p3[1]) * (p1[2] - p3[2])
+_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) = contains(GI.trait(pointlist), GI.trait(point), pointlist, point)
 
 # Implementation of a point-in-polygon algorithm
 # from Luxor.jl.  This is the Hormann-Agathos (2001) algorithm.
 
 # For the source, see https://github.com/JuliaGraphics/Luxor.jl/blob/66d60fb51f6b1bb38690fe8dcc6c0084eeb80710/src/polygons.jl#L190-L229.
-function contains(ls::GeometryBasics.LineString{2, T1}, point::Point{2, T2}) where {T1, T2}
-    pointlist = decompose(Point{2, promote_type(T1, T2)}, ls)
+function contains(::GI.LineStringTrait, ::GI.PointTrait, pointlist, point)
+    n = GI.npoint(pointlist)
     c = false
-    @inbounds for counter in eachindex(pointlist)
-        q1 = pointlist[counter]
+    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 == length(pointlist)
-            q2 = pointlist[1]
+        if counter == (n-1)
+            q2 = GI.getpoint(pointlist, 1)
         else
-            q2 = pointlist[counter + 1]
+            q2 = current_point
         end
-        if q1 == point
+        if GI.x(q1) == GI.x(point) && GI.x(q1) == GI.y(point)
             # allowonedge || error("isinside(): VertexException a")
             continue
         end
-        if q2[2] == point[2]
-            if q2[1] == point[1]
+        if GI.y(q2) == GI.y(point)
+            if GI.x(q2) == GI.x(point)
                 # allowonedge || error("isinside(): VertexException b")
                 continue
-            elseif (q1[2] == point[2]) && ((q2[1] > point[1]) == (q1[1] < point[1]))
+            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 (q1[2] < point[2]) != (q2[2] < point[2]) # crossing
-            if q1[1] >= point[1]
-                if q2[1] > point[1]
+        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) == (q2[2] > q1[2]))
+                elseif ((_cross(q1, q2, point) > 0) == (GI.y(q2) > GI.y(q1)))
                     c = !c
                 end
-            elseif q2[1] > point[1]
-                if ((_cross(q1, q2, point) > 0) == (q2[2] > q1[2]))
+            elseif GI.x(q2) > GI.x(point)
+                if ((_cross(q1, q2, point) > 0) == (GI.y(q2) > GI.y(q1)))
                     c = !c
                 end
             end
diff --git a/src/methods/signed_area.jl b/src/methods/signed_area.jl
index 0113e362c..7e202701f 100644
--- a/src/methods/signed_area.jl
+++ b/src/methods/signed_area.jl
@@ -1,4 +1,5 @@
 # # Signed area
+
 export signed_area
 
 # ## What is signed area?
@@ -11,7 +12,10 @@ export signed_area
 
 # To provide an example, consider this rectangle:
 # ```@example rect
-# using GeometryBasics
+# using GeometryOps
+# using GeometryOps.GeometryBasics
+# using Makie
+# 
 # rect = Polygon([Point(0,0), Point(0,1), Point(1,1), Point(1,0), Point(0, 0)])
 # f, a, p = poly(rect; axis = (; aspect = DataAspect()))
 # ```
diff --git a/src/methods/signed_distance.jl b/src/methods/signed_distance.jl
index 52db4e6fd..c0ccb0e20 100644
--- a/src/methods/signed_distance.jl
+++ b/src/methods/signed_distance.jl
@@ -1,5 +1,7 @@
 # # Signed distance
+
 export signed_distance
+
 # TODO: clean this up.  It already supports GeoInterface.
 
 Base.@propagate_inbounds euclid_distance(p1, p2) = sqrt((GeoInterface.x(p2)-GeoInterface.x(p1))^2 + (GeoInterface.y(p2)-GeoInterface.y(p1))^2)
diff --git a/src/primitives.jl b/src/primitives.jl
index 634f509bb..431e7fa76 100644
--- a/src/primitives.jl
+++ b/src/primitives.jl
@@ -1,3 +1,7 @@
+# # Primitive functions
+
+# This file mainly defines the [`apply`](@ref) function.
+
 """
     apply(f, target::Type{<:AbstractTrait}, obj; crs)
 
diff --git a/src/transformations/flip.jl b/src/transformations/flip.jl
index e8b3e6cdd..194119095 100644
--- a/src/transformations/flip.jl
+++ b/src/transformations/flip.jl
@@ -1,3 +1,8 @@
+# # Coordinate flipping
+
+# This is a simple example of how to use the `apply` functionality in a function,
+# by flipping the x and y coordinates of a geometry.
+
 """
     flip(obj)
 
diff --git a/src/transformations/reproject.jl b/src/transformations/reproject.jl
index 5d8e9bb1f..5e73df925 100644
--- a/src/transformations/reproject.jl
+++ b/src/transformations/reproject.jl
@@ -1,4 +1,12 @@
+# # Geometry reprojection
 
+export reproject
+
+# This file is pretty simple - it simply reprojects a geometry pointwise from one CRS
+# to another. It uses the `Proj` package for the transformation, but this could be 
+# moved to an extension if needed.
+
+# This works using the [`apply`](@ref) functionality.
 
 """
     reproject(geometry; source_crs, target_crs, transform, always_xy, time)
diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl
index 0ffbdb36a..d9b4b4b11 100644
--- a/src/transformations/simplify.jl
+++ b/src/transformations/simplify.jl
@@ -1,3 +1,23 @@
+# # Geometry simplification
+
+# This file holds implementations for the Douglas-Peucker and Visvalingam-Whyatt
+# algorithms for simplifying geometries (specifically polygons and lines).
+
+export simplify, VisvalingamWhyatt, DouglasPeucker
+
+
+"""
+    abstract type SimplifyAlg
+
+Abstract type for simplification algorithms.
+
+## API
+
+For now, the algorithm must hold the `number`, `ratio` and `tol` properties.  
+
+Simplification algorithm types can hook into the interface by implementing 
+the `_simplify(trait, alg, geom)` methods for whichever traits are necessary.
+"""
 abstract type SimplifyAlg end
 
 const SIMPLIFY_ALG_KEYWORDS = """
@@ -83,27 +103,27 @@ simplify(data; kw...) = _simplify(DouglasPeucker(; kw...), data)
 simplify(alg::SimplifyAlg, data) = _simplify(alg, data)
 
 function _simplify(alg::SimplifyAlg, data)
-    # Apply simplication to all curves, multipoints, and points,
-    # reconstructing everything else around them.
+    ## Apply simplication to all curves, multipoints, and points,
+    ## reconstructing everything else around them.
     simplifier(geom) = _simplify(trait(geom), alg, geom)
     apply(simplifier, Union{PolygonTrait,AbstractCurveTrait,MultiPoint,PointTrait}, data)
 end
-# For Point and MultiPoint traits we do nothing
+## For Point and MultiPoint traits we do nothing
 _simplify(::PointTrait, alg, geom) = geom
 _simplify(::MultiPointTrait, alg, geom) = geom
 function _simplify(::PolygonTrait, alg, geom)
-    # Force treating children as LinearRing
+    ## Force treating children as LinearRing
     rebuilder(g) = rebuild(g, _simplify(LinearRingTrait(), alg, g))
     lrs = map(rebuilder, GI.getgeom(geom))
     return rebuild(geom, lrs)
 end
-# For curves and rings we simplify
+## For curves and rings we simplify
 _simplify(::AbstractCurveTrait, alg, geom) = rebuild(geom, simplify(alg, tuple_points(geom)))
 function _simplify(::LinearRingTrait, alg, geom)
-    # Make a vector of points 
+    ## Make a vector of points 
     points = tuple_points(geom)
 
-    # Simplify it once
+    ## Simplify it once
     simple = _simplify(alg, points)
 
     return rebuild(geom, simple)
@@ -138,9 +158,9 @@ function _simplify(alg::RadialDistance, points::Vector)
         distances[i] = _squared_dist(point, previous)
         previous = point
     end
-    # Never remove the end points
+    ## Never remove the end points
     distances[begin] = distances[end] = Inf
-    # This avoids taking the square root of each distance above
+    ## This avoids taking the square root of each distance above
     if !isnothing(alg.tol)
         alg = settol(alg, (alg.tol::Float64)^2)
     end
@@ -180,8 +200,8 @@ settol(alg::DouglasPeucker, tol) = DouglasPeucker(alg.number, alg.ratio, tol, al
 
 function _simplify(alg::DouglasPeucker, points::Vector)
     length(points) <= MIN_POINTS && return points
-    # TODO do we need this?
-    # points = alg.prefilter ? simplify(RadialDistance(alg.tol), points) : points
+    ## TODO do we need this?
+    ## points = alg.prefilter ? simplify(RadialDistance(alg.tol), points) : points
 
     distances = _build_tolerances(_squared_segdist, points)
     return _get_points(alg, points, distances)
@@ -239,19 +259,19 @@ function _simplify(alg::VisvalingamWhyatt, points::Vector)
     length(points) <= MIN_POINTS && return points
     areas = _build_tolerances(_triangle_double_area, points)
 
-    # This avoids diving everything by two
+    ## This avoids diving everything by two
     if !isnothing(alg.tol)
         alg = settol(alg, (alg.tol::Float64)*2)
     end
     return _get_points(alg, points, areas)
 end
 
-# calculates the area of a triangle given its vertices
+## calculates the area of a triangle given its vertices
 _triangle_double_area(p1, p2, p3) =
     abs(p1[1] * (p2[2] - p3[2]) + p2[1] * (p3[2] - p1[2]) + p3[1] * (p1[2] - p2[2]))
 
 
-# Shared utils
+# ### Shared utils
 
 function _build_tolerances(f, points)
     nmax = length(points)
@@ -317,7 +337,11 @@ function tuple_points(geom)
 end
 
 function _get_points(alg, points, tolerances)
-    (; tol, number, ratio) = alg
+    ## This assumes that `alg` has the properties
+    ## `tol`, `number`, and `ratio` available...
+    tol = alg.tol
+    number = alg.number
+    ratio = alg.ratio
     bit_indices = if !isnothing(tol) 
         _tol_indices(alg.tol::Float64, points, tolerances)
     elseif !isnothing(number) 
@@ -336,8 +360,8 @@ function _number_indices(n, points, tolerances)
     tol = partialsort(tolerances, length(points) - n + 1)
     bit_indices = _tol_indices(tol, points, tolerances)
     nselected = sum(bit_indices)
-    # If there are multiple values exactly at `tol` we will get 
-    # the wrong output length. So we need to remove some.
+    ## If there are multiple values exactly at `tol` we will get 
+    ## the wrong output length. So we need to remove some.
     while nselected > n
         min_tol = Inf
         min_i = 0
diff --git a/src/utils.jl b/src/utils.jl
index 0fdaa8872..b59940a78 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,3 +1,4 @@
+# # Utility functions
 
 _is3d(geom) = _is3d(GI.trait(geom), geom)
 _is3d(::GI.AbstractGeometryTrait, geom) = GI.is3d(geom)
diff --git a/test/methods/barycentric.jl b/test/methods/barycentric.jl
new file mode 100644
index 000000000..27f378d93
--- /dev/null
+++ b/test/methods/barycentric.jl
@@ -0,0 +1,26 @@
+@testset "Mean value coordinates" begin
+    @testset "Preserving return type" begin
+        @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2{BigFloat}(1,1)) isa Vector{BigFloat}
+        @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2f(1,1)) isa Vector{BigFloat} # keep the most precise type
+        @test barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) isa Vector{Float64}
+        @test barycentric_coordinates(MeanValue(), Point2{Float32}[(0,0), (1,0), (0,1)], Point2{Float32}(1,1)) isa Vector{Float32}
+    end
+    @testset "Triangle coordinates" begin
+        # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (0.5,0.5) are (0.5,0.25,0.25)
+        @test all(barycentric_coordinates(MeanValue(), Point2f[(0,0), (1,0), (0,1)], Point2f(0.5,0.5)) .== (0.5,0.25,0.25))
+        # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (1,1) are (-1,1,1)
+        @test all(barycentric_coordinates(MeanValue(), Point2f[(0,0), (1,0), (0,1)], Point2f(1,1)) .≈ (-1,1,1))
+        # Test that calculations with different number types (in this case Float64) are also accurate
+        @test all(barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) .≈ (-1,1,1))
+    end
+    @testset "Tests for helper methods" begin
+        @testset "`t_value`" begin
+            @test GeometryOps.t_value(Point2f(0,0), Point2f(1,0), 1, 1) == 0
+            @test GeometryOps.t_value(Point2f(0, 1), Point2f(1, 0), 1, 2) == -0.5f0
+        end
+        @testset "`_det`" begin
+            @test GeometryOps._det((1,0), (0,1)) == 1f0
+            @test GeometryOps._det(Point2f(1, 2), Point2f(3, 4)) == -2.0f0
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index bb1917a16..dd6f4d964 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@ using GeometryOps
 using Test
 
 using GeometryOps.GeoInterface
+using GeometryOps.GeometryBasics
 using ArchGDAL
 
 const GI = GeoInterface
@@ -10,7 +11,8 @@ const AG = ArchGDAL
 @testset "GeometryOps.jl" begin
     @testset "Primitives" begin include("primitives.jl") end
     @testset "Signed Area" begin include("methods/signed_area.jl") end
-    @testset "reproject" begin include("transformations/reproject.jl") end
-    @testset "flip" begin include("transformations/flip.jl") end
-    @testset "simplify" begin include("transformations/simplify.jl") end
+    @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
+    @testset "Reproject" begin include("transformations/reproject.jl") end
+    @testset "Flip" begin include("transformations/flip.jl") end
+    @testset "Simplify" begin include("transformations/simplify.jl") end
 end