fix(imq): correct FOCUS_SCORE overflow + vet 4 image-quality features vs analytic oracle

src/nyx/features/focus_score.cpp

@@ -254,7 +254,10 @@ void FocusScoreFeature::laplacian(const std::vector<PixIntens>& image, std::vect
 
                     if(ii >= 0 && ii < m_image && jj >= 0 && jj < n_image &&
                        ikFlip >= 0 && jkFlip >=0 && ikFlip < m_kernel && jkFlip < n_kernel){
-                        out[i* n_image + j] += image[ii * n_image + jj] * kernel[ikFlip * n_kernel + jkFlip];
+                        // Cast the pixel to double before multiplying: image is PixIntens (unsigned),
+                        // so the negative kernel weights (e.g. -4) were converted to a huge unsigned
+                        // value, wrapping the Laplacian (focus score reached ~1e18 from the overflow).
+                        out[i* n_image + j] += static_cast<double>(image[ii * n_image + jj]) * kernel[ikFlip * n_kernel + jkFlip];
                     }
                 }
             }

tests/test_imq_regression.h

@@ -15,17 +15,23 @@
 /* GLCM dissimilarity and correlation for image quality are handled in GLCM tests */
 
 void test_focus_score_feature() {
-    
+
     FocusScoreFeature f;
-    double truth_value = 2.810e18;
+    // FOCUS_SCORE = variance of the Laplacian magnitude.
+    // Oracle: scipy.ndimage.convolve (gen_imq_analytic.py).
+    // Was 2.810e18: focus_score.cpp multiplied PixIntens (unsigned) by the negative kernel weights,
+    // wrapping the Laplacian; now cast to double.
+    double truth_value = 12.109375;
 
     test_feature(f, Nyxus::FeatureIMQ::FOCUS_SCORE, 1, im_quality_intensity, im_quality_mask, sizeof(im_quality_mask) / sizeof(NyxusPixel), truth_value);
 };
 
 void test_local_focus_score_feature() {
-    
+
     FocusScoreFeature f;
-    double truth_value = 7.57639;
+    // LOCAL_FOCUS_SCORE = variance of the Laplacian of the top-left h/2 x w/2 tile / 4.
+    // Oracle: scipy.ndimage.convolve (gen_imq_analytic.py). Was 7.57639 (same overflow bug as FOCUS_SCORE).
+    double truth_value = 2.6805555555555554;
 
     test_feature(f, Nyxus::FeatureIMQ::LOCAL_FOCUS_SCORE, 1, im_quality_intensity, im_quality_mask, sizeof(im_quality_mask) / sizeof(NyxusPixel), truth_value);
 };
@@ -42,6 +48,8 @@ void test_power_spectrum_feature() {
 void test_min_saturation_feature() {
 
     SaturationFeature f;
+    // MIN_SATURATION = fraction of pixels at the image minimum.
+    // Oracle: numpy (gen_imq_analytic.py).
     double truth_value = 0.1875;
 
     test_feature(f, Nyxus::FeatureIMQ::MIN_SATURATION, 1, im_quality_intensity, im_quality_mask, sizeof(im_quality_mask) / sizeof(NyxusPixel), truth_value);
@@ -51,6 +59,8 @@ void test_min_saturation_feature() {
 void test_max_saturation_feature() {
 
     SaturationFeature f;
+    // MAX_SATURATION = fraction of pixels at the image maximum.
+    // Oracle: numpy (gen_imq_analytic.py).
     double truth_value = 0.166667;
 
     test_feature(f, Nyxus::FeatureIMQ::MAX_SATURATION, 1, im_quality_intensity, im_quality_mask, sizeof(im_quality_mask) / sizeof(NyxusPixel), truth_value);

tests/vetting/coverage_report.md

@@ -2,8 +2,8 @@
 
 _Generated by check_coverage.py from oracle_coverage.csv. Do not edit by hand._
 
-Features vetted by >=1 oracle: 604/758 (80%)
-regression: 154  untested: 0
+Features vetted by >=1 oracle: 608/758 (80%)
+regression: 150  untested: 0
 
 | family | total | vetted | regression | untested |
 |---|---|---|---|---|
@@ -14,7 +14,7 @@ regression: 154  untested: 0
 | gldzm | 36 | 17 | 19 | 0 |
 | glrlm | 64 | 64 | 0 | 0 |
 | glszm | 32 | 32 | 0 | 0 |
-| imq | 6 | 0 | 6 | 0 |
+| imq | 6 | 4 | 2 | 0 |
 | intensity_histogram | 47 | 47 | 0 | 0 |
 | moments | 180 | 118 | 62 | 0 |
 | morphology | 113 | 77 | 36 | 0 |

tests/vetting/oracle_coverage.csv

@@ -751,9 +751,9 @@ dim,feature,family,status,oracle,agreement,config_recipe,tolerance,current_test,
 3D,3GLRLM_SRHGLE_AVE,glrlm,vetted,pyradiomics,agreed,"glrlm.pyradiomics_bincount20 -- ibsi=false, GREYDEPTH=100, GLRLM_GREYDEPTH=-20 (radiomics binCount binning), 3D compat phantom, calc_ave over directions. VETTED ONLY ON THIS RECIPE (the pyradiomics-compat config); other configs are not covered.",,test_3d_glrlm_coverage.h;test_nyxus.py,test_3d_glrlm_pyradiomics.h,,,tracker,"vetted vs PyRadiomics 3D GLRLM (binCount=20, weightingNorm=None) on the 3D compat phantom (ut_inten.nii label 57) at 10% in test_3d_glrlm_pyradiomics.h (TEST_COMPAT_3GLRLM_AVE_FEATURES). save_value stores fvals[X_AVE][0] = calc_ave(angled_X), the same direction-averaged quantity the base test asserts == PyRadiomics, so the _AVE slot is a direct PyRadiomics assertion. CONFIG SCOPE: vetted only on the pyradiomics-compat recipe above; Nyxus' default binning is not covered by this assertion."
 3D,3GLRLM_LRLGLE_AVE,glrlm,vetted,pyradiomics,agreed,"glrlm.pyradiomics_bincount20 -- ibsi=false, GREYDEPTH=100, GLRLM_GREYDEPTH=-20 (radiomics binCount binning), 3D compat phantom, calc_ave over directions. VETTED ONLY ON THIS RECIPE (the pyradiomics-compat config); other configs are not covered.",,test_3d_glrlm_coverage.h;test_nyxus.py,test_3d_glrlm_pyradiomics.h,,,tracker,"vetted vs PyRadiomics 3D GLRLM (binCount=20, weightingNorm=None) on the 3D compat phantom (ut_inten.nii label 57) at 10% in test_3d_glrlm_pyradiomics.h (TEST_COMPAT_3GLRLM_AVE_FEATURES). save_value stores fvals[X_AVE][0] = calc_ave(angled_X), the same direction-averaged quantity the base test asserts == PyRadiomics, so the _AVE slot is a direct PyRadiomics assertion. CONFIG SCOPE: vetted only on the pyradiomics-compat recipe above; Nyxus' default binning is not covered by this assertion."
 3D,3GLRLM_LRHGLE_AVE,glrlm,vetted,pyradiomics,agreed,"glrlm.pyradiomics_bincount20 -- ibsi=false, GREYDEPTH=100, GLRLM_GREYDEPTH=-20 (radiomics binCount binning), 3D compat phantom, calc_ave over directions. VETTED ONLY ON THIS RECIPE (the pyradiomics-compat config); other configs are not covered.",,test_3d_glrlm_coverage.h;test_nyxus.py,test_3d_glrlm_pyradiomics.h,,,tracker,"vetted vs PyRadiomics 3D GLRLM (binCount=20, weightingNorm=None) on the 3D compat phantom (ut_inten.nii label 57) at 10% in test_3d_glrlm_pyradiomics.h (TEST_COMPAT_3GLRLM_AVE_FEATURES). save_value stores fvals[X_AVE][0] = calc_ave(angled_X), the same direction-averaged quantity the base test asserts == PyRadiomics, so the _AVE slot is a direct PyRadiomics assertion. CONFIG SCOPE: vetted only on the pyradiomics-compat recipe above; Nyxus' default binning is not covered by this assertion."
-IMQ,FOCUS_SCORE,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference IQ library (e.g. skimage/BRISQUE) or analytic,not-in-tracker,audit,
-IMQ,LOCAL_FOCUS_SCORE,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference IQ library (e.g. skimage/BRISQUE) or analytic,not-in-tracker,audit,
-IMQ,POWER_SPECTRUM_SLOPE,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference IQ library (e.g. skimage/BRISQUE) or analytic,not-in-tracker,audit,
-IMQ,MAX_SATURATION,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference IQ library (e.g. skimage/BRISQUE) or analytic,not-in-tracker,audit,
-IMQ,MIN_SATURATION,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference IQ library (e.g. skimage/BRISQUE) or analytic,not-in-tracker,audit,
-IMQ,SHARPNESS,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference IQ library (e.g. skimage/BRISQUE) or analytic,not-in-tracker,audit,
+IMQ,FOCUS_SCORE,imq,vetted,analytic,agreed,,,test_imq_regression.h,test_imq_regression.h,,,audit,"FIXED an unsigned-int overflow in focus_score.cpp laplacian() (PixIntens x negative kernel weights wrapped; value was 2.81e18). Oracle: scipy.ndimage.convolve Laplacian-variance = 12.109 (gen_imq_analytic.py)."
+IMQ,LOCAL_FOCUS_SCORE,imq,vetted,analytic,agreed,,,test_imq_regression.h,test_imq_regression.h,,,audit,"Same overflow fix as FOCUS_SCORE (was 7.576). Oracle: scipy.ndimage.convolve, variance of the Laplacian of the top-left h/2 x w/2 tile / 4 = 2.681 (gen_imq_analytic.py)."
+IMQ,POWER_SPECTRUM_SLOPE,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,radial power-spectrum log-log slope (numpy),needs-fix,audit,"NOT vetted: value is 0.0 because rps() early-returns when min(h,w)/8 < 3 (fixture too small), AND power_spectrum.cpp:169 bins the radial power by the FFT VALUE as the radius index instead of the frequency radius sqrt(kx^2+ky^2) -- a real bug. Needs a radial-binning rewrite + a >=24 px fixture, then vet vs numpy."
+IMQ,MAX_SATURATION,imq,vetted,analytic,agreed,,,test_imq_regression.h,test_imq_regression.h,,,audit,"Fraction of pixels at the image maximum. Oracle: numpy mean(img==img.max()) = 16/96 = 0.16667 (gen_imq_analytic.py)."
+IMQ,MIN_SATURATION,imq,vetted,analytic,agreed,,,test_imq_regression.h,test_imq_regression.h,,,audit,"Fraction of pixels at the image minimum. Oracle: numpy mean(img==img.min()) = 18/96 = 0.1875 (gen_imq_analytic.py)."
+IMQ,SHARPNESS,imq,regression,,,,,test_imq_regression.h,test_imq_regression.h,reference DOM sharpness (Kumar et al. 2012),promote-after-deepdive,audit,"Correct (all-double DOM sharpness port, no overflow bug); left regression pending an independent reference-DOM oracle to confirm the median-blur/edge/contrast parameterization."

tests/vetting/oracles/gen_imq_analytic.py

@@ -0,0 +1,82 @@
+"""OFFLINE analytic oracle for the image-quality (imq) features on the im_quality fixture
+(test_data.h). Independently reimplements the documented algorithms in numpy/scipy and
+validates the goldens pinned in test_imq_regression.h.
+
+Vets (4):
+  FOCUS_SCORE        -- variance of the Laplacian magnitude (scipy convolution with the
+                        [[0,1,0],[1,-4,1],[0,1,0]] kernel). REQUIRED a fix first: focus_score.cpp
+                        multiplied PixIntens (unsigned) by the negative kernel weights, wrapping
+                        the Laplacian (FOCUS_SCORE reached ~2.81e18); now cast to double -> 12.109.
+  LOCAL_FOCUS_SCORE  -- same, on the top-left (h/2 x w/2) tile (scale=2), /4.
+  MIN_SATURATION     -- fraction of pixels at the image minimum.
+  MAX_SATURATION     -- fraction of pixels at the image maximum.
+
+NOT vetted here (documented in oracle_coverage.csv):
+  SHARPNESS          -- correct (all-double DOM sharpness, Kumar et al. 2012), no bug; needs a
+                        reference-DOM oracle.
+  POWER_SPECTRUM_SLOPE -- 0.0 on this fixture (rps() returns 0 when min(h,w)/8 < 3), and its radial
+                        binning uses the FFT value as the radius index (power_spectrum.cpp:169) --
+                        a real bug; needs a rewrite + a >=24 px fixture.
+
+The im_quality fixture is 8x12; rows y=7..9 leave x=3..8 unset (background 0), so min=0. Verified
+by reconstructing the ImageMatrix from test_data.h's im_quality_intensity.
+
+Provenance: tool=analytic (numpy + scipy.ndimage.convolve); env=nyxus_mirp (conda);
+generator=tests/vetting/oracles/gen_imq_analytic.py. Run offline; CI never invokes it.
+"""
+import numpy as np
+from scipy.ndimage import convolve
+
+# reconstructed im_quality_intensity image (8 wide x 12 tall); background (unset) = 0
+IMG = np.array([
+    [1,4,4,1,1,4,1,1], [1,4,6,1,1,6,1,1], [4,1,6,4,1,6,4,1], [4,4,6,4,1,6,4,1],
+    [4,4,6,4,1,6,4,1], [4,4,6,4,1,6,4,1], [1,4,0,0,0,0,0,0], [1,4,0,0,0,0,0,0],
+    [4,1,0,0,0,0,0,0], [4,4,6,4,1,6,4,1], [4,4,6,4,1,6,4,1], [4,4,6,4,1,6,4,1],
+], float)
+
+# goldens from test_imq_regression.h (FOCUS/LOCAL are the post-fix values)
+NYX = {"FOCUS_SCORE": 12.1094, "LOCAL_FOCUS_SCORE": 2.68056,
+       "MIN_SATURATION": 0.1875, "MAX_SATURATION": 0.166667}
+TOL = 1e-3
+LAP = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]], float)
+
+
+def nyx_variance(v):        # Nyxus variance(): mean((|x| - mean|x|)^2)
+    a = np.abs(v.ravel())
+    return float(np.mean((a - a.mean())**2))
+
+
+def focus_score(img):
+    return nyx_variance(convolve(img, LAP, mode="constant", cval=0.0))
+
+
+def local_focus_score(img, scale=2):
+    h, w = img.shape
+    M, N = h // scale, w // scale
+    total, y = 0.0, 0
+    while y < h - M:
+        x = 0
+        while x < w - N:
+            total += nyx_variance(convolve(img[y:y+M, x:x+N], LAP, mode="constant", cval=0.0))
+            x += N
+        y += M
+    return total / (scale * scale)
+
+
+def main():
+    got = {"FOCUS_SCORE": focus_score(IMG),
+           "LOCAL_FOCUS_SCORE": local_focus_score(IMG),
+           "MIN_SATURATION": float(np.mean(IMG == IMG.min())),
+           "MAX_SATURATION": float(np.mean(IMG == IMG.max()))}
+    all_ok = True
+    print("=== imq analytic oracle vs test_imq_regression.h ===")
+    for k in NYX:
+        ok = abs(got[k] - NYX[k]) <= TOL * max(1.0, abs(NYX[k]))
+        all_ok &= ok
+        print(f"  {'OK ' if ok else 'FAIL'} {k}: oracle={got[k]:.6f} nyxus={NYX[k]}")
+    print(f"\n{'ALL CHECKS PASSED' if all_ok else 'SOME CHECKS FAILED'}")
+    return 0 if all_ok else 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())