fix(texture): square zone sizes, run lengths and grey levels in double

src/nyx/features/3d_glrlm.cpp

@@ -342,9 +342,10 @@ double D3_GLRLM_feature::calc_SRE (const SimpleMatrix<int> &P, const double sum_
 		}
 	}
 
+	// Run lengths are squared in double: j^2 passes INT_MAX at a run of 46,341.
 	for (int j = 1; j <= Nr; ++j) 
 	{
-		f += rj[j] / double(j * j);
+		f += rj[j] / (double(j) * double(j));
 	}
 
 	double retval = f / sum_p;
@@ -366,7 +367,7 @@ double D3_GLRLM_feature::calc_LRE (const SimpleMatrix<int> &P, const double sum_
 	{
 		for (int j = 1; j <= Nr; j++)
 		{
-			f += P.matlab(i, j) * double(j * j);
+			f += P.matlab(i, j) * (double(j) * double(j));	// j^2 in double: an int overflows at j >= 46,341
 		}
 	}
 
@@ -578,10 +579,12 @@ double D3_GLRLM_feature::calc_LGLRE (const SimpleMatrix<int> &P, const std::vect
 	double f = 0.0;
 	for (int i = 1; i <= Ng; i++)
 	{
-		auto inten = I[i - 1];
+		// Grey levels are taken in double: a level squared passes UINT_MAX at 65,536, and so does a
+		// level times a run once their product reaches 65,536.
+		double inten = (double) I[i - 1];
 		for (int j = 1; j <= Nr; j++)
 		{
-			f += P.matlab(i, j) / double(inten * inten);
+			f += P.matlab(i, j) / (inten * inten);
 		}
 	}
 	double retval = f / double(sum_p);
@@ -601,10 +604,10 @@ double D3_GLRLM_feature::calc_HGLRE (const SimpleMatrix<int> &P, const std::vect
 	double f = 0.0;
 	for (int i = 1; i <= Ng; i++)
 	{
-		auto inten = I[i - 1];
+		double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 		for (int j = 1; j <= Nr; j++)
 		{
-			f += P.matlab(i, j) * double(inten) * double(inten);
+			f += P.matlab(i, j) * (inten * inten);
 		}
 	}
 	double retval = f / double(sum_p);
@@ -624,10 +627,10 @@ double D3_GLRLM_feature::calc_SRLGLE (const SimpleMatrix<int> &P, const std::vec
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I [i-1];
+			double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) / double(inten * inten * j * j);
+				f += P.matlab(i, j) / (inten * inten * (double(j) * double(j)));	// j^2 in double: an int overflows at j >= 46,341
 			}
 		}
 		double retval = f / sum_p;
@@ -647,10 +650,10 @@ double D3_GLRLM_feature::calc_SRHGLE (const SimpleMatrix<int> &P, const std::vec
 	double f = 0.0;
 	for (int i = 1; i <= Ng; i++)
 	{
-		auto inten = I[i - 1];
+		double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 		for (int j = 1; j <= Nr; j++)
 		{
-			f += P.matlab(i, j) * double(inten * inten) / double(j * j);
+			f += P.matlab(i, j) * (inten * inten) / (double(j) * double(j));	// j^2 in double: an int overflows at j >= 46,341
 		}
 	}
 	double retval = f / sum_p;
@@ -670,10 +673,10 @@ double D3_GLRLM_feature::calc_LRLGLE (const SimpleMatrix<int> &P, const std::vec
 	double f = 0.0;
 	for (int i = 1; i <= Ng; i++)
 	{
-		auto inten = I[i - 1];
+		double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 		for (int j = 1; j <= Nr; j++)
 		{
-			f += P.matlab(i, j) * double(j * j) / double(inten * inten);
+			f += P.matlab(i, j) * (double(j) * double(j)) / (inten * inten);	// j^2 in double: an int overflows at j >= 46,341
 		}
 	}
 	double retval = f / sum_p;
@@ -693,10 +696,10 @@ double D3_GLRLM_feature::calc_LRHGLE (const SimpleMatrix<int> &P, const std::vec
 	double f = 0.0;
 	for (int i = 1; i <= Ng; i++)
 	{
-		auto inten = I[i - 1];
+		double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 		for (int j = 1; j <= Nr; j++)
 		{
-			f += P.matlab(i, j) * double(inten * inten * j * j);
+			f += P.matlab(i, j) * (inten * inten * (double(j) * double(j)));	// j^2 in double: an int overflows at j >= 46,341
 		}
 	}
 	double retval = f / sum_p;

src/nyx/features/3d_glszm.cpp

@@ -820,48 +820,23 @@ void D3_GLSZM_feature::save_value(std::vector<std::vector<double>>& fvals)
 // 1. Small Area Emphasis
 double D3_GLSZM_feature::calc_SAE()
 {
-#if 0
 	// Calculate feature. 'sj' is expected to have been initialized in calc_sums_of_P()
+	// Zone sizes are squared in double: j^2 passes INT_MAX at a zone of 46,341.
 	double f = 0.0;
-	for (int j = 1; j <= /*Ns*/P.width(); j++)
-	{
-		f += sj[j] / (j * j);
-	}
-	double retval = f / sum_p;
-	return retval;
-#endif
-
-	double tot = 0;
-	for (int i=1; i<=Ng; i++)
-	{
-		for (int j=1; j<=P.width(); j++)
-			tot += double(P.matlab(i,j)) / double(j*j);
-	}
-	double retval = tot / Nz;
+	for (int j = 1; j <= Ns; j++)
+		f += sj[j] / (double(j) * double(j));
+	double retval = f / Nz;
 	return retval;
 }
 
 // 2. Large Area Emphasis
 double D3_GLSZM_feature::calc_LAE()
 {
-#if 0
 	// Calculate feature. 'sj' is expected to have been initialized in calc_sums_of_P()
 	double f = 0.0;
-	for (int j = 1; j <= /*Ns*/P.width(); j++)
-	{
-		f += sj[j] * (j * j);
-	}
-	double retval = f / sum_p;
-	return retval;
-#endif
-
-	double tot = 0;
-	for (int i = 1; i <= Ng; i++)
-	{
-		for (int j = 1; j <= P.width(); j++)
-			tot += double(P.matlab(i,j)) * double(j*j);
-	}
-	double retval = tot / Nz;
+	for (int j = 1; j <= Ns; j++)
+		f += sj[j] * (double(j) * double(j));	// j^2 in double: an int overflows at j >= 46,341
+	double retval = f / Nz;
 	return retval;
 }
 

src/nyx/features/glrlm.cpp

@@ -378,8 +378,9 @@ void GLRLMFeature::calc_SRE(AngledFtrs& af)
 			}
 		}
 
+		// Run lengths are squared in double: j^2 passes INT_MAX at a run of 46,341.
 		for (int j = 1; j <= Nr; ++j) {
-			f += rj[j] / (j * j);
+			f += rj[j] / (double(j) * double(j));
 		}
 
 		double retval = f / double(sum_p[ai]);
@@ -410,7 +411,7 @@ void GLRLMFeature::calc_LRE(AngledFtrs& af)
 		{
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) * j * j;
+				f += P.matlab(i, j) * (double(j) * double(j));	// j^2 in double: an int overflows at j >= 46,341
 			}
 		}
 
@@ -711,10 +712,12 @@ void GLRLMFeature::calc_LGLRE (AngledFtrs& af)
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I[i-1];
+			// Grey levels are taken in double: a level squared passes UINT_MAX at 65,536, and so does a
+			// level times a run once their product reaches 65,536.
+			double inten = (double) I[i-1];
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) / double(inten*inten);
+				f += P.matlab(i, j) / (inten * inten);
 			}
 		}
 		double retval = f / double(sum_p[ai]);
@@ -743,10 +746,10 @@ void GLRLMFeature::calc_HGLRE(AngledFtrs& af)
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I[i - 1];
+			double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) * double(inten * inten);
+				f += P.matlab(i, j) * (inten * inten);
 			}
 		}
 		double retval = f / double(sum_p[ai]);
@@ -775,10 +778,10 @@ void GLRLMFeature::calc_SRLGLE(AngledFtrs& af)
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I[i-1];
+			double inten = (double) I[i-1];	// level in double: its square passes UINT_MAX at 65,536
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) / double(inten * inten * j * j);
+				f += P.matlab(i, j) / (inten * inten * (double(j) * double(j)));	// j^2 in double: an int overflows at j >= 46,341
 			}
 		}
 		double retval = f / double(sum_p[ai]);
@@ -807,10 +810,10 @@ void GLRLMFeature::calc_SRHGLE(AngledFtrs& af)
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I[i-1];
+			double inten = (double) I[i-1];	// level in double: its square passes UINT_MAX at 65,536
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) * double(inten * inten) / double(j * j);
+				f += P.matlab(i, j) * (inten * inten) / (double(j) * double(j));	// j^2 in double: an int overflows at j >= 46,341
 			}
 		}
 		double retval = f / double(sum_p[ai]);
@@ -838,10 +841,10 @@ void GLRLMFeature::calc_LRLGLE(AngledFtrs& af)
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I[i-1];
+			double inten = (double) I[i-1];	// level in double: its square passes UINT_MAX at 65,536
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) * double(j * j) / double(inten*inten);
+				f += P.matlab(i, j) * (double(j) * double(j)) / (inten * inten);	// j^2 in double: an int overflows at j >= 46,341
 			}
 		}
 		double retval = f / double(sum_p[ai]);
@@ -870,10 +873,10 @@ void GLRLMFeature::calc_LRHGLE(AngledFtrs& af)
 		double f = 0.0;
 		for (int i = 1; i <= Ng; i++)
 		{
-			auto inten = I[i - 1];
+			double inten = (double) I[i - 1];	// level in double: its square passes UINT_MAX at 65,536
 			for (int j = 1; j <= Nr; j++)
 			{
-				f += P.matlab(i, j) * double(inten * inten * j * j);
+				f += P.matlab(i, j) * (inten * inten * (double(j) * double(j)));	// j^2 in double: an int overflows at j >= 46,341
 			}
 		}
 		double retval = f / double(sum_p[ai]);

src/nyx/features/glszm.cpp

@@ -421,9 +421,10 @@ double GLSZMFeature::calc_SAE()
 {
 	// Calculate feature. 'sj' is expected to have been initialized in calc_sums_of_P()
 	double f = 0.0;
+	// Zone sizes are squared in double: j^2 passes INT_MAX at a zone of 46,341.
 	for (int j = 1; j <= Ns; j++)
 	{
-		f += sj[j] / (j * j);
+		f += sj[j] / (double(j) * double(j));
 	}
 	double retval = f / sum_p;
 	return retval;
@@ -436,7 +437,7 @@ double GLSZMFeature::calc_LAE()
 	double f = 0.0;
 	for (int j = 1; j <= Ns; j++)
 	{
-		f += sj[j] * (j * j);
+		f += sj[j] * (double(j) * double(j));	// j^2 in double: an int overflows at j >= 46,341
 	}
 	double retval = f / sum_p;
 	return retval;

src/nyx/gpu/glszm.cu

@@ -117,8 +117,9 @@ namespace NyxusGpu
 			return;
 
 		float s = SJ[j];
-		atomicAdd(&state[f_SAE], s / (j * j));
-		atomicAdd(&state[f_LAE], s * (j * j));
+		// Zone sizes are squared in floating point: j^2 passes INT_MAX at a zone of 46,341.
+		atomicAdd(&state[f_SAE], s / (float(j) * float(j)));
+		atomicAdd(&state[f_LAE], s * (float(j) * float(j)));
 		atomicAdd(&state[f_SZN], s * s);
 		atomicAdd(&state[f_SZNN], s * s);
 	}

tests/test_2d_glrlm_analytic.h

@@ -0,0 +1,108 @@
+#pragma once
+
+// 2D GLRLM emphases against their closed form, on ROIs whose run lengths and grey levels are large
+// enough that a squared term no longer fits in 32 bits (SPEC 4, oracle `analytic`).
+
+// Nothing of its own: gtest, <vector>, GLRLMFeature, Feature2D, the IBSI settings recipe and the
+// ROI loader all arrive with the common header.
+#include "test_2d_glrlm_common.h"
+
+// Featurises one fully masked ROI given as (x, y, intensity) pixels, at the IBSI recipe, which keeps
+// the grey levels as they are.
+static void run_2d_glrlm_on_pixels (LR& roidata, const std::vector<NyxusPixel>& intensity)
+{
+    std::vector<NyxusPixel> mask (intensity);
+    for (auto& px : mask)
+        px.intensity = 1u;
+    Nyxus::load_masked_test_roi_data (roidata, intensity.data(), mask.data(), intensity.size());
+
+    GLRLMFeature f;
+    f.calculate (roidata, make_2d_glrlm_ibsi_settings());
+    roidata.initialize_fvals();
+    f.save_value (roidata.fvals);
+}
+
+// rel=1e-12 rather than SPEC 7's exact tier (abs 1e-9): these values run from below 1e-9 to past 1e9,
+// so an absolute 1e-9 would accept a small emphasis of 0 and lies below a large one's last bit.
+static void assert_2d_glrlm_closed_form (const LR& roidata, Feature2D fcode, int angle, double expected)
+{
+    const double actual = roidata.fvals[(int)fcode][angle];
+    ASSERT_NEAR (actual, expected, 1e-12 * expected) << "feature " << (int)fcode << " angle " << angle;
+}
+
+// bench_glrlm_2d_two_long_runs -- one row of 2L pixels, L = 65,536, the left half at intensity 1 and
+// the right half at intensity 2. At 0 degrees that is two runs, one of length L at each level, and
+// L^2 = 2^32 wraps to exactly 0 in 32 bits whether the product is read signed or unsigned, so every
+// run-length term is right here only if it is taken in double. A run between 46,341 and 65,535 also
+// overflows an int, but a compiler that reads the overflowed product as unsigned still gets its
+// square exactly. Over the Nr = 2 runs:
+//
+//   SRE    = 1/L^2                 LRE    = L^2
+//   SRLGLE = (1 + 1/4) / 2 / L^2   SRHGLE = (1 + 4) / 2 / L^2
+//   LRLGLE = (1 + 1/4) / 2 * L^2   LRHGLE = (1 + 4) / 2 * L^2
+//
+// At 45, 90 and 135 degrees a one-pixel-high row has no neighbour, so every pixel is a run of
+// length 1: 2L runs, half at each level, giving SRE = LRE = 1, SRLGLE = LRLGLE = (1 + 1/4) / 2 and
+// SRHGLE = LRHGLE = (1 + 4) / 2. Each _AVE is the mean of the four directions.
+void test_2d_glrlm_long_run_emphasis_analytic()
+{
+    const size_t L = 65536;
+    std::vector<NyxusPixel> intensity;
+    for (size_t x = 0; x < 2 * L; x++)
+        intensity.push_back ({ x, 0, x < L ? 1u : 2u });
+
+    LR roidata;
+    ASSERT_NO_FATAL_FAILURE (run_2d_glrlm_on_pixels (roidata, intensity));
+
+    const double L2 = double(L) * double(L);
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRE, 0, 1.0 / L2));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRE, 0, L2));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRLGLE, 0, 0.625 / L2));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRHGLE, 0, 2.5 / L2));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRLGLE, 0, 0.625 * L2));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRHGLE, 0, 2.5 * L2));
+
+    // every other direction, and the mean over all 4
+    for (int angle = 1; angle < 4; angle++)
+    {
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRE, angle, 1.0));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRE, angle, 1.0));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRLGLE, angle, 0.625));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRHGLE, angle, 2.5));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRLGLE, angle, 0.625));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRHGLE, angle, 2.5));
+    }
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRE_AVE, 0, (1.0 / L2 + 3 * 1.0) / 4));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRE_AVE, 0, (L2 + 3 * 1.0) / 4));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRLGLE_AVE, 0, (0.625 / L2 + 3 * 0.625) / 4));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRHGLE_AVE, 0, (2.5 / L2 + 3 * 2.5) / 4));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRLGLE_AVE, 0, (0.625 * L2 + 3 * 0.625) / 4));
+    ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRHGLE_AVE, 0, (2.5 * L2 + 3 * 2.5) / 4));
+}
+
+// bench_glrlm_2d_high_grey_pair -- two adjacent pixels at intensities 1 and G = 70,000. Every
+// direction sees two runs of length 1, and G^2 = 4.9e9 is past UINT_MAX, so the grey-level terms are
+// right here only if the level is squared in double. Over the Nr = 2 runs, in all four directions:
+//
+//   LGLRE = SRLGLE = LRLGLE = (1 + 1/G^2) / 2
+//   HGLRE = SRHGLE = LRHGLE = (1 + G^2) / 2
+void test_2d_glrlm_high_grey_level_emphasis_analytic()
+{
+    const double G = 70000.0;
+    std::vector<NyxusPixel> intensity = { { 0, 0, 1u }, { 1, 0, 70000u } };
+
+    LR roidata;
+    ASSERT_NO_FATAL_FAILURE (run_2d_glrlm_on_pixels (roidata, intensity));
+
+    const double low = (1.0 + 1.0 / (G * G)) / 2.0,
+        high = (1.0 + G * G) / 2.0;
+    for (int angle = 0; angle < 4; angle++)
+    {
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LGLRE, angle, low));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRLGLE, angle, low));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRLGLE, angle, low));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_HGLRE, angle, high));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_SRHGLE, angle, high));
+        ASSERT_NO_FATAL_FAILURE (assert_2d_glrlm_closed_form (roidata, Feature2D::GLRLM_LRHGLE, angle, high));
+    }
+}