diff --git a/src/nyx/features/focus_score.cpp b/src/nyx/features/focus_score.cpp index cb43d592..5e76dbb9 100644 --- a/src/nyx/features/focus_score.cpp +++ b/src/nyx/features/focus_score.cpp @@ -254,7 +254,10 @@ void FocusScoreFeature::laplacian(const std::vector& 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(image[ii * n_image + jj]) * kernel[ikFlip * n_kernel + jkFlip]; } } } diff --git a/tests/test_imq_regression.h b/tests/test_imq_regression.h index 90512618..47305ada 100644 --- a/tests/test_imq_regression.h +++ b/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); diff --git a/tests/vetting/coverage_report.md b/tests/vetting/coverage_report.md index 4e31cf53..b663d752 100644 --- a/tests/vetting/coverage_report.md +++ b/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 | diff --git a/tests/vetting/oracle_coverage.csv b/tests/vetting/oracle_coverage.csv index 8882b744..eae2a8ff 100644 --- a/tests/vetting/oracle_coverage.csv +++ b/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." diff --git a/tests/vetting/oracles/gen_imq_analytic.py b/tests/vetting/oracles/gen_imq_analytic.py new file mode 100644 index 00000000..9148e81c --- /dev/null +++ b/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())