From 4a2e9ff021c5d844ebb1b2aeaf859ff04682203b Mon Sep 17 00:00:00 2001
From: Matthew Alexander <1710128+quiqueg@users.noreply.github.com>
Date: Sun, 6 Dec 2020 16:55:04 -0800
Subject: [PATCH 1/4] Implement dual marching squares (linear) smoothing, add
 tests

---
 src/contours.js       |  38 +++-
 test/contours-test.js | 495 +++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 531 insertions(+), 2 deletions(-)

diff --git a/src/contours.js b/src/contours.js
index 614e9ec..3f1527b 100644
--- a/src/contours.js
+++ b/src/contours.js
@@ -182,6 +182,42 @@ export default function() {
     });
   }
 
+  /**
+   * Applies smoothing to an iso-ring according to the Dual Marching Squares algorithm.
+   * @param {[number, number][]} ring The sorted list of (x,y) coordinates of points for a given iso-ring.
+   * @param {number[]} values The underlying grid values.
+   * @param {number} value The iso-value for this iso-ring.
+   */
+  function smoothLinearDual(ring, values, value) {
+    // The first step in Dual Marching Squares smoothing is linear interpolation.
+    smoothLinear(ring, values, value);
+
+    // Check if the first point is the same as the last point
+    var isCompleteRing = ring[0][0] === ring[ring.length - 1][0] &&
+      ring[0][1] === ring[ring.length - 1][1];
+
+    ring.forEach(function(point, i) {
+      var x = point[0],
+          y = point[1],
+          x1,
+          y1;
+
+      if (i < ring.length - 1) {
+        // Next point
+        x1 = ring[i + 1][0];
+        y1 = ring[i + 1][1];
+
+        // Set the current point to the midpoint between it and the next point.
+        point[0] = x + (x1 - x) / 2;
+        point[1] = y + (y1 - y) / 2;
+      } else if (isCompleteRing) {
+        // This is the last point, complete the ring by matching the first point
+        point[0] = ring[0][0];
+        point[1] = ring[0][1];
+      }
+    });
+  }
+
   contours.contour = contour;
 
   contours.size = function(_) {
@@ -196,7 +232,7 @@ export default function() {
   };
 
   contours.smooth = function(_) {
-    return arguments.length ? (smooth = _ ? smoothLinear : noop, contours) : smooth === smoothLinear;
+    return arguments.length ? (smooth = _ === 'linearDual' ? smoothLinearDual : _ ? smoothLinear : noop, contours) : smooth === smoothLinear || smooth === smoothLinearDual;
   };
 
   return contours;
diff --git a/test/contours-test.js b/test/contours-test.js
index 650c5e7..8cc92fa 100644
--- a/test/contours-test.js
+++ b/test/contours-test.js
@@ -24,6 +24,31 @@ tape("contours(values) returns the expected result for an empty polygon", functi
   test.end();
 });
 
+tape("contours.smooth('linearDual')(values) returns the expected result for an empty polygon", function(test) {
+  var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": []
+    }
+  ]);
+  test.end();
+});
+
+
+
 tape("contours(values) returns the expected result for a simple polygon", function(test) {
   var contours = d3.contours().size([10, 10]).thresholds([0.5]);
   test.deepEqual(contours([
@@ -53,6 +78,105 @@ tape("contours(values) returns the expected result for a simple polygon", functi
   test.end();
 });
 
+tape("contours.smooth(false)(values) returns the expected result for a simple polygon", function(test) {
+  var contours = d3.contours().smooth(false).size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [6, 7.5], [6, 6.5], [6, 5.5], [6, 4.5], [6, 3.5], // Right edge
+            [5.5, 3], [4.5, 3], [3.5, 3], // Top edge
+            [3, 3.5], [3, 4.5], [3, 5.5], [3, 6.5], [3, 7.5], // Left edge
+            [3.5, 8], [4.5, 8], [5.5, 8], // Bottom edge
+            [6, 7.5]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth(true)(values) returns the expected result for a simple polygon", function(test) {
+  var contours = d3.contours().smooth(true).size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [6, 7.5], [6, 6.5], [6, 5.5], [6, 4.5], [6, 3.5], // Right edge
+            [5.5, 3], [4.5, 3], [3.5, 3], // Top edge
+            [3, 3.5], [3, 4.5], [3, 5.5], [3, 6.5], [3, 7.5], // Left edge
+            [3.5, 8], [4.5, 8], [5.5, 8], // Bottom edge
+            [6, 7.5]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth('linearDual')(values) returns the expected result for a simple polygon", function(test) {
+  var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [6, 7], [6, 6], [6, 5], [6, 4], [5.75, 3.25], // Right edge
+            [5, 3], [4, 3], [3.25, 3.25], // Top edge
+            [3, 4], [3, 5], [3, 6], [3, 7], [3.25, 7.75], // Left edge
+            [4, 8], [5, 8], [5.75, 7.75], // Bottom edge
+            [6, 7]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
 tape("contours(values).contour(value) returns the expected result for a simple polygon", function(test) {
   var contours = d3.contours().size([10, 10]);
   test.deepEqual(contours.contour([
@@ -80,7 +204,106 @@ tape("contours(values).contour(value) returns the expected result for a simple p
   test.end();
 });
 
-tape("contours.smooth(false)(values) returns the expected result for a simple polygon", function(test) {
+tape("contours(values) returns the expected result for a corner polygon", function(test) {
+  var contours = d3.contours().size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+    1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [2, 1.5], [2, 0.5], // Right edge
+            [1.5, 0], [0.5, 0], // Top edge
+            [0, 0.5], [0, 1.5], // Left edge
+            [0.5, 2], [1.5, 2], // Bottom edge
+            [2, 1.5]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth(false)(values) returns the expected result for a polygon in the corner", function(test) {
+  var contours = d3.contours().smooth(false).size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+    1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [2, 1.5], [2, 0.5], // Right edge
+            [1.5, 0], [0.5, 0], // Top edge
+            [0, 0.5], [0, 1.5], // Left edge
+            [0.5, 2], [1.5, 2], // Bottom edge
+            [2, 1.5]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth('linearDual')(values) returns the expected result for polygon in the corner", function(test) {
+  var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+    1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [2, 1], [1.75, 0.25], // Right edge
+            [1, 0], [0.25, 0.25], // Top edge
+            [0, 1], [0.25, 1.75], // Left edge
+            [1, 2], [1.75, 1.75], // Bottom edge
+            [2, 1]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth(false)(values) returns the expected result for a complex polygon", function(test) {
   var contours = d3.contours().smooth(false).size([10, 10]).thresholds([0.5]);
   test.deepEqual(contours([
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
@@ -109,6 +332,72 @@ tape("contours.smooth(false)(values) returns the expected result for a simple po
   test.end();
 });
 
+tape("contours.smooth(true)(values) returns the expected result for a complex polygon", function(test) {
+  var contours = d3.contours().smooth(true).size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 2, 1, 2, 0, 0, 0, 0,
+    0, 0, 0, 2, 2, 2, 0, 0, 0, 0,
+    0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
+    0, 0, 0, 2, 2, 2, 0, 0, 0, 0,
+    0, 0, 0, 2, 1, 2, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [6.25, 7.5], [6.25, 6.5], [6, 5.5], [6.25, 4.5], [6.25, 3.5], // Right edge
+            [5.5, 2.75], [4.5, 3], [3.5, 2.75], // Top edge
+            [2.75, 3.5], [2.75, 4.5], [3, 5.5], [2.75, 6.5], [2.75, 7.5], // Left edge
+            [3.5,  8.25], [4.5,  8], [5.5, 8.25], // Bottom edge
+            [6.25, 7.5]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth('linearDual')(values) returns the expected result for a complex polygon", function(test) {
+  var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 2, 1, 2, 0, 0, 0, 0,
+    0, 0, 0, 2, 2, 2, 0, 0, 0, 0,
+    0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
+    0, 0, 0, 2, 2, 2, 0, 0, 0, 0,
+    0, 0, 0, 2, 1, 2, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            [6.25, 7], [6.125, 6], [6.125, 5], [6.25, 4], [5.875, 3.125], // Right edge
+            [5, 2.875], [4, 2.875], [3.125, 3.125], // Top edge
+            [2.75, 4], [2.875, 5], [2.875, 6], [2.75, 7], [3.125, 7.875], // Left edge
+            [4, 8.125], [5, 8.125], [5.875, 7.875], // Bottom edge
+            [6.25, 7]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
 tape("contours(values) returns the expected result for a polygon with a hole", function(test) {
   var contours = d3.contours().size([10, 10]).thresholds([0.5]);
   test.deepEqual(contours([
@@ -140,6 +429,90 @@ tape("contours(values) returns the expected result for a polygon with a hole", f
   test.end();
 });
 
+tape("contours.smooth(true)(values) returns the expected result for a polygon with a hole", function(test) {
+  var contours = d3.contours().smooth(true).size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            // Outer polygon, counter-clockwise path
+            [6, 7.5], [6, 6.5], [6, 5.5], [6, 4.5], [6, 3.5], // Right edge
+            [5.5, 3], [4.5, 3], [3.5, 3], // Top edge
+            [3, 3.5], [3, 4.5], [3, 5.5], [3, 6.5], [3, 7.5], // Left edge
+            [3.5, 8], [4.5, 8], [5.5, 8], // Bottom edge
+            [6, 7.5]
+          ],
+          [
+            // Inner polygon, clockwise path
+            [4.5, 7], // Bottom point
+            [4, 6.5], [4, 5.5], [4, 4.5], // Left edge
+            [4.5, 4], // Top point
+            [5, 4.5], [5, 5.5], [5, 6.5], // Right edge
+            [4.5, 7]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth('linearDual')(values) returns the expected result for a polygon with a hole", function(test) {
+  var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 0, 1, 0, 0, 0, 0,
+    0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            // Outer polygon, counter-clockwise path
+            [6, 7], [6, 6], [6, 5], [6, 4], [5.75, 3.25], // Right edge
+            [5, 3], [4, 3], [3.25, 3.25], // Top edge
+            [3, 4], [3, 5], [3, 6], [3, 7], [3.25, 7.75], // Left edge
+            [4, 8], [5, 8], [5.75, 7.75], // Bottom edge
+            [6, 7]
+          ],
+          [
+            // Inner polygon, clockwise path
+            [4.25, 6.75], // Bottom edge
+            [4, 6], [4, 5], // Left edge
+            [4.25, 4.25], [4.75, 4.25], // Top edge
+            [5, 5], [5, 6], // Right edge
+            [4.75, 6.75], [4.25, 6.75] // Bottom edge
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
 tape("contours(values) returns the expected result for a multipolygon", function(test) {
   var contours = d3.contours().size([10, 10]).thresholds([0.5]);
   test.deepEqual(contours([
@@ -207,6 +580,126 @@ tape("contours(values) returns the expected result for a multipolygon with holes
   test.end();
 });
 
+tape("contours.smooth(false)(values) returns the expected result for a multipolygon with holes", function(test) {
+  var contours = d3.contours().smooth(false).size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 1, 1, 1, 0, 1, 1, 1, 0, 0,
+    0, 1, 0, 1, 0, 1, 0, 1, 0, 0,
+    0, 1, 1, 1, 0, 1, 1, 1, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            // Outer left polygon, counter-clockwise path
+            [4, 5.5], [4, 4.5], [4, 3.5], // Right edge
+            [3.5, 3], [2.5, 3], [1.5, 3], // Top edge
+            [1, 3.5], [1, 4.5], [1, 5.5], // Left edge
+            [1.5, 6], [2.5, 6], [3.5, 6], // Bottom edge
+            [4,5.5]
+          ],
+          [
+            // Inner left polygon, clockwise path
+            [2.5, 5], // Bottom
+            [2, 4.5], // Left
+            [2.5, 4], // Top
+            [3, 4.5], // Right
+            [2.5, 5]
+          ]
+        ],
+        [
+          [
+            // Outer right polygon, counter-clockwise path
+            [8, 5.5], [8, 4.5], [8, 3.5], // Right edge
+            [7.5, 3], [6.5, 3], [5.5, 3], // Top edge
+            [5, 3.5], [5, 4.5], [5, 5.5], // Left edge
+            [5.5, 6], [6.5, 6], [7.5, 6], // Bottom edge
+            [8, 5.5]
+          ],
+          [
+            // Inner right polygon, clockwise path
+            [6.5, 5], // Bottom
+            [6, 4.5], // Left
+            [6.5, 4], // Top
+            [7, 4.5], // Right
+            [6.5, 5]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
+tape("contours.smooth('linearDual')(values) returns the expected result for a multipolygon with holes", function(test) {
+  var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
+  test.deepEqual(contours([
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 1, 1, 1, 0, 1, 1, 1, 0, 0,
+    0, 1, 0, 1, 0, 1, 0, 1, 0, 0,
+    0, 1, 1, 1, 0, 1, 1, 1, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+    0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+  ]), [
+    {
+      "type": "MultiPolygon",
+      "value": 0.5,
+      "coordinates": [
+        [
+          [
+            // Outer left polygon, counter-clockwise path
+            [4, 5], [4, 4], // Right edge
+            [3.75, 3.25], [3, 3], [2, 3], [1.25, 3.25], // Top edge
+            [1, 4], [1, 5], // Left edge
+            [1.25, 5.75], [2, 6], [3, 6], [3.75, 5.75], // Bottom edge
+            [4, 5]
+          ],
+          [
+            // Inner left polygon, clockwise path
+            [2.25, 4.75], // Bottom-left
+            [2.25, 4.25], // Top-left
+            [2.75, 4.25], // Top-right
+            [2.75, 4.75], // Bottom-right
+            [2.25, 4.75]
+          ]
+        ],
+        [
+          [
+            // Outer right polygon, counter-clockwise path
+            [8, 5], [8, 4], // Right edge
+            [7.75, 3.25], [7, 3], [6, 3], [5.25, 3.25], // Top edge
+            [5, 4], [5, 5], // Left edge
+            [5.25, 5.75], [6, 6], [7, 6], [7.75, 5.75], // Bottom edge
+            [8, 5]
+          ],
+          [
+            // Inner right polygon, clockwise path
+            [6.25, 4.75], // Bottom-left
+            [6.25, 4.25], // Top-left
+            [6.75, 4.25], // Top-right
+            [6.75, 4.75], // Bottom-right
+            [6.25, 4.75]
+          ]
+        ]
+      ]
+    }
+  ]);
+  test.end();
+});
+
 tape("contours.size(…) validates the specified size", function(test) {
   test.deepEqual(d3.contours().size([1, 2]).size(), [1, 2]);
   test.deepEqual(d3.contours().size([0, 0]).size(), [0, 0]);

From 07e58d394166aee028632f293d574ee989740386 Mon Sep 17 00:00:00 2001
From: Matthew Alexander <1710128+quiqueg@users.noreply.github.com>
Date: Sun, 6 Dec 2020 16:55:43 -0800
Subject: [PATCH 2/4] Add VS Code settings and build task

---
 .vscode/settings.json |  3 +++
 .vscode/tasks.json    | 16 ++++++++++++++++
 2 files changed, 19 insertions(+)
 create mode 100644 .vscode/settings.json
 create mode 100644 .vscode/tasks.json

diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 0000000..ad94bc7
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1,3 @@
+{
+  "npm.packageManager": "yarn"
+}
diff --git a/.vscode/tasks.json b/.vscode/tasks.json
new file mode 100644
index 0000000..bfd6d50
--- /dev/null
+++ b/.vscode/tasks.json
@@ -0,0 +1,16 @@
+{
+  "version": "2.0.0",
+  "tasks": [
+    {
+      "type": "npm",
+      "script": "test",
+      "group": {
+        "kind": "build",
+        "isDefault": true
+      },
+      "problemMatcher": [],
+      "label": "yarn: test",
+      "detail": "tape 'test/**/*-test.js' && eslint src"
+    }
+  ]
+}

From 1031d2c377ad9607c864f1376b623eebaece1f7f Mon Sep 17 00:00:00 2001
From: Matthew Alexander <1710128+quiqueg@users.noreply.github.com>
Date: Thu, 10 Dec 2020 23:01:46 -0800
Subject: [PATCH 3/4] Add documentation for new contour smoothing options

---
 README.md | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 54 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 5abb632..738fbf0 100644
--- a/README.md
+++ b/README.md
@@ -98,7 +98,60 @@ If *size* is specified, sets the expected size of the input *values* grid to the
 
 <a name="contours_smooth" href="#contours_smooth">#</a> <i>contours</i>.<b>smooth</b>([<i>smooth</i>]) · [Source](https://github.com/d3/d3-contour/blob/master/src/contours.js), [Examples](https://observablehq.com/@d3/contours-smooth)
 
-If *smooth* is specified, sets whether or not the generated contour polygons are smoothed using linear interpolation. If *smooth* is not specified, returns the current smoothing flag, which defaults to true.
+If *smooth* is specified, sets the smoothing method to use when generating the contour polygons. If *smooth* is not specified, returns the current smoothing method, which defaults to true.
+
+The available *smooth* options are:
+
+- `false`: Smoothing is disabled.
+- `true` or `"linear"` (Default): Linear interpolation smoothing, as described in the original [Marching Squares algorithm](https://en.wikipedia.org/wiki/Marching_squares).
+- `"linearDual"`: Dual linear interpolation smoothing, as described in the [Dual Marching Squares algorithm](https://ieeexplore.ieee.org/document/7459173).
+
+In general, the quality of each smoothing method is inversely proportional to its runtime performance.
+
+Additionally, the density of the source data impacts contour smoothness: the differences between the smoothing methods are more noticeable with low-density data than they are with high-density data. With _very_ high-density data, there is no discernible quality difference between the three methods.
+
+<table>
+  <thead>
+    <tr>
+      <th colspan="2" scope="colgroup"></th>
+      <th colspan="2" scope="colgroup">Low-density data</th>
+      <th colspan="2" scope="colgroup">High-density data</th>
+    </tr>
+    <tr>
+      <th colspan="2" scope="col">Smoothing method</th>
+      <th scope="col">Contour quality</th>
+      <th scope="col">Performance</th>
+      <th scope="col">Contour quality</th>
+      <th scope="col">Performance</th>
+    </tr>
+  </thead>
+  <tbody>
+    <tr>
+      <th scope="row">None</th>
+      <th scope="row"><code>false</code></th>
+      <td>Poor</td>
+      <td>Best</td>
+      <td>Good</td>
+      <td>Best</td>
+    </tr>
+    <tr>
+      <th scope="row">Linear</th>
+      <th scope="row"><code>true</code> or<br /><code>"linear"</code></th>
+      <td>Good</td>
+      <td>Good</td>
+      <td>Best</td>
+      <td>Poor</td>
+    </tr>
+    <tr>
+      <th scope="row">Dual Linear</th>
+      <th scope="row"><code>"linearDual"</code></th>
+      <td>Best</td>
+      <td>Good</td>
+      <td>Best</td>
+      <td>Poor</td>
+    </tr>
+  </tbody>
+</table>
 
 <a name="contours_thresholds" href="#contours_thresholds">#</a> <i>contours</i>.<b>thresholds</b>([<i>thresholds</i>]) · [Source](https://github.com/d3/d3-contour/blob/master/src/contours.js), [Examples](https://observablehq.com/@d3/volcano-contours)
 

From fa89f82b3ebeb9076ae4d11097e18c3029f23076 Mon Sep 17 00:00:00 2001
From: Matthew Alexander <1710128+quiqueg@users.noreply.github.com>
Date: Sat, 12 Dec 2020 23:06:10 -0800
Subject: [PATCH 4/4] Improve performance of dual smoothing, update README with
 performance numbers

---
 README.md       | 30 ++++++++++++++++--------------
 src/contours.js | 19 ++++++++-----------
 2 files changed, 24 insertions(+), 25 deletions(-)

diff --git a/README.md b/README.md
index 738fbf0..04bc6df 100644
--- a/README.md
+++ b/README.md
@@ -114,15 +114,15 @@ Additionally, the density of the source data impacts contour smoothness: the dif
   <thead>
     <tr>
       <th colspan="2" scope="colgroup"></th>
-      <th colspan="2" scope="colgroup">Low-density data</th>
-      <th colspan="2" scope="colgroup">High-density data</th>
+      <th colspan="2" scope="colgroup" style="text-align: center;">Low-density data</th>
+      <th colspan="2" scope="colgroup" style="text-align: center;">High-density data</th>
     </tr>
     <tr>
       <th colspan="2" scope="col">Smoothing method</th>
-      <th scope="col">Contour quality</th>
-      <th scope="col">Performance</th>
-      <th scope="col">Contour quality</th>
-      <th scope="col">Performance</th>
+      <th scope="col" style="text-align: center;">Contour<br/>quality</th>
+      <th scope="col" style="text-align: center;">Performance<br/>cost*</th>
+      <th scope="col" style="text-align: center;">Contour<br/>quality</th>
+      <th scope="col" style="text-align: center;">Performance<br/>cost*</th>
     </tr>
   </thead>
   <tbody>
@@ -130,29 +130,31 @@ Additionally, the density of the source data impacts contour smoothness: the dif
       <th scope="row">None</th>
       <th scope="row"><code>false</code></th>
       <td>Poor</td>
-      <td>Best</td>
+      <td>+0%</td>
       <td>Good</td>
-      <td>Best</td>
+      <td>+0%</td>
     </tr>
     <tr>
-      <th scope="row">Linear</th>
+      <th scope="row">Linear (Default)</th>
       <th scope="row"><code>true</code> or<br /><code>"linear"</code></th>
       <td>Good</td>
-      <td>Good</td>
+      <td>+2.40%</td>
       <td>Best</td>
-      <td>Poor</td>
+      <td>+1.55%</td>
     </tr>
     <tr>
-      <th scope="row">Dual Linear</th>
+      <th scope="row">Dual linear</th>
       <th scope="row"><code>"linearDual"</code></th>
       <td>Best</td>
-      <td>Good</td>
+      <td>+2.70%</td>
       <td>Best</td>
-      <td>Poor</td>
+      <td>+2.00%</td>
     </tr>
   </tbody>
 </table>
 
+\* Estimated performance cost based on 2 sample datasets.
+
 <a name="contours_thresholds" href="#contours_thresholds">#</a> <i>contours</i>.<b>thresholds</b>([<i>thresholds</i>]) · [Source](https://github.com/d3/d3-contour/blob/master/src/contours.js), [Examples](https://observablehq.com/@d3/volcano-contours)
 
 If *thresholds* is specified, sets the threshold generator to the specified function or array and returns this contour generator. If *thresholds* is not specified, returns the current threshold generator, which by default implements [Sturges’ formula](https://github.com/d3/d3-array/blob/master/README.md#thresholdSturges).
diff --git a/src/contours.js b/src/contours.js
index 3f1527b..6b7b6fb 100644
--- a/src/contours.js
+++ b/src/contours.js
@@ -189,18 +189,15 @@ export default function() {
    * @param {number} value The iso-value for this iso-ring.
    */
   function smoothLinearDual(ring, values, value) {
+    var point, x, y, x1, y1;
+
     // The first step in Dual Marching Squares smoothing is linear interpolation.
     smoothLinear(ring, values, value);
 
-    // Check if the first point is the same as the last point
-    var isCompleteRing = ring[0][0] === ring[ring.length - 1][0] &&
-      ring[0][1] === ring[ring.length - 1][1];
-
-    ring.forEach(function(point, i) {
-      var x = point[0],
-          y = point[1],
-          x1,
-          y1;
+    for (var i = 0; i < ring.length; i++) {
+      point = ring[i];
+      x = point[0];
+      y = point[1];
 
       if (i < ring.length - 1) {
         // Next point
@@ -210,12 +207,12 @@ export default function() {
         // Set the current point to the midpoint between it and the next point.
         point[0] = x + (x1 - x) / 2;
         point[1] = y + (y1 - y) / 2;
-      } else if (isCompleteRing) {
+      } else  {
         // This is the last point, complete the ring by matching the first point
         point[0] = ring[0][0];
         point[1] = ring[0][1];
       }
-    });
+    }
   }
 
   contours.contour = contour;