re-implement tri_delaunay

src/experimental/_impl/_tri_delaunay_impl.scad

@@ -1,48 +1,229 @@
 use <experimental/tri_circumcircle.scad>;
+use <util/map/hashmap.scad>;
+use <util/map/hashmap_get.scad>;
+use <util/map/hashmap_keys.scad>;
+use <util/map/hashmap_entries.scad>;
+use <util/some.scad>;
+use <util/has.scad>;
+use <util/slice.scad>;
+use <util/find_index.scad>;
 
-// [max, min, min, max]      
+function cc_center(cc) = cc[0];
-function _tri_delaunay_boundIndices(points, leng, indices, i = 0) =
+function cc_rr(cc) = cc[2];
-    let(
-        bounds = i < leng - 1 ? _tri_delaunay_boundIndices(points, leng, indices, i + 1)
-                              : [undef, undef, undef, undef],
-        p0 = points[indices[0]],
-        p1 = points[indices[1]],
-        p2 = points[i],
-        c_circle = tri_circumcircle([p0, p1, p2]),
-        r = c_circle[1],
-        cside = c_circle == [] ? undef : cross(p1 - p0, c_circle[0] - p1),
-        pside = cross(p1 - p0, p2 - p1),
-        pdot = (p2 - p0) * (p2 - p1)
-    )
-    pdot < 0 && 0 == pside  ? [-1, 1, 1, -1] :
-    c_circle == [] ? bounds :
-    pside > 0 && cside > 0  ? [bounds[0], bounds[1], bounds[2], bounds[3] ? min(bounds[3], r) : r]  :
-    pside > 0 && cside <= 0 ? [bounds[0], bounds[1], bounds[2] ? max(bounds[2], r) : r, bounds[3]] :
-    pside <= 0 && cside > 0 ? [bounds[0], bounds[1] ? max(bounds[1],r) : r, bounds[2], bounds[3]] :
-                              [bounds[0] ? min(bounds[0], r) : r, bounds[1], bounds[2], bounds[3]]; // is_undef(cside) also returns this
 
-function _tri_delaunay_try_triangle(points, indices_boundIndices, i) = 
+function delaunay_init(points) =
     let(
-      indices = indices_boundIndices[0],
+		xs = [for(p = points) p[0]],
-      bounds = indices_boundIndices[1],
+		ys = [for(p = points) p[1]],
-      p0 = points[indices[0]],
+		max_x = max(xs),
-      p1 = points[indices[1]],
+		min_x = min(xs),
-      p2 = points[i],
+		max_y = max(ys),
-      c_circle = tri_circumcircle([p0, p1, p2]),
+		min_y = min(ys),
-      r = c_circle[1],
+		center = [max_x + min_x, max_y + min_y] / 2,
-      cside = c_circle == [] ? undef : cross(p1 - p0, c_circle[0] - p1)
+		halfW = abs(max_x - center[0]) * 2,
+		halfH = abs(max_y - center[1]) * 2,
+		coords = [
+		    center + [-halfW, -halfH],
+			center + [-halfW,  halfH],
+			center + [ halfW,  halfH],
+			center + [ halfW, -halfH],
+		],
+		t1 = [0, 1, 3], // indices
+		t2 = [2, 3, 1],
+		triangles = hashmap([
+		        [t1, [t2, undef, undef]],
+				[t2, [t1, undef, undef]]
+		    ]
+		),
+		circles = hashmap([
+		        [t1, tri_circumcircle([for(i = t1) coords[i]])],
+				[t2, tri_circumcircle([for(i = t2) coords[i]])]
+		    ]
 		)
-   c_circle == [] ? undef : 
+	)
-   (cside > 1  && (!bounds[2] && (!bounds[1] || r>=bounds[1]) && (!bounds[3] || r<=bounds[3]))) ||
+	[coords, triangles, circles];
-   (cside <= 1 && (!bounds[1] && (!bounds[2] || r>=bounds[2]) && (!bounds[0] || r<=bounds[0]))) ? concat(indices, [i]) :
-   undef;
 	
-function _tri_delaunay_triangleIndices(points, leng, indices_boundIndices, triangles, i = 0) =
+function delaunay_coords(d) = d[0];
+function delaunay_triangles(d) = d[1];
+function delaunay_circles(d) = d[2];
+
+function delaunay_addpoint(d, p) =
     let(
-      indices = indices_boundIndices[0],
+	    idx = len(delaunay_coords(d)),
-      newtriangle = _tri_delaunay_try_triangle(points, indices_boundIndices, i),
+		ndelaunay = delaunayAddCoords(d, p),
-      newtriangles = newtriangle ? concat(triangles, [newtriangle]) : triangles
+		badTriangles = delaunayBadTriangles(ndelaunay, p),
+		boundaries = delaunayBoundaries(ndelaunay, badTriangles),
+		ndelaunay2 = delBadTriangles(ndelaunay, badTriangles),
+		newTriangles = [
+		    for(b = boundaries) [
+			    [idx, b[0][0], b[0][1]], // t
+				b[0],                    // edge
+				b[1]                     // delaunayTri
+			]
+		]
 	)
-   i <= indices[1] ? _tri_delaunay_triangleIndices(points, leng, indices_boundIndices, [], indices[1] + 1) : 
+	adjustNeighbors(ndelaunay2, newTriangles);
-   i< leng - 1 ? _tri_delaunay_triangleIndices(points, leng, indices_boundIndices, newtriangles,i + 1) :
+
-   newtriangles;
+function adjustNeighbors(d, newTriangles) = 
+    let(
+	    coords = delaunay_coords(d),
+		nts = [
+			for(nt = newTriangles)
+			[nt[0], [nt[2], undef, undef]]
+		],
+	    ncs = [
+		    for(nt = newTriangles)
+			[nt[0], tri_circumcircle([for(i = nt[0]) coords[i]])]
+		],
+		nd = [
+		    coords, 
+			hashmap(concat(hashmap_entries(delaunay_triangles(d)), nts)),
+			hashmap(concat(hashmap_entries(delaunay_circles(d)), ncs))
+		],
+		leng = len(newTriangles),
+		aDtrid = _adjustNeighborsDtri(nd, newTriangles, leng)
+	)
+	// aDtrid;
+	_adjustNeighborsOtri(aDtrid, newTriangles, leng);
+
+function _adjustNeighborsOtri(d, newTriangles, leng, i = 0) = 
+    i == leng ? d :
+	let(
+	    t = newTriangles[i][0],
+		nbr1 = newTriangles[(i + 1) % leng][0],
+		nbr2 = newTriangles[i > 0 ? i - 1 : leng + i - 1][0],
+		triangles = delaunay_triangles(d),
+        entries = hashmap_entries(triangles),
+		nTriangles = hashmap([for(entry = entries) entry[0] == t ? [t, [entry[1][0], nbr1, nbr2]] : entry]),
+		nd = [delaunay_coords(d), nTriangles, delaunay_circles(d)]
+	)
+	_adjustNeighborsOtri(nd, newTriangles, leng, i + 1);
+ 
+
+function _adjustNeighborsDtri(d, newTriangles, leng, i = 0) =
+    i == leng ? d :
+    let(
+	    t = newTriangles[i][0],
+		edge = newTriangles[i][1],
+		delaunayTri = newTriangles[i][2]
+	)
+	delaunayTri != undef ? 
+	let(
+	    neighbors = hashmap_get(delaunay_triangles(d), delaunayTri),
+		nbri = find_index(neighbors, function(nbr) nbr != undef && has(nbr, edge[1]) && has(nbr, edge[0])),
+		nd = nbri == -1 ? d : updateNbrs(d, delaunayTri, concat(slice(neighbors, 0, nbri), [t], slice(neighbors, nbri + 1)))
+	)
+    _adjustNeighborsDtri(nd, newTriangles, leng, i + 1) :
+	_adjustNeighborsDtri(d, newTriangles, leng, i + 1);
+	
+function updateNbrs(d, delaunayTri, neighbors) =
+    let(
+	    coords = delaunay_coords(d),
+	    triangles = delaunay_triangles(d),
+		circles = delaunay_circles(d),
+		tri_entries = hashmap_entries(triangles),
+		nTriangles = hashmap([
+		    for(entry = tri_entries)
+			entry[0] == delaunayTri ? [delaunayTri, neighbors] : entry
+		])
+	)
+	[coords, nTriangles, circles];
+	
+function delaunayAddCoords(d, p) = 
+    [
+	    concat(delaunay_coords(d), [p]), 
+		delaunay_triangles(d), 
+		delaunay_circles(d)
+	];
+
+function delaunayBadTriangles(d, p) = 
+    let(
+	     triangles = delaunay_triangles(d),
+		 circles = delaunay_circles(d)
+	) 
+	[
+        for(t = hashmap_keys(triangles))
+	        if(inCircumcircle(t, p, circles))
+	            t
+    ];
+
+/* 
+    is p in t?
+    t: triangle indices
+	circles: a hashmap
+*/
+function inCircumcircle(t, p, circles) = 
+    let(
+	    c = hashmap_get(circles, t),
+		v = cc_center(c) - p,
+		rr = v[0] ^ 2 + v[1] ^ 2
+	)
+	rr <= cc_rr(c);
+
+function delaunayBoundaries(d, badTriangles) = 
+    let(
+	    boundaries = [],
+		t = badTriangles[0],
+		vi = 0
+	)
+	_delaunayBoundaries(d, badTriangles, boundaries, t, vi);
+
+function _delaunayBoundaries(d, badTriangles, boundaries, t, vi) = 
+    let(
+	    triangles = delaunay_triangles(d),
+	    opTri = hashmap_get(triangles, t)[vi]
+	)
+	some(badTriangles, function(tri) tri == opTri) ?
+		let(
+		    i = find_index(hashmap_get(triangles, opTri), function(tri) tri == t),
+			nvi = (i + 1) % 3,
+			nt = opTri
+		)
+		_delaunayBoundaries(d, badTriangles, boundaries, nt, nvi) : 
+	    let(
+		    nboundaries = concat(boundaries, [[
+			    [t[(vi + 1) % 3], t[vi > 0 ? vi - 1 : (vi + 2)]], // edge
+				opTri                                             // delaunayTri
+			]]),
+			nvi = (vi + 1) % 3,
+			v1 = nboundaries[0][0][0],
+			v2 = nboundaries[len(nboundaries) - 1][0][1]
+		)
+		v1 == v2 ? nboundaries : _delaunayBoundaries(d, badTriangles, nboundaries, t, nvi);
+
+function delBadTriangles(d, badTriangles) = 
+    let(
+	    triangles = delaunay_triangles(d),
+		circles = delaunay_circles(d),
+		nTriangles = hashmap([
+		     for(t = hashmap_keys(triangles))
+			 if(!has(badTriangles, t))
+			 [t, hashmap_get(triangles, t)]
+		]),
+		nCircles = hashmap([
+		     for(t = hashmap_keys(circles))
+			 if(!has(badTriangles, t))
+			 [t, hashmap_get(circles, t)]
+		])
+	)
+	[delaunay_coords(d), nTriangles, nCircles];
+
+
+	
+function _tri_delaunay(d, points, leng, i = 0) =
+    i == leng ? d :
+	_tri_delaunay(delaunay_addpoint(d, points[i]), points, leng, i + 1);
+
+function tri_delaunay_shapes(d) = 
+    let(coords = delaunay_coords(d))
+	[
+		for(tri = hashmap_keys(delaunay_triangles(d))) 
+		if(tri[0] > 3 && tri[1] > 3 && tri[2] > 3)
+		[coords[tri[0]], coords[tri[1]], coords[tri[2]]]
+	];
+
+function tri_delaunay_indices(d) =	[
+	for(tri = hashmap_keys(delaunay_triangles(d))) 
+	if(tri[0] > 3 && tri[1] > 3 && tri[2] > 3)
+	[tri[0] - 4, tri[1] - 4, tri[2] - 4]
+];

src/experimental/demo/tri_delaunay_demo.scad

@@ -1,16 +1,19 @@
-use <hull_polyline2d.scad>;
 use <experimental/tri_delaunay.scad>;
+use <hull_polyline2d.scad>;
 
-points=[for(i = [0:50]) rands(-20, 20, 2)];
+points = [for(i = [0:20]) rands(-100, 100, 2)]; 
 
-for(tri = tri_delaunay(points)) {
+drawTris(tri_delaunay(points));
-    hull_polyline2d(
+module drawTris(pointsOfTriangles) {
-        [points[tri[0]], points[tri[1]], points[tri[2]], points[tri[0]]], width = .2
+	#for(t = pointsOfTriangles) {
-    );
+	    hull_polyline2d(concat(t, [t[0]]));
+	}	
 }
 
-color("red")
+drawTris2(points, tri_delaunay(points, ret = "INDICES"));
-for(point = points) {
+module drawTris2(points, indices) {
-   translate(point) 
+    pointsOfTriangles = [for(i = indices) [points[i[0]], points[i[1]], points[i[2]]]];
-       circle(.5);
+	%for(t = pointsOfTriangles) {
+	    hull_polyline2d(concat(t, [t[0]]), 2);
+	}	
 }

src/experimental/tri_delaunay.scad

@@ -1,18 +1,8 @@
-use <experimental/_impl/_tri_delaunay_impl.scad>;
+use<experimental/_impl/_tri_delaunay_impl.scad>;
-use <util/flat.scad>;
 
-function tri_delaunay(points)=
+// ret: "SHAPES", "INDICES", "DELAUNAY"
-   let(
+function tri_delaunay(points, ret = "SHAPES") = 
-      leng = len(points),
+    let(d = _tri_delaunay(delaunay_init(points), points, len(points)))
-      indices_lt = [for(i=[0:leng - 3]) for(j = [i + 1:leng - 2]) [i, j]],
+    ret == "SHAPES" ?  tri_delaunay_shapes(d) : 
-      indices_boundIndices_lt = [
+    ret == "INDICES" ? tri_delaunay_indices(d) :
-          for(indices = indices_lt) 
+    d; // [coords(list), triangles(hashmap), circles(hashmap)]
-              [indices, _tri_delaunay_boundIndices(points, leng, indices)]
-      ]
-   )
-   flat(
-       [
-           for(indices_boundIndices = indices_boundIndices_lt) 
-               _tri_delaunay_triangleIndices(points, leng, indices_boundIndices, [])
-       ]
-   );