diff --git a/Makefile b/Makefile
index 57227d4c..467393ea 100644
--- a/Makefile
+++ b/Makefile
@@ -156,6 +156,7 @@ test: \
 	test_append_vector \
 	test_attribute_store \
 	test_config_validator \
+	test_correct \
 	test_deque_map \
 	test_geom \
 	test_helpers \
@@ -193,6 +194,10 @@ src/config_schema.h: resources/config-schema.json
 
 src/config_validator.o: src/config_schema.h
 
+test_correct: \
+	test/correct.test.o
+	$(CXX) $(CXXFLAGS) -o test.correct $^ $(INC) $(LIB) $(LDFLAGS) && ./test.correct
+
 test_deque_map: \
 	test/deque_map.test.o
 	$(CXX) $(CXXFLAGS) -o test.deque_map $^ $(INC) $(LIB) $(LDFLAGS) && ./test.deque_map
diff --git a/include/geometry/correct.hpp b/include/geometry/correct.hpp
index fba96c69..a385d86c 100644
--- a/include/geometry/correct.hpp
+++ b/include/geometry/correct.hpp
@@ -23,13 +23,18 @@ namespace geometry {
 
 namespace impl {
 
-template<typename C, typename T>
-static inline void result_combine(C &result, T &&new_element)
+// boxes[i] is the envelope of result[i]. Polygons whose envelopes are disjoint
+// cannot intersect, so the box test skips most of the costly intersects() calls
+// (otherwise quadratic in the number of rings, e.g. a lake with 14000 islands).
+template<typename C, typename B, typename T>
+static inline void result_combine(C &result, std::vector<B> &boxes, T &&new_element)
 {
     result.push_back(std::forward<T>(new_element));
+    boxes.push_back(boost::geometry::return_envelope<B>(result.back()));
 
    	for(std::size_t i = 0; i < result.size() - 1; ) {
-        if(!boost::geometry::intersects(result[i], result.back())) {
+        if(boost::geometry::disjoint(boxes[i], boxes.back()) ||
+           !boost::geometry::intersects(result[i], result.back())) {
             ++i;
             continue;
         }
@@ -43,8 +48,10 @@ static inline void result_combine(C &result, T &&new_element)
 		}
 
        	result.back() = std::move(union_result[0]);
+       	boxes.back() = boost::geometry::return_envelope<B>(result.back());
        	result.erase(result.begin() + i);
-    } 
+       	boxes.erase(boxes.begin() + i);
+    }
 }
 
 struct pseudo_vertice_key
@@ -302,13 +309,15 @@ template<
 	>
 struct combine_non_zero_winding
 {
-	inline void operator()(multi_polygon_t &combined_outers, multi_polygon_t &combined_inners, polygon_t &poly) 
+	using box_t = boost::geometry::model::box<point_t>;
+
+	inline void operator()(multi_polygon_t &combined_outers, std::vector<box_t> &outer_boxes, multi_polygon_t &combined_inners, std::vector<box_t> &inner_boxes, polygon_t &poly)
 	{
 		if(boost::geometry::area(poly) > 0)
-			result_combine(combined_outers, std::move(poly));
+			result_combine(combined_outers, outer_boxes, std::move(poly));
 		else {
 			std::reverse(poly.outer().begin(), poly.outer().end());
-			result_combine(combined_inners, std::move(poly));
+			result_combine(combined_inners, inner_boxes, std::move(poly));
 		}
 	}
 };
@@ -320,7 +329,9 @@ template<
 	>
 struct combine_odd_even
 {
-	inline void operator()(multi_polygon_t &combined_outers, multi_polygon_t &combined_inners, polygon_t &poly) 
+	using box_t = boost::geometry::model::box<point_t>;
+
+	inline void operator()(multi_polygon_t &combined_outers, std::vector<box_t> &, multi_polygon_t &combined_inners, std::vector<box_t> &, polygon_t &poly)
 	{
 		if(boost::geometry::area(poly) < 0)
 			std::reverse(poly.outer().begin(), poly.outer().end());
@@ -346,11 +357,12 @@ static inline void correct(polygon_t const &input, multi_polygon_t &output, doub
 	// Calculate all outers and combine them if possible
 	multi_polygon_t combined_outers;
 	multi_polygon_t combined_inners;
+	std::vector<boost::geometry::model::box<point_t>> outer_boxes, inner_boxes;
 
 	for(auto &ring: outer_rings) {
 		polygon_t poly;
 		poly.outer() = std::move(ring);
-		combine(combined_outers, combined_inners, poly);
+		combine(combined_outers, outer_boxes, combined_inners, inner_boxes, poly);
 	}
 
 	// Calculate all inners and combine them if possible
@@ -362,8 +374,8 @@ static inline void correct(polygon_t const &input, multi_polygon_t &output, doub
 		correct(poly, new_inners, remove_spike_min_area, combine);
 
 		for(auto &poly: new_inners) {
-			result_combine(combined_inners, std::move(poly));
-		} 
+			result_combine(combined_inners, inner_boxes, std::move(poly));
+		}
 	}
 
 	// Cut out all inners from all the outers
@@ -379,13 +391,18 @@ template<
 	>
 static inline void correct(multi_polygon_t const &input, multi_polygon_t &output, double remove_spike_min_area, combine_function_t combine)
 {
+	using box_t = boost::geometry::model::box<point_t>;
+	std::vector<box_t> output_boxes;
+	for(auto const &polygon: output)
+		output_boxes.push_back(boost::geometry::return_envelope<box_t>(polygon));
+
 	for(auto const &polygon: input)
 	{
 		multi_polygon_t new_polygons;
 		correct(polygon, new_polygons, remove_spike_min_area, combine);
 
-		for(auto &new_polygon: new_polygons) 
-			result_combine(output, std::move(new_polygon));
+		for(auto &new_polygon: new_polygons)
+			result_combine(output, output_boxes, std::move(new_polygon));
 	}
 }
 
diff --git a/test/correct.test.cpp b/test/correct.test.cpp
new file mode 100644
index 00000000..7178d173
--- /dev/null
+++ b/test/correct.test.cpp
@@ -0,0 +1,72 @@
+#include <iostream>
+#include <cmath>
+#include <string>
+#include "external/minunit.h"
+#include "geometry/correct.hpp"
+
+typedef boost::geometry::model::d2::point_xy<double> point_t;
+typedef boost::geometry::model::polygon<point_t> polygon_t;
+typedef boost::geometry::model::multi_polygon<polygon_t> multi_polygon_t;
+
+static bool near(double a, double b) { return std::abs(a - b) < 1e-9; }
+
+MU_TEST(test_correct_merges_overlapping_holes) {
+	// Two overlapping holes are merged into one; a disjoint hole is kept as is
+	polygon_t input;
+	boost::geometry::read_wkt(
+		"POLYGON((0 0,0 10,10 10,10 0,0 0),"
+		"(1 1,3 1,3 3,1 3,1 1),(2 2,4 2,4 4,2 4,2 2),(6 6,7 6,7 7,6 7,6 6))", input);
+	multi_polygon_t output;
+	geometry::correct(input, output);
+
+	mu_check(output.size() == 1);
+	mu_check(output[0].inners().size() == 2);
+	mu_check(near(boost::geometry::area(output), 100 - 7 - 1));
+	mu_check(boost::geometry::is_valid(output));
+}
+
+MU_TEST(test_correct_merges_overlapping_outers) {
+	multi_polygon_t input;
+	boost::geometry::read_wkt(
+		"MULTIPOLYGON(((0 0,0 2,2 2,2 0,0 0)),((1 1,1 3,3 3,3 1,1 1)),((5 5,5 6,6 6,6 5,5 5)))", input);
+	multi_polygon_t output;
+	geometry::correct(input, output);
+
+	mu_check(output.size() == 2);
+	mu_check(near(boost::geometry::area(output), 4 + 4 - 1 + 1));
+	mu_check(boost::geometry::is_valid(output));
+}
+
+MU_TEST(test_correct_many_disjoint_holes) {
+	// 400 disjoint holes, e.g. a lake with many islands: all are kept
+	std::string wkt = "POLYGON((0 0,0 100,100 100,100 0,0 0)";
+	for (int x = 0; x < 20; x++)
+		for (int y = 0; y < 20; y++) {
+			double x0 = 1 + x * 5, y0 = 1 + y * 5;
+			wkt += ",(" + std::to_string(x0) + " " + std::to_string(y0) + "," + std::to_string(x0 + 2) + " " + std::to_string(y0) + "," +
+				std::to_string(x0 + 2) + " " + std::to_string(y0 + 2) + "," + std::to_string(x0) + " " + std::to_string(y0 + 2) + "," +
+				std::to_string(x0) + " " + std::to_string(y0) + ")";
+		}
+	wkt += ")";
+	polygon_t input;
+	boost::geometry::read_wkt(wkt, input);
+	multi_polygon_t output;
+	geometry::correct(input, output);
+
+	mu_check(output.size() == 1);
+	mu_check(output[0].inners().size() == 400);
+	mu_check(near(boost::geometry::area(output), 100 * 100 - 400 * 4));
+	mu_check(boost::geometry::is_valid(output));
+}
+
+MU_TEST_SUITE(test_suite_correct) {
+	MU_RUN_TEST(test_correct_merges_overlapping_holes);
+	MU_RUN_TEST(test_correct_merges_overlapping_outers);
+	MU_RUN_TEST(test_correct_many_disjoint_holes);
+}
+
+int main() {
+	MU_RUN_SUITE(test_suite_correct);
+	MU_REPORT();
+	return MU_EXIT_CODE;
+}