diff --git a/engine/CMakeLists.txt b/engine/CMakeLists.txt
index f09040a..e46bdf7 100644
--- a/engine/CMakeLists.txt
+++ b/engine/CMakeLists.txt
@@ -29,9 +29,12 @@ add_library(engine
     src/evaluator/evaluator.cpp
     src/evaluator/knn.cpp
     src/query/knn.cpp
+    src/query/frechet_helper.cpp
     src/evaluator/where_at.cpp
     src/query/range_query.cpp
     src/evaluator/range.cpp
+    src/query/join.cpp
+    src/evaluator/join.cpp
 )
 
 target_include_directories(engine PUBLIC include)
@@ -62,6 +65,7 @@ add_executable(engine_tests
     tests/parallel_test.cpp
     tests/knn_test.cpp
     tests/config_test.cpp
+    tests/join_test.cpp
 )
 target_link_libraries(engine_tests PRIVATE engine GTest::gtest_main Arrow::arrow_shared Parquet::parquet_shared)
 gtest_discover_tests(engine_tests)
diff --git a/engine/configs/geolife.json b/engine/configs/geolife.json
index cabd7c5..29e6d8a 100644
--- a/engine/configs/geolife.json
+++ b/engine/configs/geolife.json
@@ -29,5 +29,9 @@
         "hit_radius_m": 100.0,
         "min_tolerance_m": 25.0,
         "max_tolerance_m": 500.0
+    },
+    "join": {
+        "query_count": 0,
+        "eps_m": 500.0
     }
 }
diff --git a/engine/include/backend/engine/evaluator/evaluator.hpp b/engine/include/backend/engine/evaluator/evaluator.hpp
index 72c4f2b..9d1c15b 100644
--- a/engine/include/backend/engine/evaluator/evaluator.hpp
+++ b/engine/include/backend/engine/evaluator/evaluator.hpp
@@ -1,8 +1,10 @@
 #pragma once
 
 #include "backend/engine/core/trajectory_db.hpp"
+#include "backend/engine/evaluator/join.hpp"
 #include "backend/engine/evaluator/knn.hpp"
 #include "backend/engine/evaluator/where_at.hpp"
+#include "backend/engine/query/join.hpp"
 #include "backend/engine/query/knn.hpp"
 #include "backend/engine/query/range_query.hpp"
 
@@ -15,6 +17,7 @@ struct Workloads {
     const query::KnnWorkload *knn = nullptr; // a pointer, as it also carries k and the query db
     std::span<const query::WhereAtQuery> where_at = {};
     std::span<const query::RangeQuery> range = {};
+    const query::JoinWorkload *join = nullptr; // a pointer, as it also carries eps and the queries
 };
 
 // Kendall's tau is in [-1, 1] and the others in [0, 1], or NaN if it could not be computed
@@ -22,6 +25,7 @@ struct Scoretable {
     KnnScores knn;
     double range_recall;
     WhereAtScores where_at;
+    JoinScores join;
 };
 
 struct EvaluationResult {
@@ -41,6 +45,7 @@ class Evaluator {
     KnnScores get_knn_scores() const;
     double get_range_recall() const;
     WhereAtScores get_where_at_scores() const;
+    JoinScores get_join_scores() const;
 
   private:
     const core::TrajectoryDb &m_simplified;
diff --git a/engine/include/backend/engine/evaluator/join.hpp b/engine/include/backend/engine/evaluator/join.hpp
new file mode 100644
index 0000000..99a3548
--- /dev/null
+++ b/engine/include/backend/engine/evaluator/join.hpp
@@ -0,0 +1,23 @@
+#pragma once
+
+#include "backend/engine/core/trajectory_db.hpp"
+#include "backend/engine/query/join.hpp"
+
+#include <cstddef>
+#include <limits>
+
+namespace evaluator {
+
+// Averages over queries. Precision skips queries with an empty simplified result, recall those
+// with an empty reference, F1 those where both are empty. NaN if no query counts.
+struct JoinScores {
+    double f1 = std::numeric_limits<double>::quiet_NaN();
+    double precision = std::numeric_limits<double>::quiet_NaN();
+    double recall = std::numeric_limits<double>::quiet_NaN();
+    double eps_m = std::numeric_limits<double>::quiet_NaN();
+    std::size_t query_count = 0;
+};
+
+JoinScores evaluate_join(const core::TrajectoryDb &compressed, const query::JoinWorkload &workload);
+
+} // namespace evaluator
diff --git a/engine/include/backend/engine/pipeline/config.hpp b/engine/include/backend/engine/pipeline/config.hpp
index 7b7c1c6..dff0a74 100644
--- a/engine/include/backend/engine/pipeline/config.hpp
+++ b/engine/include/backend/engine/pipeline/config.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include "backend/engine/evaluator/where_at.hpp"
+#include "backend/engine/query/join.hpp"
 #include "backend/engine/query/knn.hpp"
 #include "backend/engine/query/range_query.hpp"
 #include "backend/engine/query/where_at_query.hpp"
@@ -26,6 +27,7 @@ struct Config {
     query::RangeConfig range;
     query::WhereAtConfig where_at;
     evaluator::WhereAtScoringConfig where_at_scoring;
+    query::JoinConfig join;
 };
 
 // Reads a JSON config. Only data_path is required, a missing key keeps its default.
diff --git a/engine/include/backend/engine/query/frechet_helper.hpp b/engine/include/backend/engine/query/frechet_helper.hpp
new file mode 100644
index 0000000..53e6c03
--- /dev/null
+++ b/engine/include/backend/engine/query/frechet_helper.hpp
@@ -0,0 +1,34 @@
+#pragma once
+
+#include "backend/engine/core/trajectory_db.hpp"
+
+#include <vector>
+
+namespace query {
+
+// Squared instead of std::hypot, which is slow. Fréchet only takes max and min of distances,
+// Square still has order
+inline double squared_point_dist(double a_x, double a_y, double b_x, double b_y) {
+    const double dx = a_x - b_x;
+    const double dy = a_y - b_y;
+    return dx * dx + dy * dy;
+}
+
+// Squared discrete Fréchet distance. Values below max_squared_dist are exact, anything else is
+// at or above it.
+double get_squared_frechet_distance(const core::TrajectoryView &query,
+                                    const core::TrajectoryView &reference, double max_squared_dist,
+                                    std::vector<double> &prev, std::vector<double> &curr);
+
+double squared_endpoint_lower_bound(const core::TrajectoryView &query,
+                                    const core::TrajectoryView &reference);
+
+double frechet_distance(const core::TrajectoryView &a, const core::TrajectoryView &b);
+
+// Whether the discrete Fréchet distance is at most max_dist, up to a relative 1e-12 for
+// rounding. Stops as soon as it is not.
+bool frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b, double max_dist);
+
+double endpoint_lower_bound(const core::TrajectoryView &a, const core::TrajectoryView &b);
+
+} // namespace query
diff --git a/engine/include/backend/engine/query/join.hpp b/engine/include/backend/engine/query/join.hpp
new file mode 100644
index 0000000..7e25e45
--- /dev/null
+++ b/engine/include/backend/engine/query/join.hpp
@@ -0,0 +1,37 @@
+#pragma once
+
+#include "backend/engine/core/trajectory_db.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <vector>
+
+namespace query {
+
+struct JoinConfig {
+    std::size_t query_count = 10; // 0 = use every trajectory in the query pool
+    double eps_m = 500.0;         // join threshold in metres, must be positive
+};
+
+struct JoinQuery {
+    std::size_t query_index;
+    std::vector<std::size_t> reference;
+    std::vector<std::size_t> candidates;
+};
+
+struct JoinWorkload {
+    core::TrajectoryDb queries;
+    std::vector<JoinQuery> entries;
+    double eps_m = 0.0;
+};
+
+// Samples query_count trajectories from query_pool (the held-out split) and stores their join
+// result in original as the reference answer
+JoinWorkload generate_join_queries(const core::TrajectoryDb &original,
+                                   const core::TrajectoryDb &query_pool,
+                                   const JoinConfig &config = {}, std::uint64_t seed = 42);
+
+std::vector<std::size_t> join(const core::TrajectoryDb &trajectories,
+                              const core::TrajectoryView &query, double eps);
+
+} // namespace query
diff --git a/engine/src/evaluator/evaluator.cpp b/engine/src/evaluator/evaluator.cpp
index c1ed296..e34f47d 100644
--- a/engine/src/evaluator/evaluator.cpp
+++ b/engine/src/evaluator/evaluator.cpp
@@ -25,7 +25,8 @@ EvaluationResult Evaluator::evaluate() const {
     auto range = evaluate_range(m_simplified, m_workloads.range);
     return {.scores = {.knn = knn.scores,
                        .range_recall = range.recall,
-                       .where_at = get_where_at_scores()},
+                       .where_at = get_where_at_scores(),
+                       .join = get_join_scores()},
             .knn = std::move(knn.answers),
             .range = std::move(range.answers)};
 }
@@ -42,4 +43,8 @@ WhereAtScores Evaluator::get_where_at_scores() const {
     return evaluate_where_at(m_simplified, m_workloads.where_at, m_where_at_config);
 }
 
+JoinScores Evaluator::get_join_scores() const {
+    return m_workloads.join ? evaluate_join(m_simplified, *m_workloads.join) : JoinScores{};
+}
+
 } // namespace evaluator
diff --git a/engine/src/evaluator/join.cpp b/engine/src/evaluator/join.cpp
new file mode 100644
index 0000000..13648b9
--- /dev/null
+++ b/engine/src/evaluator/join.cpp
@@ -0,0 +1,73 @@
+#include "backend/engine/evaluator/join.hpp"
+#include "backend/engine/core/parallel.hpp"
+#include "backend/engine/query/frechet_helper.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <vector>
+
+namespace evaluator {
+
+namespace {
+
+struct Counts {
+    std::size_t found = 0;
+    std::size_t hits = 0;
+};
+
+double average(const std::vector<double> &values) {
+    double sum = 0.0;
+    std::size_t counted = 0;
+    for (const double value : values) {
+        if (!std::isnan(value)) {
+            sum += value;
+            ++counted;
+        }
+    }
+    return counted == 0 ? std::numeric_limits<double>::quiet_NaN()
+                        : sum / static_cast<double>(counted);
+}
+
+} // namespace
+
+JoinScores evaluate_join(const core::TrajectoryDb &compressed,
+                         const query::JoinWorkload &workload) {
+    const std::size_t n = workload.entries.size();
+    const double nan = std::numeric_limits<double>::quiet_NaN();
+    std::vector<double> f1(n, nan), precision(n, nan), recall(n, nan);
+
+    core::parallel_for(n, [&](std::size_t i) {
+        const auto &entry = workload.entries[i];
+        const auto query = workload.queries.trajectory(entry.query_index);
+
+        // Simplification keeps endpoints, so the original candidates cover every possible match
+        Counts counts;
+        for (const auto j : entry.candidates) {
+            if (!query::frechet_within(query, compressed.trajectory(j), workload.eps_m))
+                continue;
+            ++counts.found;
+            if (std::ranges::binary_search(entry.reference, j))
+                ++counts.hits;
+        }
+
+        const auto hits = static_cast<double>(counts.hits);
+        const auto found = static_cast<double>(counts.found);
+        const auto truth = static_cast<double>(entry.reference.size());
+        if (counts.found != 0)
+            precision[i] = hits / found;
+        if (!entry.reference.empty())
+            recall[i] = hits / truth;
+        if (counts.found != 0 || !entry.reference.empty())
+            f1[i] = 2.0 * hits / (found + truth);
+    });
+
+    return {.f1 = average(f1),
+            .precision = average(precision),
+            .recall = average(recall),
+            .eps_m = workload.eps_m,
+            .query_count = n};
+}
+
+} // namespace evaluator
diff --git a/engine/src/pipeline/config.cpp b/engine/src/pipeline/config.cpp
index ab667c3..9a43674 100644
--- a/engine/src/pipeline/config.cpp
+++ b/engine/src/pipeline/config.cpp
@@ -84,10 +84,16 @@ void parse_where_at_scoring_config(Config &config, const json &scoring) {
         scoring.value("max_tolerance_m", config.where_at_scoring.max_tolerance_m);
 }
 
+void parse_join_config(Config &config, const json &join) {
+    validate_keys(join, {"query_count", "eps_m"}, "join");
+    config.join.query_count = join.value("query_count", config.join.query_count);
+    config.join.eps_m = join.value("eps_m", config.join.eps_m);
+}
+
 Config parse_config(const json &data) {
     validate_keys(data,
                   {"data_path", "export_results", "export_dir", "seed", "budget_shares",
-                   "algorithms", "split", "knn", "range", "where_at", "where_at_scoring"},
+                   "algorithms", "split", "knn", "range", "where_at", "where_at_scoring", "join"},
                   "root");
     Config config;
 
@@ -117,6 +123,7 @@ Config parse_config(const json &data) {
     parse_range_config(config, data.value("range", json::object()));
     parse_where_at_config(config, data.value("where_at", json::object()));
     parse_where_at_scoring_config(config, data.value("where_at_scoring", json::object()));
+    parse_join_config(config, data.value("join", json::object()));
 
     return config;
 }
diff --git a/engine/src/pipeline/pipeline.cpp b/engine/src/pipeline/pipeline.cpp
index 197e44a..d45222d 100644
--- a/engine/src/pipeline/pipeline.cpp
+++ b/engine/src/pipeline/pipeline.cpp
@@ -5,6 +5,7 @@
 #include "backend/engine/evaluator/evaluator.hpp"
 #include "backend/engine/io/parquet_loader.hpp"
 #include "backend/engine/io/workload_export.hpp"
+#include "backend/engine/query/join.hpp"
 #include "backend/engine/query/knn.hpp"
 #include "backend/engine/query/range_query.hpp"
 #include "backend/engine/query/where_at_query.hpp"
@@ -54,6 +55,8 @@ void print_scores(const evaluator::Scoretable &scores) {
     std::cout << "  where-at: within radius=" << scores.where_at.within_radius
               << ", average hit=" << scores.where_at.average_hit << '\n';
     std::cout << "  range recall: " << scores.range_recall << '\n';
+    std::cout << "  join F1: " << scores.join.f1 << ", precision: " << scores.join.precision
+              << ", recall: " << scores.join.recall << '\n';
 }
 
 } // namespace
@@ -71,6 +74,12 @@ void run(const Config &config) {
               << where_at.skipped_trajectories << " trips without eligible time\n";
     const auto range = query::generate_range_queries(database, config.range, config.seed);
     std::cout << "range: " << range.queries.size() << " queries\n";
+    const auto join = query::generate_join_queries(database, queries, config.join, config.seed);
+    std::size_t join_pairs = 0;
+    for (const auto &entry : join.entries)
+        join_pairs += entry.reference.size();
+    std::cout << "join: " << join.entries.size() << " queries, eps=" << join.eps_m << " m, "
+              << join_pairs << " matching pairs\n";
 
     if (config.export_results) {
         io::export_workloads(config.export_dir, database, range.queries, knn_workload);
@@ -87,10 +96,12 @@ void run(const Config &config) {
             std::cout << share * 100 << "% " << magic_enum::enum_name(algorithm) << ": ";
             print_size("compressed", compressed);
 
-            const evaluator::Evaluator evaluator(
-                database, compressed,
-                {.knn = &knn_workload, .where_at = where_at.queries, .range = range.queries},
-                config.where_at_scoring);
+            const evaluator::Evaluator evaluator(database, compressed,
+                                                 {.knn = &knn_workload,
+                                                  .where_at = where_at.queries,
+                                                  .range = range.queries,
+                                                  .join = &join},
+                                                 config.where_at_scoring);
             const auto evaluation = evaluator.evaluate();
             print_scores(evaluation.scores);
 
diff --git a/engine/src/query/frechet_helper.cpp b/engine/src/query/frechet_helper.cpp
new file mode 100644
index 0000000..ea76914
--- /dev/null
+++ b/engine/src/query/frechet_helper.cpp
@@ -0,0 +1,124 @@
+#include "backend/engine/query/frechet_helper.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <utility>
+#include <vector>
+
+namespace query {
+
+double get_squared_frechet_distance(const core::TrajectoryView &query,
+                                    const core::TrajectoryView &reference, double max_squared_dist,
+                                    std::vector<double> &prev, std::vector<double> &curr) {
+    const std::size_t n = query.size();
+    const std::size_t m = reference.size();
+    if (n == 0 || m == 0) {
+        return std::numeric_limits<double>::infinity();
+    }
+
+    prev.resize(m);
+    curr.resize(m);
+    const double inf = std::numeric_limits<double>::infinity();
+
+    // Only cells below the current k-th distance can improve the answer. Track the
+    // reachable span in each row instead of evaluating the entire n-by-m grid.
+    double value = 0.0;
+    std::size_t last = 0;
+    for (std::size_t j = 0; j < m; ++j) {
+        value = std::max(
+            value, squared_point_dist(query.x[0], query.y[0], reference.x[j], reference.y[j]));
+        if (value >= max_squared_dist)
+            break;
+        prev[j] = value;
+        last = j + 1;
+    }
+    if (last == 0)
+        return inf;
+
+    // Cells at or above the cutoff keep their value instead of becoming infinity. They can only
+    // lead to values at or above the cutoff, so every value below it stays exact.
+    std::size_t first = 0;
+    for (std::size_t i = 1; i < n; ++i) {
+        const double qx = query.x[i];
+        const double qy = query.y[i];
+        const auto dist = [&](std::size_t j) {
+            return squared_point_dist(qx, qy, reference.x[j], reference.y[j]);
+        };
+        std::size_t next_first = m;
+        std::size_t next_last = 0;
+
+        // Within the previous row's reachable span [first, last), the cell above always exists,
+        // and the diagonal does for every cell but the first
+        double left = std::max(prev[first], dist(first));
+        curr[first] = left;
+        if (left < max_squared_dist) {
+            next_first = first;
+            next_last = first + 1;
+        }
+        for (std::size_t j = first + 1; j < last; ++j) {
+            left = std::max(std::min({prev[j], prev[j - 1], left}), dist(j));
+            curr[j] = left;
+            if (left < max_squared_dist) {
+                next_first = std::min(next_first, j);
+                next_last = j + 1;
+            }
+        }
+
+        // Past the span only the cell on the left continues a path, plus the diagonal at j == last
+        for (std::size_t j = last; j < m; ++j) {
+            const double predecessor = j == last ? std::min(prev[last - 1], left) : left;
+            if (predecessor >= max_squared_dist)
+                break;
+            left = std::max(predecessor, dist(j));
+            if (left >= max_squared_dist)
+                break;
+            curr[j] = left;
+            next_first = std::min(next_first, j);
+            next_last = j + 1;
+        }
+
+        if (next_last == 0)
+            return inf;
+        first = next_first;
+        last = next_last;
+        std::swap(prev, curr);
+    }
+    return last == m ? prev[m - 1] : inf;
+}
+
+double squared_endpoint_lower_bound(const core::TrajectoryView &query,
+                                    const core::TrajectoryView &reference) {
+    if (query.size() == 0 || reference.size() == 0) {
+        return std::numeric_limits<double>::infinity();
+    }
+    const std::size_t q_last = query.size() - 1;
+    const std::size_t r_last = reference.size() - 1;
+    return std::max(squared_point_dist(query.x[0], query.y[0], reference.x[0], reference.y[0]),
+                    squared_point_dist(query.x[q_last], query.y[q_last], reference.x[r_last],
+                                       reference.y[r_last]));
+}
+
+double frechet_distance(const core::TrajectoryView &a, const core::TrajectoryView &b) {
+    std::vector<double> prev;
+    std::vector<double> curr;
+    return std::sqrt(
+        get_squared_frechet_distance(a, b, std::numeric_limits<double>::infinity(), prev, curr));
+}
+
+bool frechet_within(const core::TrajectoryView &a, const core::TrajectoryView &b, double max_dist) {
+    // Values below the cutoff are exact, so anything at most max_dist^2 is found. The relative
+    // slack absorbs rounding when max_dist is itself a distance, as sqrt(d)^2 can fall below d.
+    const double cutoff = std::nextafter(max_dist * max_dist * (1.0 + 1e-12),
+                                         std::numeric_limits<double>::infinity());
+    std::vector<double> prev;
+    std::vector<double> curr;
+    return get_squared_frechet_distance(a, b, cutoff, prev, curr) < cutoff;
+}
+
+double endpoint_lower_bound(const core::TrajectoryView &a, const core::TrajectoryView &b) {
+    return std::sqrt(squared_endpoint_lower_bound(a, b));
+}
+
+} // namespace query
diff --git a/engine/src/query/join.cpp b/engine/src/query/join.cpp
new file mode 100644
index 0000000..1fb4813
--- /dev/null
+++ b/engine/src/query/join.cpp
@@ -0,0 +1,72 @@
+#include "backend/engine/query/join.hpp"
+#include "backend/engine/core/parallel.hpp"
+#include "backend/engine/query/frechet_helper.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+namespace query {
+
+namespace {
+
+std::vector<std::size_t> endpoint_candidates(const core::TrajectoryDb &trajectories,
+                                             const core::TrajectoryView &query, double eps) {
+    std::vector<std::size_t> candidates;
+    for (std::size_t i = 0; i < trajectories.size(); ++i) {
+        if (endpoint_lower_bound(query, trajectories.trajectory(i)) <= eps)
+            candidates.push_back(i);
+    }
+    return candidates;
+}
+
+} // namespace
+
+std::vector<std::size_t> join(const core::TrajectoryDb &trajectories,
+                              const core::TrajectoryView &query, double eps) {
+    std::vector<std::size_t> result;
+    for (const auto i : endpoint_candidates(trajectories, query, eps)) {
+        if (frechet_within(query, trajectories.trajectory(i), eps))
+            result.push_back(i);
+    }
+    return result;
+}
+
+JoinWorkload generate_join_queries(const core::TrajectoryDb &original,
+                                   const core::TrajectoryDb &query_pool, const JoinConfig &config,
+                                   std::uint64_t seed) {
+    if (!std::isfinite(config.eps_m) || config.eps_m <= 0.0)
+        throw std::invalid_argument("join: eps_m must be finite and positive");
+
+    std::vector<std::size_t> picked(query_pool.size());
+    std::iota(picked.begin(), picked.end(), 0);
+    if (config.query_count != 0 && config.query_count < picked.size()) {
+        std::mt19937_64 rng(seed);
+        std::ranges::shuffle(picked, rng);
+        picked.resize(config.query_count);
+    }
+
+    JoinWorkload results{core::subset(query_pool, picked), {}, config.eps_m};
+    const std::size_t n = results.queries.size();
+    if (n == 0 || original.size() == 0)
+        return results;
+
+    results.entries.resize(n);
+    core::parallel_for(n, [&](std::size_t i) {
+        const auto query = results.queries.trajectory(i);
+        auto &entry = results.entries[i];
+        entry.query_index = i;
+        entry.candidates = endpoint_candidates(original, query, results.eps_m);
+        for (const auto j : entry.candidates) {
+            if (frechet_within(query, original.trajectory(j), results.eps_m))
+                entry.reference.push_back(j);
+        }
+    });
+    return results;
+}
+
+} // namespace query
diff --git a/engine/src/query/knn.cpp b/engine/src/query/knn.cpp
index 07dbfbb..8dd4a23 100644
--- a/engine/src/query/knn.cpp
+++ b/engine/src/query/knn.cpp
@@ -1,6 +1,7 @@
 #include "backend/engine/query/knn.hpp"
 #include "backend/engine/core/parallel.hpp"
 #include "backend/engine/core/trajectory_db.hpp"
+#include "backend/engine/query/frechet_helper.hpp"
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
@@ -16,105 +17,6 @@ namespace query {
 
 namespace {
 
-// Squared instead of std::hypot, which is slow. Fréchet only takes max and min of distances,
-// and squaring keeps their order, so the whole computation can stay squared.
-double squared_point_dist(double a_x, double a_y, double b_x, double b_y) {
-    const double dx = a_x - b_x;
-    const double dy = a_y - b_y;
-    return dx * dx + dy * dy;
-}
-
-double get_squared_frechet_distance(const core::TrajectoryView &query,
-                                    const core::TrajectoryView &reference, double max_squared_dist,
-                                    std::vector<double> &prev, std::vector<double> &curr) {
-    const std::size_t n = query.size();
-    const std::size_t m = reference.size();
-    if (n == 0 || m == 0) {
-        return std::numeric_limits<double>::infinity();
-    }
-
-    prev.resize(m);
-    curr.resize(m);
-    const double inf = std::numeric_limits<double>::infinity();
-
-    // Only cells below the current k-th distance can improve the answer. Track the
-    // reachable span in each row instead of evaluating the entire n-by-m grid.
-    double value = 0.0;
-    std::size_t last = 0;
-    for (std::size_t j = 0; j < m; ++j) {
-        value = std::max(
-            value, squared_point_dist(query.x[0], query.y[0], reference.x[j], reference.y[j]));
-        if (value >= max_squared_dist)
-            break;
-        prev[j] = value;
-        last = j + 1;
-    }
-    if (last == 0)
-        return inf;
-
-    // Cells at or above the cutoff keep their value instead of becoming infinity. They can only
-    // lead to values at or above the cutoff, so every value below it stays exact.
-    std::size_t first = 0;
-    for (std::size_t i = 1; i < n; ++i) {
-        const double qx = query.x[i];
-        const double qy = query.y[i];
-        const auto dist = [&](std::size_t j) {
-            return squared_point_dist(qx, qy, reference.x[j], reference.y[j]);
-        };
-        std::size_t next_first = m;
-        std::size_t next_last = 0;
-
-        // Within the previous row's reachable span [first, last), the cell above always exists,
-        // and the diagonal does for every cell but the first
-        double left = std::max(prev[first], dist(first));
-        curr[first] = left;
-        if (left < max_squared_dist) {
-            next_first = first;
-            next_last = first + 1;
-        }
-        for (std::size_t j = first + 1; j < last; ++j) {
-            left = std::max(std::min({prev[j], prev[j - 1], left}), dist(j));
-            curr[j] = left;
-            if (left < max_squared_dist) {
-                next_first = std::min(next_first, j);
-                next_last = j + 1;
-            }
-        }
-
-        // Past the span only the cell on the left continues a path, plus the diagonal at j == last
-        for (std::size_t j = last; j < m; ++j) {
-            const double predecessor = j == last ? std::min(prev[last - 1], left) : left;
-            if (predecessor >= max_squared_dist)
-                break;
-            left = std::max(predecessor, dist(j));
-            if (left >= max_squared_dist)
-                break;
-            curr[j] = left;
-            next_first = std::min(next_first, j);
-            next_last = j + 1;
-        }
-
-        if (next_last == 0)
-            return inf;
-        first = next_first;
-        last = next_last;
-        std::swap(prev, curr);
-    }
-    return last == m ? prev[m - 1] : inf;
-}
-
-double squared_endpoint_lower_bound(const core::TrajectoryView &query,
-                                    const core::TrajectoryView &reference) {
-    if (query.size() == 0 || reference.size() == 0) {
-        return std::numeric_limits<double>::infinity();
-    }
-    const std::size_t q_last = query.size() - 1;
-    const std::size_t r_last = reference.size() - 1;
-    return std::max(squared_point_dist(query.x[0], query.y[0], reference.x[0], reference.y[0]),
-                    squared_point_dist(query.x[q_last], query.y[q_last], reference.x[r_last],
-                                       reference.y[r_last]));
-}
-
 // Point distance to bounding rectangle of a given trajectory bounding box
 double squared_dist_to_box(double x, double y, const core::BoundingBox &box) {
     const double dx = std::max({box.x_min - x, 0.0, x - box.x_max});
diff --git a/engine/tests/config_test.cpp b/engine/tests/config_test.cpp
index 1aaf259..b9b12a5 100644
--- a/engine/tests/config_test.cpp
+++ b/engine/tests/config_test.cpp
@@ -32,6 +32,17 @@ TEST(Config, MissingKeysKeepDefaults) {
     EXPECT_EQ(config.range.query_count, defaults.range.query_count);
     EXPECT_EQ(config.where_at.probes_per_trajectory, defaults.where_at.probes_per_trajectory);
     EXPECT_DOUBLE_EQ(config.where_at_scoring.hit_radius_m, defaults.where_at_scoring.hit_radius_m);
+    EXPECT_EQ(config.join.query_count, defaults.join.query_count);
+    EXPECT_DOUBLE_EQ(config.join.eps_m, defaults.join.eps_m);
+}
+
+TEST(Config, ReadsJoinKeys) {
+    const auto config = pipeline::load_config(write_config("join.json", R"({
+        "data_path": "data.parquet",
+        "join": {"query_count": 0, "eps_m": 250}
+    })"));
+    EXPECT_EQ(config.join.query_count, 0);
+    EXPECT_DOUBLE_EQ(config.join.eps_m, 250);
 }
 
 TEST(Config, ReadsTopLevelAndNestedKeys) {
@@ -111,7 +122,7 @@ TEST(Config, RejectsUnknownKeysAndNestedSeeds) {
     EXPECT_THROW(
         pipeline::load_config(write_config("unknown_key.json", R"({"data_path": "d", "sead": 7})")),
         std::invalid_argument);
-    for (const auto &section : {"split", "knn", "range", "where_at"}) {
+    for (const auto &section : {"split", "knn", "range", "where_at", "join"}) {
         const auto contents =
             std::string(R"({"data_path": "d", ")") + section + R"(": {"seed": 7}})";
         EXPECT_THROW(pipeline::load_config(write_config("nested_seed.json", contents)),
diff --git a/engine/tests/join_test.cpp b/engine/tests/join_test.cpp
new file mode 100644
index 0000000..4404411
--- /dev/null
+++ b/engine/tests/join_test.cpp
@@ -0,0 +1,198 @@
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <iterator>
+#include <random>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include "backend/engine/core/simplification.hpp"
+#include "backend/engine/core/trajectory_db.hpp"
+#include "backend/engine/evaluator/join.hpp"
+#include "backend/engine/query/frechet_helper.hpp"
+#include "backend/engine/query/join.hpp"
+#include "backend/engine/selection/selector.hpp"
+
+namespace {
+
+using Track = std::vector<std::pair<double, double>>; // (lat, lon) per point
+
+core::TrajectoryDb make_db(const std::vector<Track> &tracks) {
+    std::vector<std::uint64_t> offsets{0};
+    std::vector<std::int32_t> ids;
+    std::vector<std::int64_t> t;
+    std::vector<double> lat;
+    std::vector<double> lon;
+    for (std::size_t i = 0; i < tracks.size(); ++i) {
+        ids.push_back(static_cast<std::int32_t>(i));
+        for (std::size_t k = 0; k < tracks[i].size(); ++k) {
+            t.push_back(static_cast<std::int64_t>(k));
+            lat.push_back(tracks[i][k].first);
+            lon.push_back(tracks[i][k].second);
+        }
+        offsets.push_back(t.size());
+    }
+    return {std::move(offsets), std::move(ids), std::move(t), std::move(lat), std::move(lon)};
+}
+
+// Random walks starting in a small area, so many are within a few hundred metres of each other
+std::vector<Track> random_tracks(std::mt19937_64 &rng, std::size_t count) {
+    std::uniform_real_distribution<double> start_lat(39.90, 39.92);
+    std::uniform_real_distribution<double> start_lon(116.30, 116.32);
+    std::uniform_int_distribution<std::size_t> length(2, 40);
+    std::normal_distribution<double> step(0.0, 0.0005);
+    std::vector<Track> tracks(count);
+    for (auto &track : tracks) {
+        double la = start_lat(rng);
+        double lo = start_lon(rng);
+        const std::size_t n = length(rng);
+        for (std::size_t k = 0; k < n; ++k) {
+            track.emplace_back(la, lo);
+            la += step(rng);
+            lo += step(rng);
+        }
+    }
+    return tracks;
+}
+
+double frechet(const core::TrajectoryView &a, const core::TrajectoryView &b) {
+    std::vector<std::vector<double>> dp(a.size(), std::vector<double>(b.size()));
+    for (std::size_t i = 0; i < a.size(); ++i) {
+        for (std::size_t j = 0; j < b.size(); ++j) {
+            const double d = std::hypot(a.x[i] - b.x[j], a.y[i] - b.y[j]);
+            if (i == 0 && j == 0)
+                dp[i][j] = d;
+            else if (i == 0)
+                dp[i][j] = std::max(d, dp[i][j - 1]);
+            else if (j == 0)
+                dp[i][j] = std::max(d, dp[i - 1][j]);
+            else
+                dp[i][j] = std::max(d, std::min({dp[i - 1][j], dp[i][j - 1], dp[i - 1][j - 1]}));
+        }
+    }
+    return dp.back().back();
+}
+
+} // namespace
+
+TEST(Join, FrechetHelpersMatchBruteForce) {
+    std::mt19937_64 rng(5);
+    const auto db = make_db(random_tracks(rng, 60));
+    for (std::size_t a = 0; a < db.size(); a += 7) {
+        for (std::size_t b = 0; b < db.size(); ++b) {
+            const double expected = frechet(db.trajectory(a), db.trajectory(b));
+            EXPECT_NEAR(query::frechet_distance(db.trajectory(a), db.trajectory(b)), expected,
+                        1e-9);
+            EXPECT_TRUE(query::frechet_within(db.trajectory(a), db.trajectory(b), expected));
+            if (expected > 0.0) {
+                EXPECT_FALSE(
+                    query::frechet_within(db.trajectory(a), db.trajectory(b), expected * 0.999));
+            }
+            EXPECT_LE(query::endpoint_lower_bound(db.trajectory(a), db.trajectory(b)),
+                      expected + 1e-9);
+        }
+    }
+}
+
+TEST(Join, MatchesBruteForce) {
+    std::mt19937_64 rng(11);
+    const auto db = make_db(random_tracks(rng, 200));
+    const auto queries = make_db(random_tracks(rng, 20));
+    for (const double eps : {100.0, 300.0, 800.0}) {
+        for (std::size_t q = 0; q < queries.size(); ++q) {
+            std::vector<std::size_t> expected;
+            for (std::size_t i = 0; i < db.size(); ++i) {
+                if (frechet(queries.trajectory(q), db.trajectory(i)) <= eps)
+                    expected.push_back(i);
+            }
+            EXPECT_EQ(query::join(db, queries.trajectory(q), eps), expected)
+                << "query " << q << ", eps " << eps;
+        }
+    }
+}
+
+TEST(Join, WorkloadStoresTheJoinOnTheOriginalData) {
+    std::mt19937_64 rng(13);
+    const auto db = make_db(random_tracks(rng, 150));
+    const auto pool = make_db(random_tracks(rng, 30));
+    const auto workload =
+        query::generate_join_queries(db, pool, {.query_count = 0, .eps_m = 300.0});
+    ASSERT_EQ(workload.entries.size(), pool.size());
+    std::size_t pairs = 0;
+    for (const auto &entry : workload.entries) {
+        const auto query = workload.queries.trajectory(entry.query_index);
+        EXPECT_EQ(entry.reference, query::join(db, query, workload.eps_m));
+        EXPECT_TRUE(std::ranges::includes(entry.candidates, entry.reference));
+        pairs += entry.reference.size();
+    }
+    EXPECT_GT(pairs, 0u);
+}
+
+TEST(Join, FixedEpsIsUsedAsGiven) {
+    std::mt19937_64 rng(17);
+    const auto db = make_db(random_tracks(rng, 50));
+    const auto pool = make_db(random_tracks(rng, 10));
+    const auto workload =
+        query::generate_join_queries(db, pool, {.query_count = 4, .eps_m = 400.0});
+    EXPECT_EQ(workload.entries.size(), 4u);
+    EXPECT_DOUBLE_EQ(workload.eps_m, 400.0);
+}
+
+TEST(Join, RejectsInvalidConfig) {
+    const auto db = make_db({{{39.9, 116.3}, {39.91, 116.3}}});
+    EXPECT_THROW(query::generate_join_queries(db, db, {.eps_m = -1.0}), std::invalid_argument);
+    EXPECT_THROW(query::generate_join_queries(db, db, {.eps_m = 0.0}), std::invalid_argument);
+}
+
+TEST(JoinEvaluator, UnsimplifiedDatabaseScoresPerfectly) {
+    std::mt19937_64 rng(19);
+    const auto db = make_db(random_tracks(rng, 120));
+    const auto pool = make_db(random_tracks(rng, 25));
+    const auto workload = query::generate_join_queries(db, pool, {.query_count = 0});
+    const auto scores = evaluator::evaluate_join(db, workload);
+    EXPECT_DOUBLE_EQ(scores.f1, 1.0);
+    EXPECT_DOUBLE_EQ(scores.precision, 1.0);
+    EXPECT_DOUBLE_EQ(scores.recall, 1.0);
+    EXPECT_EQ(scores.query_count, pool.size());
+}
+
+TEST(JoinEvaluator, SimplifiedScoresMatchBruteForceJoin) {
+    std::mt19937_64 rng(23);
+    const auto db = make_db(random_tracks(rng, 150));
+    const auto pool = make_db(random_tracks(rng, 20));
+    const auto workload = query::generate_join_queries(db, pool, {.query_count = 0});
+    const std::vector<std::size_t> budgets(db.size(), 3);
+    const auto compressed = core::materialize(
+        db, selection::simplify_database(db, selection::AlgorithmType::SED, budgets));
+
+    // F1 per query against a join run on the whole simplified database, averaged by hand
+    double f1_sum = 0.0;
+    std::size_t counted = 0;
+    for (const auto &entry : workload.entries) {
+        const auto found =
+            query::join(compressed, workload.queries.trajectory(entry.query_index), workload.eps_m);
+        std::vector<std::size_t> hits;
+        std::ranges::set_intersection(entry.reference, found, std::back_inserter(hits));
+        if (found.empty() && entry.reference.empty())
+            continue;
+        f1_sum += 2.0 * static_cast<double>(hits.size()) /
+                  static_cast<double>(found.size() + entry.reference.size());
+        ++counted;
+    }
+    const auto scores = evaluator::evaluate_join(compressed, workload);
+    ASSERT_GT(counted, 0u);
+    EXPECT_NEAR(scores.f1, f1_sum / static_cast<double>(counted), 1e-12);
+    EXPECT_LT(scores.f1, 1.0);
+}
+
+TEST(JoinEvaluator, EmptyWorkloadIsNaN) {
+    const auto db = make_db({{{39.9, 116.3}, {39.91, 116.3}}});
+    const query::JoinWorkload workload{core::subset(db, {}), {}, 100.0};
+    const auto scores = evaluator::evaluate_join(db, workload);
+    EXPECT_TRUE(std::isnan(scores.f1));
+    EXPECT_EQ(scores.query_count, 0u);
+}