Introduce a new arcwelder

kikit/substrate.py

@@ -11,7 +11,7 @@
 from enum import IntEnum
 from itertools import product
 
-from typing import List, Tuple, Union
+from typing import Iterable, List, Tuple, Union
 
 from kikit.common import *
 from kikit.units import deg
@@ -337,19 +337,160 @@ def substratesFrom(polygons):
         substrates.append(Polygon(polygon.exterior, holes))
     return substrates
 
-def commonCircle(a, b, c):
-    """
-    Given three 2D points return (x, y, r) of the circle they lie on or None if
-    they lie in a line
-    """
-    arc = pcbnew.PCB_SHAPE()
-    arc.SetShape(STROKE_T.S_ARC)
-    arc.SetArcGeometry(toKiCADPoint(a), toKiCADPoint(b), toKiCADPoint(c))
-    center = [int(x) for x in arc.GetCenter()]
-    mid = [(a[i] + b[i]) // 2 for i in range(2)]
-    if center == mid:
-        return None
-    return roundPoint(center, -8)
+class CircleFitCandidates:
+    def __init__(self, tolerance: int = fromMm(0.01)):
+        self.tolerance = tolerance
+
+        self._xs = []
+        self._xSum: float = 0
+        self._ys = []
+        self._ySum: float = 0
+
+        self.foundCircle: Optional[Tuple[np.array, float]] = None
+
+    def __len__(self) -> int:
+        return len(self._xs)
+
+    def addPoint(self, point: VECTOR2I) -> bool:
+        """
+        Try adding a candidate point. Returns true, if the point could lie on
+        already found a circle. If the point doesn't lie on an already found
+        circle, it is not added to the collection.
+        """
+        self._xs.append(point[0])
+        self._xSum += point[0]
+
+        self._ys.append(point[1])
+        self._ySum += point[1]
+
+        if len(self) < 5:
+            return True
+
+        newCenter, newRadius = self._fitCircle()
+
+        needsRevalidation = False
+        if self.foundCircle is None:
+            needsRevalidation = True
+        else:
+            c, r = self.foundCircle
+            needsRevalidation = np.linalg.norm(newCenter - c) > self.tolerance or np.abs(newRadius - r) > self.tolerance
+
+        if needsRevalidation:
+            points = [np.array([x, y]) for x, y in zip(self._xs, self._ys)]
+            toRevalidate = list(zip(points, points[1:]))
+        else:
+            toRevalidate = [(np.array([self._xs[-2], self._ys[-2]]), np.array([self._xs[-1], self._ys[-1]]))]
+
+        if self._doLinesFitCircle(toRevalidate, newCenter, newRadius):
+            self.foundCircle = newCenter, newRadius
+            return True
+
+        self._popLast()
+        return False
+
+    def _doLinesFitCircle(self, lines: Iterable[Tuple[np.array, np.array]],
+                          c: np.array, r: float) -> bool:
+        for start, end in lines:
+            # We don't use list of candidates in this as it is in the hot path
+            # and constructing the list of candidates adds a significant
+            # overhead (both technically, and for early return)
+            #
+            # The extreme occurs either in one of the endpoints or in the
+            # projection of center of the circle to the line (if it lies on the
+            # segment).
+            if np.abs(np.linalg.norm(start - c) - r) > self.tolerance:
+                return False
+            if np.abs(np.linalg.norm(end - c) - r) > self.tolerance:
+                return False
+
+            # Project center to the line; if it doesn't fit line, continue
+            ap = c - start
+            ab = end - start
+            t = np.dot(ap, ab) / np.dot(ab, ab)
+            if t < 0 or t > 1:
+                continue
+            projection = start + t * ab
+            if np.abs(np.linalg.norm(projection - c) - r) > self.tolerance:
+                return False
+        return True
+
+    @property
+    def start(self) -> np.array:
+        return np.array((self._xs[0], self._ys[0]))
+
+    @property
+    def end(self) -> np.array:
+        return np.array((self._xs[-1], self._ys[-1]))
+
+    @property
+    def mid(self) -> np.array:
+        idx = len(self) // 2
+        return np.array((self._xs[idx], self._ys[idx]))
+
+    def _popLast(self):
+        self._xSum -= self._xs[-1]
+        self._xs.pop()
+        self._ySum -= self._ys[-1]
+        self._ys.pop()
+
+    def _fitCircle(self, maxIter = 10) -> Tuple[np.array, float]:
+        """
+        Implements Kenichi Kanatani, Prasanna Rangarajan, "Hyper least squares fitting of circles and ellipses"
+        Computational Statistics & Data Analysis, Vol. 55, pages 2197-2208, (2011)
+
+        Implementation is based on https://github.com/AlliedToasters/circle-fit/blob/master/src/circle_fit/circle_fit.py
+
+        Returns center and radius of the circle fit.
+        """
+        n = len(self._xs)
+
+        xMean = self._xSum / n
+        yMean = self._ySum / n
+
+        Xi = np.array(self._xs) - xMean
+        Yi = np.array(self._ys) - yMean
+        Zi = Xi * Xi + Yi * Yi
+
+        # compute moments
+        Mxy = (Xi * Yi).sum() / n
+        Mxx = (Xi * Xi).sum() / n
+        Myy = (Yi * Yi).sum() / n
+        Mxz = (Xi * Zi).sum() / n
+        Myz = (Yi * Zi).sum() / n
+        Mzz = (Zi * Zi).sum() / n
+
+        # computing the coefficients of characteristic polynomial
+        Mz = Mxx + Myy
+        Cov_xy = Mxx * Myy - Mxy * Mxy
+        Var_z = Mzz - Mz * Mz
+
+        A2 = 4 * Cov_xy - 3 * Mz * Mz - Mzz
+        A1 = Var_z * Mz + 4. * Cov_xy * Mz - Mxz * Mxz - Myz * Myz
+        A0 = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy
+        A22 = A2 + A2
+
+        # finding the root of the characteristic polynomial
+        Y = A0
+        X = 0.0
+        for i in range(maxIter):
+            Dy = A1 + X * (A22 + 16. * (X ** 2))
+            xnew = X - Y / Dy
+            if xnew == X or not np.isfinite(xnew):
+                break
+            ynew = A0 + xnew * (A1 + xnew * (A2 + 4. * xnew * xnew))
+            if abs(ynew) >= abs(Y):
+                break
+            X, Y = xnew, ynew
+
+        det = X ** 2 - X * Mz + Cov_xy
+        Xcenter = (Mxz * (Myy - X) - Myz * Mxy) / det / 2.
+        Ycenter = (Myz * (Mxx - X) - Mxz * Mxy) / det / 2.
+
+        xc: float = Xcenter + xMean
+        yc: float = Ycenter + yMean
+        r = np.sqrt(abs(Xcenter ** 2 + Ycenter ** 2 + Mz))
+        return np.array((xc, yc)), r
+
 
 
 def liesOnSegment(start, end, point, tolerance=fromMm(0.01)):
@@ -529,32 +670,73 @@ def serialize(self, reconstructArcs=False):
         return items
 
     def _serializeRing(self, ring, reconstructArcs):
+        TOLERANCE = fromMm(0.01)
         coords = ring.coords
-        segments = []
         if coords[0] != coords[-1]:
             raise RuntimeError("Ring is incomplete")
+
+        # Always start with the longes semgent (so we do not start in the middle
+        # of an arc if possible)
+        coords = np.array(coords[:-1])  # Exclude the last point since it's a duplicate of the first
+        distances = np.sqrt(np.sum(np.diff(coords, axis=0, append=coords[:1,:])**2, axis=1))
+        max_dist_index = np.argmax(distances)
+        rearranged = np.roll(coords, -max_dist_index, axis=0)
+        coords = np.vstack((rearranged, rearranged[0]))
+
+        segments = []
         i = 0
         while i < len(coords):
             j = i # in the case the following cycle never happens
+            candidateCircle = CircleFitCandidates(tolerance=TOLERANCE)
             if reconstructArcs:
-                for j in range(i, len(coords) - 3): # Just walk until there is an arc
-                    cc1 = commonCircle(coords[j], coords[j + 1], coords[j + 2])
-                    cc2 = commonCircle(coords[j + 1], coords[j + 2], coords[j + 3])
-                    if cc1 is None or cc2 is None or cc1 != cc2:
+                for j in range(i, len(coords)):
+                    # Just walk edge segments until there is an arc
+                    if not candidateCircle.addPoint(coords[j]):
                         break
-            if j - i > 10:
-                j += 1
-                # Yield a circle
-                a = coords[i]
-                b = coords[(i + j) // 2]
-                c = coords[j]
-                segments.append(self._constructArc(a, b, c))
-                i = j
+
+            if candidateCircle.foundCircle is not None:
+                center, radius = candidateCircle.foundCircle
+                start, end, mid = candidateCircle.start, candidateCircle.end, candidateCircle.mid
+
+                # We prefer to preserve arc start and end points, adjust center
+                # so it is true:
+                dir = end - start
+                chordLength = np.linalg.norm(dir)
+                if chordLength == 0:
+                    segments.append(self._constructCircle(center, radius))
+                else:
+                    dir /= chordLength
+                    normal = np.array([dir[1], -dir[0]])
+                    chordMidpoint = (start + end) / 2
+
+                    # Due to numerical errors, the height might be invalid, if
+                    # so, assume it is zero
+                    heightSquared = radius ** 2 - chordLength ** 2 / 4
+                    if heightSquared < 0:
+                        height = 0
+                    else:
+                        height = np.sqrt(heightSquared)
+
+                    centerCandidates = [chordMidpoint + normal * height, chordMidpoint - normal * height]
+                    center = min(centerCandidates, key=lambda c: np.linalg.norm(c - center))
+
+                    centerToMidpoint = chordMidpoint - center
+                    centerToMidpointLength = np.linalg.norm(centerToMidpoint)
+                    if centerToMidpointLength == 0:
+                        centerToMidpoint = normal
+                    else:
+                        centerToMidpoint /= np.linalg.norm(centerToMidpoint)
+                    middleCandidates = [center + centerToMidpoint * radius, center - centerToMidpoint * radius]
+                    arcMiddle = min(middleCandidates, key=lambda c: np.linalg.norm(c - mid))
+
+                    segments.append(self._constructArc(toKiCADPoint(start), toKiCADPoint(arcMiddle), toKiCADPoint(end)))
+
+                i += len(candidateCircle) - 1
             else:
                 # Yield a line
                 a = coords[i]
                 b = coords[(i + 1) % len(coords)]
-                if a != b:
+                if np.linalg.norm(np.array(a) - np.array(b)) > SHP_EPSILON:
                     segments.append(self._constructEdgeSegment(a, b))
                 i += 1
         return segments
@@ -577,6 +759,14 @@ def _constructArc(self, a, b, c):
         arc.SetArcGeometry(toKiCADPoint(a), toKiCADPoint(b), toKiCADPoint(c))
         return arc
 
+    def _constructCircle(self, c, r):
+        circle = pcbnew.PCB_SHAPE()
+        circle.SetShape(STROKE_T.S_CIRCLE)
+        circle.SetLayer(Layer.Edge_Cuts)
+        circle.SetCenter(toKiCADPoint(c))
+        circle.SetRadius(int(r))
+        return circle
+
     def boundingBox(self):
         """
         Return bounding box as BOX2I
@@ -674,6 +864,7 @@ def tab(self, origin, direction, width, partitionLine=None,
                         # numerical instability on small slopes that yields
                         # artifacts on substrate union
                         offsetTabFace = [(p[0] - SHP_EPSILON * direction[0], p[1] - SHP_EPSILON * direction[1]) for p in tabFace.coords]
+                        partitionFaceCoord = [(p[0] + SHP_EPSILON * direction[0], p[1] + SHP_EPSILON * direction[1]) for p in partitionFaceCoord]
                         tab = Polygon(offsetTabFace + partitionFaceCoord)
                         return self._makeTabFillet(tab, tabFace, fillet)
                 return None, None
@@ -739,7 +930,7 @@ def millFillets(self, millRadius):
         Add fillets to inner corners which will be produced by a mill with
         given radius.
         """
-        EPS = fromMm(0.01)
+        EPS = 0.01 # This number is intentionally below KiCAD's resolution of 1nm to not enclose narrow slots, but to preserve radius
         RES = 32
         if millRadius < EPS:
             return

setup.py

@@ -59,7 +59,7 @@
     install_requires=[
         "numpy", # Required for MacOS
         "pcbnewTransition >= 0.4, <=0.5",
-        "shapely>=1.7",
+        "shapely>=2.0.3",
         "click>=7.1",
         "markdown2>=2.4",
         "pybars3>=0.9",