Skip to content

test(moments): assert the normalized raw moments against scikit-image #430

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
darkclad wants to merge 3 commits into
base: main
Choose a base branch
from
Open

test(moments): assert the normalized raw moments against scikit-image #430

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1cd0338
test(moments): assert the normalized raw moments against scikit-image
Aug 14, 2026
b746ac1
test(vetting): re-verify every pinned intensity-histogram golden, not…
Aug 14, 2026
efd2818
test(vetting): re-verify every pinned GLRLM golden, not a subset
Aug 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 5 additions & 6 deletions tests/test_2d_moments_common.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@

#include <algorithm>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>

Expand Down Expand Up @@ -53,9 +52,9 @@ static void load_geomoment_fixture(LR& roidata)
{
const auto intensity = static_cast<PixIntens>(geomoment_fixture_intensity(x, y));
if (roidata.aux_area == 0)
init_label_record_3(roidata, x, y, intensity);
Nyxus::init_label_record_3(roidata, x, y, intensity);
else
update_label_record_3(roidata, x, y, intensity);
Nyxus::update_label_record_3(roidata, x, y, intensity);
}
}

Expand Down Expand Up @@ -100,7 +99,7 @@ static void calculate_2d_geomoment_feature_values(std::vector<std::vector<double

// Thin right wedge: pixels (x,y) with 0<=x<40, 0<=y<8, 5*y <= x. Elongated AND skewed, so the
// odd-order normalized central moments are large (eta30 ~ -0.23) and the Hu invariants h5/h6
// discriminate the fixed calcHu_imp formulas from the historic defective ones by ~440x / ~1900x
// discriminate a correct h5/h6 from the two classic formula slips by ~440x / ~1900x
// the assertion tolerance -- the symmetric 48x40 rectangle fixture above cannot catch an h5/h6
// regression at all (its odd etas vanish). Goldens: tests/vetting/oracles/gen_moments_skimage.py.
static void load_wedge_fixture(LR& roidata)
Expand All @@ -114,9 +113,9 @@ static void load_wedge_fixture(LR& roidata)
if (5 * y > x)
continue;
if (roidata.aux_area == 0)
init_label_record_3(roidata, x, y, 1);
Nyxus::init_label_record_3(roidata, x, y, 1);
else
update_label_record_3(roidata, x, y, 1);
Nyxus::update_label_record_3(roidata, x, y, 1);
}

for (int y = 0; y < height; y++)
Expand Down
8 changes: 4 additions & 4 deletions tests/test_2d_moments_regression.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,8 @@ static const ref_vals_list<GeomomentGoldenValue> moments_2d_regression_shape_ref
{Nyxus::Feature2D::WEIGHTED_HU_M2, "WEIGHTED_HU_M2", 0.00024306185106698786},
{Nyxus::Feature2D::WEIGHTED_HU_M3, "WEIGHTED_HU_M3", 1.1262214820995351e-05},
{Nyxus::Feature2D::WEIGHTED_HU_M4, "WEIGHTED_HU_M4", 7.674925625409561e-05},
{Nyxus::Feature2D::WEIGHTED_HU_M5, "WEIGHTED_HU_M5", -1.8387946856938181e-09}, // FIX(h5): snapshot after the calcHu_imp h5 fix (skimage moments_hu on the pinned WT_NORM_CTR_MOM_*; gen_moments_skimage.py)
{Nyxus::Feature2D::WEIGHTED_HU_M6, "WEIGHTED_HU_M6", 2.1453724299933224e-07}, // FIX(h6): snapshot after the calcHu_imp h6 fix; old pin 0.004652... == WT_NORM_CTR_MOM_03, i.e. the stray "+eta03" of the precedence bug (gen_moments_skimage.py)
{Nyxus::Feature2D::WEIGHTED_HU_M5, "WEIGHTED_HU_M5", -1.8387946856938181e-09}, // skimage moments_hu over the pinned WT_NORM_CTR_MOM_* (gen_moments_skimage.py)
{Nyxus::Feature2D::WEIGHTED_HU_M6, "WEIGHTED_HU_M6", 2.1453724299933224e-07}, // distinct from WT_NORM_CTR_MOM_03; equality there means h6 lost the (eta21+eta03) grouping
{Nyxus::Feature2D::WEIGHTED_HU_M7, "WEIGHTED_HU_M7", -1.3078000499389243e-09},
};

Expand Down Expand Up @@ -98,8 +98,8 @@ static const ref_vals_list<GeomomentGoldenValue> moments_2d_regression_intensity
{Nyxus::Feature2D::IMOM_WHU2, "IMOM_WHU2", 5.6161730111679807e-09},
{Nyxus::Feature2D::IMOM_WHU3, "IMOM_WHU3", 1.1923046864432925e-12},
{Nyxus::Feature2D::IMOM_WHU4, "IMOM_WHU4", 3.1287168309169019e-12},
{Nyxus::Feature2D::IMOM_WHU5, "IMOM_WHU5", -3.578057915124686e-25}, // FIX(h5): snapshot after the calcHu_imp h5 fix (skimage moments_hu on the pinned IMOM_WNCM_*; gen_moments_skimage.py)
{Nyxus::Feature2D::IMOM_WHU6, "IMOM_WHU6", 3.297985694277869e-17}, // FIX(h6): snapshot after the calcHu_imp h6 fix; old pin 4.6265447915851515e-07 == IMOM_WNCM_03, i.e. the stray "+eta03" of the precedence bug (gen_moments_skimage.py)
{Nyxus::Feature2D::IMOM_WHU5, "IMOM_WHU5", -3.578057915124686e-25}, // skimage moments_hu over the pinned IMOM_WNCM_* (gen_moments_skimage.py)
{Nyxus::Feature2D::IMOM_WHU6, "IMOM_WHU6", 3.297985694277869e-17}, // distinct from IMOM_WNCM_03; equality there means h6 lost the (eta21+eta03) grouping
{Nyxus::Feature2D::IMOM_WHU7, "IMOM_WHU7", -5.6202519491182231e-24},
};

Expand Down
75 changes: 68 additions & 7 deletions tests/test_2d_moments_skimage.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ static const ref_vals_list<GeomomentGoldenValue> moments_2d_skimage_shape_ref_va
{Nyxus::Feature2D::HU_M2, "HU_M2", 0.0009336419753086421},
{Nyxus::Feature2D::HU_M3, "HU_M3", 0},
{Nyxus::Feature2D::HU_M4, "HU_M4", 0},
{Nyxus::Feature2D::HU_M5, "HU_M5", 0}, // FIX(h5): oracle-exact 0 on this symmetric fixture (odd etas vanish); the old pin 4.598e-10 was summation noise of the defective h5 formula (gen_moments_skimage.py)
{Nyxus::Feature2D::HU_M5, "HU_M5", 0}, // exactly 0: the odd etas vanish on this symmetric fixture
{Nyxus::Feature2D::HU_M6, "HU_M6", 0},
{Nyxus::Feature2D::HU_M7, "HU_M7", -2.3604559630908652e-10},
};
Expand Down Expand Up @@ -90,8 +90,8 @@ static const ref_vals_list<GeomomentGoldenValue> moments_2d_skimage_intensity_re
{Nyxus::Feature2D::IMOM_HU2, "IMOM_HU2", 3.7112807351259333e-08},
{Nyxus::Feature2D::IMOM_HU3, "IMOM_HU3", 2.6788774435431954e-11},
{Nyxus::Feature2D::IMOM_HU4, "IMOM_HU4", 2.8991603200922661e-12},
{Nyxus::Feature2D::IMOM_HU5, "IMOM_HU5", -1.9898126265836865e-22}, // FIX(h5): skimage-oracle value after the calcHu_imp h5 fix; old pin encoded the 9x-bracket defect (gen_moments_skimage.py)
{Nyxus::Feature2D::IMOM_HU6, "IMOM_HU6", -6.66827346717698e-17}, // FIX(h6): skimage-oracle value after the calcHu_imp h6 fix; old pin -2.3976723312959013e-06 == IMOM_NCM_03, i.e. the stray "+eta03" of the precedence bug (gen_moments_skimage.py)
{Nyxus::Feature2D::IMOM_HU5, "IMOM_HU5", -1.9898126265836865e-22},
{Nyxus::Feature2D::IMOM_HU6, "IMOM_HU6", -6.66827346717698e-17}, // distinct from IMOM_NCM_03; a value equal to it means h6 lost the (eta21+eta03) grouping
{Nyxus::Feature2D::IMOM_HU7, "IMOM_HU7", 2.3835594243218353e-23},
};

Expand Down Expand Up @@ -121,13 +121,15 @@ static const ref_vals_list<GeomomentGoldenValue> moments_2d_skimage_wedge_hu_ref
{Nyxus::Feature2D::HU_M2, "HU_M2", 0.2921768785246152},
{Nyxus::Feature2D::HU_M3, "HU_M3", 0.06341348298776381},
{Nyxus::Feature2D::HU_M4, "HU_M4", 0.05042335332144146},
{Nyxus::Feature2D::HU_M5, "HU_M5", 0.002851264767850148}, // defective h5 gives 0.0032935269006466807 (442x tolerance away)
{Nyxus::Feature2D::HU_M6, "HU_M6", 0.027252861297278195}, // defective h6 gives 0.02916606310377036 (1913x tolerance away)
{Nyxus::Feature2D::HU_M5, "HU_M5", 0.002851264767850148}, // a 3*(eta30+eta12)^2 -> (3*(eta30+eta12))^2 slip gives 0.0032935269006466807, 442x the tolerance
{Nyxus::Feature2D::HU_M6, "HU_M6", 0.027252861297278195}, // dropping the (eta21+eta03) grouping gives 0.02916606310377036, 1913x the tolerance
{Nyxus::Feature2D::HU_M7, "HU_M7", 5.616461607732461e-06},
};

// Vets Hu h5/h6 (and the full moment chain) against scikit-image on an asymmetric fixture --
// the historic calcHu_imp h5 9x-bracket / h6 "+eta03"-precedence defects fail this assertion.
// Vets Hu h5/h6 (and the full moment chain) against scikit-image on an asymmetric fixture, where
// the odd-order etas are large enough that a mis-bracketed h5 or an h6 missing the (eta21+eta03)
// grouping misses by hundreds of times the tolerance. The symmetric rectangle cannot: its odd etas
// vanish, so both formulas agree there.
void test_2d_moments_hu_wedge_skimage()
{
std::vector<std::vector<double>> fvals;
Expand All @@ -148,3 +150,62 @@ void test_2d_moments_intensity_skimage()
calculate_2d_geomoment_feature_values(fvals);
assert_2d_geomoment_features(fvals, moments_2d_skimage_intensity_ref_vals, "VERIFIABLE_WITH_3P_BUILTIN_ORACLE__");
}

// Normalized RAW moments: NORM_SPAT_MOMENT_pq = m_pq / m_00^((p+q)/2 + 1), and IMOM_NRM_pq the
// same over the intensity-weighted raw moments. skimage has no native function for this - its
// moments_normalized() applies that exponent to the CENTRAL moments, which is a different
// quantity Nyxus exposes separately as NORM_CENTRAL_MOMENT_pq / IMOM_NCM_pq. So the oracle here
// is skimage's raw moments carried through skimage's own normalization exponent; the generator
// asserts the two are distinct so a future edit cannot quietly swap one for the other.
// Goldens: tests/vetting/oracles/gen_moments_skimage.py section D.
static const ref_vals_list<GeomomentGoldenValue> moments_2d_skimage_normraw_shape_ref_vals{
{Nyxus::Feature2D::NORM_SPAT_MOMENT_00, "NORM_SPAT_MOMENT_00", 1.0},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_01, "NORM_SPAT_MOMENT_01", 0.4450245779729475},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_02, "NORM_SPAT_MOMENT_02", 0.2674479166666667},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_03, "NORM_SPAT_MOMENT_03", 0.1807912348015099},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_10, "NORM_SPAT_MOMENT_10", 0.5363116708904752},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_11, "NORM_SPAT_MOMENT_11", 0.238671875},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_12, "NORM_SPAT_MOMENT_12", 0.14343543906367653},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_13, "NORM_SPAT_MOMENT_13", 0.09696044921875},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_20, "NORM_SPAT_MOMENT_20", 0.3875868055555556},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_21, "NORM_SPAT_MOMENT_21", 0.17248565457024395},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_22, "NORM_SPAT_MOMENT_22", 0.10365928367332176},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_23, "NORM_SPAT_MOMENT_23", 0.07007229716916162},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_30, "NORM_SPAT_MOMENT_30", 0.31508310664815414},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_31, "NORM_SPAT_MOMENT_31", 0.1402197265625},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_32, "NORM_SPAT_MOMENT_32", 0.08426832044990999},
{Nyxus::Feature2D::NORM_SPAT_MOMENT_33, "NORM_SPAT_MOMENT_33", 0.056964263916015626},
};

static const ref_vals_list<GeomomentGoldenValue> moments_2d_skimage_normraw_intensity_ref_vals{
{Nyxus::Feature2D::IMOM_NRM_00, "IMOM_NRM_00", 1.0},
{Nyxus::Feature2D::IMOM_NRM_01, "IMOM_NRM_01", 0.039397166363954135},
{Nyxus::Feature2D::IMOM_NRM_02, "IMOM_NRM_02", 0.001897046009773198},
{Nyxus::Feature2D::IMOM_NRM_03, "IMOM_NRM_03", 9.951746245055554e-05},
{Nyxus::Feature2D::IMOM_NRM_10, "IMOM_NRM_10", 0.04534163016617505},
{Nyxus::Feature2D::IMOM_NRM_11, "IMOM_NRM_11", 0.0017522720110346945},
{Nyxus::Feature2D::IMOM_NRM_12, "IMOM_NRM_12", 8.375967849944927e-05},
{Nyxus::Feature2D::IMOM_NRM_13, "IMOM_NRM_13", 4.3766925283775395e-06},
{Nyxus::Feature2D::IMOM_NRM_20, "IMOM_NRM_20", 0.002579984204506706},
{Nyxus::Feature2D::IMOM_NRM_21, "IMOM_NRM_21", 9.892567767206171e-05},
{Nyxus::Feature2D::IMOM_NRM_22, "IMOM_NRM_22", 4.7143264687233126e-06},
{Nyxus::Feature2D::IMOM_NRM_23, "IMOM_NRM_23", 2.459309328437938e-07},
{Nyxus::Feature2D::IMOM_NRM_30, "IMOM_NRM_30", 0.0001612294112519097},
{Nyxus::Feature2D::IMOM_NRM_31, "IMOM_NRM_31", 6.155292790140729e-06},
{Nyxus::Feature2D::IMOM_NRM_32, "IMOM_NRM_32", 2.928338360295916e-07},
{Nyxus::Feature2D::IMOM_NRM_33, "IMOM_NRM_33", 1.5262053460086997e-08},
};

void test_2d_moments_normraw_shape_skimage()
{
std::vector<std::vector<double>> fvals;
calculate_2d_geomoment_feature_values(fvals);
assert_2d_geomoment_features(fvals, moments_2d_skimage_normraw_shape_ref_vals, "VERIFIABLE_WITH_3P_BUILTIN_ORACLE__");
}

void test_2d_moments_normraw_intensity_skimage()
{
std::vector<std::vector<double>> fvals;
calculate_2d_geomoment_feature_values(fvals);
assert_2d_geomoment_features(fvals, moments_2d_skimage_normraw_intensity_ref_vals, "VERIFIABLE_WITH_3P_BUILTIN_ORACLE__");
}
10 changes: 10 additions & 0 deletions tests/test_all.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -992,6 +992,16 @@ TEST(TEST_NYXUS, TEST_2D_MOMENTS_SHAPE_REGRESSION)
ASSERT_NO_THROW(test_2d_moments_shape_regression());
}

TEST(TEST_NYXUS, TEST_2D_MOMENTS_NORMRAW_SHAPE_SKIMAGE)
{
ASSERT_NO_THROW(test_2d_moments_normraw_shape_skimage());
}

TEST(TEST_NYXUS, TEST_2D_MOMENTS_NORMRAW_INTENSITY_SKIMAGE)
{
ASSERT_NO_THROW(test_2d_moments_normraw_intensity_skimage());
}

TEST(TEST_NYXUS, TEST_2D_MOMENTS_INTENSITY_SKIMAGE)
{
ASSERT_NO_THROW(test_2d_moments_intensity_skimage());
Expand Down
7 changes: 7 additions & 0 deletions tests/vetting/TOOLS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,13 @@ research pass per tool; see per-tool detail below and the setup matrix first.
method 5 — useful when checking whether a Nyxus percentile matches a tool's *native* percentile
rather than a reimplementation of Nyxus' own. See
`audit/intensity_histogram_2d_analytic_vetting_report.md` for why that distinction matters.
- **`skimage.measure.moments_normalized` answers the CENTRAL-moment question only.** Nyxus also
exposes normalized *raw* moments (`NORM_SPAT_MOMENT_pq`, `IMOM_NRM_pq`) = `m_pq /
m_00^((p+q)/2+1)`, for which skimage has no native function. Do not vet one against the other:
on the 48x40 moments fixture they are 0.3876 vs 0.0999 at p=2,q=0. Build the normalized raw value
from `moments()` and skimage's own exponent instead; see
`audit/moments_2d_skimage_vetting_report.md`. Related trap: Nyxus `m_pq` has p on x and q on y, so
index arrays `A[x, y]` and skimage's `moments(A)[i,j]` maps to `i==p, j==q` with no transpose.
- **`radiomicsj` 2.1.2 is a phantom version** — no GitHub releases/tags; Maven Central publishes
2.1.3 … 2.1.18 (no 2.1.0/2.1.1/2.1.2). Pin **2.1.3** (earliest, nearest the requested tag) or
**2.1.18** (latest). Ships a thin jar → must shade an uber-jar; has a built-in commons-cli `main()`
Expand Down
Loading
Loading