Drop sketch optimization, add a custom example

Sketch optimization by 'minimizing quad energy' does not work in
real-life as it is too unstable numerically and slow.

RoundSquare sketch works because it is convex and smoothing is
sufficient.

Also smoothing is now done by averaging distances to neighbour corners
AND quad centers. Not sure if it helps.

Sketches with concave corners will not profit from smoothing.
Now to solve that; options are numerous, one more difficult than the
other.

src/classy_blocks/construct/flat/map.py

@@ -1,7 +1,6 @@
-from typing import Dict, List, Set
+from typing import Dict, List, Optional, Set
 
 import numpy as np
-import scipy.optimize
 
 from classy_blocks.construct.flat.quad import Quad
 from classy_blocks.types import NPPointListType, QuadIndexType
@@ -11,13 +10,25 @@
 class QuadMap:
     def __init__(self, positions: NPPointListType, indexes: List[QuadIndexType]):
         self.indexes = indexes
-        self.positions = positions
 
-        self.quads = [Quad(self.positions, quad_indexes) for quad_indexes in indexes]
+        self.quads = [Quad(positions, quad_indexes) for quad_indexes in indexes]
+
+    def update(self, positions: Optional[NPPointListType] = None) -> None:
+        if positions is None:
+            positions = self.positions
 
-    def update(self) -> None:
         for quad in self.quads:
-            quad.update()
+            quad.update(positions)
+
+    @property
+    def positions(self) -> NPPointListType:
+        """Reconstructs positions back from faces so they are always up-to-date,
+        even after transforms"""
+        indexes = list(np.array(self.indexes).flatten())
+        max_index = max(indexes)
+        all_points = f.flatten_2d_list([quad.face.point_array.tolist() for quad in self.quads])
+
+        return np.array([all_points[indexes.index(i)] for i in range(max_index + 1)])
 
     @property
     def connections(self) -> List[Set[int]]:
@@ -39,7 +50,7 @@ def neighbours(self) -> Dict[int, Set[int]]:
         return neighbours
 
     @property
-    def fixed_points(self) -> Set[int]:
+    def boundary_points(self) -> Set[int]:
         """Returns indexes of points that can be smoothed"""
         connections = self.connections
         fixed_points: Set[int] = set()
@@ -50,65 +61,42 @@ def fixed_points(self) -> Set[int]:
 
         return fixed_points
 
-    def smooth_laplacian(self, iterations: int = 5) -> None:
-        """Smooth the points using laplacian smoothing;
-        each point is moved to the average of its neighbours"""
-        neighbours = self.neighbours
-        fixed_points = self.fixed_points
-
-        for _ in range(iterations):
-            for point_index, point_neighbours in neighbours.items():
-                if point_index in fixed_points:
-                    continue
-
-                nei_positions = [self.positions[i] for i in point_neighbours]
-
-                self.positions[point_index] = np.average(nei_positions, axis=0)
-
-        self.update()
-
-    def get_nearby_quads(self) -> Dict[int, List[Quad]]:
+    def get_nearby_quads(self, index: int) -> Set[int]:
         """Returns a list of quads that contain each movable point"""
-        corner_points: Dict[int, List[Quad]] = {}
-
-        fixed_points = self.fixed_points
-
-        for i, point in enumerate(self.positions):
-            if i in fixed_points:
-                continue
-
-            corner_points[i] = []
+        indexes = set()
 
-            for quad in self.quads:
-                if quad.contains(point):
-                    corner_points[i].append(quad)
+        for i, quad in enumerate(self.quads):
+            if index in quad.indexes:
+                indexes.add(i)
 
-        return corner_points
+        return indexes
 
-    def optimize_energy(self):
-        """Replace quad edges by springs and moves their positions
-        so that all springs are in the most relaxed state possible,
-        minimizing the energy of the system"""
-        # get quads that are defined by each point
-        fixed_points = self.fixed_points
+    def smooth_laplacian(
+        self,
+        fix_points: Optional[Set[int]] = None,
+    ) -> None:
+        """Smooth the points using laplacian smoothing;
+        each point is moved to the average of its neighbours"""
+        if fix_points is None:
+            fix_points = set()
 
-        e1 = self.quads[0].e1
-        e2 = self.quads[0].e2
+        fixed_points = self.boundary_points.union(fix_points)
+        positions = self.positions
 
-        def energy(i, x) -> float:
-            e = 0
+        for point_index, point_neighbours in self.neighbours.items():
+            if point_index in fixed_points:
+                continue
 
-            self.positions[i] = x[0] * e1 + x[1] * e2
+            nei_positions = [positions[i] for i in point_neighbours]
+            corner_center = np.average(nei_positions, axis=0)
 
-            for quad in self.quads:
-                quad.update()
-                e += quad.energy
+            nei_quads = [self.quads[i] for i in self.get_nearby_quads(point_index)]
+            quad_center = np.average([quad.center for quad in nei_quads], axis=0)
 
-            return e
+            # distances = np.array([f.norm(pos - center) for pos in nei_positions])
+            # ratios = distances / np.average(distances)
+            # weights = ratios
 
-        for i in range(len(self.positions)):
-            if i in fixed_points:
-                continue
+            positions[point_index] = (corner_center + quad_center) / 2
 
-            initial = [1, 1]
-            scipy.optimize.minimize(lambda x, i=i: energy(i, x), initial)
+        self.update(positions)

src/classy_blocks/construct/flat/quad.py

@@ -13,12 +13,11 @@ class Quad:
 
     def __init__(self, positions: NPPointListType, indexes: QuadIndexType):
         self.indexes = indexes
-        self.positions = positions
-        self.face = Face([self.positions[i] for i in self.indexes])
+        self.face = Face([positions[i] for i in self.indexes])
 
-    def update(self) -> None:
+    def update(self, positions: NPPointListType) -> None:
         """Update Face position"""
-        self.face.update([self.positions[i] for i in self.indexes])
+        self.face.update([positions[i] for i in self.indexes])
 
     def contains(self, point: NPPointType) -> bool:
         """Returns True if the given point is a part of this quad"""
@@ -45,18 +44,17 @@ def center(self):
         return np.average(self.points, axis=0)
 
     @property
-    def energy(self):
-        e = 0
-
-        ideal_side = self.perimeter / 4
-        ideal_diagonal = (ideal_side / 2) * 2**0.5
+    def area(self):
         center = self.center
+        sum_area = 0
 
         for i in range(4):
-            e += (f.norm(self.points[i] - self.points[(i + 1) % 4]) - ideal_side) ** 2
-            e += (f.norm(center - self.points[i]) - ideal_diagonal) ** 2
+            side_1 = self.points[i] - center
+            side_2 = self.points[(i + 1) % 4] - center
+
+            sum_area += 0.5 * f.norm(np.cross(side_1, side_2))
 
-        return e / 8
+        return sum_area
 
     @property
     def e1(self):

src/classy_blocks/construct/flat/sketches/mapped.py

@@ -37,9 +37,5 @@ def center(self) -> NPPointType:
 
     def smooth(self, n_iter: int = 5) -> None:
         """Smooth the internal points using laplacian smoothing"""
-        self.quad_map.smooth_laplacian(n_iter)
-
-    def optimize(self, n_iter: int = 3) -> None:
-        """Optimize internal vertices using the spring analogy"""
         for _ in range(n_iter):
-            self.quad_map.optimize_energy()
+            self.quad_map.smooth_laplacian()

tests/test_construct/test_flat/test_map.py

@@ -20,7 +20,7 @@ def test_quad_update(self):
         new_position = [-1, -1, 0]
 
         self.positions[0] = new_position
-        quad.update()
+        quad.update(self.positions)
 
         np.testing.assert_equal(quad.points, self.positions)
 
@@ -38,21 +38,13 @@ def test_perimeter_rectangle(self):
         self.positions[1] += delta
         self.positions[2] += delta
 
-        self.quad.update()
+        self.quad.update(self.positions)
 
         self.assertEqual(self.quad.perimeter, 6)
 
     def test_center(self):
         np.testing.assert_equal(self.quad.center, [0.5, 0.5, 0])
 
-    def test_energy_zero(self):
-        self.assertAlmostEqual(self.quad.energy, 0)
-
-    def test_energy_nonzero(self):
-        self.positions[0] += f.vector(-1, 0, 0)
-
-        self.assertGreater(self.quad.energy, 0)
-
     def test_e1(self):
         np.testing.assert_almost_equal(self.quad.e1, [1, 0, 0])
 
@@ -93,6 +85,9 @@ def quads(self):
     def quad_map(self):
         return QuadMap(self.positions, self.quads)
 
+    def test_positions(self):
+        np.testing.assert_equal(self.quad_map.positions, self.positions)
+
     def test_find_neighbours(self):
         # A random blocking (quadding)
         positions = np.zeros((9, 3))
@@ -128,25 +123,26 @@ def test_fixed_points(self):
 
         quad_map = QuadMap(positions, indexes)
 
-        fixed_points = quad_map.fixed_points
+        fixed_points = quad_map.boundary_points
 
         self.assertSetEqual(
             fixed_points,
             {4, 5, 6, 7},
         )
 
-    def test_smooth(self):
-        # a grid of vertices 3x3
-        quad_map = QuadMap(self.positions, self.quads)
-        quad_map.smooth_laplacian(10)
-
-        np.testing.assert_almost_equal(quad_map.positions[4], [1, 1, 0])
+    def test_update(self):
+        positions = self.positions
+        positions[0] = [0.5, 0.5, 0]
 
-    def test_optimize(self):
-        quad_map = QuadMap(self.positions, self.quads)
+        quad_map = self.quad_map
+        quad_map.update(positions)
 
-        self.assertGreater(quad_map.quads[0].energy, 0)
+        np.testing.assert_equal(quad_map.quads[0].points[0], [0.5, 0.5, 0])
 
-        quad_map.optimize_energy()
+    def test_smooth(self):
+        # a grid of vertices 3x3
+        quad_map = QuadMap(self.positions, self.quads)
+        for _ in range(10):
+            quad_map.smooth_laplacian()
 
-        self.assertAlmostEqual(quad_map.quads[0].energy, 0)
+        np.testing.assert_almost_equal(quad_map.positions[4], [1, 1, 0], decimal=5)