[core] Introduce distance expression

Add distance unit choices

Fix cmake and add license

Add support for LineString Features

Add template to geometry helper function

Only support line and point

Rename geometry_within.cpp hpp file to geometry_util.cpp .hpp

Remove incorrect indexFilter, fix pointSetsDistance

Fix distance expression

8 files changed, 789 insertions, 175 deletions

diff --git a/src/mbgl/style/expression/distance.cpp b/src/mbgl/style/expression/distance.cpp
new file mode 100644
index 0000000000..67a460b7d7
--- /dev/null
+++ b/src/mbgl/style/expression/distance.cpp
@@ -0,0 +1,601 @@
+#include <mbgl/style/expression/distance.hpp>
+
+#include <mapbox/cheap_ruler.hpp>
+#include <mapbox/geojson.hpp>
+#include <mapbox/geometry.hpp>
+
+#include <mbgl/style/conversion/json.hpp>
+#include <mbgl/tile/geometry_tile_data.hpp>
+#include <mbgl/util/geometry_util.hpp>
+#include <mbgl/util/logging.hpp>
+#include <mbgl/util/string.hpp>
+
+#include <rapidjson/document.h>
+
+#include <algorithm>
+#include <deque>
+#include <limits>
+#include <queue>
+#include <tuple>
+
+namespace mbgl {
+namespace {
+
+const std::size_t MinPointsSize = 500;
+const std::size_t MinLinePointsSize = 200;
+
+using BBox = std::array<double, 4>;
+const BBox DefaultBBox = BBox{std::numeric_limits<double>::infinity(),
+                              std::numeric_limits<double>::infinity(),
+                              -std::numeric_limits<double>::infinity(),
+                              -std::numeric_limits<double>::infinity()};
+
+// bbox[minX, minY, maxX, maxY]
+void updateBBox(BBox& bbox, const mapbox::geometry::point<double>& p) {
+    bbox[0] = std::min(p.x, bbox[0]);
+    bbox[1] = std::min(p.y, bbox[1]);
+    bbox[2] = std::max(p.x, bbox[2]);
+    bbox[3] = std::max(p.y, bbox[3]);
+}
+
+// Inclusive index range for multipoint or linestring container
+using IndexRange = std::pair<std::size_t, std::size_t>;
+
+std::size_t getRangeSize(const IndexRange& range) {
+    return range.second - range.first + 1;
+}
+
+BBox getBBox(const mapbox::geometry::multi_point<double>& points, const IndexRange& range) {
+    BBox bbox = DefaultBBox;
+    for (std::size_t i = range.first; i <= range.second; ++i) {
+        updateBBox(bbox, points[i]);
+    }
+    return bbox;
+}
+
+BBox getBBox(const mapbox::geometry::line_string<double>& line, const IndexRange& range) {
+    BBox bbox = DefaultBBox;
+    for (std::size_t i = range.first; i <= range.second; ++i) {
+        updateBBox(bbox, line[i]);
+    }
+    return bbox;
+}
+
+// Calculate the distance between two bounding boxes.
+// Calculate the delta in x and y direction, and use two fake points {0, 0} and {dx, dy} to calculate the distance.
+// Distance will be 0 if bounding box are overlapping.
+double bboxToBBoxDistance(const BBox& bbox1, const BBox& bbox2, mapbox::cheap_ruler::CheapRuler& ruler) {
+    double dx = 0.;
+    double dy = 0.;
+    // bbox1 in left side
+    if (bbox1[2] < bbox2[0]) {
+        dx = bbox2[0] - bbox1[2];
+    }
+    // bbox1 in right side
+    if (bbox1[0] > bbox2[2]) {
+        dx = bbox1[0] - bbox2[2];
+    }
+    // bbox1 in above side
+    if (bbox1[1] > bbox2[3]) {
+        dy = bbox1[1] - bbox2[3];
+    }
+    // bbox1 in down side
+    if (bbox1[3] < bbox2[1]) {
+        dy = bbox2[1] - bbox1[3];
+    }
+    return ruler.distance(mapbox::geometry::point<double>{0., 0.}, mapbox::geometry::point<double>{dx, dy});
+}
+
+double pointToLineDistance(const mapbox::geometry::point<double>& point,
+                           const mapbox::geometry::line_string<double>& line,
+                           mapbox::cheap_ruler::CheapRuler& ruler) {
+    const auto nearestPoint = std::get<0>(ruler.pointOnLine(line, point));
+    return ruler.distance(point, nearestPoint);
+}
+
+double lineToLineDistance(const mapbox::geometry::line_string<double>& line1,
+                          IndexRange& range1,
+                          const mapbox::geometry::line_string<double>& line2,
+                          IndexRange& range2,
+                          mapbox::cheap_ruler::CheapRuler& ruler) {
+    double dist = std::numeric_limits<double>::infinity();
+    for (std::size_t i = range1.first; i < range1.second; ++i) {
+        const auto& p1 = line1[i];
+        const auto& p2 = line1[i + 1];
+        for (std::size_t j = range2.first; j < range2.second; ++j) {
+            const auto& q1 = line2[j];
+            const auto& q2 = line2[j + 1];
+            if (GeometryUtil<double>::segmentIntersectSegment(p1, p2, q1, q2)) return 0.;
+            auto dist1 = std::min(pointToLineDistance(p1, mapbox::geometry::line_string<double>{q1, q2}, ruler),
+                                  pointToLineDistance(p2, mapbox::geometry::line_string<double>{q1, q2}, ruler));
+            auto dist2 = std::min(pointToLineDistance(q1, mapbox::geometry::line_string<double>{p1, p2}, ruler),
+                                  pointToLineDistance(q2, mapbox::geometry::line_string<double>{p1, p2}, ruler));
+            dist = std::min(dist, std::min(dist1, dist2));
+        }
+    }
+    return dist;
+}
+
+double pointsToPointsDistance(const mapbox::geometry::multi_point<double>& points1,
+                              IndexRange& range1,
+                              const mapbox::geometry::multi_point<double>& points2,
+                              IndexRange& range2,
+                              mapbox::cheap_ruler::CheapRuler& ruler) {
+    double dist = std::numeric_limits<double>::infinity();
+    for (std::size_t i = range1.first; i <= range1.second; ++i) {
+        for (std::size_t j = range2.first; j <= range2.second; ++j) {
+            dist = std::min(dist, ruler.distance(points1[i], points2[j]));
+            if (dist == 0.) return dist;
+        }
+    }
+    return dist;
+}
+
+std::pair<mbgl::optional<IndexRange>, mbgl::optional<IndexRange>> splitRange(const IndexRange& range, bool isLine) {
+    auto size = getRangeSize(range);
+    if (isLine) {
+        if (size == 2) {
+            return std::make_pair(range, nullopt);
+        }
+        auto size1 = size / 2;
+        IndexRange range1(range.first, range.first + size1);
+        IndexRange range2(range.first + size1, range.second);
+        return std::make_pair(std::move(range1), std::move(range2));
+    } else {
+        if (size == 1) {
+            return std::make_pair(range, nullopt);
+        }
+        auto size1 = size / 2 - 1;
+        IndexRange range1(range.first, range.first + size1);
+        IndexRange range2(range.first + size1 + 1, range.second);
+        return std::make_pair(std::move(range1), std::move(range2));
+    }
+}
+
+// <distance, range1, range2>
+using DistPair = std::tuple<double, IndexRange, IndexRange>;
+class Comparator {
+public:
+    bool operator()(DistPair& left, DistPair& right) { return std::get<0>(left) < std::get<0>(right); }
+};
+// The priority queue will ensure the top element would always be the pair that has the biggest distance
+using DistQueue = std::priority_queue<DistPair, std::deque<DistPair>, Comparator>;
+
+// Divide and conqure, the time complexity is O(n*lgn), faster than Brute force O(n*n)
+// Use index for in-place processing.
+double lineToLineDistance(const mapbox::geometry::line_string<double>& line1,
+                          const mapbox::geometry::multi_point<double>& line2,
+                          mapbox::cheap_ruler::CheapRuler& ruler) {
+    auto miniDist = ruler.distance(line1[0], line2[0]);
+    DistQueue distQueue;
+    distQueue.push(std::forward_as_tuple(0, IndexRange(0, line1.size() - 1), IndexRange(0, line2.size() - 1)));
+
+    while (!distQueue.empty()) {
+        auto distPair = distQueue.top();
+        distQueue.pop();
+        if (std::get<0>(distPair) > miniDist) continue;
+        auto& rangeA = std::get<1>(distPair);
+        auto& rangeB = std::get<2>(distPair);
+
+        // In case the set size are relatively small, we could use brute-force directly
+        if (getRangeSize(rangeA) <= MinLinePointsSize && getRangeSize(rangeB) <= MinLinePointsSize) {
+            miniDist = std::min(miniDist, lineToLineDistance(line1, rangeA, line2, rangeB, ruler));
+            if (miniDist == 0.) return 0.;
+        } else {
+            auto newRangesA = splitRange(rangeA, true /*isLine*/);
+            auto newRangesB = splitRange(rangeB, true /*isLine*/);
+            const auto updateQueue = [&distQueue, &miniDist, &ruler, &line1, &line2](
+                                         mbgl::optional<IndexRange>& range1, mbgl::optional<IndexRange>& range2) {
+                if (!range1 || !range2) return;
+                auto tempDist = bboxToBBoxDistance(getBBox(line1, *range1), getBBox(line2, *range2), ruler);
+                // Insert new pair to the queue if the bbox distance is less or equal to miniDist,
+                // The pair with biggest distance will be at the top
+                if (tempDist <= miniDist)
+                    distQueue.push(std::make_tuple(tempDist, std::move(*range1), std::move(*range2)));
+            };
+            updateQueue(newRangesA.first, newRangesB.first);
+            updateQueue(newRangesA.first, newRangesB.second);
+            updateQueue(newRangesA.second, newRangesB.first);
+            updateQueue(newRangesA.second, newRangesB.second);
+        }
+    }
+    return miniDist;
+}
+
+// Divide and conqure, the time complexity is O(n*lgn), faster than Brute force O(n*n)
+// Use index for in-place processing.
+double pointsToPointsDistance(const mapbox::geometry::multi_point<double>& pointSet1,
+                              const mapbox::geometry::multi_point<double>& pointSet2,
+                              mapbox::cheap_ruler::CheapRuler& ruler) {
+    auto miniDist = ruler.distance(pointSet1[0], pointSet2[0]);
+    DistQueue distQueue;
+    distQueue.push(std::forward_as_tuple(0, IndexRange(0, pointSet1.size() - 1), IndexRange(0, pointSet2.size() - 1)));
+
+    while (!distQueue.empty()) {
+        auto distPair = distQueue.top();
+        distQueue.pop();
+        if (std::get<0>(distPair) > miniDist) {
+            continue;
+        }
+        auto& rangeA = std::get<1>(distPair);
+        auto& rangeB = std::get<2>(distPair);
+
+        // In case the set size are relatively small, we could use brute-force directly
+        if (getRangeSize(rangeA) <= MinPointsSize && getRangeSize(rangeB) <= MinPointsSize) {
+            miniDist = std::min(miniDist, pointsToPointsDistance(pointSet1, rangeA, pointSet2, rangeB, ruler));
+            if (miniDist == 0.) return 0.;
+        } else {
+            auto newRangesA = splitRange(rangeA, false /*isLine*/);
+            auto newRangesB = splitRange(rangeB, false /*isLine*/);
+            const auto updateQueue = [&distQueue, &miniDist, &ruler, &pointSet1, &pointSet2](
+                                         mbgl::optional<IndexRange>& range1, mbgl::optional<IndexRange>& range2) {
+                if (!range1 || !range2) return;
+                auto tempDist = bboxToBBoxDistance(getBBox(pointSet1, *range1), getBBox(pointSet2, *range2), ruler);
+                // Insert new pair to the queue if the bbox distance is less or equal to miniDist,
+                // The pair with biggest distance will be at the top
+                if (tempDist <= miniDist)
+                    distQueue.push(std::make_tuple(tempDist, std::move(*range1), std::move(*range2)));
+            };
+            updateQueue(newRangesA.first, newRangesB.first);
+            updateQueue(newRangesA.first, newRangesB.second);
+            updateQueue(newRangesA.second, newRangesB.first);
+            updateQueue(newRangesA.second, newRangesB.second);
+        }
+    }
+    return miniDist;
+}
+
+// Divide and conqure, the time complexity is O(n*lgn), faster than Brute force O(n*n)
+// Most of the time, use index for in-place processing.
+double pointsToLineDistance(const mapbox::geometry::multi_point<double>& points,
+                            const mapbox::geometry::line_string<double>& line,
+                            mapbox::cheap_ruler::CheapRuler& ruler) {
+    auto miniDist = ruler.distance(points[0], line[0]);
+    DistQueue distQueue;
+    distQueue.push(std::forward_as_tuple(0, IndexRange(0, points.size() - 1), IndexRange(0, line.size() - 1)));
+
+    while (!distQueue.empty()) {
+        auto distPair = distQueue.top();
+        distQueue.pop();
+        if (std::get<0>(distPair) > miniDist) continue;
+        auto& rangeA = std::get<1>(distPair);
+        auto& rangeB = std::get<2>(distPair);
+
+        // In case the set size are relatively small, we could use brute-force directly
+        if (getRangeSize(rangeA) <= MinPointsSize && getRangeSize(rangeB) <= MinLinePointsSize) {
+            auto subLine =
+                mapbox::geometry::multi_point<double>(line.begin() + rangeB.first, line.begin() + rangeB.second + 1);
+            for (std::size_t i = rangeA.first; i <= rangeA.second; ++i) {
+                miniDist = std::min(miniDist, pointToLineDistance(points[i], subLine, ruler));
+                if (miniDist == 0.) return 0.;
+            }
+        } else {
+            auto newRangesA = splitRange(rangeA, false /*isLine*/);
+            auto newRangesB = splitRange(rangeB, true /*isLine*/);
+            const auto updateQueue = [&distQueue, &miniDist, &ruler, &points, &line](
+                                         mbgl::optional<IndexRange>& range1, mbgl::optional<IndexRange>& range2) {
+                if (!range1 || !range2) return;
+                auto tempDist = bboxToBBoxDistance(getBBox(points, *range1), getBBox(line, *range2), ruler);
+                // Insert new pair to the queue if the bbox distance is less or equal to miniDist,
+                // The pair with biggest distance will be at the top
+                if (tempDist <= miniDist)
+                    distQueue.push(std::make_tuple(tempDist, std::move(*range1), std::move(*range2)));
+            };
+            updateQueue(newRangesA.first, newRangesB.first);
+            updateQueue(newRangesA.first, newRangesB.second);
+            updateQueue(newRangesA.second, newRangesB.first);
+            updateQueue(newRangesA.second, newRangesB.second);
+        }
+    }
+    return miniDist;
+}
+
+double pointsToLinesDistance(const mapbox::geometry::multi_point<double>& points,
+                             const mapbox::geometry::multi_line_string<double>& lines,
+                             mapbox::cheap_ruler::CheapRuler& ruler) {
+    double dist = std::numeric_limits<double>::infinity();
+    for (const auto& line : lines) {
+        dist = std::min(dist, pointsToLineDistance(points, line, ruler));
+        if (dist == 0.) return dist;
+    }
+    return dist;
+}
+
+double lineToLinesDistance(const mapbox::geometry::line_string<double>& line,
+                           const mapbox::geometry::multi_line_string<double>& lines,
+                           mapbox::cheap_ruler::CheapRuler& ruler) {
+    double dist = std::numeric_limits<double>::infinity();
+    for (const auto& l : lines) {
+        dist = std::min(dist, lineToLineDistance(line, l, ruler));
+        if (dist == 0.) return dist;
+    }
+    return dist;
+}
+
+double pointsToGeometryDistance(const mapbox::geometry::multi_point<double>& points,
+                                const Feature::geometry_type& geoSet,
+                                mapbox::cheap_ruler::CheapRuler::Unit unit) {
+    mapbox::cheap_ruler::CheapRuler ruler(points.front().y, unit);
+    return geoSet.match(
+        [&points, &ruler](const mapbox::geometry::point<double>& p) {
+            return pointsToPointsDistance(mapbox::geometry::multi_point<double>{p}, points, ruler);
+        },
+        [&points, &ruler](const mapbox::geometry::multi_point<double>& points1) {
+            return pointsToPointsDistance(points, points1, ruler);
+        },
+        [&points, &ruler](const mapbox::geometry::line_string<double>& line) {
+            return pointsToLineDistance(points, line, ruler);
+        },
+        [&points, &ruler](const mapbox::geometry::multi_line_string<double>& lines) {
+            return pointsToLinesDistance(points, lines, ruler);
+        },
+        [](const auto&) { return std::numeric_limits<double>::quiet_NaN(); });
+}
+
+double lineToGeometryDistance(const mapbox::geometry::line_string<double>& line,
+                              const Feature::geometry_type& geoSet,
+                              mapbox::cheap_ruler::CheapRuler::Unit unit) {
+    assert(!line.empty());
+    mapbox::cheap_ruler::CheapRuler ruler(line.front().y, unit);
+    return geoSet.match(
+        [&line, &ruler](const mapbox::geometry::point<double>& p) {
+            return pointsToLineDistance(mapbox::geometry::multi_point<double>{p}, line, ruler);
+        },
+        [&line, &ruler](const mapbox::geometry::multi_point<double>& points) {
+            return pointsToLineDistance(points, line, ruler);
+        },
+        [&line, &ruler](const mapbox::geometry::line_string<double>& line1) {
+            return lineToLineDistance(line, line1, ruler);
+        },
+        [&line, &ruler](const mapbox::geometry::multi_line_string<double>& lines) {
+            return lineToLinesDistance(line, lines, ruler);
+        },
+        [](const auto&) { return std::numeric_limits<double>::quiet_NaN(); });
+}
+
+double calculateDistance(const GeometryTileFeature& feature,
+                         const CanonicalTileID& canonical,
+                         const Feature::geometry_type& geoSet,
+                         mapbox::cheap_ruler::CheapRuler::Unit unit) {
+    return convertGeometry(feature, canonical)
+        .match(
+            [&geoSet, &unit](const mapbox::geometry::point<double>& point) -> double {
+                return pointsToGeometryDistance(mapbox::geometry::multi_point<double>{point}, geoSet, unit);
+            },
+            [&geoSet, &unit](const mapbox::geometry::multi_point<double>& points) -> double {
+                return pointsToGeometryDistance(points, geoSet, unit);
+            },
+            [&geoSet, &unit](const mapbox::geometry::line_string<double>& line) -> double {
+                return lineToGeometryDistance(line, geoSet, unit);
+            },
+            [&geoSet, &unit](const mapbox::geometry::multi_line_string<double>& lines) -> double {
+                double dist = std::numeric_limits<double>::infinity();
+                for (const auto& line : lines) {
+                    dist = std::min(dist, lineToGeometryDistance(line, geoSet, unit));
+                    if (dist == 0.) return dist;
+                }
+                return dist;
+            },
+            [](const auto&) -> double { return std::numeric_limits<double>::quiet_NaN(); });
+}
+
+struct Arguments {
+    Arguments(GeoJSON& geojson_, mapbox::cheap_ruler::CheapRuler::Unit unit_)
+        : geojson(std::move(geojson_)), unit(unit_) {}
+
+    GeoJSON geojson;
+    mapbox::cheap_ruler::CheapRuler::Unit unit;
+};
+
+optional<Arguments> parseValue(const style::conversion::Convertible& value, style::expression::ParsingContext& ctx) {
+    if (isArray(value)) {
+        // object value, quoted with ["Distance", GeoJSONObj, "unit(optional)"]
+        auto length = arrayLength(value);
+        if (length != 2 && length != 3) {
+            ctx.error("'distance' expression requires either one argument or two arguments, but found " +
+                      util::toString(arrayLength(value) - 1) + " instead.");
+            return nullopt;
+        }
+
+        // Parse Unit info for distance calculation, "Meters" by default
+        mapbox::cheap_ruler::CheapRuler::Unit unit = mapbox::cheap_ruler::CheapRuler::Unit::Meters;
+        if (length == 3) {
+            auto input = toString(arrayMember(value, 2));
+            if (input == nullopt) {
+                ctx.error("Failed to parse unit argument from 'distance' expression");
+                return nullopt;
+            }
+            if (*input == "Meters" || *input == "Metres") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::Meters;
+            } else if (*input == "Kilometers") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::Kilometers;
+            } else if (*input == "Miles") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::Miles;
+            } else if (*input == "NauticalMiles") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::NauticalMiles;
+            } else if (*input == "Yards") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::Kilometers;
+            } else if (*input == "Feet") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::Miles;
+            } else if (*input == "Inches") {
+                unit = mapbox::cheap_ruler::CheapRuler::Unit::Inches;
+            } else {
+                ctx.error(
+                    "'distance' expression only accepts following Units:  'Kilometers', 'Miles', 'NauticalMiles', "
+                    "'Meters', 'Metres', 'Yards', 'Feet', 'Inches'.");
+                return nullopt;
+            }
+        }
+        // Parse geometry info
+        const auto& argument1 = arrayMember(value, 1);
+        if (isObject(argument1)) {
+            style::conversion::Error error;
+            auto geojson = toGeoJSON(argument1, error);
+            if (geojson && error.message.empty()) {
+                return Arguments(*geojson, unit);
+            }
+            ctx.error(error.message);
+        }
+    }
+    ctx.error(
+        "'distance' expression needs to be an array with format [\"Distance\", GeoJSONObj, \"unit(optional, Meters by "
+        "default)\"].");
+    return nullopt;
+}
+
+optional<Feature::geometry_type> getGeometry(const Feature& feature, mbgl::style::expression::ParsingContext& ctx) {
+    const auto type = apply_visitor(ToFeatureType(), feature.geometry);
+    if (type == FeatureType::Point || type == FeatureType::LineString) {
+        return feature.geometry;
+    }
+    ctx.error("'distance' expression requires valid geojson object with valid geometry type: Point/LineString.");
+    return nullopt;
+}
+} // namespace
+
+namespace style {
+namespace expression {
+
+Distance::Distance(GeoJSON geojson, Feature::geometry_type geometries_, mapbox::cheap_ruler::CheapRuler::Unit unit_)
+    : Expression(Kind::Distance, type::Number),
+      geoJSONSource(std::move(geojson)),
+      geometries(std::move(geometries_)),
+      unit(unit_) {}
+
+Distance::~Distance() = default;
+
+using namespace mbgl::style::conversion;
+
+EvaluationResult Distance::evaluate(const EvaluationContext& params) const {
+    if (!params.feature || !params.canonical) {
+        return EvaluationError{"distance expression requirs valid feature and canonical information."};
+    }
+    auto geometryType = params.feature->getType();
+    if (geometryType == FeatureType::Point || geometryType == FeatureType::LineString) {
+        auto distance = calculateDistance(*params.feature, *params.canonical, geometries, unit);
+        if (!std::isnan(distance)) {
+            return distance;
+        }
+    }
+    return EvaluationError{"distance expression currently only evaluates Point/LineString geometries."};
+}
+
+ParseResult Distance::parse(const Convertible& value, ParsingContext& ctx) {
+    auto parsedValue = parseValue(value, ctx);
+    if (!parsedValue) {
+        return ParseResult();
+    }
+
+    return parsedValue->geojson.match(
+        [&parsedValue, &ctx](const mapbox::geometry::geometry<double>& geometrySet) {
+            if (auto ret = getGeometry(mbgl::Feature(geometrySet), ctx)) {
+                return ParseResult(
+                    std::make_unique<Distance>(parsedValue->geojson, std::move(*ret), parsedValue->unit));
+            }
+            return ParseResult();
+        },
+        [&parsedValue, &ctx](const mapbox::feature::feature<double>& feature) {
+            if (auto ret = getGeometry(mbgl::Feature(feature), ctx)) {
+                return ParseResult(
+                    std::make_unique<Distance>(parsedValue->geojson, std::move(*ret), parsedValue->unit));
+            }
+            return ParseResult();
+        },
+        [&parsedValue, &ctx](const mapbox::feature::feature_collection<double>& features) {
+            for (const auto& feature : features) {
+                if (auto ret = getGeometry(mbgl::Feature(feature), ctx)) {
+                    return ParseResult(
+                        std::make_unique<Distance>(parsedValue->geojson, std::move(*ret), parsedValue->unit));
+                }
+            }
+            return ParseResult();
+        },
+        [&ctx](const auto&) {
+            ctx.error("'distance' expression requires valid geojson that contains LineString/Point geometries.");
+            return ParseResult();
+        });
+
+    return ParseResult();
+}
+
+Value convertValue(const mapbox::geojson::rapidjson_value& v) {
+    if (v.IsDouble()) {
+        return v.GetDouble();
+    }
+    if (v.IsString()) {
+        return std::string(v.GetString());
+    }
+    if (v.IsArray()) {
+        std::vector<Value> result;
+        result.reserve(v.Size());
+        for (const auto& m : v.GetArray()) {
+            result.push_back(convertValue(m));
+        }
+        return result;
+    }
+    if (v.IsObject()) {
+        std::unordered_map<std::string, Value> result;
+        for (const auto& m : v.GetObject()) {
+            result.emplace(m.name.GetString(), convertValue(m.value));
+        }
+        return result;
+    }
+    // Ignore other types as valid geojson only contains above types.
+    return Null;
+}
+
+std::string getUnits(const mapbox::cheap_ruler::CheapRuler::Unit& unit) {
+    switch (unit) {
+        case mapbox::cheap_ruler::CheapRuler::Kilometers:
+            return "Kilometers";
+        case mapbox::cheap_ruler::CheapRuler::Miles:
+            return "Miles";
+        case mapbox::cheap_ruler::CheapRuler::NauticalMiles:
+            return "NauticalMiles";
+        case mapbox::cheap_ruler::CheapRuler::Meters:
+            return "Meters";
+        case mapbox::cheap_ruler::CheapRuler::Yards:
+            return "Yards";
+        case mapbox::cheap_ruler::CheapRuler::Feet:
+            return "Feet";
+        case mapbox::cheap_ruler::CheapRuler::Inches:
+            return "Inches";
+        default:
+            return "Error";
+    }
+}
+
+mbgl::Value Distance::serialize() const {
+    std::unordered_map<std::string, Value> serialized;
+    rapidjson::CrtAllocator allocator;
+    const mapbox::geojson::rapidjson_value value = mapbox::geojson::convert(geoJSONSource, allocator);
+    if (value.IsObject()) {
+        for (const auto& m : value.GetObject()) {
+            serialized.emplace(m.name.GetString(), convertValue(m.value));
+        }
+    } else {
+        mbgl::Log::Error(mbgl::Event::General,
+                         "Failed to serialize 'distance' expression, converted rapidJSON is not an object");
+    }
+    return std::vector<mbgl::Value>{{getOperator(), *fromExpressionValue<mbgl::Value>(serialized), getUnits(unit)}};
+}
+
+bool Distance::operator==(const Expression& e) const {
+    if (e.getKind() == Kind::Distance) {
+        auto rhs = static_cast<const Distance*>(&e);
+        return geoJSONSource == rhs->geoJSONSource && geometries == rhs->geometries && unit == rhs->unit;
+    }
+    return false;
+}
+
+std::vector<optional<Value>> Distance::possibleOutputs() const {
+    return {nullopt};
+}
+
+std::string Distance::getOperator() const {
+    return "distance";
+}
+
+} // namespace expression
+} // namespace style
+} // namespace mbgl
diff --git a/src/mbgl/style/expression/is_constant.cpp b/src/mbgl/style/expression/is_constant.cpp
index e6bd05418f..7b14e79eda 100644
--- a/src/mbgl/style/expression/is_constant.cpp
+++ b/src/mbgl/style/expression/is_constant.cpp
@@ -33,6 +33,10 @@ bool isFeatureConstant(const Expression& expression) {
         return false;
     }
 
+    if (expression.getKind() == Kind::Distance) {
+        return false;
+    }
+
     if (expression.getKind() == Kind::CollatorExpression) {
         // Although the results of a Collator expression with fixed arguments
         // generally shouldn't change between executions, we can't serialize them
diff --git a/src/mbgl/style/expression/parsing_context.cpp b/src/mbgl/style/expression/parsing_context.cpp
index 56556d770c..f3842613d1 100644
--- a/src/mbgl/style/expression/parsing_context.cpp
+++ b/src/mbgl/style/expression/parsing_context.cpp
@@ -12,6 +12,7 @@
 #include <mbgl/style/expression/coercion.hpp>
 #include <mbgl/style/expression/comparison.hpp>
 #include <mbgl/style/expression/compound_expression.hpp>
+#include <mbgl/style/expression/distance.hpp>
 #include <mbgl/style/expression/expression.hpp>
 #include <mbgl/style/expression/format_expression.hpp>
 #include <mbgl/style/expression/image_expression.hpp>
@@ -121,6 +122,7 @@ MAPBOX_ETERNAL_CONSTEXPR const auto expressionRegistry =
         {"case", Case::parse},
         {"coalesce", Coalesce::parse},
         {"collator", CollatorExpression::parse},
+        {"distance", Distance::parse},
         {"format", FormatExpression::parse},
         {"image", ImageExpression::parse},
         {"interpolate", parseInterpolate},
diff --git a/src/mbgl/style/expression/within.cpp b/src/mbgl/style/expression/within.cpp
index 99449a50bd..a3e2500319 100644
--- a/src/mbgl/style/expression/within.cpp
+++ b/src/mbgl/style/expression/within.cpp
@@ -5,6 +5,7 @@
 #include <mbgl/style/conversion/json.hpp>
 #include <mbgl/tile/geometry_tile_data.hpp>
 
+#include <mbgl/util/geometry_util.hpp>
 #include <mbgl/util/logging.hpp>
 #include <mbgl/util/string.hpp>
 
@@ -14,6 +15,29 @@
 namespace mbgl {
 namespace {
 
+// contains minX, minY, maxX, maxY
+using WithinBBox = std::array<int64_t, 4>;
+const WithinBBox DefaultBBox = WithinBBox{std::numeric_limits<int64_t>::max(),
+                                          std::numeric_limits<int64_t>::max(),
+                                          std::numeric_limits<int64_t>::min(),
+                                          std::numeric_limits<int64_t>::min()};
+
+// check if bbox1 is within bbox2
+void updateBBox(WithinBBox& bbox, const Point<int64_t>& p) {
+    bbox[0] = std::min(p.x, bbox[0]);
+    bbox[1] = std::min(p.y, bbox[1]);
+    bbox[2] = std::max(p.x, bbox[2]);
+    bbox[3] = std::max(p.y, bbox[3]);
+}
+
+bool boxWithinBox(const WithinBBox& bbox1, const WithinBBox& bbox2) {
+    if (bbox1[0] <= bbox2[0]) return false;
+    if (bbox1[2] >= bbox2[2]) return false;
+    if (bbox1[1] <= bbox2[1]) return false;
+    if (bbox1[3] >= bbox2[3]) return false;
+    return true;
+}
+
 Point<int64_t> latLonToTileCoodinates(const Point<double>& point, const mbgl::CanonicalTileID& canonical) {
     const double size = util::EXTENT * std::pow(2, canonical.z);
 
@@ -149,8 +173,9 @@ bool featureWithinPolygons(const GeometryTileFeature& feature,
             MultiPoint<int64_t> points = getTilePoints(geometries.at(0), canonical, pointBBox, polyBBox);
             if (!boxWithinBox(pointBBox, polyBBox)) return false;
 
-            return std::all_of(
-                points.begin(), points.end(), [&polygons](const auto& p) { return pointWithinPolygons(p, polygons); });
+            return std::all_of(points.begin(), points.end(), [&polygons](const auto& p) {
+                return GeometryUtil<int64_t>::pointWithinPolygons(p, polygons);
+            });
         }
         case FeatureType::LineString: {
             WithinBBox lineBBox = DefaultBBox;
@@ -158,7 +183,7 @@ bool featureWithinPolygons(const GeometryTileFeature& feature,
             if (!boxWithinBox(lineBBox, polyBBox)) return false;
 
             return std::all_of(multiLineString.begin(), multiLineString.end(), [&polygons](const auto& line) {
-                return lineStringWithinPolygons(line, polygons);
+                return GeometryUtil<int64_t>::lineStringWithinPolygons(line, polygons);
             });
         }
         default:
diff --git a/src/mbgl/util/geometry_util.cpp b/src/mbgl/util/geometry_util.cpp
new file mode 100644
index 0000000000..df4119ce6a
--- /dev/null
+++ b/src/mbgl/util/geometry_util.cpp
@@ -0,0 +1,129 @@
+#include <mbgl/util/geometry_util.hpp>
+
+#include <algorithm>
+
+namespace mbgl {
+
+template <typename T>
+bool GeometryUtil<T>::rayIntersect(const Point<T>& p, const Point<T>& p1, const Point<T>& p2) {
+    return ((p1.y > p.y) != (p2.y > p.y)) && (p.x < (p2.x - p1.x) * (p.y - p1.y) / (p2.y - p1.y) + p1.x);
+}
+
+// check if point p is on line segment with end points p1 and p2
+template <typename T>
+bool GeometryUtil<T>::pointOnBoundary(const Point<T>& p, const Point<T>& p1, const Point<T>& p2) {
+    // requirements of point p on line segment:
+    // 1. colinear: cross product of vector p->p1(x1, y1) and vector p->p2(x2, y2) equals to 0
+    // 2. p is between p1 and p2
+    const auto x1 = p.x - p1.x;
+    const auto y1 = p.y - p1.y;
+    const auto x2 = p.x - p2.x;
+    const auto y2 = p.y - p2.y;
+    return (x1 * y2 - x2 * y1 == 0) && (x1 * x2 <= 0) && (y1 * y2 <= 0);
+}
+
+template <typename T>
+bool GeometryUtil<T>::segmentIntersectSegment(const Point<T>& a,
+                                              const Point<T>& b,
+                                              const Point<T>& c,
+                                              const Point<T>& d) {
+    // a, b are end points for line segment1, c and d are end points for line segment2
+    const auto perp = [](const Point<T>& v1, const Point<T>& v2) { return (v1.x * v2.y - v1.y * v2.x); };
+
+    // check if two segments are parallel or not
+    // precondition is end point a, b is inside polygon, if line a->b is
+    // parallel to polygon edge c->d, then a->b won't intersect with c->d
+    auto vectorP = Point<T>(b.x - a.x, b.y - a.y);
+    auto vectorQ = Point<T>(d.x - c.x, d.y - c.y);
+    if (perp(vectorQ, vectorP) == 0) return false;
+
+    // check if p1 and p2 are in different sides of line segment q1->q2
+    const auto twoSided = [](const Point<T>& p1, const Point<T>& p2, const Point<T>& q1, const Point<T>& q2) {
+        // q1->p1 (x1, y1), q1->p2 (x2, y2), q1->q2 (x3, y3)
+        T x1 = p1.x - q1.x;
+        T y1 = p1.y - q1.y;
+        T x2 = p2.x - q1.x;
+        T y2 = p2.y - q1.y;
+        T x3 = q2.x - q1.x;
+        T y3 = q2.y - q1.y;
+        auto ret1 = (x1 * y3 - x3 * y1);
+        auto ret2 = (x2 * y3 - x3 * y2);
+        return (ret1 > 0 && ret2 < 0) || (ret1 < 0 && ret2 > 0);
+    };
+
+    // If lines are intersecting with each other, the relative location should be:
+    // a and b lie in different sides of segment c->d
+    // c and d lie in different sides of segment a->b
+    return twoSided(a, b, c, d) && twoSided(c, d, a, b);
+}
+
+template <typename T>
+bool GeometryUtil<T>::lineIntersectPolygon(const Point<T>& p1, const Point<T>& p2, const Polygon<T>& polygon) {
+    for (auto ring : polygon) {
+        auto length = ring.size();
+        // loop through every edge of the ring
+        for (std::size_t i = 0; i < length - 1; ++i) {
+            if (segmentIntersectSegment(p1, p2, ring[i], ring[i + 1])) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+// ray casting algorithm for detecting if point is in polygon
+template <typename T>
+bool GeometryUtil<T>::pointWithinPolygon(const Point<T>& point, const Polygon<T>& polygon, bool trueOnBoundary) {
+    bool within = false;
+    for (const auto& ring : polygon) {
+        const auto length = ring.size();
+        // loop through every edge of the ring
+        for (std::size_t i = 0; i < length - 1; ++i) {
+            if (pointOnBoundary(point, ring[i], ring[i + 1])) return trueOnBoundary;
+            if (rayIntersect(point, ring[i], ring[i + 1])) {
+                within = !within;
+            }
+        }
+    }
+    return within;
+}
+
+template <typename T>
+bool GeometryUtil<T>::pointWithinPolygons(const Point<T>& point, const MultiPolygon<T>& polygons, bool trueOnBoundary) {
+    for (const auto& polygon : polygons) {
+        if (pointWithinPolygon(point, polygon, trueOnBoundary)) return true;
+    }
+    return false;
+}
+
+template <typename T>
+bool GeometryUtil<T>::lineStringWithinPolygon(const LineString<T>& line, const Polygon<T>& polygon) {
+    const auto length = line.size();
+    // First, check if geometry points of line segments are all inside polygon
+    for (std::size_t i = 0; i < length; ++i) {
+        if (!pointWithinPolygon(line[i], polygon)) {
+            return false;
+        }
+    }
+
+    // Second, check if there is line segment intersecting polygon edge
+    for (std::size_t i = 0; i < length - 1; ++i) {
+        if (lineIntersectPolygon(line[i], line[i + 1], polygon)) {
+            return false;
+        }
+    }
+    return true;
+}
+
+template <typename T>
+bool GeometryUtil<T>::lineStringWithinPolygons(const LineString<T>& line, const MultiPolygon<T>& polygons) {
+    for (const auto& polygon : polygons) {
+        if (lineStringWithinPolygon(line, polygon)) return true;
+    }
+    return false;
+}
+
+template struct GeometryUtil<double>;
+template struct GeometryUtil<int64_t>;
+
+} // namespace mbgl
diff --git a/src/mbgl/util/geometry_util.hpp b/src/mbgl/util/geometry_util.hpp
new file mode 100644
index 0000000000..ccf7bfc504
--- /dev/null
+++ b/src/mbgl/util/geometry_util.hpp
@@ -0,0 +1,25 @@
+#pragma once
+
+#include <mbgl/util/geometry.hpp>
+
+namespace mbgl {
+
+template <typename T>
+struct GeometryUtil {
+    using type = T;
+    static bool segmentIntersectSegment(const Point<type>& a,
+                                        const Point<type>& b,
+                                        const Point<type>& c,
+                                        const Point<type>& d);
+    static bool rayIntersect(const Point<type>& p, const Point<type>& p1, const Point<type>& p2);
+    static bool pointOnBoundary(const Point<type>& p, const Point<type>& p1, const Point<type>& p2);
+    static bool lineIntersectPolygon(const Point<type>& p1, const Point<type>& p2, const Polygon<type>& polygon);
+    static bool pointWithinPolygon(const Point<type>& point, const Polygon<type>& polygon, bool trueOnBoundary = false);
+    static bool pointWithinPolygons(const Point<type>& point,
+                                    const MultiPolygon<type>& polygons,
+                                    bool trueOnBoundary = false);
+    static bool lineStringWithinPolygon(const LineString<type>& line, const Polygon<type>& polygon);
+    static bool lineStringWithinPolygons(const LineString<type>& line, const MultiPolygon<type>& polygons);
+};
+
+} // namespace mbgl
diff --git a/src/mbgl/util/geometry_within.cpp b/src/mbgl/util/geometry_within.cpp
deleted file mode 100644
index 618567a8e1..0000000000
--- a/src/mbgl/util/geometry_within.cpp
+++ /dev/null
@@ -1,143 +0,0 @@
-#include <mbgl/util/geometry_within.hpp>
-
-#include <algorithm>
-
-namespace mbgl {
-
-namespace {
-bool rayIntersect(const Point<int64_t>& p, const Point<int64_t>& p1, const Point<int64_t>& p2) {
-    return ((p1.y > p.y) != (p2.y > p.y)) && (p.x < (p2.x - p1.x) * (p.y - p1.y) / (p2.y - p1.y) + p1.x);
-}
-
-// check if point p in on line segment with end points p1 and p2
-bool onBoundary(const Point<int64_t>& p, const Point<int64_t>& p1, const Point<int64_t>& p2) {
-    // requirements of point p on line segment:
-    // 1. colinear: cross product of vector p->p1(x1, y1) and vector p->p2(x2, y2) equals to 0
-    // 2. p is between p1 and p2
-    const auto x1 = p.x - p1.x;
-    const auto y1 = p.y - p1.y;
-    const auto x2 = p.x - p2.x;
-    const auto y2 = p.y - p2.y;
-    return (x1 * y2 - x2 * y1 == 0) && (x1 * x2 <= 0) && (y1 * y2 <= 0);
-}
-
-// a, b are end points for line segment1, c and d are end points for line segment2
-bool lineIntersectLine(const Point<int64_t>& a,
-                       const Point<int64_t>& b,
-                       const Point<int64_t>& c,
-                       const Point<int64_t>& d) {
-    const auto perp = [](const Point<int64_t>& v1, const Point<int64_t>& v2) { return (v1.x * v2.y - v1.y * v2.x); };
-
-    // check if two segments are parallel or not
-    // precondition is end point a, b is inside polygon, if line a->b is
-    // parallel to polygon edge c->d, then a->b won't intersect with c->d
-    auto vectorP = Point<int64_t>(b.x - a.x, b.y - a.y);
-    auto vectorQ = Point<int64_t>(d.x - c.x, d.y - c.y);
-    if (perp(vectorQ, vectorP) == 0) return false;
-
-    // check if p1 and p2 are in different sides of line segment q1->q2
-    const auto twoSided =
-        [](const Point<int64_t>& p1, const Point<int64_t>& p2, const Point<int64_t>& q1, const Point<int64_t>& q2) {
-            int64_t x1;
-            int64_t y1;
-            int64_t x2;
-            int64_t y2;
-            int64_t x3;
-            int64_t y3;
-
-            // q1->p1 (x1, y1), q1->p2 (x2, y2), q1->q2 (x3, y3)
-            x1 = p1.x - q1.x;
-            y1 = p1.y - q1.y;
-            x2 = p2.x - q1.x;
-            y2 = p2.y - q1.y;
-            x3 = q2.x - q1.x;
-            y3 = q2.y - q1.y;
-            auto ret1 = (x1 * y3 - x3 * y1);
-            auto ret2 = (x2 * y3 - x3 * y2);
-            return (ret1 > 0 && ret2 < 0) || (ret1 < 0 && ret2 > 0);
-        };
-
-    // If lines are intersecting with each other, the relative location should be:
-    // a and b lie in different sides of segment c->d
-    // c and d lie in different sides of segment a->b
-    return twoSided(a, b, c, d) && twoSided(c, d, a, b);
-}
-
-bool lineIntersectPolygon(const Point<int64_t>& p1, const Point<int64_t>& p2, const Polygon<int64_t>& polygon) {
-    for (auto ring : polygon) {
-        auto length = ring.size();
-        // loop through every edge of the ring
-        for (std::size_t i = 0; i < length - 1; ++i) {
-            if (lineIntersectLine(p1, p2, ring[i], ring[i + 1])) {
-                return true;
-            }
-        }
-    }
-    return false;
-}
-
-} // namespace
-
-void updateBBox(WithinBBox& bbox, const Point<int64_t>& p) {
-    bbox[0] = std::min(p.x, bbox[0]);
-    bbox[1] = std::min(p.y, bbox[1]);
-    bbox[2] = std::max(p.x, bbox[2]);
-    bbox[3] = std::max(p.y, bbox[3]);
-}
-
-bool boxWithinBox(const WithinBBox& bbox1, const WithinBBox& bbox2) {
-    if (bbox1[0] <= bbox2[0]) return false;
-    if (bbox1[2] >= bbox2[2]) return false;
-    if (bbox1[1] <= bbox2[1]) return false;
-    if (bbox1[3] >= bbox2[3]) return false;
-    return true;
-}
-
-// ray casting algorithm for detecting if point is in polygon
-bool pointWithinPolygon(const Point<int64_t>& point, const Polygon<int64_t>& polygon) {
-    bool within = false;
-    for (const auto& ring : polygon) {
-        const auto length = ring.size();
-        // loop through every edge of the ring
-        for (std::size_t i = 0; i < length - 1; ++i) {
-            if (onBoundary(point, ring[i], ring[i + 1])) return false;
-            if (rayIntersect(point, ring[i], ring[i + 1])) {
-                within = !within;
-            }
-        }
-    }
-    return within;
-}
-
-bool pointWithinPolygons(const Point<int64_t>& point, const MultiPolygon<int64_t>& polygons) {
-    for (const auto& polygon : polygons) {
-        if (pointWithinPolygon(point, polygon)) return true;
-    }
-    return false;
-}
-
-bool lineStringWithinPolygon(const LineString<int64_t>& line, const Polygon<int64_t>& polygon) {
-    const auto length = line.size();
-    // First, check if geometry points of line segments are all inside polygon
-    for (std::size_t i = 0; i < length; ++i) {
-        if (!pointWithinPolygon(line[i], polygon)) {
-            return false;
-        }
-    }
-
-    // Second, check if there is line segment intersecting polygon edge
-    for (std::size_t i = 0; i < length - 1; ++i) {
-        if (lineIntersectPolygon(line[i], line[i + 1], polygon)) {
-            return false;
-        }
-    }
-    return true;
-}
-
-bool lineStringWithinPolygons(const LineString<int64_t>& line, const MultiPolygon<int64_t>& polygons) {
-    for (const auto& polygon : polygons) {
-        if (lineStringWithinPolygon(line, polygon)) return true;
-    }
-    return false;
-}
-} // namespace mbgl
diff --git a/src/mbgl/util/geometry_within.hpp b/src/mbgl/util/geometry_within.hpp
deleted file mode 100644
index b078783501..0000000000
--- a/src/mbgl/util/geometry_within.hpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#pragma once
-
-#include <array>
-#include <limits>
-#include <mbgl/util/geometry.hpp>
-
-namespace mbgl {
-
-// contains minX, minY, maxX, maxY
-using WithinBBox = std::array<int64_t, 4>;
-const WithinBBox DefaultBBox = WithinBBox{std::numeric_limits<int64_t>::max(),
-                                          std::numeric_limits<int64_t>::max(),
-                                          std::numeric_limits<int64_t>::min(),
-                                          std::numeric_limits<int64_t>::min()};
-
-// check if bbox1 is within bbox2
-bool boxWithinBox(const WithinBBox& bbox1, const WithinBBox& bbox2);
-
-void updateBBox(WithinBBox& bbox, const Point<int64_t>& p);
-
-bool pointWithinPolygon(const Point<int64_t>& point, const Polygon<int64_t>& polygon);
-
-bool pointWithinPolygons(const Point<int64_t>& point, const MultiPolygon<int64_t>& polygons);
-
-bool lineStringWithinPolygon(const LineString<int64_t>& lineString, const Polygon<int64_t>& polygon);
-
-bool lineStringWithinPolygons(const LineString<int64_t>& line, const MultiPolygon<int64_t>& polygons);
-
-} // namespace mbgl