fix(contour): trace every ROI border in order, in absolute coordinates

docs/examples/fractal_worked_example.py

@@ -78,18 +78,27 @@ def disk(r):
 check(fd(disk(40)), [1.686, 1.032], "disk 40")
 
 # One-pixel-wide straight lines: box counting reaches 1 only at power-of-two lengths, and the
-# divider walk reads above 1 because the closing chord spans the whole line
-check(fd(np.ones((1, 16))), [1.0, 1.190], "line 16")
-check(fd(np.ones((1, 10))), [0.7551, 1.205], "line 10")
+# divider walk reads above 1 because the walk runs along the line and back, and the chords that
+# span an end fold back
+check(fd(np.ones((1, 16))), [1.0, 1.035], "line 16")
+check(fd(np.ones((1, 10))), [0.7551, 1.101], "line 10")
 check(fd(np.ones((1, 100)))[:1], [0.9257], "line 100")
 
-# Two separate squares form one ROI with two contours; the divider walk joins them
+# Two separate squares form one ROI with two borders, each walked on its own: the divider value is
+# that of one square
 two = np.zeros((40, 80), np.uint32)
 two[5:30, 5:30] = 1
 two[5:30, 45:70] = 1
-check(fd(two), [1.423, 1.107], "two squares")
+check(fd(two), [1.423, 1.0], "two squares")
 check(fd(square(25)), [1.842, 1.0], "one square")
 
+# A hole's border, shorter than the large strides, is measured whole
+holed = np.zeros((30, 30), np.uint32)
+holed[3:27, 3:27] = 1
+check(fd(holed)[1:], [1.005], "24 x 24 square")
+holed[10:20, 10:20] = 0
+check(fd(holed)[1:], [1.045], "24 x 24 square with a 10 x 10 hole")
+
 # A 40 x 40 square and one separate pixel. In memory an intensity of 4294967295 makes the box
 # count's marker overflow its 32-bit conversion, which is undefined, so only 100 is checked here
 lone = np.zeros((80, 80), np.uint32)
@@ -143,23 +152,37 @@ def from_files(m, intensity=None, single_roi=False, codes=F + ["PERIMETER"], **k
 
 check(from_files(mask, ram_limit=0)[0], [1.3966, 1.1311, 26.0], "example out of core")
 # pins known defect https://github.com/PolusAI/nyxus/issues/509
-check(from_files(two, ram_limit=0)[0][1:2], [1.0], "two squares out of core")
+check(from_files(holed, ram_limit=0)[0][1:2], [1.005], "square with a hole out of core")
 for value in (100, 4294967295):
     inten = lone * 100
     inten[70, 70] = value
     check(from_files(lone, inten, ram_limit=0)[0][:1], [1.4156], f"lone pixel {value} out of core")
 check(from_files(mask)[0][:2], [1.3732, 1.1311], "example, no anisotropy")
-check(from_files(mask, anisotropy_x=2)[0][:2], [1.2459, 1.1731], "anisotropy_x=2")
-check(from_files(mask, anisotropy_y=3)[0][:2], [1.4458, 1.2546], "anisotropy_y=3")
+check(from_files(mask, anisotropy_x=2)[0][:2], [1.2459, 1.1566], "anisotropy_x=2")
+check(from_files(mask, anisotropy_y=3)[0][:2], [1.4458, 1.1635], "anisotropy_y=3")
 # With anisotropy the out-of-core contour differs, so a single-contour disk differs too
-check(from_files(disk(8), anisotropy_x=2)[0][1:2], [1.1431], "disk 8, anisotropy_x=2")
+check(from_files(disk(8), anisotropy_x=2)[0][1:2], [1.0917], "disk 8, anisotropy_x=2")
 # pins known defect https://github.com/PolusAI/nyxus/issues/509
 check(from_files(disk(8), anisotropy_x=2, ram_limit=0)[0][1:2], [1.0482], "disk 8, anisotropy_x=2, out of core")
+# Without anisotropy a single contour can still differ out of core, which walks a one-pixel-wide
+# part once where the in-memory walk goes along and back: a 1 x 16 line and an 8 x 8 square with a
+# one-pixel-wide branch, at the image origin
+line16 = np.zeros((4, 19), np.uint32)
+line16[0, 0:16] = 1
+check(from_files(line16)[0][1:2], [1.035], "line 16")
+# pins known defect https://github.com/PolusAI/nyxus/issues/509
+check(from_files(line16, ram_limit=0)[0][1:2], [1.190], "line 16, out of core")
+branch = np.zeros((18, 18), np.uint32)
+branch[0:8, 0:8] = 1
+branch[3, 8:16] = 1
+check(from_files(branch)[0][1:2], [1.190], "square with a branch")
+# pins known defect https://github.com/PolusAI/nyxus/issues/509
+check(from_files(branch, ram_limit=0)[0][1:2], [1.102], "square with a branch, out of core")
 labels = two.copy()
 labels[5:30, 45:70] = 2
 for row in from_files(labels, (labels > 0) * 100):
     check(row[:2], [1.8416, 1.0], "two labels")
-check(from_files(labels, (labels > 0) * 100, mergerois=True)[0][:2], [1.4231, 1.107], "merged")
+check(from_files(labels, (labels > 0) * 100, mergerois=True)[0][:2], [1.4231, 1.0], "merged")
 # pins known defect https://github.com/PolusAI/nyxus/issues/493: the box is 201 x 201, which rounds
 # up to the same 256 canvas as 200 x 200
 check(from_files(np.ones((200, 200)), single_roi=True)[0][:2], [1.8703, 1.0], "whole slide")

docs/examples/moments_worked_example.py

@@ -54,8 +54,8 @@ def same(got, expected, rel=5e-4):
 assert same(features["SPAT_MOMENT_10"][0] / features["SPAT_MOMENT_00"][0], 3.275)
 assert same(features["SPAT_MOMENT_01"][0] / features["SPAT_MOMENT_00"][0], 2.925)
 
-# 48 x 40 rectangle: Nyxus, and the same moment with exact distances to the contour as traced
-# (one pixel right of and below the ROI's border) and to the ROI's own border pixels
+# 48 x 40 rectangle: Nyxus, and the same moment with exact distances to the contour, the ROI's
+# border pixels
 H, W = 40, 48
 rect = np.ones((H, W), dtype=np.uint32)
 w00 = Nyxus(["WEIGHTED_SPAT_MOMENT_00"]).featurize(rect, rect)["WEIGHTED_SPAT_MOMENT_00"][0]
@@ -70,13 +70,12 @@ def weight_sum(cx, cy):
     return float(np.sum(np.log(d + 0.001)))
 
 
-assert same(w00, 2498.7), w00
-assert same(weight_sum(bx + 1, by + 1), 2422.7)
+assert same(w00, 1895.0), w00
 assert same(weight_sum(bx, by), 1821.3)
 
-# Digital disks: the sum of the weights changes sign between radius 5 and radius 8, and the
+# Digital disks: the sum of the weights changes sign between radius 5 and radius 12, and the
 # odd-order normalized weighted moments follow it
-for r, m00, ncm30 in ((5, -68.94, 0.0), (8, 6.806, 7881.0), (20, 1549.1, 0.007713)):
+for r, m00, ncm30 in ((5, -153.9, 0.0), (12, 2.195, -33355.0), (20, 1210.6, -0.001704)):
     n = 2 * r + 3
     yy, xx = np.mgrid[0:n, 0:n]
     disk = (((yy - r - 1) ** 2 + (xx - r - 1) ** 2) <= r * r).astype(np.uint32)
@@ -105,9 +104,9 @@ def weight_sum(cx, cy):
 for code in ("CENTRAL_MOMENT_01", "CENTRAL_MOMENT_10", "IMOM_CM_01", "IMOM_CM_10"):
     assert abs(features[code][0]) < 1e-13, (code, features[code][0])
 
-# The 48 x 40 rectangle has 172 border pixels; 85 ROI pixels lie on the traced (shifted) contour
-on_traced = (bx + 1 < W) & (by + 1 < H)
-assert len(bx) == 172 and on_traced.sum() == 85
+# The 48 x 40 rectangle's contour is its 172 border pixels: with intensity 1 they sum to 172
+edge = Nyxus(["EDGE_INTEGRATED_INTENSITY"]).featurize(rect, rect)["EDGE_INTEGRATED_INTENSITY"][0]
+assert len(bx) == 172 and edge == 172
 
 # Out of core (ram_limit=0, file input): the weighted moments of a 245-pixel ellipse depend on where
 # its bounding box starts; the 118 unweighted moments do not
@@ -129,23 +128,34 @@ def ellipse(n, cx, cy):
 weighted = [c for c in codes if c.startswith(("WEIGHTED_", "WT_NORM_", "IMOM_W"))]
 assert len(weighted) == 62
 # pins known defect https://github.com/PolusAI/nyxus/issues/509
-for n, cx, cy, corner, ooc in ((30, 12, 10, (2, 2), 7.484), (26, 13, 13, (3, 5), 58.50),
-                               (40, 20, 20, (10, 12), 338.2)):
+for n, cx, cy, corner, ooc in ((30, 12, 10, (2, 2), 35.62), (26, 13, 13, (3, 5), 78.42),
+                               (40, 20, 20, (10, 12), 425.0)):
     e = ellipse(n, cx, cy)
     ys, xs = np.nonzero(e)
     assert e.sum() == 245 and (xs.min(), ys.min()) == corner
     mem, out = from_files(e), from_files(e, ram_limit=0)
-    assert same(mem["WEIGHTED_SPAT_MOMENT_00"][0], 49.81)
+    assert same(mem["WEIGHTED_SPAT_MOMENT_00"][0], -93.36)
     assert same(out["WEIGHTED_SPAT_MOMENT_00"][0], ooc), (corner, out["WEIGHTED_SPAT_MOMENT_00"][0])
     for code in codes:
         if code not in weighted:
             assert math.isclose(out[code][0], mem[code][0], rel_tol=1e-9, abs_tol=1e-9), code
-# At the image origin: the same ellipse still differs, the rectangle does not
+# At the image origin: the ellipse keeps WEIGHTED_SPAT_MOMENT_00 but not its weighted first moment,
+# the rectangle gets the in-memory values
 e = ellipse(30, 12, 10)[2:, 2:]
+mem, out = from_files(e), from_files(e, ram_limit=0)
+assert same(out["WEIGHTED_SPAT_MOMENT_00"][0], -93.36)
+assert same(mem["WEIGHTED_SPAT_MOMENT_10"][0], -1158.2)
 # pins known defect https://github.com/PolusAI/nyxus/issues/509
-assert same(from_files(e, ram_limit=0)["WEIGHTED_SPAT_MOMENT_00"][0], 43.75)
+assert same(out["WEIGHTED_SPAT_MOMENT_10"][0], -708.9)
 corner_rect = np.zeros((50, 60), dtype=np.uint16)
 corner_rect[:40, :48] = 1
 mem, out = from_files(corner_rect), from_files(corner_rect, ram_limit=0)
 for code in codes:
     assert math.isclose(out[code][0], mem[code][0], rel_tol=1e-9, abs_tol=1e-9), code
+# A ROI with a hole at the origin differs too: the out-of-core contour leaves the hole out
+# pins known defect https://github.com/PolusAI/nyxus/issues/509
+holed = np.zeros((20, 20), dtype=np.uint16)
+holed[:9, :9] = 1
+holed[3:6, 3:6] = 0
+assert same(from_files(holed)["WEIGHTED_SPAT_MOMENT_00"][0], -274.2)
+assert same(from_files(holed, ram_limit=0)["WEIGHTED_SPAT_MOMENT_00"][0], -194.7)

docs/examples/morphology_worked_example.py

@@ -68,12 +68,11 @@ def run(image, labels, **kwargs):
 # The worked example as the page walks through it
 check(one, {"POLYGONALITY_AVE": -1, "HEXAGONALITY_AVE": -1, "EDGE_INTEGRATED_INTENSITY": 228}, "ROI 1")
 assert round(one["EDGE_INTEGRATED_INTENSITY"][0] / one["EDGE_MEAN_INTENSITY"][0]) == 44   # contour pixels
-# pins known defect https://github.com/PolusAI/nyxus/issues/497
 check(two, {"AREA_PIXELS_COUNT": 15, "PERIMETER": 12, "CIRCULARITY": 1.144, "ASPECT_RATIO": 0.6,
             "ORIENTATION": 90, "CONVEX_HULL_AREA": 15, "PERCENT_TOUCHING": 41.67,
             "CLOSEST_NEIGHBOR1_ANG": 174.4, "DIAMETER_MIN_ENCLOSING_CIRCLE": 4.4723,
-            "DIAMETER_CIRCUMSCRIBING_CIRCLE": 10, "DIAMETER_INSCRIBING_CIRCLE": 2,
-            "ROI_RADIUS_MEAN": 0.6943, "ROI_RADIUS_MAX": 1.414, "ROI_RADIUS_MEDIAN": 1,
+            "DIAMETER_CIRCUMSCRIBING_CIRCLE": 4.472, "DIAMETER_INSCRIBING_CIRCLE": 2,
+            "ROI_RADIUS_MEAN": 0.2, "ROI_RADIUS_MAX": 1, "ROI_RADIUS_MEDIAN": 0,
             "MAXCHORDS_MIN": 2, "MAXCHORDS_MEAN": 3.722}, "ROI 2")
 
 # Settings that change nothing in memory: IBSI mode, the grey-level depth, the anisotropy keywords
@@ -107,23 +106,35 @@ def run(image, labels, **kwargs):
 check(run(hole, mask), {"AREA_PIXELS_COUNT": 116, "MAXCHORDS_MEAN": 8.8, "ALLCHORDS_MEAN": 5.197},
       "zero pixel")
 
-# PERIMETER of a ROI with a hole, and of a label in two pieces
+# PERIMETER of a ROI with a hole, and of a label in two pieces: the sum of their closed walks. The
+# circles of the holed square are measured to its outer border only.
 sq = np.zeros((13, 13), dtype=np.uint32)
 sq[2:11, 2:11] = 1
 sq[5:8, 5:8] = 0
-check(run(sq * 4, sq), {"PERIMETER": 50.01, "EULER_NUMBER": 0}, "holed square")
+check(run(sq * 4, sq), {"PERIMETER": 45.66, "EULER_NUMBER": 0, "DIAMETER_CIRCUMSCRIBING_CIRCLE": 11.31,
+                        "DIAMETER_INSCRIBING_CIRCLE": 8}, "holed square")
 pieces = np.zeros((10, 16), dtype=np.uint32)
 pieces[2:7, 2:6] = 1
 pieces[2:7, 10:14] = 1
-check(run(pieces * 4, pieces), {"PERIMETER": 42.12, "EULER_NUMBER": 2}, "two pieces")
+check(run(pieces * 4, pieces), {"PERIMETER": 28, "EULER_NUMBER": 2}, "two pieces")
 one_piece = np.zeros((10, 10), dtype=np.uint32)
 one_piece[2:7, 2:6] = 1
 check(run(one_piece * 4, one_piece), {"PERIMETER": 14}, "one piece")
 
-# A 1 x 13 bar: hull of two vertices, Nassenstein 0, Martin only rounding residue below 1e-6
+# A 1 x 13 bar: walked on both sides, so PERIMETER is 24; hull of two vertices, Nassenstein 0,
+# Martin only rounding residue below 1e-6
 line = np.zeros((7, 20), dtype=np.uint32)
 line[3, 3:16] = 1
 flat = run(line * 100, line)
+check(flat, {"PERIMETER": 24}, "1 x 13 bar")
+
+# A single pixel is its own contour; two pixels are walked there and back
+dot = run(np.pad(np.full((1, 1), 7, np.uint32), 2), np.pad(np.ones((1, 1), np.uint32), 2))
+check(dot, {"PERIMETER": 0, "EDGE_MEAN_INTENSITY": 7, "EDGE_INTEGRATED_INTENSITY": 7, "CIRCULARITY": 0,
+            "DIAMETER_MIN_ENCLOSING_CIRCLE": 2e-4, "DIAMETER_CIRCUMSCRIBING_CIRCLE": 0,
+            "DIAMETER_INSCRIBING_CIRCLE": 0, "ROI_RADIUS_MAX": 0, "CONVEX_HULL_AREA": 1}, "one pixel")
+pair = run(np.pad(np.full((1, 2), 7, np.uint32), 2), np.pad(np.ones((1, 2), np.uint32), 2))
+check(pair, {"PERIMETER": 2, "CIRCULARITY": 2.507}, "two pixels")
 for stat in ("MIN", "MAX", "MEAN", "MEDIAN", "STDDEV", "MODE"):
     assert flat[f"STAT_NASSENSTEIN_DIAM_{stat}"][0] == 0, stat
     assert 0 <= flat[f"STAT_MARTIN_DIAM_{stat}"][0] < 1e-6, stat
@@ -174,7 +185,7 @@ def alone(codes, image=img, labels=mask):
 circle_codes = ["DIAMETER_MIN_ENCLOSING_CIRCLE", "DIAMETER_CIRCUMSCRIBING_CIRCLE", "DIAMETER_INSCRIBING_CIRCLE"]
 assert (alone(circle_codes)[circle_codes].to_numpy() == 0).all()
 check(alone(circle_codes + ["PERIMETER"]),
-      {"DIAMETER_MIN_ENCLOSING_CIRCLE": 14.17, "DIAMETER_CIRCUMSCRIBING_CIRCLE": 35.61}, "circles with PERIMETER")
+      {"DIAMETER_MIN_ENCLOSING_CIRCLE": 14.17, "DIAMETER_CIRCUMSCRIBING_CIRCLE": 30}, "circles with PERIMETER")
 hex_codes = ["HEXAGONALITY_AVE", "HEXAGONALITY_STDDEV"]
 assert (alone(hex_codes, dimg, disk)[hex_codes].to_numpy() == 0).all()