From 59519a8c6ad663a6b4ce2152e0d9b07954e3394a Mon Sep 17 00:00:00 2001
From: Justin Lin <caterpillar@openhome.cc>
Date: Fri, 20 Mar 2020 12:58:13 +0800
Subject: [PATCH] add marching_squares

---
 .../_impl/_marching_squares_impl.scad         | 88 +++++++++++++++++++
 .../demo/marching_squares_demo.scad           | 19 ++++
 src/experimental/marching_squares.scad        | 20 +++++
 3 files changed, 127 insertions(+)
 create mode 100644 src/experimental/_impl/_marching_squares_impl.scad
 create mode 100644 src/experimental/demo/marching_squares_demo.scad
 create mode 100644 src/experimental/marching_squares.scad

diff --git a/src/experimental/_impl/_marching_squares_impl.scad b/src/experimental/_impl/_marching_squares_impl.scad
new file mode 100644
index 00000000..10415410
--- /dev/null
+++ b/src/experimental/_impl/_marching_squares_impl.scad
@@ -0,0 +1,88 @@
+function pn_label(pts, sigma) =
+    [
+        for(row = pts)
+        [for(p = row) [p[0], p[1], p[2], p[2] - sigma >= 0]]
+    ]; 
+    
+function corner_value(cell_pts) =
+    let(
+        c0 = cell_pts[0][3] ? 1 : 0,
+        c1 = cell_pts[1][3] ? 2 : 0,
+        c2 = cell_pts[2][3] ? 4 : 0,
+        c3 = cell_pts[3][3] ? 8 : 0
+    )
+    c0 + c1 + c2 + c3;
+
+function interpolated_pt(p0, p1, sigma) =
+    let(
+        x0 = p0[0],
+        y0 = p0[1],
+        z0 = p0[2],
+        v = p1 - p0,
+        t = (sigma - z0) / v[2]
+    )
+    [x0 + v[0] * t, y0 + v[1] * t, sigma];
+
+function case1_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[0], cell_pts[1], sigma), interpolated_pt(cell_pts[0], cell_pts[3], sigma)]
+];
+
+function case2_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[0], cell_pts[3], sigma), interpolated_pt(cell_pts[2], cell_pts[3], sigma)]
+];
+
+function case3_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[0], cell_pts[1], sigma), interpolated_pt(cell_pts[2], cell_pts[3], sigma)]
+];
+
+function case4_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[1], cell_pts[2], sigma), interpolated_pt(cell_pts[2], cell_pts[3], sigma)]
+];
+
+function case5_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[0], cell_pts[1], sigma), interpolated_pt(cell_pts[1], cell_pts[2], sigma)],
+    [interpolated_pt(cell_pts[0], cell_pts[3], sigma), interpolated_pt(cell_pts[2], cell_pts[3], sigma)]
+];
+
+function case6_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[1], cell_pts[2], sigma), interpolated_pt(cell_pts[0], cell_pts[3], sigma)]
+];
+
+function case7_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[0], cell_pts[1], sigma), interpolated_pt(cell_pts[1], cell_pts[2], sigma)]
+];
+    
+function case8_contours(cell_pts, sigma) = case7_contours(cell_pts, sigma);
+
+function case9_contours(cell_pts, sigma) = case6_contours(cell_pts, sigma);
+
+function case10_contours(cell_pts, sigma) = [
+    [interpolated_pt(cell_pts[0], cell_pts[1], sigma), interpolated_pt(cell_pts[0], cell_pts[3], sigma)],
+    [interpolated_pt(cell_pts[1], cell_pts[2], sigma), interpolated_pt(cell_pts[2], cell_pts[3], sigma)]
+];
+    
+function case11_contours(cell_pts, sigma) = case4_contours(cell_pts, sigma);
+
+function case12_contours(cell_pts, sigma) = case3_contours(cell_pts, sigma);
+
+function case13_contours(cell_pts, sigma) = case2_contours(cell_pts, sigma);
+
+function case14_contours(cell_pts, sigma) = case1_contours(cell_pts, sigma);
+
+function contours_of(cell_pts, sigma) =
+    let(cv = corner_value(cell_pts))
+    cv == 0  ? [] :
+    cv == 1  ? case1_contours(cell_pts, sigma) :
+    cv == 8  ? case2_contours(cell_pts, sigma) :
+    cv == 9  ? case3_contours(cell_pts, sigma) :
+    cv == 4  ? case4_contours(cell_pts, sigma) :
+    cv == 5  ? case5_contours(cell_pts, sigma) :
+    cv == 12 ? case6_contours(cell_pts, sigma) :
+    cv == 13 ? case7_contours(cell_pts, sigma) :
+    cv == 2  ? case8_contours(cell_pts, sigma) :
+    cv == 3  ? case9_contours(cell_pts, sigma) :
+    cv == 10 ? case10_contours(cell_pts, sigma) :
+    cv == 11 ? case11_contours(cell_pts, sigma) :
+    cv == 6  ? case12_contours(cell_pts, sigma) :
+    cv == 7  ? case13_contours(cell_pts, sigma) :
+    cv == 14 ? case14_contours(cell_pts, sigma) : [];
diff --git a/src/experimental/demo/marching_squares_demo.scad b/src/experimental/demo/marching_squares_demo.scad
new file mode 100644
index 00000000..d86bc043
--- /dev/null
+++ b/src/experimental/demo/marching_squares_demo.scad
@@ -0,0 +1,19 @@
+use <hull_polyline2d.scad>;
+use <util/rand.scad>;
+use <experimental/pnoise2.scad>;
+use <experimental/marching_squares.scad>;
+
+seed = rand(0, 256);
+points = [
+    for(y = [0:.2:10]) [
+        for(x = [0:.2:10]) [x, y, pnoise2(x, y, seed)]
+    ]
+];
+
+for(row = marching_squares(points, 0.1)) {
+    for(line = row) {
+        p0 = [line[0][0], line[0][1]];
+        p1 = [line[1][0], line[1][1]];
+        hull_polyline2d([p0, p1], width = .1);
+    }
+}    
\ No newline at end of file
diff --git a/src/experimental/marching_squares.scad b/src/experimental/marching_squares.scad
new file mode 100644
index 00000000..fe533f4c
--- /dev/null
+++ b/src/experimental/marching_squares.scad
@@ -0,0 +1,20 @@
+use <experimental/_impl/_marching_squares_impl.scad>;
+
+function marching_squares(pts, sigma) = 
+    let(labeled_pts = pn_label(pts, sigma))
+    [
+        for(y = [0:len(labeled_pts) - 2])
+            [
+                 for(x = [0:len(labeled_pts[0]) - 2])
+                 let(
+                    p0 = labeled_pts[y][x],
+                    p1 = labeled_pts[y + 1][x],
+                    p2 = labeled_pts[y + 1][x + 1],
+                    p3 = labeled_pts[y][x + 1],
+                    cell_pts = [p0, p1, p2, p3],
+                    contours_lt = contours_of(cell_pts, sigma)
+                 )
+                 if(contours_lt != [])
+                 each contours_lt
+           ]
+    ];
\ No newline at end of file