diff --git a/.github/workflows/build_cuda11_wheels.yml b/.github/workflows/build_cuda11_wheels.yml
index 1614823c..609d1a28 100644
--- a/.github/workflows/build_cuda11_wheels.yml
+++ b/.github/workflows/build_cuda11_wheels.yml
@@ -12,7 +12,7 @@ jobs:
       matrix:
         os: [windows-2019, ubuntu-20.04]
         cibw_archs: ["auto64"]
-        cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
diff --git a/.github/workflows/build_cuda12_wheels.yml b/.github/workflows/build_cuda12_wheels.yml
index b58d3ad9..91e23747 100644
--- a/.github/workflows/build_cuda12_wheels.yml
+++ b/.github/workflows/build_cuda12_wheels.yml
@@ -12,7 +12,7 @@ jobs:
       matrix:
         os: [windows-2019, ubuntu-20.04]
         cibw_archs: ["auto64"]
-        cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
diff --git a/.github/workflows/build_wheels.yml b/.github/workflows/build_wheels.yml
index b0fb06aa..f62b56af 100644
--- a/.github/workflows/build_wheels.yml
+++ b/.github/workflows/build_wheels.yml
@@ -14,7 +14,7 @@ jobs:
       matrix:
         os: [ubuntu-20.04, macos-13, windows-latest]
         cibw_archs: ["auto64"]
-        cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
diff --git a/.github/workflows/publish_cuda11_pypi.yml b/.github/workflows/publish_cuda11_pypi.yml
index 5d844615..711a11f7 100644
--- a/.github/workflows/publish_cuda11_pypi.yml
+++ b/.github/workflows/publish_cuda11_pypi.yml
@@ -14,7 +14,7 @@ jobs:
       matrix:
         os: [windows-2019, ubuntu-20.04]
         cibw_archs: ["auto64"]
-        cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
diff --git a/.github/workflows/publish_cuda12_pypi.yml b/.github/workflows/publish_cuda12_pypi.yml
index 8a104c00..20db04d0 100644
--- a/.github/workflows/publish_cuda12_pypi.yml
+++ b/.github/workflows/publish_cuda12_pypi.yml
@@ -14,7 +14,7 @@ jobs:
       matrix:
         os: [windows-2019, ubuntu-20.04]
         cibw_archs: ["auto64"]
-        cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
diff --git a/.github/workflows/publish_pypi.yml b/.github/workflows/publish_pypi.yml
index 22d3f60a..385d0780 100644
--- a/.github/workflows/publish_pypi.yml
+++ b/.github/workflows/publish_pypi.yml
@@ -16,7 +16,7 @@ jobs:
       matrix:
         os: [ubuntu-20.04, macos-13, windows-latest]
         cibw_archs: ["auto64"]
-        cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
@@ -86,7 +86,7 @@ jobs:
       matrix:
         os: [macos-13-xlarge]
         cibw_archs: ["arm64"]
-        cibw_build: ["cp39-*", "cp310-*", "cp311-*",  "cp312-*"]
+        cibw_build: ["cp39-*", "cp310-*", "cp311-*", "cp312-*"]
 
     steps:
       - uses: actions/checkout@v3
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c566e906..c4d795ca 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -168,6 +168,8 @@ set(SOURCE
 	src/nyx/features/3d_glcm.cpp
 	src/nyx/features/3d_glcm_nontriv.cpp
 	src/nyx/features/3d_glszm.cpp
+	src/nyx/features/3d_glrlm.cpp
+	src/nyx/features/3d_glrlm_nontriv.cpp
 	src/nyx/features/intensity.cpp
 	src/nyx/features/neighbors.cpp
 	src/nyx/features/ngldm.cpp
diff --git a/ci-utils/install_prereq_linux.sh b/ci-utils/install_prereq_linux.sh
index 2361a962..cb1564e0 100755
--- a/ci-utils/install_prereq_linux.sh
+++ b/ci-utils/install_prereq_linux.sh
@@ -69,7 +69,7 @@ cd ../../
 if [[ $BUILD_BOOST_DEP -eq 1 ]]; then
     for i in {1..5}
     do
-        curl -L https://boostorg.jfrog.io/artifactory/main/release/1.79.0/source/boost_1_79_0.tar.bz2 -o boost_1_79_0.tar.bz2 
+        curl -L https://archives.boost.io/release/1.79.0/source/boost_1_79_0.tar.bz2 -o boost_1_79_0.tar.bz2 
         if [ -f "boost_1_79_0.tar.bz2" ] ; then
             break
         fi
diff --git a/ci-utils/install_prereq_win.bat b/ci-utils/install_prereq_win.bat
index c1d8486c..8b8965e5 100644
--- a/ci-utils/install_prereq_win.bat
+++ b/ci-utils/install_prereq_win.bat
@@ -42,7 +42,7 @@ popd
 
 if "%BUILD_Z5_DEP%" == "1" (
     for /l %%x in (1, 1, 5) do (
-        curl -L https://boostorg.jfrog.io/artifactory/main/release/1.79.0/source/boost_1_79_0.zip -o boost_1_79_0.zip
+        curl -L https://archives.boost.io/release/1.79.0/source/boost_1_79_0.zip -o boost_1_79_0.zip
         if  exist boost_1_79_0.zip (
             goto :continue_boost
         )
diff --git a/src/nyx/env_features.cpp b/src/nyx/env_features.cpp
index c61e80e8..86d15a65 100644
--- a/src/nyx/env_features.cpp
+++ b/src/nyx/env_features.cpp
@@ -451,6 +451,76 @@ bool Environment::expand_3D_featuregroup (const std::string& s)
 		return true;
 	}
 
+	if (s == Nyxus::theFeatureSet.findGroupNameByCode(Fgroup3D::FG3_GLCM))
+	{
+		theFeatureSet.enableAll (false);
+
+		auto F =
+		{
+			Nyxus::Feature3D::GLCM_ACOR,	
+			Nyxus::Feature3D::GLCM_ASM,
+			Nyxus::Feature3D::GLCM_CLUPROM,
+			Nyxus::Feature3D::GLCM_CLUSHADE,
+			Nyxus::Feature3D::GLCM_CLUTEND,
+			Nyxus::Feature3D::GLCM_CONTRAST,	
+			Nyxus::Feature3D::GLCM_CORRELATION,
+			Nyxus::Feature3D::GLCM_DIFAVE,
+			Nyxus::Feature3D::GLCM_DIFENTRO,
+			Nyxus::Feature3D::GLCM_DIFVAR,
+			Nyxus::Feature3D::GLCM_DIS,
+			Nyxus::Feature3D::GLCM_ENERGY,
+			Nyxus::Feature3D::GLCM_ENTROPY,
+			Nyxus::Feature3D::GLCM_HOM1,
+			Nyxus::Feature3D::GLCM_HOM2,
+			Nyxus::Feature3D::GLCM_ID,
+			Nyxus::Feature3D::GLCM_IDN,
+			Nyxus::Feature3D::GLCM_IDM,
+			Nyxus::Feature3D::GLCM_IDMN,	
+			Nyxus::Feature3D::GLCM_INFOMEAS1,
+			Nyxus::Feature3D::GLCM_INFOMEAS2,
+			Nyxus::Feature3D::GLCM_IV,
+			Nyxus::Feature3D::GLCM_JAVE,
+			Nyxus::Feature3D::GLCM_JE,
+			Nyxus::Feature3D::GLCM_JMAX,
+			Nyxus::Feature3D::GLCM_JVAR,
+			Nyxus::Feature3D::GLCM_SUMAVERAGE,	
+			Nyxus::Feature3D::GLCM_SUMENTROPY,	
+			Nyxus::Feature3D::GLCM_SUMVARIANCE,
+			Nyxus::Feature3D::GLCM_VARIANCE,
+			Nyxus::Feature3D::GLCM_ASM_AVE,
+			Nyxus::Feature3D::GLCM_ACOR_AVE,
+			Nyxus::Feature3D::GLCM_CLUPROM_AVE,
+			Nyxus::Feature3D::GLCM_CLUSHADE_AVE,
+			Nyxus::Feature3D::GLCM_CLUTEND_AVE,
+			Nyxus::Feature3D::GLCM_CONTRAST_AVE,
+			Nyxus::Feature3D::GLCM_CORRELATION_AVE,
+			Nyxus::Feature3D::GLCM_DIFAVE_AVE,
+			Nyxus::Feature3D::GLCM_DIFENTRO_AVE,
+			Nyxus::Feature3D::GLCM_DIFVAR_AVE,
+			Nyxus::Feature3D::GLCM_DIS_AVE,
+			Nyxus::Feature3D::GLCM_ENERGY_AVE,
+			Nyxus::Feature3D::GLCM_ENTROPY_AVE,
+			Nyxus::Feature3D::GLCM_HOM1_AVE,
+			Nyxus::Feature3D::GLCM_ID_AVE,
+			Nyxus::Feature3D::GLCM_IDN_AVE,
+			Nyxus::Feature3D::GLCM_IDM_AVE,
+			Nyxus::Feature3D::GLCM_IDMN_AVE,
+			Nyxus::Feature3D::GLCM_IV_AVE,
+			Nyxus::Feature3D::GLCM_JAVE_AVE,
+			Nyxus::Feature3D::GLCM_JE_AVE,
+			Nyxus::Feature3D::GLCM_INFOMEAS1_AVE,
+			Nyxus::Feature3D::GLCM_INFOMEAS2_AVE,
+			Nyxus::Feature3D::GLCM_VARIANCE_AVE,
+			Nyxus::Feature3D::GLCM_JMAX_AVE,
+			Nyxus::Feature3D::GLCM_JVAR_AVE,
+			Nyxus::Feature3D::GLCM_SUMAVERAGE_AVE,
+			Nyxus::Feature3D::GLCM_SUMENTROPY_AVE,
+			Nyxus::Feature3D::GLCM_SUMVARIANCE_AVE
+		};
+
+		theFeatureSet.enableFeatures(F);
+		return true;
+	}
 	if (s == Nyxus::theFeatureSet.findGroupNameByCode(Fgroup3D::FG3_GLDZM))
 	{
 		theFeatureSet.enableAll (false);
@@ -465,6 +535,51 @@ bool Environment::expand_3D_featuregroup (const std::string& s)
 		return true;
 	}	
 	
+	if (s == Nyxus::theFeatureSet.findGroupNameByCode(Fgroup3D::FG3_GLRLM))
+	{
+		theFeatureSet.enableAll(false);
+
+		auto F =
+		{
+			Nyxus::Feature3D::GLRLM_SRE,	
+			Nyxus::Feature3D::GLRLM_LRE, 
+			Nyxus::Feature3D::GLRLM_GLN, 
+			Nyxus::Feature3D::GLRLM_GLNN, 
+			Nyxus::Feature3D::GLRLM_RLN,
+			Nyxus::Feature3D::GLRLM_RLNN, 
+			Nyxus::Feature3D::GLRLM_RP,
+			Nyxus::Feature3D::GLRLM_GLV, 
+			Nyxus::Feature3D::GLRLM_RV, 
+			Nyxus::Feature3D::GLRLM_RE, 
+			Nyxus::Feature3D::GLRLM_LGLRE, 
+			Nyxus::Feature3D::GLRLM_HGLRE, 
+			Nyxus::Feature3D::GLRLM_SRLGLE, 
+			Nyxus::Feature3D::GLRLM_SRHGLE, 
+			Nyxus::Feature3D::GLRLM_LRLGLE, 
+			Nyxus::Feature3D::GLRLM_LRHGLE, 
+			
+			Nyxus::Feature3D::GLRLM_SRE_AVE,
+			Nyxus::Feature3D::GLRLM_LRE_AVE,
+			Nyxus::Feature3D::GLRLM_GLN_AVE,
+			Nyxus::Feature3D::GLRLM_GLNN_AVE,
+			Nyxus::Feature3D::GLRLM_RLN_AVE,
+			Nyxus::Feature3D::GLRLM_RLNN_AVE,
+			Nyxus::Feature3D::GLRLM_RP_AVE,
+			Nyxus::Feature3D::GLRLM_GLV_AVE,
+			Nyxus::Feature3D::GLRLM_RV_AVE,
+			Nyxus::Feature3D::GLRLM_RE_AVE,
+			Nyxus::Feature3D::GLRLM_LGLRE_AVE,
+			Nyxus::Feature3D::GLRLM_HGLRE_AVE,
+			Nyxus::Feature3D::GLRLM_SRLGLE_AVE,
+			Nyxus::Feature3D::GLRLM_SRHGLE_AVE,
+			Nyxus::Feature3D::GLRLM_LRLGLE_AVE,
+			Nyxus::Feature3D::GLRLM_LRHGLE_AVE
+		};
+
+		theFeatureSet.enableFeatures(F);
+		return true;
+	}
+
 	return false;
 }
 
diff --git a/src/nyx/feature_mgr_init.cpp b/src/nyx/feature_mgr_init.cpp
index 9e51f8d4..8406e7c1 100644
--- a/src/nyx/feature_mgr_init.cpp
+++ b/src/nyx/feature_mgr_init.cpp
@@ -32,6 +32,7 @@
 #include "features/3d_glcm.h"
 #include "features/3d_gldzm.h"
 #include "features/3d_glszm.h"
+#include "features/3d_glrlm.h"
 
 #include "features/focus_score.h"
 #include "features/power_spectrum.h"
@@ -76,6 +77,7 @@ FeatureManager::FeatureManager()
 	register_feature (new D3_GLCM_feature());
 	register_feature (new D3_GLDZM_feature());
 	register_feature (new D3_GLSZM_feature());
+	register_feature (new D3_GLRLM_feature());
 
 	// image quality
 	register_feature (new FocusScoreFeature());
diff --git a/src/nyx/features/3d_gldzm.h b/src/nyx/features/3d_gldzm.h
index f465f1f9..697e64b6 100644
--- a/src/nyx/features/3d_gldzm.h
+++ b/src/nyx/features/3d_gldzm.h
@@ -1,13 +1,11 @@
 #pragma once
 
-#pragma once
-
 #include "../feature_method.h"
 #include "texture_feature.h"
 
 /// @brief Grey Level Distance Zone (GLDZM) features
 /// 
-/// Grey Level Dsitance Zone (GLDZM) quantifies distances zones of same intensity to the ROI border
+/// GLDZM quantifies distance zones of same intensity to the ROI border
 
 class D3_GLDZM_feature : public FeatureMethod, public TextureFeature
 {
diff --git a/src/nyx/features/3d_glrlm.cpp b/src/nyx/features/3d_glrlm.cpp
new file mode 100644
index 00000000..ee5ac2b3
--- /dev/null
+++ b/src/nyx/features/3d_glrlm.cpp
@@ -0,0 +1,953 @@
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <numeric>
+#include <sstream>
+#include <unordered_set>
+#include "3d_glrlm.h"
+#include "../environment.h"
+
+using namespace Nyxus;
+
+int D3_GLRLM_feature::n_levels = 0;
+
+
+D3_GLRLM_feature::D3_GLRLM_feature() : FeatureMethod("D3_GLRLM_feature")
+{
+	provide_features(D3_GLRLM_feature::featureset);
+}
+
+void D3_GLRLM_feature::calculate(LR& r)
+{
+	//==== Clear the feature values buffers
+	clear_buffers();
+
+	auto minI = r.aux_min,
+		maxI = r.aux_max;
+
+	// intercept blank ROIs
+	if (minI == maxI)
+	{
+		// insert a non-NAN value for all 4 angles to make the output expecting 4-angled values happy
+		auto w = theEnvironment.nan_substitute;	// safe NAN
+		angled_SRE.resize(4, w);
+		angled_LRE.resize(4, w);
+		angled_GLN.resize(4, w);
+		angled_GLNN.resize(4, w);
+		angled_RLN.resize(4, w);
+		angled_RLNN.resize(4, w);
+		angled_RP.resize(4, w);
+		angled_GLV.resize(4, w);
+		angled_RV.resize(4, w);
+		angled_RE.resize(4, w);
+		angled_LGLRE.resize(4, w);
+		angled_HGLRE.resize(4, w);
+		angled_SRLGLE.resize(4, w);
+		angled_SRHGLE.resize(4, w);
+		angled_LRLGLE.resize(4, w);
+		angled_LRHGLE.resize(4, w);
+
+		bad_roi_data = true;
+		return;
+	}
+
+	//==== Make a list of intensity clusters (zones)
+	using ACluster = std::pair<PixIntens, int>;
+	using AngleZones = std::vector<ACluster>;
+
+	//==== While scanning clusters, learn unique intensities 
+	using AngleUniqInte = std::unordered_set<PixIntens>;
+
+	//==== Iterate angles 
+	for (int angleIdx = 0; angleIdx < 8; angleIdx++)
+	{
+		// Clusters at angle 'angleIdx'
+		AngleZones Z;
+
+		// Unique intensities at angle 'angleIdx'
+		AngleUniqInte U;
+
+		// We need it to estimate the x-dimension of matrix P
+		int maxZoneArea = 0;
+
+		// Copy the image matrix. We'll use it to maintain state of cluster scanning 
+		int w = r.aux_image_cube.width(),
+			h = r.aux_image_cube.height(),
+			d = r.aux_image_cube.depth();
+		SimpleCube<PixIntens> D;
+		D.allocate (w, h, d);
+
+		// Squeeze the intensity range
+		auto greyInfo = theEnvironment.get_coarse_gray_depth();
+		auto greyInfo_localFeature = D3_GLRLM_feature::n_levels;
+		if (greyInfo_localFeature != 0 && greyInfo != greyInfo_localFeature)
+			greyInfo = greyInfo_localFeature;
+		if (Nyxus::theEnvironment.ibsi_compliance)
+			greyInfo = 0;
+
+		bin_intensities_3d (D, r.aux_image_cube, r.aux_min, r.aux_max, greyInfo);
+
+		// allocate intensities matrix
+		if (ibsi_grey_binning(greyInfo))
+		{
+			auto n_ibsi_levels = *std::max_element(D.begin(), D.end());
+			I.resize(n_ibsi_levels);
+			for (int i = 0; i < n_ibsi_levels; i++)
+				I[i] = i + 1;
+		}
+		else // radiomics and matlab
+		{
+			std::unordered_set<PixIntens> U(D.begin(), D.end());
+			U.erase(0);	// discard intensity '0'
+			I.assign(U.begin(), U.end());
+			std::sort(I.begin(), I.end());
+		}
+
+		// Find zones
+		const int VISITED = -1;
+
+		// --- Scan the image and check non-blank pixels' clusters
+		for (int zslice = 0; zslice < d; zslice++)
+		{
+			for (int row = 0; row < h; row++)
+			{
+				for (int col = 0; col < w; col++)
+				{
+					// Find a non-blank pixel
+					auto pi = D.zyx (zslice, row, col);
+					if (pi == 0 || int(pi) == VISITED)
+						continue;
+
+					// Found a non-blank pixel. Find same-intensity neighbourhood of it.
+					std::vector<std::tuple<int, int, int>> history;
+					int x = col, y = row, z = zslice;
+					int zoneArea = 1;
+					D.zyx (z, y, x) = VISITED;
+
+					// State machine scanning the rest of the cluster
+					for (;;)
+					{
+						int dx, dy, dz;
+
+						//********** ang XY,Z = 0
+						// ang X,Y = 0
+						if (angleIdx == 0 && D.safe(z, y, x+1) && D.zyx(z, y, x+1) == pi)
+						{
+							D.zyx(z, y, x + 1) = VISITED;
+							zoneArea++;
+
+							// Remember this pixel
+							history.push_back({ x,y,z });
+							// Advance 
+							x = x + 1;
+							// Proceed
+							continue;
+						}
+
+						// ang X,Y = 45
+						if (angleIdx == 1 && D.safe(z, y+1, x+1) && D.zyx(z, y+1, x+1) == pi)
+						{
+							D.zyx (z, y+1, x+1) = VISITED;
+							zoneArea++;
+
+							history.push_back({ x,y,z });
+							x = x + 1;
+							y = y + 1;
+							continue;
+						}
+
+						// ang X,Y = 90
+						if (angleIdx == 2 && D.safe(z, y+1, x) && D.zyx(z, y+1, x) == pi)
+						{
+							D.zyx(z, y+1, x) = VISITED;
+							zoneArea++;
+
+							history.push_back({ x,y,z });
+							y = y + 1;
+							continue;
+						}
+
+						// ang X,Y = 135
+						if (angleIdx == 3 && D.safe(z, y+1, x-1) && D.zyx(z, y+1, x-1) == pi)
+						{
+							D.zyx(z, y+1, x-1) = VISITED;
+							zoneArea++;
+
+							history.push_back({ x,y,z });
+							x = x - 1;
+							y = y + 1;
+							continue;
+						}
+
+						//********** ang XY,Z != 0
+						// ang X,Y = 0
+						if (angleIdx == 4 && D.safe(z+1, y, x) && D.zyx(z+1, y, x) == pi)
+						{
+							D.zyx(z+1, y, x) = VISITED;
+							zoneArea++;
+
+							// remember this pixel
+							history.push_back({ x,y,z });
+							// advance 
+							z = z + 1;
+							// proceed
+							continue;
+						}
+
+						// ang X,Y = 45
+						if (angleIdx == 5 && D.safe(z+1, y+1, x) && D.zyx(z+1, y+1, x) == pi)
+						{
+							D.zyx(z+1, y+1, x) = VISITED;
+							zoneArea++;
+
+							history.push_back({ x,y,z });
+							z = z + 1;
+							y = y + 1;
+							continue;
+						}
+
+						// ang X,Y = 90
+						if (angleIdx == 6 && D.safe(z+1, y, x+1) && D.zyx(z+1, y, x+1) == pi)
+						{
+							D.zyx(z+1, y, x+1) = VISITED;
+							zoneArea++;
+
+							history.push_back({ x,y,z });
+							z = z + 1;
+							x = x + 1;
+							continue;
+						}
+
+						// ang X,Y = 135
+						if (angleIdx == 7 && D.safe(z+1, y+1, x-1) && D.zyx(z+1, y+1, x-1) == pi)
+						{
+							D.zyx(z+1, y+1, x-1) = VISITED;
+							zoneArea++;
+
+							history.push_back({ x,y,z });
+							z = z + 1;
+							y = y + 1;
+							x = x - 1;
+							continue;
+						}
+
+						// Return from the branch
+						if (history.size() > 0)
+						{
+							// Recollect the coordinate where we diverted from
+							std::tuple<int, int, int> prev = history [history.size() - 1];
+							history.pop_back();
+						}
+
+						// We are done exploring this cluster
+						break;
+					}
+
+					// --2
+					maxZoneArea = std::max(maxZoneArea, zoneArea);
+
+					// --3
+					ACluster clu = { pi, zoneArea };
+					Z.push_back(clu);
+				}
+			}
+		}
+
+		// count non-zero pixels
+		int count = 0;
+
+		for (const auto& px : r.aux_image_cube)
+		{
+			if (px != 0)
+				count++;
+		}
+
+		//==== Create a zone matrix
+
+		int Ng = Environment::ibsi_compliance ? *std::max_element(I.begin(), I.end()) : I.size();
+		int Nr = maxZoneArea;
+		int Nz = (int)Z.size();
+		int Np = count;
+
+		// --allocate the matrix
+		P_matrix P;
+		P.allocate(Nr, Ng);	// w = Nr, h = card(I) aka Ng
+
+		// --iterate zones and fill the matrix
+		for (auto& z : Z)
+		{
+			auto inten = z.first;
+			// row (grey level)
+			int row = -1;
+			if (Environment::ibsi_compliance)
+				row = inten - 1;
+			else
+			{
+				auto lower = std::lower_bound(I.begin(), I.end(), inten);	// enjoy sorted vector 'I'
+				row = int(lower - I.begin());	// intensity index in array of unique intensities 'I'
+			}
+			// col
+			int col = z.second - 1;	// 0-based => -1
+			// update the matrix
+			auto& k = P.xy(col, row);
+			k++;
+		}
+
+		// --save this angle's results
+		angles_P.push_back(P);
+		angles_Ng.push_back(Ng);
+		angles_Nr.push_back(Nr);
+		angles_Np.push_back(Np);
+
+		double sum = 0;
+		for (int i = 1; i <= P.height(); ++i)
+			for (int j = 1; j <= P.width(); ++j)
+				sum += P.matlab(i, j);
+		sum = 0;
+		for (auto p : P)
+			sum += p;
+
+		sum_p.push_back(sum);
+	} //- angles
+
+	calc_SRE(angled_SRE);
+	calc_LRE(angled_LRE);
+	calc_GLN(angled_GLN);
+	calc_GLNN(angled_GLNN);
+	calc_RLN(angled_RLN);
+	calc_RLNN(angled_RLNN);
+	calc_RP(angled_RP);
+	calc_GLV(angled_GLV);
+	calc_RV(angled_RV);
+	calc_RE(angled_RE);
+	calc_LGLRE(angled_LGLRE);
+	calc_HGLRE(angled_HGLRE);
+	calc_SRLGLE(angled_SRLGLE);
+	calc_SRHGLE(angled_SRHGLE);
+	calc_LRLGLE(angled_LRLGLE);
+	calc_LRHGLE(angled_LRHGLE);
+}
+
+void D3_GLRLM_feature::clear_buffers()
+{
+	bad_roi_data = false;
+	angles_Ng.clear();
+	angles_Nr.clear();
+	angles_Np.clear();
+	angles_P.clear();
+	sum_p.clear();
+	I.clear();
+
+	angled_SRE.clear();
+	angled_LRE.clear();
+	angled_GLN.clear();
+	angled_GLNN.clear();
+	angled_RLN.clear();
+	angled_RLNN.clear();
+	angled_RP.clear();
+	angled_GLV.clear();
+	angled_RV.clear();
+	angled_RE.clear();
+	angled_LGLRE.clear();
+	angled_HGLRE.clear();
+	angled_SRLGLE.clear();
+	angled_SRHGLE.clear();
+	angled_LRLGLE.clear();
+	angled_LRHGLE.clear();
+}
+
+// Not supporting the online mode
+void D3_GLRLM_feature::osized_add_online_pixel(size_t x, size_t y, uint32_t intensity) {} // Not supporting
+
+void D3_GLRLM_feature::save_value(std::vector<std::vector<double>>& fvals)
+{
+	fvals[(int)Feature3D::GLRLM_SRE] = angled_SRE;
+	fvals[(int)Feature3D::GLRLM_LRE] = angled_LRE;
+	fvals[(int)Feature3D::GLRLM_GLN] = angled_GLN;
+	fvals[(int)Feature3D::GLRLM_GLNN] = angled_GLNN;
+	fvals[(int)Feature3D::GLRLM_RLN] = angled_RLN;
+	fvals[(int)Feature3D::GLRLM_RLNN] = angled_RLNN;
+	fvals[(int)Feature3D::GLRLM_RP] = angled_RP;
+	fvals[(int)Feature3D::GLRLM_GLV] = angled_GLV;
+	fvals[(int)Feature3D::GLRLM_RV] = angled_RV;
+	fvals[(int)Feature3D::GLRLM_RE] = angled_RE;
+	fvals[(int)Feature3D::GLRLM_LGLRE] = angled_LGLRE;
+	fvals[(int)Feature3D::GLRLM_HGLRE] = angled_HGLRE;
+	fvals[(int)Feature3D::GLRLM_SRLGLE] = angled_SRLGLE;
+	fvals[(int)Feature3D::GLRLM_SRHGLE] = angled_SRHGLE;
+	fvals[(int)Feature3D::GLRLM_LRLGLE] = angled_LRLGLE;
+	fvals[(int)Feature3D::GLRLM_LRHGLE] = angled_LRHGLE;
+
+	// -- averages --
+	fvals[(int)Feature3D::GLRLM_SRE_AVE][0] = calc_ave(angled_SRE);
+	fvals[(int)Feature3D::GLRLM_LRE_AVE][0] = calc_ave(angled_LRE);
+	fvals[(int)Feature3D::GLRLM_GLN_AVE][0] = calc_ave(angled_GLN);
+	fvals[(int)Feature3D::GLRLM_GLNN_AVE][0] = calc_ave(angled_GLNN);
+	fvals[(int)Feature3D::GLRLM_RLN_AVE][0] = calc_ave(angled_RLN);
+	fvals[(int)Feature3D::GLRLM_RLNN_AVE][0] = calc_ave(angled_RLNN);
+	fvals[(int)Feature3D::GLRLM_RP_AVE][0] = calc_ave(angled_RP);
+	fvals[(int)Feature3D::GLRLM_GLV_AVE][0] = calc_ave(angled_GLV);
+	fvals[(int)Feature3D::GLRLM_RV_AVE][0] = calc_ave(angled_RV);
+	fvals[(int)Feature3D::GLRLM_RE_AVE][0] = calc_ave(angled_RE);
+	fvals[(int)Feature3D::GLRLM_LGLRE_AVE][0] = calc_ave(angled_LGLRE);
+	fvals[(int)Feature3D::GLRLM_HGLRE_AVE][0] = calc_ave(angled_HGLRE);
+	fvals[(int)Feature3D::GLRLM_SRLGLE_AVE][0] = calc_ave(angled_SRLGLE);
+	fvals[(int)Feature3D::GLRLM_SRHGLE_AVE][0] = calc_ave(angled_SRHGLE);
+	fvals[(int)Feature3D::GLRLM_LRLGLE_AVE][0] = calc_ave(angled_LRLGLE);
+	fvals[(int)Feature3D::GLRLM_LRHGLE_AVE][0] = calc_ave(angled_LRHGLE);
+}
+
+
+// 1. Short Run Emphasis 
+// ai - angle index
+void D3_GLRLM_feature::calc_SRE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		double f = 0.;
+		std::vector<double> rj(Nr + 1, 0.);
+		for (int i = 1; i <= Ng; ++i) {
+			for (int j = 1; j <= Nr; ++j) {
+				rj[j] += P.matlab(i, j);
+			}
+		}
+
+		for (int j = 1; j <= Nr; ++j) {
+			f += rj[j] / (j * j);
+		}
+
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 2. Long Run Emphasis 
+void D3_GLRLM_feature::calc_LRE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) * j * j;
+			}
+		}
+
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 3. Gray Level Non-Uniformity 
+void D3_GLRLM_feature::calc_GLN(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			double sum = 0.0;
+			for (int j = 1; j <= Nr; j++)
+			{
+				sum += P.matlab(i, j);
+			}
+			f += sum * sum;
+		}
+
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 4. Gray Level Non-Uniformity Normalized 
+void D3_GLRLM_feature::calc_GLNN(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			double sum = 0.0;
+			for (int j = 1; j <= Nr; j++)
+			{
+				sum += P.matlab(i, j);
+			}
+			f += sum * sum;
+		}
+
+		double retval = f / double(sum_p[ai] * sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 5. Run Length Non-Uniformity
+void D3_GLRLM_feature::calc_RLN(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int j = 1; j <= Nr; j++)
+		{
+			double sum = 0.0;
+			for (int i = 1; i <= Ng; i++)
+			{
+				sum += P.matlab(i, j);
+			}
+			f += sum * sum;
+		}
+
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 6. Run Length Non-Uniformity Normalized 
+void D3_GLRLM_feature::calc_RLNN(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int j = 1; j <= Nr; j++)
+		{
+			double sum = 0.0;
+			for (int i = 1; i <= Ng; i++)
+			{
+				sum += P.matlab(i, j);
+			}
+			f += sum * sum;
+		}
+
+		double retval = f / double(sum_p[ai] * sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 7. Run Percentage
+void D3_GLRLM_feature::calc_RP(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Np = angles_Np[ai];
+
+		double retval = double(sum_p[ai] / Np);
+		af.push_back(retval);
+	}
+}
+
+// 8. Gray Level Variance 
+void D3_GLRLM_feature::calc_GLV(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0)
+		{
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double mu = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				mu += P.matlab(i, j) / sum_p[ai] * inten;
+			}
+		}
+
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				double mu2 = (inten - mu) * (inten - mu);
+				f += P.matlab(i, j) / sum_p[ai] * mu2;
+			}
+		}
+		af.push_back(f);
+	}
+}
+
+// 9. Run Variance 
+void D3_GLRLM_feature::calc_RV(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double mu = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			for (int j = 1; j <= Nr; j++)
+			{
+				mu += P.matlab(i, j) / sum_p[ai] * j;
+			}
+		}
+
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			for (int j = 1; j <= Nr; j++)
+			{
+				double mu2 = (j - mu) * (j - mu);
+				f += P.matlab(i, j) / sum_p[ai] * mu2;
+			}
+		}
+		af.push_back(f);
+	}
+}
+
+// 10. Run Entropy 
+void D3_GLRLM_feature::calc_RE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			for (int j = 1; j <= Nr; j++)
+			{
+				double entrTerm = fast_log10(P.matlab(i, j) / sum_p[ai] + EPS) / LOG10_2;
+				f += P.matlab(i, j) / sum_p[ai] * entrTerm;
+			}
+		}
+		double retval = -f;
+		af.push_back(retval);
+	}
+}
+
+// 11. Low Gray Level Run Emphasis 
+void D3_GLRLM_feature::calc_LGLRE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) / double(inten * inten);
+			}
+		}
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 12. High Gray Level Run Emphasis 
+void D3_GLRLM_feature::calc_HGLRE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) * double(inten * inten);
+			}
+		}
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 13. Short Run Low Gray Level Emphasis 
+void D3_GLRLM_feature::calc_SRLGLE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) / double(inten * inten * j * j);
+			}
+		}
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 14. Short Run High Gray Level Emphasis 
+void D3_GLRLM_feature::calc_SRHGLE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) * double(inten * inten) / double(j * j);
+			}
+		}
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 15. Long Run Low Gray Level Emphasis 
+void D3_GLRLM_feature::calc_LRLGLE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0)
+		{
+			af.push_back(0.0);
+			continue;
+		}
+
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) * double(j * j) / double(inten * inten);
+			}
+		}
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+// 16. Long Run High Gray Level Emphasis 
+void D3_GLRLM_feature::calc_LRHGLE(AngledFtrs& af)
+{
+	af.clear();
+
+	for (int ai = 0; ai < 4; ai++)
+	{
+		if (sum_p[ai] == 0) {
+			af.push_back(0.0);
+			continue;
+		}
+
+		// Get ahold of the requested angle's matrix and its related N parameters 
+		int Ng = angles_Ng[ai],
+			Nr = angles_Nr[ai];
+		const SimpleMatrix<int>& P = angles_P[ai];
+
+		// Calculate
+		double f = 0.0;
+		for (int i = 1; i <= Ng; i++)
+		{
+			auto inten = I[i - 1];
+			for (int j = 1; j <= Nr; j++)
+			{
+				f += P.matlab(i, j) * double(inten * inten * j * j);
+			}
+		}
+		double retval = f / double(sum_p[ai]);
+		af.push_back(retval);
+	}
+}
+
+void D3_GLRLM_feature::reduce (size_t start, size_t end, std::vector<int>* ptrLabels, std::unordered_map <int, LR>* ptrLabelData)
+{
+	for (auto i = start; i < end; i++)
+	{
+		int lab = (*ptrLabels)[i];
+		LR& r = (*ptrLabelData)[lab];
+
+		D3_GLRLM_feature glrlm;
+		glrlm.calculate(r);
+		glrlm.save_value(r.fvals);
+	}
+}
+
+// 'afv' is angled feature values
+double D3_GLRLM_feature::calc_ave(const std::vector<double>& afv)
+{
+	if (afv.empty())
+		return 0;
+
+	double n = static_cast<double> (afv.size()),
+		ave = std::reduce(afv.begin(), afv.end()) / n;
+
+	return ave;
+}
\ No newline at end of file
diff --git a/src/nyx/features/3d_glrlm.h b/src/nyx/features/3d_glrlm.h
new file mode 100644
index 00000000..e657c6c7
--- /dev/null
+++ b/src/nyx/features/3d_glrlm.h
@@ -0,0 +1,167 @@
+#pragma once
+
+#include <unordered_map>
+#include "../roi_cache.h"
+#include "../feature_method.h"
+#include "image_matrix.h"
+#include "texture_feature.h"
+
+/// @brief Gray Level Run Length Matrix(GLRLM) features
+/// Gray Level Run Length Matrix(GLRLM) quantifies gray level runs, which are defined as the length in number of
+/// pixels, of consecutive pixels that have the same gray level value.In a gray level run length matrix
+/// 	: math:`\textbf{ P }(i, j | \theta)`, the :math:`(i, j)^ {\text{ th }}` element describes the number of runs with gray level
+/// 	: math:`i`and length :math:`j` occur in the image(ROI) along angle : math:`\theta`.
+/// 
+
+class D3_GLRLM_feature : public FeatureMethod, public TextureFeature
+{
+public:
+
+	// Codes of features implemented by this class. Used in feature manager's mechanisms, 
+	// in the feature group nickname expansion, and in the feature value output 
+	const constexpr static std::initializer_list<Nyxus::Feature3D> featureset =
+	{
+		Nyxus::Feature3D::GLRLM_SRE,		// Short Run Emphasis 
+		Nyxus::Feature3D::GLRLM_LRE,		// Long Run Emphasis 
+		Nyxus::Feature3D::GLRLM_GLN,		// Gray Level Non-Uniformity 
+		Nyxus::Feature3D::GLRLM_GLNN,		// Gray Level Non-Uniformity Normalized 
+		Nyxus::Feature3D::GLRLM_RLN,		// Run Length Non-Uniformity
+		Nyxus::Feature3D::GLRLM_RLNN,		// Run Length Non-Uniformity Normalized 
+		Nyxus::Feature3D::GLRLM_RP,		// Run Percentage
+		Nyxus::Feature3D::GLRLM_GLV,		// Gray Level Variance 
+		Nyxus::Feature3D::GLRLM_RV,		// Run Variance 
+		Nyxus::Feature3D::GLRLM_RE,		// Run Entropy 
+		Nyxus::Feature3D::GLRLM_LGLRE,	// Low Gray Level Run Emphasis 
+		Nyxus::Feature3D::GLRLM_HGLRE,	// High Gray Level Run Emphasis 
+		Nyxus::Feature3D::GLRLM_SRLGLE,	// Short Run Low Gray Level Emphasis 
+		Nyxus::Feature3D::GLRLM_SRHGLE,	// Short Run High Gray Level Emphasis 
+		Nyxus::Feature3D::GLRLM_LRLGLE,	// Long Run Low Gray Level Emphasis 
+		Nyxus::Feature3D::GLRLM_LRHGLE,	// Long Run High Gray Level Emphasis 
+		// averaged features:
+		Nyxus::Feature3D::GLRLM_SRE_AVE,
+		Nyxus::Feature3D::GLRLM_LRE_AVE,
+		Nyxus::Feature3D::GLRLM_GLN_AVE,
+		Nyxus::Feature3D::GLRLM_GLNN_AVE,
+		Nyxus::Feature3D::GLRLM_RLN_AVE,
+		Nyxus::Feature3D::GLRLM_RLNN_AVE,
+		Nyxus::Feature3D::GLRLM_RP_AVE,
+		Nyxus::Feature3D::GLRLM_GLV_AVE,
+		Nyxus::Feature3D::GLRLM_RV_AVE,
+		Nyxus::Feature3D::GLRLM_RE_AVE,
+		Nyxus::Feature3D::GLRLM_LGLRE_AVE,
+		Nyxus::Feature3D::GLRLM_HGLRE_AVE,
+		Nyxus::Feature3D::GLRLM_SRLGLE_AVE,
+		Nyxus::Feature3D::GLRLM_SRHGLE_AVE,
+		Nyxus::Feature3D::GLRLM_LRLGLE_AVE,
+		Nyxus::Feature3D::GLRLM_LRHGLE_AVE
+	};
+
+	// Features implemented by this class that do not require vector-like angled output. Instead, they are output as a single values
+	const constexpr static std::initializer_list<Nyxus::Feature3D> nonAngledFeatures =
+	{
+		Nyxus::Feature3D::GLRLM_SRE_AVE,
+		Nyxus::Feature3D::GLRLM_LRE_AVE,
+		Nyxus::Feature3D::GLRLM_GLN_AVE,
+		Nyxus::Feature3D::GLRLM_GLNN_AVE,
+		Nyxus::Feature3D::GLRLM_RLN_AVE,
+		Nyxus::Feature3D::GLRLM_RLNN_AVE,
+		Nyxus::Feature3D::GLRLM_RP_AVE,
+		Nyxus::Feature3D::GLRLM_GLV_AVE,
+		Nyxus::Feature3D::GLRLM_RV_AVE,
+		Nyxus::Feature3D::GLRLM_RE_AVE,
+		Nyxus::Feature3D::GLRLM_LGLRE_AVE,
+		Nyxus::Feature3D::GLRLM_HGLRE_AVE,
+		Nyxus::Feature3D::GLRLM_SRLGLE_AVE,
+		Nyxus::Feature3D::GLRLM_SRHGLE_AVE,
+		Nyxus::Feature3D::GLRLM_LRLGLE_AVE,
+		Nyxus::Feature3D::GLRLM_LRHGLE_AVE
+	};
+
+	D3_GLRLM_feature();
+
+	void calculate(LR& r);
+	void osized_add_online_pixel(size_t x, size_t y, uint32_t intensity);
+	void osized_calculate(LR& r, ImageLoader& imloader);
+	void save_value(std::vector<std::vector<double>>& feature_vals);
+	static void reduce(size_t start, size_t end, std::vector<int>* ptrLabels, std::unordered_map <int, LR>* ptrLabelData);
+
+	// Compatibility with the manual reduce
+	static int required(const FeatureSet& fs)
+	{
+		return fs.anyEnabled(D3_GLRLM_feature::featureset);
+	}
+
+	static int n_levels;	// default value: 0
+	using P_matrix = SimpleMatrix<int>;
+	using AngledFtrs = std::vector<double>;
+
+	// 1. Short Run Emphasis 
+	void calc_SRE(AngledFtrs& af);
+	// 2. Long Run Emphasis 
+	void calc_LRE(AngledFtrs& af);
+	// 3. Gray Level Non-Uniformity 
+	void calc_GLN(AngledFtrs& af);
+	// 4. Gray Level Non-Uniformity Normalized 
+	void calc_GLNN(AngledFtrs& af);
+	// 5. Run Length Non-Uniformity
+	void calc_RLN(AngledFtrs& af);
+	// 6. Run Length Non-Uniformity Normalized 
+	void calc_RLNN(AngledFtrs& af);
+	// 7. Run Percentage
+	void calc_RP(AngledFtrs& af);
+	// 8. Gray Level Variance 
+	void calc_GLV(AngledFtrs& af);
+	// 9. Run Variance 
+	void calc_RV(AngledFtrs& af);
+	// 10. Run Entropy 
+	void calc_RE(AngledFtrs& af);
+	// 11. Low Gray Level Run Emphasis 
+	void calc_LGLRE(AngledFtrs& af);
+	// 12. High Gray Level Run Emphasis 
+	void calc_HGLRE(AngledFtrs& af);
+	// 13. Short Run Low Gray Level Emphasis 
+	void calc_SRLGLE(AngledFtrs& af);
+	// 14. Short Run High Gray Level Emphasis 
+	void calc_SRHGLE(AngledFtrs& af);
+	// 15. Long Run Low Gray Level Emphasis 
+	void calc_LRLGLE(AngledFtrs& af);
+	// 16. Long Run High Gray Level Emphasis 
+	void calc_LRHGLE(AngledFtrs& af);
+
+	constexpr static int rotAngles[] = { 0, 45, 90, 135 };
+
+private:
+
+	std::vector<double> angled_SRE,
+		angled_LRE,
+		angled_GLN,
+		angled_GLNN,
+		angled_RLN,
+		angled_RLNN,
+		angled_RP,
+		angled_GLV,
+		angled_RV,
+		angled_RE,
+		angled_LGLRE,
+		angled_HGLRE,
+		angled_SRLGLE,
+		angled_SRHGLE,
+		angled_LRLGLE,
+		angled_LRHGLE;
+
+	double calc_ave(const std::vector<double>& angled_feature_vals);
+
+	bool bad_roi_data = false;	// used to prevent calculation of degenerate ROIs
+	std::vector<int> angles_Ng;	// number of discrete intensity values in the image
+	std::vector<int> angles_Nr; // number of discrete run lengths in the image
+	std::vector<int> angles_Np; // number of voxels in the image
+	std::vector<P_matrix> angles_P;
+	std::vector<double> sum_p;
+	std::vector<PixIntens> I;	// sorted unique intensities
+
+	void clear_buffers();
+
+	const double EPS = 2.2e-16;
+	const double BAD_ROI_FVAL = 0.0;
+	const double LOG10_2 = 0.30102999566;	// precalculated log 2 base 10
+};
\ No newline at end of file
diff --git a/src/nyx/features/3d_glrlm_nontriv.cpp b/src/nyx/features/3d_glrlm_nontriv.cpp
new file mode 100644
index 00000000..d4831526
--- /dev/null
+++ b/src/nyx/features/3d_glrlm_nontriv.cpp
@@ -0,0 +1,11 @@
+#include "../environment.h"
+#include "3d_glrlm.h"
+#include "image_matrix_nontriv.h"
+
+using namespace Nyxus;
+
+void D3_GLRLM_feature::osized_calculate (LR& r, ImageLoader& imloader)
+{
+	calculate(r);
+}
+
diff --git a/src/nyx/featureset.cpp b/src/nyx/featureset.cpp
index 670fa45b..e3ccdd61 100644
--- a/src/nyx/featureset.cpp
+++ b/src/nyx/featureset.cpp
@@ -796,7 +796,41 @@ namespace Nyxus
 		{ "GLSZM_SALGLE", Feature3D::GLSZM_SALGLE },
 		{ "GLSZM_SAHGLE", Feature3D::GLSZM_SAHGLE },
 		{ "GLSZM_LALGLE", Feature3D::GLSZM_LALGLE },
-		{ "GLSZM_LAHGLE", Feature3D::GLSZM_LAHGLE }
+		{ "GLSZM_LAHGLE", Feature3D::GLSZM_LAHGLE },
+
+		{ "3GLRLM_SRE",					Feature3D::GLRLM_SRE },
+		{ "3GLRLM_LRE",					Feature3D::GLRLM_LRE },
+		{ "3GLRLM_GLN",					Feature3D::GLRLM_GLN },
+		{ "3GLRLM_GLNN",				Feature3D::GLRLM_GLNN },
+		{ "3GLRLM_RLN",					Feature3D::GLRLM_RLN },
+		{ "3GLRLM_RLNN",				Feature3D::GLRLM_RLNN },
+		{ "3GLRLM_RP",					Feature3D::GLRLM_RP },
+		{ "3GLRLM_GLV",					Feature3D::GLRLM_GLV },
+		{ "3GLRLM_RV",					Feature3D::GLRLM_RV },
+		{ "3GLRLM_RE",					Feature3D::GLRLM_RE },
+		{ "3GLRLM_LGLRE",				Feature3D::GLRLM_LGLRE },
+		{ "3GLRLM_HGLRE",				Feature3D::GLRLM_HGLRE },
+		{ "3GLRLM_SRLGLE",			Feature3D::GLRLM_SRLGLE },
+		{ "3GLRLM_SRHGLE",			Feature3D::GLRLM_SRHGLE },
+		{ "3GLRLM_LRLGLE",			Feature3D::GLRLM_LRLGLE },
+		{ "3GLRLM_LRHGLE",			Feature3D::GLRLM_LRHGLE },
+		{ "3GLRLM_SRE_AVE",			Feature3D::GLRLM_SRE_AVE },
+		{ "3GLRLM_LRE_AVE",			Feature3D::GLRLM_LRE_AVE },
+		{ "3GLRLM_GLN_AVE",			Feature3D::GLRLM_GLN_AVE },
+		{ "3GLRLM_GLNN_AVE",		Feature3D::GLRLM_GLNN_AVE },
+		{ "3GLRLM_RLN_AVE",			Feature3D::GLRLM_RLN_AVE },
+		{ "3GLRLM_RLNN_AVE",		Feature3D::GLRLM_RLNN_AVE },
+		{ "3GLRLM_RP_AVE",			Feature3D::GLRLM_RP_AVE },
+		{ "3GLRLM_GLV_AVE",			Feature3D::GLRLM_GLV_AVE },
+		{ "3GLRLM_RV_AVE",			Feature3D::GLRLM_RV_AVE },
+		{ "3GLRLM_RE_AVE",			Feature3D::GLRLM_RE_AVE },
+		{ "3GLRLM_LGLRE_AVE",		Feature3D::GLRLM_LGLRE_AVE },
+		{ "3GLRLM_HGLRE_AVE",		Feature3D::GLRLM_HGLRE_AVE },
+		{ "3GLRLM_SRLGLE_AVE",	Feature3D::GLRLM_SRLGLE_AVE },
+		{ "3GLRLM_SRHGLE_AVE",	Feature3D::GLRLM_SRHGLE_AVE },
+		{ "3GLRLM_LRLGLE_AVE",	Feature3D::GLRLM_LRLGLE_AVE },
+		{ "3GLRLM_LRHGLE_AVE",	Feature3D::GLRLM_LRHGLE_AVE }
+
 	};
 
 	std::map <std::string, Fgroup3D> UserFacing3dFeaturegroupNames =
@@ -808,6 +842,7 @@ namespace Nyxus
 		{ "*3D_GLCM*", Fgroup3D::FG3_GLCM },
 		{ "*3D_GLDZM*", Fgroup3D::FG3_GLDZM },
 		{ "*3D_GLSZM*", Fgroup3D::FG3_GLSZM },
+		{ "*3D_GLRLM*", Fgroup3D::FG3_GLRLM },
 		{ "*3D_ALL_NEIGHBOR*", Fgroup3D::FG3_NEIG },
 		{ "*3D_MOMENTS*", Fgroup3D::FG3_MOMENTS },
 	};
diff --git a/src/nyx/featureset.h b/src/nyx/featureset.h
index 3442efa0..3dfbee65 100644
--- a/src/nyx/featureset.h
+++ b/src/nyx/featureset.h
@@ -757,6 +757,41 @@ namespace Nyxus
 		GLSZM_LALGLE,	// Large Area Low Gray Level Emphasis
 		GLSZM_LAHGLE,	// Large Area High Gray Level Emphasis
 
+		// GLRLM:
+		GLRLM_SRE,	// Short Run Emphasis 
+		GLRLM_LRE,	// Long Run Emphasis 
+		GLRLM_GLN,	// Gray Level Non-Uniformity 
+		GLRLM_GLNN,	// Gray Level Non-Uniformity Normalized 
+		GLRLM_RLN,	// Run Length Non-Uniformity
+		GLRLM_RLNN,	// Run Length Non-Uniformity Normalized 
+		GLRLM_RP,	// Run Percentage
+		GLRLM_GLV,	// Gray Level Variance 
+		GLRLM_RV,	// Run Variance 
+		GLRLM_RE,	// Run Entropy 
+		GLRLM_LGLRE,	// Low Gray Level Run Emphasis 
+		GLRLM_HGLRE,	// High Gray Level Run Emphasis 
+		GLRLM_SRLGLE,	// Short Run Low Gray Level Emphasis 
+		GLRLM_SRHGLE,	// Short Run High Gray Level Emphasis 
+		GLRLM_LRLGLE,	// Long Run Low Gray Level Emphasis 
+		GLRLM_LRHGLE,	// Long Run High Gray Level Emphasis 
+		// -- averages --
+		GLRLM_SRE_AVE,
+		GLRLM_LRE_AVE,
+		GLRLM_GLN_AVE,
+		GLRLM_GLNN_AVE,
+		GLRLM_RLN_AVE,
+		GLRLM_RLNN_AVE,
+		GLRLM_RP_AVE,
+		GLRLM_GLV_AVE,
+		GLRLM_RV_AVE,
+		GLRLM_RE_AVE,
+		GLRLM_LGLRE_AVE,
+		GLRLM_HGLRE_AVE,
+		GLRLM_SRLGLE_AVE,
+		GLRLM_SRHGLE_AVE,
+		GLRLM_LRLGLE_AVE,
+		GLRLM_LRHGLE_AVE,
+
 		_COUNT_
 	};
 
@@ -804,6 +839,7 @@ namespace Nyxus
 		FG3_GLCM,
 		FG3_GLDZM,
 		FG3_GLSZM,
+		FG3_GLRLM,
 		FG3_NEIG,
 		FG3_MOMENTS,
 		_COUNT_
diff --git a/src/nyx/reduce_trivial_rois.cpp b/src/nyx/reduce_trivial_rois.cpp
index 5b8a8873..2b06dcdc 100644
--- a/src/nyx/reduce_trivial_rois.cpp
+++ b/src/nyx/reduce_trivial_rois.cpp
@@ -45,6 +45,7 @@
 #include "features/3d_glcm.h"
 #include "features/3d_gldzm.h"
 #include "features/3d_glszm.h"
+#include "features/3d_glrlm.h"
 //--future-- #include "features/3d_surface.h"
 
 #include "features/focus_score.h"
@@ -389,6 +390,11 @@ namespace Nyxus
 		{
 			STOPWATCH("3D GLSZM/3DGLSZM/3DGLSZM/#FFFF00", "\t=");
 			runParallel (D3_GLSZM_feature::reduce, n_threads, work_per_thread, job_size, &L, &roiData);
+		}		
+		if (D3_GLRLM_feature::required(theFeatureSet))
+		{
+			STOPWATCH("3D GLRLM/3DGLRLM/3DGLRLM/#FFFF00", "\t=");
+			runParallel(D3_GLRLM_feature::reduce, n_threads, work_per_thread, job_size, &L, &roiData);
 		}
 
 		//==== morphology/surface