From 80fdbd755df698244f42f1a006909a9ebe5ae6b4 Mon Sep 17 00:00:00 2001 From: Nejc Jurkovic Date: Fri, 24 May 2024 22:26:02 +0200 Subject: [PATCH] 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 | 104 ++++++++---------- src/classy_blocks/construct/flat/quad.py | 22 ++-- .../construct/flat/sketches/mapped.py | 6 +- tests/test_construct/test_flat/test_map.py | 40 +++---- 4 files changed, 75 insertions(+), 97 deletions(-) diff --git a/src/classy_blocks/construct/flat/map.py b/src/classy_blocks/construct/flat/map.py index 0f1d951..da685bf 100644 --- a/src/classy_blocks/construct/flat/map.py +++ b/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) diff --git a/src/classy_blocks/construct/flat/quad.py b/src/classy_blocks/construct/flat/quad.py index 16f46e1..30a0d69 100644 --- a/src/classy_blocks/construct/flat/quad.py +++ b/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): diff --git a/src/classy_blocks/construct/flat/sketches/mapped.py b/src/classy_blocks/construct/flat/sketches/mapped.py index 66fee22..143e157 100644 --- a/src/classy_blocks/construct/flat/sketches/mapped.py +++ b/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() diff --git a/tests/test_construct/test_flat/test_map.py b/tests/test_construct/test_flat/test_map.py index 70e60ec..b39cec5 100644 --- a/tests/test_construct/test_flat/test_map.py +++ b/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)