Merge works.

gerbonara/gerber/gerber_statements.py

@@ -21,6 +21,14 @@ Gerber (RS-274X) Statements
 
 """
 
+def convert(value, src, dst):
+        if src == dst or src is None or dst is None or value is None:
+            return value
+        elif dst == 'mm':
+            return value * 25.4
+        else:
+            return value / 25.4
+
 class Statement:
     pass
 
@@ -88,6 +96,9 @@ class ApertureDefStmt(ParamStmt):
     def __str__(self):
         return f'<AD aperture def for {str(self.aperture).strip("<>")}>'
 
+    def __repr__(self):
+        return f'ApertureDefStmt({self.number}, {repr(self.aperture)})'
+
 
 class ApertureMacroStmt(ParamStmt):
     """ AM - Aperture Macro Statement """
@@ -117,13 +128,14 @@ class ImagePolarityStmt(ParamStmt):
 class CoordStmt(Statement):
     """ D01 - D03 operation statements """
 
-    def __init__(self, x, y, i=None, j=None):
+    def __init__(self, x, y, i=None, j=None, unit=None):
         self.x, self.y, self.i, self.j = x, y, i, j
+        self.unit = unit
 
     def to_gerber(self, settings=None):
         ret = ''
         for var in 'xyij':
-            val = getattr(self, var)
+            val = convert(getattr(self, var), self.unit, settings.unit)
             if val is not None:
                 ret += var.upper() + settings.write_gerber_value(val)
         return ret + self.code + '*'

gerbonara/gerber/graphic_objects.py

@@ -9,6 +9,7 @@ from .gerber_statements import *
 class GerberObject:
     _ : KW_ONLY
     polarity_dark : bool = True
+    unit : str = None
 
     def to_primitives(self):
         raise NotImplementedError()
@@ -31,12 +32,12 @@ class Flash(GerberObject):
     def to_statements(self, gs):
         yield from gs.set_polarity(self.polarity_dark)
         yield from gs.set_aperture(self.aperture)
-        yield FlashStmt(self.x, self.y)
-        gs.update_point(self.x, self.y)
+        yield FlashStmt(self.x, self.y, unit=self.unit)
+        gs.update_point(self.x, self.y, unit=self.unit)
 
 class Region(GerberObject):
-    def __init__(self, outline=None, arc_centers=None, *, polarity_dark):
-        super().__init__(polarity_dark=polarity_dark)
+    def __init__(self, outline=None, arc_centers=None, *, unit, polarity_dark):
+        super().__init__(unit=unit, polarity_dark=polarity_dark)
         outline = [] if outline is None else outline
         arc_centers = [] if arc_centers is None else arc_centers
         self.poly = gp.ArcPoly(outline, arc_centers)
@@ -50,7 +51,8 @@ class Region(GerberObject):
     def with_offset(self, dx, dy):
         return Region([ (x+dx, y+dy) for x, y in self.poly.outline ],
                 self.poly.arc_centers,
-                polarity_dark=self.polarity_dark)
+                polarity_dark=self.polarity_dark,
+                unit=self.unit)
 
     def rotate(self, angle, cx=0, cy=0):
         self.poly.outline = [ gp.rotate_point(x, y, angle, cx, cy) for x, y in self.poly.outline ]
@@ -76,18 +78,20 @@ class Region(GerberObject):
         yield from gs.set_polarity(self.polarity_dark)
         yield RegionStartStmt()
 
-        yield from gs.set_current_point(self.poly.outline[0])
+        yield from gs.set_current_point(self.poly.outline[0], unit=self.unit)
 
         for point, arc_center in zip(self.poly.outline[1:], self.poly.arc_centers):
             if arc_center is None:
                 yield from gs.set_interpolation_mode(LinearModeStmt)
-                yield InterpolateStmt(*point)
+                yield InterpolateStmt(*point, unit=self.unit)
+                gs.update_point(*point, unit=self.unit)
 
             else:
                 cx, cy = arc_center
                 x2, y2 = point
                 yield from gs.set_interpolation_mode(CircularCCWModeStmt)
-                yield InterpolateStmt(x2, y2, cx-x2, cy-y2)
+                yield InterpolateStmt(x2, y2, cx-x2, cy-y2, unit=self.unit)
+                gs.update_point(x2, y2, unit=self.unit)
 
         yield RegionEndStmt()
 
@@ -123,9 +127,9 @@ class Line(GerberObject):
         yield from gs.set_polarity(self.polarity_dark)
         yield from gs.set_aperture(self.aperture)
         yield from gs.set_interpolation_mode(LinearModeStmt)
-        yield from gs.set_current_point(self.p1)
-        yield InterpolateStmt(*self.p2)
-        gs.update_point(*self.p2)
+        yield from gs.set_current_point(self.p1, unit=self.unit)
+        yield InterpolateStmt(*self.p2, unit=self.unit)
+        gs.update_point(*self.p2, unit=self.unit)
 
 
 @dataclass
@@ -214,7 +218,8 @@ class Arc(GerberObject):
         yield from gs.set_polarity(self.polarity_dark)
         yield from gs.set_aperture(self.aperture)
         yield from gs.set_interpolation_mode(CircularCCWModeStmt)
-        yield from gs.set_current_point(self.p1)
-        yield InterpolateStmt(self.x2, self.y2, self.cx, self.cy)
+        yield from gs.set_current_point(self.p1, unit=self.unit)
+        yield InterpolateStmt(self.x2, self.y2, self.cx, self.cy, unit=self.unit)
+        gs.update_point(*self.p2, unit=self.unit)
 
 

gerbonara/gerber/rs274x.py

@@ -41,6 +41,14 @@ from . import graphic_objects as go
 from . import apertures
 
 
+def convert(self, value, src, dst):
+        if src == dst or src is None or dst is None or value is None:
+            return value
+        elif dst == 'mm':
+            return value * 25.4
+        else:
+            return value / 25.4
+
 class GerberFile(CamFile):
     """ A class representing a single gerber file
 
@@ -55,8 +63,6 @@ class GerberFile(CamFile):
 
     def merge(self, other):
         """ Merge other GerberFile into this one """
-        # FIXME unit handling
-
         self.comments += other.comments
 
         # dedup apertures
@@ -68,18 +74,19 @@ class GerberFile(CamFile):
                 new_apertures[gbr] = ap
             else:
                 replace_apertures[id(ap)] = new_apertures[gbr]
-        self.apertures = new_apertures
+        self.apertures = list(new_apertures.values())
 
         self.objects += other.objects
         for obj in self.objects:
-            if (ap := replace_apertures.get(id(getattr(obj, 'aperture')))):
+            # If object has an aperture attribute, replace that aperture.
+            if (ap := replace_apertures.get(id(getattr(obj, 'aperture', None)))):
                 obj.aperture = ap
 
         # dedup aperture macros
         macros = { m.to_gerber(): m
                 for m in [ GenericMacros.circle, GenericMacros.rect, GenericMacros.obround, GenericMacros.polygon] }
-        for ap in new_apertures:
-            if isinstance(aperture, apertures.ApertureMacroInstance):
+        for ap in new_apertures.values():
+            if isinstance(ap, apertures.ApertureMacroInstance):
                 macro_grb = ap.macro.to_gerber() # use native unit to compare macros
                 if macro_grb in macros:
                     ap.macro = macros[macro_grb]
@@ -246,6 +253,7 @@ class GraphicsState:
     image_polarity : str = 'positive' # IP image polarity; deprecated
     point : tuple = None
     aperture : apertures.Aperture = None
+    file_settings : FileSettings = None
     interpolation_mode : InterpolationModeStmt = LinearModeStmt
     multi_quadrant_mode : bool = None # used only for syntax checking
     aperture_mirroring = (False, False) # LM mirroring (x, y)
@@ -261,8 +269,9 @@ class GraphicsState:
     aperture_map = {}
 
 
-    def __init__(self, aperture_map=None):
+    def __init__(self, file_settings=None, aperture_map=None):
         self._mat = None
+        self.file_settings = file_settings
         if aperture_map is not None:
             self.aperture_map = aperture_map
 
@@ -333,7 +342,9 @@ class GraphicsState:
 
     def flash(self, x, y):
         self.update_point(x, y)
-        return go.Flash(*self.map_coord(*self.point), self.aperture, polarity_dark=self.polarity_dark)
+        return go.Flash(*self.map_coord(*self.point), self.aperture,
+                polarity_dark=self.polarity_dark,
+                unit=self.file_settings.unit)
 
     def interpolate(self, x, y, i=None, j=None, aperture=True):
         if self.point is None:
@@ -363,21 +374,24 @@ class GraphicsState:
                 return self._create_arc(old_point, self.map_coord(*self.point), (i, j), aperture)
 
     def _create_line(self, old_point, new_point, aperture=True):
-        return go.Line(*old_point, *new_point, self.aperture if aperture else None, polarity_dark=self.polarity_dark)
+        return go.Line(*old_point, *new_point, self.aperture if aperture else None,
+                polarity_dark=self.polarity_dark, unit=self.file_settings.unit)
 
     def _create_arc(self, old_point, new_point, control_point, aperture=True):
         direction = 'ccw' if self.interpolation_mode == CircularCCWModeStmt else 'cw'
         return go.Arc(*old_point, *new_point,* self.map_coord(*control_point, relative=True),
-                flipped=(direction == 'cw'), aperture=(self.aperture if aperture else None), polarity_dark=self.polarity_dark)
+                flipped=(direction == 'cw'), aperture=(self.aperture if aperture else None),
+                polarity_dark=self.polarity_dark, unit=self.file_settings.unit)
 
-    def update_point(self, x, y):
+    def update_point(self, x, y, unit=None):
         old_point = self.point
         if x is None:
             x = self.point[0]
         if y is None:
             y = self.point[1]
-        if self.point != (x, y):
-            self.point = (x, y)
+        if unit == 'inch':
+            x, y = x*25.4, y*25.4
+        self.point = (x, y)
         return old_point
 
     # Helpers for gerber generation
@@ -391,10 +405,17 @@ class GraphicsState:
             self.aperture = aperture
             yield ApertureStmt(self.aperture_map[id(aperture)])
 
-    def set_current_point(self, point):
-        if self.point != point:
-            self.point = point
-            yield MoveStmt(*point)
+    def set_current_point(self, point, unit=None):
+        # FIXME use math.isclose for point comparisons here and elsewhere due to converted coords
+        # FIXME maybe even calculate appropriate precision given file_settings.notation
+        if unit == 'inch':
+            point_mm = point[0]*25.4, point[1]*25.4
+        else:
+            point_mm = point
+
+        if self.point != point_mm:
+            self.point = point_mm
+            yield MoveStmt(*point, unit=unit)
 
     def set_interpolation_mode(self, mode):
         if self.interpolation_mode != mode:
@@ -446,7 +467,7 @@ class GerberParser:
         self.include_dir = include_dir
         self.include_stack = []
         self.file_settings = FileSettings()
-        self.graphics_state = GraphicsState()
+        self.graphics_state = GraphicsState(file_settings=self.file_settings)
         self.aperture_map = {}
         self.aperture_macros = {}
         self.current_region = None
@@ -566,7 +587,9 @@ class GerberParser:
                     # Start a new region for every outline. As gerber has no concept of fill rules or winding numbers,
                     # it does not make a graphical difference, and it makes the implementation slightly easier.
                     self.target.objects.append(self.current_region)
-                    self.current_region = go.Region(polarity_dark=self.graphics_state.polarity_dark)
+                    self.current_region = go.Region(
+                            polarity_dark=self.graphics_state.polarity_dark,
+                            unit=self.file_settings.unit)
 
             else: # D03
                 if self.current_region is None:
@@ -699,7 +722,9 @@ class GerberParser:
         self.target.comments.append(match["comment"])
 
     def _parse_region_start(self, _match):
-        self.current_region = go.Region(polarity_dark=self.graphics_state.polarity_dark)
+        self.current_region = go.Region(
+                polarity_dark=self.graphics_state.polarity_dark,
+                unit=self.file_settings.unit)
 
     def _parse_region_end(self, _match):
         if self.current_region is None:

gerbonara/gerber/tests/image_support.py

@@ -109,12 +109,61 @@ def gerber_difference(reference, actual, diff_out=None, svg_transform=None, size
         with svg_soup(act_svg.name) as soup:
             cleanup_clips(soup)
 
-        # FIXME DEBUG
-        shutil.copyfile(act_svg.name, '/tmp/test-act.svg')
-        shutil.copyfile(ref_svg.name, '/tmp/test-ref.svg')
-
         return svg_difference(ref_svg.name, act_svg.name, diff_out=diff_out)
 
+def gerber_difference_merge(ref1, ref2, actual, diff_out=None, composite_out=None, svg_transform1=None, svg_transform2=None, size=(10,10)):
+    with tempfile.NamedTemporaryFile(suffix='.svg') as act_svg,\
+        tempfile.NamedTemporaryFile(suffix='.svg') as ref1_svg,\
+        tempfile.NamedTemporaryFile(suffix='.svg') as ref2_svg:
+
+        gbr_to_svg(ref1, ref1_svg.name, size=size)
+        gbr_to_svg(ref2, ref2_svg.name, size=size)
+        gbr_to_svg(actual, act_svg.name, size=size)
+
+        with svg_soup(ref1_svg.name) as soup1:
+            if svg_transform1 is not None:
+                soup1.find('g', attrs={'id': 'surface1'})['transform'] = svg_transform1
+            cleanup_clips(soup1)
+
+            with svg_soup(ref2_svg.name) as soup2:
+                if svg_transform2 is not None:
+                    soup2.find('g', attrs={'id': 'surface1'})['transform'] = svg_transform2
+                cleanup_clips(soup2)
+
+                defs1 = soup1.find('defs')
+                if not defs1:
+                    defs1 = soup1.new_tag('defs')
+                    soup1.find('svg').insert(0, defs1)
+
+                defs2 = soup2.find('defs')
+                if defs2:
+                    defs2 = defs2.extract()
+                    # explicitly convert .contents into list here and below because else bs4 stumbles over itself
+                    # iterating because we modify the tree in the loop body.
+                    for c in list(defs2.contents):
+                        if hasattr(c, 'attrs'):
+                            c['id'] = 'gn-merge-b-' + c.attrs.get('id', str(id(c)))
+                        defs1.append(c)
+
+                for use in soup2.find_all('use', recursive=True):
+                    if (href := use.get('xlink:href', '')).startswith('#'):
+                        use['xlink:href'] = f'#gn-merge-b-{href[1:]}'
+
+                svg1 = soup1.find('svg')
+                for c in list(soup2.find('svg').contents):
+                    if hasattr(c, 'attrs'):
+                        c['id'] = 'gn-merge-b-' + c.attrs.get('id', str(id(c)))
+                    svg1.append(c)
+            # FIXME prefix all group ids with "b-"
+
+        if composite_out:
+            shutil.copyfile(ref1_svg.name, composite_out)
+
+        with svg_soup(act_svg.name) as soup:
+            cleanup_clips(soup)
+
+        return svg_difference(ref1_svg.name, act_svg.name, diff_out=diff_out)
+
 def svg_difference(reference, actual, diff_out=None):
     with tempfile.NamedTemporaryFile(suffix='-ref.png') as ref_png,\
         tempfile.NamedTemporaryFile(suffix='-act.png') as act_png:

gerbonara/gerber/tests/test_rs274x.py

@@ -17,7 +17,7 @@ import pytest
 from ..rs274x import GerberFile
 from ..cam import FileSettings
 
-from .image_support import gerber_difference
+from .image_support import gerber_difference, gerber_difference_merge
 
 
 deg_to_rad = lambda a: a/180 * math.pi
@@ -28,9 +28,10 @@ reference_path = lambda reference: Path(__file__).parent / 'resources' / referen
 @pytest.fixture
 def temp_files(request):
     with tempfile.NamedTemporaryFile(suffix='.gbr') as tmp_out_gbr,\
+         tempfile.NamedTemporaryFile(suffix='.svg') as tmp_out_svg,\
          tempfile.NamedTemporaryFile(suffix='.png') as tmp_out_png:
 
-        yield Path(tmp_out_gbr.name), Path(tmp_out_png.name)
+        yield Path(tmp_out_gbr.name), Path(tmp_out_svg.name), Path(tmp_out_png.name)
 
         if request.node.rep_call.failed:
             module, _, test_name = request.node.nodeid.rpartition('::')
@@ -39,14 +40,27 @@ def temp_files(request):
             test_name = re.sub(r'[^\w\d]', '_', test_name)
             fail_dir.mkdir(exist_ok=True)
             perm_path_gbr = fail_dir / f'failure_{test_name}.gbr'
+            perm_path_svg = fail_dir / f'failure_{test_name}.svg'
             perm_path_png = fail_dir / f'failure_{test_name}.png'
             shutil.copy(tmp_out_gbr.name, perm_path_gbr)
+            if Path(tmp_out_svg.name).is_file():
+                shutil.copy(tmp_out_svg.name, perm_path_svg)
             shutil.copy(tmp_out_png.name, perm_path_png)
             print(f'Failing output saved to {perm_path_gbr}')
-            print(f'Reference file is {reference_path(request.node.funcargs["reference"])}')
+            args = request.node.funcargs
+            if 'reference' in args:
+                print(f'Reference file is {reference_path(args["reference"])}')
+            else:
+                print(f'Reference file A is {reference_path(args["file_a"])}')
+                print(f'Reference file B is {reference_path(args["file_b"])}')
             print(f'Difference image saved to {perm_path_png}')
+            if Path(tmp_out_svg.name).is_file():
+                print(f'Sum SVG saved to {perm_path_svg}')
             print(f'gerbv command line:')
-            print(f'gerbv {perm_path_gbr} {reference_path(request.node.funcargs["reference"])}')
+            if 'reference' in args:
+                print(f'gerbv {perm_path_gbr} {reference_path(request.node.funcargs["reference"])}')
+            else:
+                print(f'gerbv {perm_path_gbr} {reference_path(args["file_a"])} {reference_path(args["file_b"])}')
 
 to_gerbv_svg_units = lambda val, unit='mm': val*72 if unit == 'inch' else val/25.4*72
 
@@ -109,11 +123,12 @@ MIN_REFERENCE_FILES = [
     'eagle_files/copper_bottom_l4.gbr'
     ]
 
+
 @pytest.mark.filterwarnings('ignore:Deprecated.*statement found.*:DeprecationWarning')
 @pytest.mark.filterwarnings('ignore::SyntaxWarning')
 @pytest.mark.parametrize('reference', REFERENCE_FILES)
 def test_round_trip(temp_files, reference):
-    tmp_gbr, tmp_png = temp_files
+    tmp_gbr, _tmp_svg, tmp_png = temp_files
     ref = reference_path(reference)
 
     GerberFile.open(ref).save(tmp_gbr)
@@ -135,7 +150,7 @@ def test_rotation(temp_files, reference, angle):
         # gerbv's rendering of this is broken, the hole is missing.
         return
 
-    tmp_gbr, tmp_png = temp_files
+    tmp_gbr, _tmp_svg, tmp_png = temp_files
     ref = reference_path(reference)
 
     f = GerberFile.open(ref)
@@ -156,7 +171,7 @@ def test_rotation_center(temp_files, reference, angle, center):
     if 'flash_rectangle' in reference and angle in (30, 1024):
         # gerbv's rendering of this is broken, the hole is missing.
         return
-    tmp_gbr, tmp_png = temp_files
+    tmp_gbr, _tmp_svg, tmp_png = temp_files
     ref = reference_path(reference)
 
     f = GerberFile.open(ref)
@@ -165,7 +180,7 @@ def test_rotation_center(temp_files, reference, angle, center):
 
     # calculate circle center in SVG coordinates 
     size = (10, 10) # inches
-    cx, cy = to_gerbv_svg_units(center[0]), to_gerbv_svg_units(10, 'inch')-to_gerbv_svg_units(center[1], 'mm')
+    cx, cy = to_gerbv_svg_units(center[0]), to_gerbv_svg_units(size[1], 'inch')-to_gerbv_svg_units(center[1], 'mm')
     mean, _max, hist = gerber_difference(ref, tmp_gbr, diff_out=tmp_png,
             svg_transform=f'rotate({angle} {cx} {cy})',
             size=size)
@@ -178,7 +193,7 @@ def test_rotation_center(temp_files, reference, angle, center):
 @pytest.mark.parametrize('reference', MIN_REFERENCE_FILES)
 @pytest.mark.parametrize('offset', TEST_OFFSETS)
 def test_offset(temp_files, reference, offset):
-    tmp_gbr, tmp_png = temp_files
+    tmp_gbr, _tmp_svg, tmp_png = temp_files
     ref = reference_path(reference)
 
     f = GerberFile.open(ref)
@@ -201,7 +216,7 @@ def test_combined(temp_files, reference, angle, center, offset):
     if 'flash_rectangle' in reference and angle in (30, 1024):
         # gerbv's rendering of this is broken, the hole is missing.
         return
-    tmp_gbr, tmp_png = temp_files
+    tmp_gbr, _tmp_svg, tmp_png = temp_files
     ref = reference_path(reference)
 
     f = GerberFile.open(ref)
@@ -210,7 +225,7 @@ def test_combined(temp_files, reference, angle, center, offset):
     f.save(tmp_gbr, settings=FileSettings(unit=f.unit, number_format=(4,7)))
 
     size = (10, 10) # inches
-    cx, cy = to_gerbv_svg_units(center[0]), to_gerbv_svg_units(10, 'inch')-to_gerbv_svg_units(center[1], 'mm')
+    cx, cy = to_gerbv_svg_units(center[0]), to_gerbv_svg_units(size[1], 'inch')-to_gerbv_svg_units(center[1], 'mm')
     dx, dy = to_gerbv_svg_units(offset[0]), -to_gerbv_svg_units(offset[1])
     mean, _max, hist = gerber_difference(ref, tmp_gbr, diff_out=tmp_png,
             svg_transform=f'translate({dx} {dy}) rotate({angle} {cx} {cy})',
@@ -219,3 +234,54 @@ def test_combined(temp_files, reference, angle, center, offset):
     assert hist[9] < 100
     assert hist[3:].sum() < 1e-3*hist.size
 
+@pytest.mark.filterwarnings('ignore:Deprecated.*statement found.*:DeprecationWarning')
+@pytest.mark.filterwarnings('ignore::SyntaxWarning')
+@pytest.mark.parametrize('file_a', MIN_REFERENCE_FILES)
+@pytest.mark.parametrize('file_b', [
+    'example_two_square_boxes.gbr',
+    'example_outline_with_arcs.gbr',
+    'example_am_exposure_modifier.gbr',
+    'bottom_silk.GBO',
+    'eagle_files/copper_bottom_l4.gbr', ])
+@pytest.mark.parametrize('angle', [0, 10, 90])
+@pytest.mark.parametrize('offset', [(0, 0, 0, 0), (100, 0, 0, 0), (0, 0, 0, 100), (100, 0, 0, 100)])
+def test_compositing(temp_files, file_a, file_b, angle, offset):
+
+    # TODO bottom_silk.GBO renders incorrectly with gerbv: the outline does not exist in svg. In GUI, the logo only
+    # renders at very high magnification. Skip, and once we have our own SVG export maybe use that instead. Or just use
+    # KiCAD's gerbview.
+    # TODO check if this and the issue with aperture holes not rendering in test_combined actually are bugs in gerbv
+    # and fix/report upstream.
+    if file_a == 'bottom_silk.GBO' or file_b == 'bottom_silk.GBO':
+        return
+
+    tmp_gbr, tmp_svg, tmp_png = temp_files
+    ref_a = reference_path(file_a)
+    ref_b = reference_path(file_b)
+
+    ax, ay, bx, by = offset
+    grb_a = GerberFile.open(ref_a)
+    grb_a.rotate(deg_to_rad(angle))
+    grb_a.offset(ax, ay)
+
+    grb_b = GerberFile.open(ref_b)
+    grb_b.offset(bx, by)
+
+    grb_a.merge(grb_b)
+    grb_a.save(tmp_gbr, settings=FileSettings(unit=grb_a.unit, number_format=(4,7)))
+
+    size = (10, 10) # inches
+    ax, ay = to_gerbv_svg_units(ax), -to_gerbv_svg_units(ay)
+    bx, by = to_gerbv_svg_units(bx), -to_gerbv_svg_units(by)
+    # note that we have to specify cx, cy even if we rotate around the origin since gerber's origin lies at (x=0
+    # y=+document size) in SVG's coordinate space because svg's y axis is flipped compared to gerber's.
+    cx, cy = 0, to_gerbv_svg_units(size[1], 'inch')
+    mean, _max, hist = gerber_difference_merge(ref_a, ref_b, tmp_gbr, composite_out=tmp_svg, diff_out=tmp_png,
+            svg_transform1=f'translate({ax} {ay}) rotate({angle} {cx} {cy})',
+            svg_transform2=f'translate({bx} {by})',
+            size=size)
+    assert mean < 1e-3
+    assert hist[9] < 100
+    assert hist[3:].sum() < 1e-3*hist.size
+
+