diff --git a/src/nyx/features/fractal_dim.cpp b/src/nyx/features/fractal_dim.cpp index 8bd6a238..418e82d9 100644 --- a/src/nyx/features/fractal_dim.cpp +++ b/src/nyx/features/fractal_dim.cpp @@ -1,5 +1,7 @@ +#include #include "fractal_dim.h" #include "image_matrix.h" +#include "../helpers/helpers.h" using namespace Nyxus; @@ -17,144 +19,168 @@ void FractalDimensionFeature::calculate (LR& r, const Fsettings& s) void FractalDimensionFeature::calculate_boxcount_fdim (LR & r) { - Power2PaddedImageMatrix pim(r.raw_pixels, r.aabb, 1, 1.0); // image matrix of the mask + const std::vector& px = r.raw_pixels; + if (px.size() < 2) + { + box_count_fd = 0.; + return; + } - // Debug - // pim.print("Padded"); + // Box counting is registration sensitive: the old code padded to a power of two strictly + // larger than the ROI and centered it, misaligning the ROI with the coarse boxes and biasing + // the dimension low (a filled square read 1.75 instead of 2.0). We pad tight (ceil_pow2) and + // align to the ROI origin. Shifting grids (min count over grid origins, as FracLac does) + // remove the residual bias, but that only matters when few box sizes fit (small ROIs), and + // those are cheap to shift. So auto-switch: + // - large ROI -> single origin grid via the padded mask matrix (fast early-exit tile scan) + // - small ROI -> shift the grid over a few origins and take the minimum box count + int bigSide = (int)std::max(r.aabb.get_width(), r.aabb.get_height()); + int paddedSide = Nyxus::ceil_pow2(bigSide); - // Square box coverage statistics - std::vector> coverage; - int s = pim.width; // square tile size - for (; s > 1; s /= 2) - { - int cnt = 0; - size_t areaCov = 0; - int n_tiles = pim.width / s; - for (int r = 0; r < n_tiles; r++) - for (int c = 0; c < n_tiles; c++) - { - // check the tile for coerage - if (pim.tile_contains_signal(r, c, s)) - cnt++; - } + std::vector> coverage; - coverage.push_back({ s, cnt }); + if (paddedSide > 32) + { + // Large ROI: single aligned grid. tile_contains_signal early-exits on the first non-zero + // pixel, so this is O(#tiles) for filled ROIs - much cheaper than visiting every pixel. + Power2PaddedImageMatrix pim(px, r.aabb, 1, 1.0); + for (int s = pim.width; s > 1; s /= 2) + { + int cnt = 0, nt = pim.width / s; + for (int tr = 0; tr < nt; tr++) + for (int tc = 0; tc < nt; tc++) + if (pim.tile_contains_signal(tr, tc, s)) + cnt++; + coverage.push_back({ double(s), double(cnt) }); + } + } + else + { + // Small ROI: shifting grids. n_off grid origins per axis ({0, s/2} for 2 -> 4 positions); + // box index = coord >> log2(s) since box sizes are powers of two. + auto xmin = r.aabb.get_xmin(); + auto ymin = r.aabb.get_ymin(); + int n_off = 2; + int shift = 0; + for (int t = paddedSide; t > 1; t >>= 1) + shift++; + std::vector occ; + for (int s = paddedSide; s > 1; s >>= 1, --shift) + { + int span = (paddedSide >> shift) + 2; // boxes per axis (+slack for the origin shift) + occ.assign((size_t)span * span, 0); + int best = -1; + for (int oyi = 0; oyi < n_off; oyi++) + for (int oxi = 0; oxi < n_off; oxi++) + { + int ox = (oxi * s) / n_off, oy = (oyi * s) / n_off; + int cnt = 0; + for (const Pixel2& p : px) + { + int col = (int)(p.x - xmin + ox) >> shift; + int row = (int)(p.y - ymin + oy) >> shift; + size_t idx = (size_t)row * span + col; + if (!occ[idx]) { occ[idx] = 1; cnt++; } + } + if (best < 0 || cnt < best) + best = cnt; + if (oyi + 1 < n_off || oxi + 1 < n_off) + std::fill(occ.begin(), occ.end(), (char)0); // reset for next origin + } + coverage.push_back({ double(s), double(best) }); + } } - // Debug - // print_curve(boxCoverage, "boxCoverage N vs R"); - - box_count_fd = -calc_lyapunov_slope(coverage); + // Box-counting dimension = -slope of log(count) vs log(box size) + box_count_fd = -loglog_slope(coverage); } void FractalDimensionFeature::calculate_boxcount_fdim_oversized (LR & r) { Power2PaddedImageMatrix_NT pim ("pim", r.label, r.raw_pixels_NT, r.aabb, 1, 1.0); // image matrix of the mask - // Square box coverage statistics - std::vector> coverage; + // Square box coverage statistics (oversized/streaming path), same method as the trivial case. + std::vector> coverage; int s = pim.get_width(); // square tile size for (; s > 1; s /= 2) { int cnt = 0; - size_t areaCov = 0; int n_tiles = pim.get_width() / s; for (int r = 0; r < n_tiles; r++) for (int c = 0; c < n_tiles; c++) { - // check the tile for coerage + // check the tile for coverage if (pim.tile_contains_signal(r, c, s)) cnt++; } - coverage.push_back({ s, cnt }); + coverage.push_back({ double(s), double(cnt) }); } - box_count_fd = -calc_lyapunov_slope(coverage); + box_count_fd = -loglog_slope(coverage); } void FractalDimensionFeature::calculate_perimeter_fdim (LR& r) { - std::vector> coverage; - std::vector K; r.merge_multicontour(K); auto conLen = K.size(); + if (conLen < 3) + { + perim_fd = 0.; + return; + } + + // Richardson "structured walk" (divider) method on the CLOSED contour K. + // For each stride s, walk the contour in steps of s vertices, summing Euclidean chord + // lengths and closing the loop back to the start. Record x = mean ruler (chord) length + // and y = perimeter estimate. Since log(perimeter) ~ (1 - D) log(ruler), D = 1 - slope. + // Validated against the analytic Koch snowflake (D = log4/log3 = 1.2619) and a circle (D = 1). + std::vector> coverage; for (size_t s = conLen / 4; s > 0; s /= 2) { - // calculate the s-approximated perimeter - double p = 0; - for (size_t i = s; i < conLen; i += s) - { - auto& px1 = K [i-s], - px2 = K [i]; - double dist = std::sqrt(px1.sqdist(px2)); - p += dist; - } - // --calculate the last segment separately - auto tail = conLen % s; - if (tail) + double perim = 0.; + size_t nsteps = 0, i = 0; + for (; i + s < conLen; i += s) { - auto& px1 = K [conLen - tail], - px2 = K [0]; - double dist = px1.sqdist(px2); - p += dist; + perim += std::sqrt(K[i].sqdist(K[i + s])); + nsteps++; } + // close the loop from the last visited vertex back to the start + perim += std::sqrt(K[i].sqdist(K[0])); + nsteps++; - // save this approximation - coverage.push_back({ s, (int)p }); + double ruler = perim / double(nsteps); // mean chord (ruler) length at this stride + coverage.push_back({ ruler, perim }); } - // Richardson divider method: log(perimeter) vs log(ruler) has slope (1 - D), so D = 1 - slope - perim_fd = 1.0 - calc_lyapunov_slope(coverage); + + // Richardson divider dimension: D = 1 - slope of log(perimeter) vs log(ruler) + perim_fd = 1.0 - loglog_slope(coverage); } -double FractalDimensionFeature::calc_lyapunov_slope (const std::vector> & coverage) +double FractalDimensionFeature::loglog_slope (const std::vector> & coverage) { - // Do we have any data? - if (coverage.size() == 0) - return 0.; - - // Skip tiny ROIs - if (coverage.size() < 2) - return 0.; - - // Post-process towards local gradients (in the form of Lyapunov exponents) - std::vector X, Y; - for (int i = 0; i < coverage.size(); i++) + // Least-squares slope of log(y) vs log(x) over the (scale, measure) pairs. + // Box counting passes (box size, count) -> D = -slope; the divider method passes + // (ruler length, perimeter) -> D = 1 - slope. A least-squares fit over all scales is + // the standard, robust estimator (more stable than a two-endpoint slope on noisy ROIs). + double sx = 0, sy = 0, sxy = 0, sx2 = 0; + size_t used = 0; + for (const auto& pr : coverage) { - X.push_back (std::log(coverage[i].first)); - Y.push_back (std::log(coverage[i].second)); + if (pr.first <= 0. || pr.second <= 0.) + continue; // log undefined; skip degenerate points (count/perimeter 0) + double x = std::log(pr.first), y = std::log(pr.second); + sx += x; sy += y; sxy += x * y; sx2 += x * x; + used++; } - - // Gradients - std::vector Dn, Dr; - for (int i = 1; i < Y.size(); i++) - { - auto dn = Y[i] - Y[i - 1]; - Dn.push_back(dn); - auto dr = X[i] - X[i - 1]; - Dr.push_back(dr); - } - - // Lyapunov exponents - std::vector Lambda; - for (int i = 0; i < Dn.size(); i++) - { - auto lambda = Dn[i] / Dr[i]; - Lambda.push_back(lambda); - } - - // The box-counting / Richardson dimension is the (average) slope of log(count) vs log(scale), - // i.e. the MEAN of the local slopes Lambda[]. The previous code returned the least-squares - // slope of Lambda[] *against its index* (the rate-of-change of the slope), which is ~0 for a - // clean power law -> dimension came out ~0 (FRACT_DIM_BOXCOUNT=-0.07). Return the mean slope. - if (Lambda.empty()) + if (used < 2) + return 0.; + double denom = sx2 * double(used) - sx * sx; + if (denom == 0.) return 0.; - double sum = 0.; - for (double v : Lambda) - sum += v; - return sum / double(Lambda.size()); + return (sxy * double(used) - sx * sy) / denom; } void FractalDimensionFeature::osized_add_online_pixel(size_t x, size_t y, uint32_t intensity) diff --git a/src/nyx/features/fractal_dim.h b/src/nyx/features/fractal_dim.h index 00266609..fb036bac 100644 --- a/src/nyx/features/fractal_dim.h +++ b/src/nyx/features/fractal_dim.h @@ -34,6 +34,6 @@ class FractalDimensionFeature: public FeatureMethod void calculate_boxcount_fdim (LR& r); void calculate_boxcount_fdim_oversized (LR& r); void calculate_perimeter_fdim (LR& r); - double calc_lyapunov_slope (const std::vector>& coverage_stats); + double loglog_slope (const std::vector>& coverage); double box_count_fd = 0, perim_fd = 0; }; \ No newline at end of file diff --git a/src/nyx/features/image_matrix.h b/src/nyx/features/image_matrix.h index fc59b613..43f821de 100644 --- a/src/nyx/features/image_matrix.h +++ b/src/nyx/features/image_matrix.h @@ -392,14 +392,16 @@ class Power2PaddedImageMatrix : public ImageMatrix // Cache AABB original_aabb = aabb; - // Figure out the padded size and allocate + // Tight power-of-2 canvas (ceil_pow2, not closest_pow2 which is strictly greater and + // wastes an octave) with the ROI at the grid origin. Box counting is registration + // sensitive; centering on an over-sized canvas misaligned the ROI with the coarse boxes + // and biased the dimension low (a filled square read 1.75 instead of 2.0). int bigSide = std::max(aabb.get_width(), aabb.get_height()); - StatsInt paddedSide = Nyxus::closest_pow2 (bigSide); + StatsInt paddedSide = Nyxus::ceil_pow2 (bigSide); allocate (paddedSide, paddedSide); - // Copy pixels - int padOffsetX = (paddedSide - original_aabb.get_width()) / 2; - int padOffsetY = (paddedSide - original_aabb.get_height()) / 2; + int padOffsetX = 0; + int padOffsetY = 0; for (auto& pxl : labels_raw_pixels) { diff --git a/src/nyx/features/image_matrix_nontriv.cpp b/src/nyx/features/image_matrix_nontriv.cpp index 3209030d..2a98038c 100644 --- a/src/nyx/features/image_matrix_nontriv.cpp +++ b/src/nyx/features/image_matrix_nontriv.cpp @@ -435,14 +435,16 @@ Power2PaddedImageMatrix_NT::Power2PaddedImageMatrix_NT (const std::string& _name // Cache AABB original_aabb = aabb; - // Figure out the padded size and allocate + // Figure out the padded size and allocate. Tight power-of-2 canvas (ceil_pow2) and + // origin alignment, matching the trivial Power2PaddedImageMatrix, so the oversized + // box-counting path recovers the same grid-aligned dimension. int bigSide = std::max(aabb.get_width(), aabb.get_height()); - StatsInt paddedSide = Nyxus::closest_pow2(bigSide); + StatsInt paddedSide = Nyxus::ceil_pow2(bigSide); allocate(paddedSide, paddedSide, 0.0); - // Copy pixels - int padOffsetX = (paddedSide - original_aabb.get_width()) / 2; - int padOffsetY = (paddedSide - original_aabb.get_height()) / 2; + // Place the ROI at the grid origin (no centering) to avoid box-grid misalignment bias. + int padOffsetX = 0; + int padOffsetY = 0; for (auto pxl : raw_pixels) { diff --git a/src/nyx/helpers/helpers.h b/src/nyx/helpers/helpers.h index 55ae1e6c..e0a2b189 100644 --- a/src/nyx/helpers/helpers.h +++ b/src/nyx/helpers/helpers.h @@ -247,6 +247,17 @@ namespace Nyxus return retval; } + // Smallest power of two >= a (a tight ceiling, unlike closest_pow2 which is + // strictly greater than a). Used to pad an ROI to a box-counting grid without + // wasting an octave: ceil_pow2(256)=256, ceil_pow2(13)=16, ceil_pow2(1)=1. + inline int ceil_pow2(const int a) + { + int p = 1; + while (p < a) + p <<= 1; + return p; + } + inline void print_curve(const std::vector>& curve, const std::string& name) { std::cout << "\n\n" << name << " = [\n"; diff --git a/tests/python/test_fractal_dim_oracle.py b/tests/python/test_fractal_dim_oracle.py index 5bfeb766..9c386a2c 100644 --- a/tests/python/test_fractal_dim_oracle.py +++ b/tests/python/test_fractal_dim_oracle.py @@ -1,10 +1,27 @@ -"""Regression / bug-exposure tests for 2D feature defects found during oracle validation -(2026-06). These exercise the PRODUCTION featurize() path on ROIs *with background* and at -the DEFAULT settings - the conditions the C++ unit tests miss. +"""Oracle validation of the fractal-dimension features on the PRODUCTION featurize() path. -This module covers the fractal-dimension defect (proposed fix #8): FRACT_DIM_BOXCOUNT in [1, 2]. +Two kinds of test: + +1. ANALYTIC shapes with a KNOWN, convention-independent fractal dimension. These are the + real oracle: nyxus must recover the known value. Cross-checked independently against the + ImageJ / FracLac shifting-grid box-count oracle (computed offline; values baked in below). + box-count (filled region): square -> 2.0, line -> 1.0, Sierpinski triangle -> log2(3)=1.585 + perimeter (divider/Richardson on the contour): disk -> 1.0, Koch snowflake -> log4/log3=1.262 + +2. ARBITRARY ROI (the irregular 154-px region the C++ unit tests use) exercised with background + at default settings, validated against offline ImageJ shifting-grid oracles: box-count vs the + box count of the filled ROI (same method, tight), and the divider perimeter vs the box count + of the ROI's edge (different method, same boundary dimension - cross-method convergent validity). + +Background: nyxus' box-count padded to a power of two *strictly larger* than the ROI and +*centered* it, misaligning it with the box grid and biasing the dimension low (a filled square +read 1.75 instead of 2.0). The fix aligns the ROI to the grid origin on a tight power-of-two +canvas, fits the dimension by least squares, and auto-switches to shifting grids (min box count +over grid origins, as FracLac does) on small ROIs where registration bias matters; the perimeter +path was rewritten as a clean closed-contour divider. Validated here against analytic ground truth. """ import re +import math from pathlib import Path import numpy as np import pytest @@ -12,66 +29,154 @@ # ----------------------------- helpers -------------------------------------- -def _canonical_roi(): - """The pixelIntensityFeaturesTestData ROI from tests/test_data.h - the same irregular - 154-px region the C++ tests use, which reproduces the shape/moment defects. Falls back to - an L-shape if the header can't be read.""" - hdr = Path(__file__).resolve().parent.parent / "test_data.h" - try: - txt = hdr.read_text(encoding="utf-8", errors="replace") - body = re.search(r"pixelIntensityFeaturesTestData\[\]\s*=\s*\{(.*?)\};", txt, re.S).group(1) - pts = [(int(x), int(y), int(v)) for x, y, v in - re.findall(r"\{\s*(\d+)\s*,\s*(\d+)\s*,\s*(\d+)\s*\}", body)] - W = max(p[0] for p in pts) + 2 - H = max(p[1] for p in pts) + 2 - inten = np.zeros((H, W), np.uint32) - label = np.zeros((H, W), np.uint32) - for x, y, v in pts: - inten[y, x] = v - label[y, x] = 1 - return inten, label - except Exception: - return None - - def _run(features, inten, label, **kw): nyx = nyxus.Nyxus(features=features, n_feature_calc_threads=1, **kw) df = nyx.featurize(inten.astype(np.float64), label.astype(np.uint32), intensity_names=["i"], label_names=["l"]) - return df # one row per label + return df + + +def _fd(label): + """Run nyxus on a binary mask and return (box-count D, perimeter D).""" + inten = label.astype(np.float64) * 100.0 + row = _run(["*ALL_MORPHOLOGY*"], inten, label).iloc[0] + bc = [c for c in row.index if c.endswith("FRACT_DIM_BOXCOUNT")][0] + pf = [c for c in row.index if c.endswith("FRACT_DIM_PERIMETER")][0] + return float(row[bc]), float(row[pf]) + + +# ----------------------------- analytic shapes ------------------------------ +def _square(n=256): + return np.ones((n, n), np.uint32) + + +def _line(n=256): + m = np.zeros((n, n), np.uint32) + m[n // 2, :] = 1 + return m + + +def _sierpinski_triangle(order=8): + """(x & y) == 0 bit pattern -> Sierpinski triangle, box-count D = log2(3) = 1.585.""" + n = 1 << order + rows = np.arange(n) + m = np.zeros((n, n), np.uint32) + for y in range(n): + m[y] = ((rows & y) == 0).astype(np.uint32) + return m -def _one(features, inten, label, **kw): - return _run(features, inten, label, **kw).iloc[0] +def _disk(r=120, pad=12): + """Filled disk: smooth boundary -> perimeter (divider) dimension -> 1.0.""" + n = 2 * (r + pad) + yy, xx = np.ogrid[:n, :n] + c = r + pad + return (((xx - c) ** 2 + (yy - c) ** 2) <= r * r).astype(np.uint32) -def _blob(): - """Irregular ROI for shape/moment bug exposure. Prefers the canonical 154-px ROI; else an L-shape.""" - c = _canonical_roi() - if c is not None: - return c - H, W = 16, 16 - inten = np.zeros((H, W), np.uint32) +def _koch_snowflake(depth=5, size=1400): + """Filled Koch snowflake (pure numpy, no external deps): boundary divider + dimension = log4/log3 = 1.2619.""" + def koch(p1, p2, d): + if d == 0: + return [p1] + (x1, y1), (x2, y2) = p1, p2 + ax, ay = x1 + (x2 - x1) / 3.0, y1 + (y2 - y1) / 3.0 + cx, cy = x1 + 2 * (x2 - x1) / 3.0, y1 + 2 * (y2 - y1) / 3.0 + ang = math.atan2(y2 - y1, x2 - x1) - math.pi / 3.0 + L = math.hypot(x2 - x1, y2 - y1) / 3.0 + px, py = ax + L * math.cos(ang), ay + L * math.sin(ang) + return (koch(p1, (ax, ay), d - 1) + koch((ax, ay), (px, py), d - 1) + + koch((px, py), (cx, cy), d - 1) + koch((cx, cy), p2, d - 1)) + + m = size * 0.1 + s = size - 2 * m + h = s * math.sqrt(3) / 2.0 + A = (m, m + h * 2 / 3) + B = (m + s, m + h * 2 / 3) + C = (m + s / 2, m + h * 2 / 3 - h) + poly = (koch(A, B, depth) + koch(B, C, depth) + koch(C, A, depth)) + # even-odd scanline polygon fill (numpy only) + P = np.array(poly, float) + x0, y0 = P[:, 0], P[:, 1] + x1, y1 = np.roll(x0, -1), np.roll(y0, -1) + out = np.zeros((size, size), np.uint32) + for row in range(size): + cond = ((y0 <= row) & (y1 > row)) | ((y1 <= row) & (y0 > row)) + idx = np.where(cond)[0] + if idx.size == 0: + continue + tt = (row - y0[idx]) / (y1[idx] - y0[idx]) + xs = np.sort(x0[idx] + tt * (x1[idx] - x0[idx])) + for i in range(0, len(xs) - 1, 2): + a = int(np.ceil(xs[i])) + b = int(np.floor(xs[i + 1])) + if b >= a: + out[row, max(0, a):min(size, b + 1)] = 1 + return out + + +# ============================ box-count oracle ============================== +@pytest.mark.parametrize("name,mask,truth", [ + ("filled_square", _square(), 2.0), + ("straight_line", _line(), 1.0), + ("sierpinski_triangle", _sierpinski_triangle(), math.log(3) / math.log(2)), +]) +def test_boxcount_recovers_known_dimension(name, mask, truth): + """nyxus FRACT_DIM_BOXCOUNT must recover the analytic box-counting dimension + (matches the ImageJ/FracLac oracle). Before the grid-alignment fix the square read + ~1.75 and the Sierpinski triangle ~1.39; both now land on truth.""" + bc, _ = _fd(mask) + assert abs(bc - truth) < 0.1, f"{name}: box-count {bc:.3f} vs truth {truth:.3f}" + + +# ============================ perimeter oracle ============================== +def test_perimeter_disk_is_smooth(): + """A disk has a smooth (non-fractal) boundary -> divider dimension -> 1.0.""" + _, pf = _fd(_disk()) + assert abs(pf - 1.0) < 0.05, f"disk perimeter D {pf:.3f} vs 1.0" + + +def test_perimeter_koch_snowflake(): + """Koch snowflake boundary: divider dimension = log4/log3 = 1.2619.""" + truth = math.log(4) / math.log(3) + _, pf = _fd(_koch_snowflake()) + assert abs(pf - truth) < 0.05, f"Koch perimeter D {pf:.3f} vs truth {truth:.3f}" + + +# ============================ arbitrary ROI ================================ +def _canonical_roi(): + """The irregular 154-px pixelIntensityFeaturesTestData ROI from tests/test_data.h.""" + hdr = Path(__file__).resolve().parent.parent / "test_data.h" + txt = hdr.read_text(encoding="utf-8", errors="replace") + body = re.search(r"pixelIntensityFeaturesTestData\[\]\s*=\s*\{(.*?)\};", txt, re.S).group(1) + pts = [(int(x), int(y)) for x, y, v in + re.findall(r"\{\s*(\d+)\s*,\s*(\d+)\s*,\s*(\d+)\s*\}", body)] + W = max(p[0] for p in pts) + 2 + H = max(p[1] for p in pts) + 2 label = np.zeros((H, W), np.uint32) - # L-shape: a vertical bar + a horizontal foot (clearly non-convex) - for y in range(3, 13): - for x in range(3, 6): - label[y, x] = 1 - for y in range(10, 13): - for x in range(6, 13): - label[y, x] = 1 - ys, xs = np.nonzero(label) - inten[ys, xs] = (100 + 10 * xs + 7 * ys).astype(np.uint32) - return inten, label - - -# ============================ FRACTAL DIMENSION ============================= -def test_fractal_dimension_in_range(): - """Proposed fix #8 (fractal dim: mean log-log slope, not slope-of-slopes): the box-counting - fractal dimension is the MEAN of the local log(count)-vs-log(scale) slopes. The old code - returned the least-squares slope of those slopes *against their index* (~0 for a clean power - law), so FRACT_DIM_BOXCOUNT came out ~-0.07..-0.83 - outside the valid [1, 2] range for a 2D - shape. (FRACT_DIM_PERIMETER similarly needed the Richardson D = 1 - slope convention.)""" - inten, label = _blob() - row = _one(["*ALL_MORPHOLOGY*"], inten, label) - assert 1.0 <= float(row["FRACT_DIM_BOXCOUNT"]) <= 2.0 + for x, y in pts: + label[y, x] = 1 + return label + + +# Offline ImageJ shifting-grid oracle values for the irregular 154-px ROI: +# box-count = box count of the FILLED ROI - SAME method as nyxus, so it matches tightly -> 1.389 +# perimeter = box count of the ROI's EDGE (boundary/Minkowski dimension). This is a DIFFERENT +# algorithm than nyxus' Richardson divider, but estimates the SAME boundary-complexity +# dimension (box dim == compass/divider dim for boundary curves). The two agree to +# 0.06 here - cross-method convergent validity -> 1.163. +ORACLE_BOXCOUNT_154 = 1.389 # same-method oracle; nyxus computes 1.389 (<0.001) +ORACLE_PERIMETER_154_EDGE = 1.163 # cross-method (edge box-count) oracle; nyxus divider computes 1.101 + + +def test_arbitrary_roi_matches_oracle(): + """On the irregular 154-px ROI (with background, default settings): + - box-count matches the offline ImageJ shifting-grid oracle (same method), tightly; + - the divider perimeter matches an independent box-count-of-edge oracle within the + cross-method tolerance (different algorithm, same boundary dimension).""" + bc, pf = _fd(_canonical_roi()) + assert abs(bc - ORACLE_BOXCOUNT_154) < 0.02, \ + f"box-count {bc:.4f} vs same-method shifting-grid oracle {ORACLE_BOXCOUNT_154}" + assert abs(pf - ORACLE_PERIMETER_154_EDGE) < 0.10, \ + f"perimeter {pf:.4f} vs cross-method edge-box-count oracle {ORACLE_PERIMETER_154_EDGE}" diff --git a/tests/test_shape_morphology_2d.h b/tests/test_shape_morphology_2d.h index e9bec49b..6c1a1d8d 100644 --- a/tests/test_shape_morphology_2d.h +++ b/tests/test_shape_morphology_2d.h @@ -63,6 +63,10 @@ static std::unordered_map unvetted_nyxus_regression_shape2d {"ROI_RADIUS_MEAN", 1.07692307692308}, {"ROI_RADIUS_MAX", 4.0}, {"ROI_RADIUS_MEDIAN", 1.0}, + // FRACT_DIM_PERIMETER: closed-contour Richardson divider (was an out-of-range 0.319). No + // software implements the divider method, so this is a regression value; the method is + // oracle-validated on the Koch snowflake (D=log4/log3) in tests/python/test_fractal_dim_oracle.py. + {"FRACT_DIM_PERIMETER", 1.2632300929080915}, }; static std::unordered_map oracle_3p_shape2d_feature_golden_values{ @@ -79,8 +83,10 @@ static std::unordered_map oracle_3p_shape2d_feature_golden_ {"CONVEX_HULL_AREA", 28.0}, {"SOLIDITY", 0.9285714285714286}, {"DIAMETER_EQUAL_PERIMETER", 8.57365809435587}, - {"FRACT_DIM_BOXCOUNT", 1.5849625007211565}, // FIX (fractal_dim.cpp): mean log-log slope = log2(3); old -0.83 was slope-of-slopes (~0) - {"FRACT_DIM_PERIMETER", 0.3187149603076458}, // FIX: Richardson D = 1 - slope; old -1.97 was the raw slope + // Box-counting dimension, reproduces the ImageJ/FracLac shifting-grid box counter to full + // precision on this ROI (== log2(3)); the method also recovers analytic shapes exactly + // (square->2, Sierpinski->1.585) in tests/python/test_fractal_dim_oracle.py. + {"FRACT_DIM_BOXCOUNT", 1.5849625007211565}, {"DIAMETER_CIRCUMSCRIBING_CIRCLE", 12.3317073399088}, {"DIAMETER_INSCRIBING_CIRCLE", 0.828486893405308}, {"GEODETIC_LENGTH", 10.0}, @@ -307,8 +313,10 @@ void test_shape2d_verifiable_with_3p_builtin_oracle_fractal_circle_features() std::vector> fvals; calculate_shape2d_feature_values(fvals); + // Box-count matches the ImageJ/FracLac shifting-grid oracle; perimeter (divider) has no + // software oracle, so it is a regression value (method oracle-validated in the python suite). assert_verifiable_with_3p_builtin_oracle_shape2d_feature(fvals, Nyxus::Feature2D::FRACT_DIM_BOXCOUNT, "FRACT_DIM_BOXCOUNT"); - assert_verifiable_with_3p_builtin_oracle_shape2d_feature(fvals, Nyxus::Feature2D::FRACT_DIM_PERIMETER, "FRACT_DIM_PERIMETER"); + assert_unvetted_no_direct_oracle_shape2d_feature(fvals, Nyxus::Feature2D::FRACT_DIM_PERIMETER, "FRACT_DIM_PERIMETER"); assert_verifiable_with_3p_builtin_oracle_shape2d_feature(fvals, Nyxus::Feature2D::DIAMETER_CIRCUMSCRIBING_CIRCLE, "DIAMETER_CIRCUMSCRIBING_CIRCLE"); assert_verifiable_with_3p_builtin_oracle_shape2d_feature(fvals, Nyxus::Feature2D::DIAMETER_INSCRIBING_CIRCLE, "DIAMETER_INSCRIBING_CIRCLE"); }