Merge pull request PolusAI#216 from JesseMckinzie/image_quality_v5

Add image quality features

CMakeLists.txt

@@ -147,6 +147,7 @@ set(SOURCE
 	src/nyx/features/erosion_pixels.cpp
 	src/nyx/features/euler_number.cpp
 	src/nyx/features/extrema.cpp
+	src/nyx/features/focus_score.cpp
 	src/nyx/features/fractal_dim.cpp
 	src/nyx/features/gabor.cpp
 	src/nyx/features/gabor_nontriv.cpp
@@ -168,9 +169,12 @@ set(SOURCE
 	src/nyx/features/ngldm.cpp
 	src/nyx/features/ngtdm.cpp
 	src/nyx/features/pixel.cpp
+	src/nyx/features/power_spectrum.cpp
 	src/nyx/features/radial_distribution.cpp
 	src/nyx/features/roi_radius.cpp
 	src/nyx/features/rotation.cpp
+	src/nyx/features/saturation.cpp
+	src/nyx/features/sharpness.cpp
 	src/nyx/features/specfunc.cpp
 	src/nyx/features/texture_feature.cpp
 	src/nyx/features/zernike.cpp
@@ -209,7 +213,6 @@ set(SOURCE
 	src/nyx/strpat.cpp
 )
 
-
 # CLI
 if(BUILD_CLI)
 	add_executable(nyxus ${SOURCE} src/nyx/main_nyxus.cpp)

ci-utils/envs/conda_cpp.txt

@@ -9,4 +9,4 @@ xsimd >=8,<9
 cmake
 dcmtk >=3.6.7
 fmjpeg2koj >=1.0.3
-libarrow ==12.0.0
+libarrow ==12.0.0

docs/source/Math.rst

@@ -21,3 +21,4 @@ Math formulas
    Math/f_ngtdm
    Math/f_morphology
    Math/f_zernike
+   Math/f_image_quality

docs/source/Math/f_image_quality.rst

@@ -0,0 +1,140 @@
+Image Quality
+=============
+
+Image quality features provide calculations to determine how blurry or noisy an image is. 
+
+
+Focus Score
+-----------
+
+Focus score first applies the Laplacian operator to the image to detect the edges. The Laplacian 
+operator for :math:`\text{ksize}>1` is a filter using the kernel 
+
+.. math::
+
+    \begin{bmatrix}
+    2 & 0 & 2\\
+    0 & -8 & 0\\
+    2 & 0 & 2 \end{bmatrix}.
+
+For the default value of :math:`\text{ksize}=1`, the kernel becomes
+
+.. math::
+    
+    \begin{bmatrix}
+    0 & 1 & 0\\
+    1 & -4 & 1\\
+    0 & 1 & 0 \end{bmatrix}.
+
+The filtering uses a reflective border condition. For example, given the matrix
+
+.. math::
+    
+    \begin{bmatrix}
+    1 & 2 & 3\\
+    4 & 5 & 6\\
+    7 & 8 & 9 \end{bmatrix},
+
+the padded matrix becomes
+
+.. math::
+    
+    \begin{bmatrix}
+    9 & 8 & 7 & 8 & 9 & 8 & 7\\
+    6 & 5 & 4 & 5 & 6 & 5 & 4\\
+    3 & 2 & 1 & 2 & 3 & 2 & 1\\
+    6 & 5 & 4 & 5 & 6 & 5 & 4\\
+    9 & 8 & 7 & 8 & 9 & 8 & 7\\
+    6 & 5 & 4 & 5 & 6 & 5 & 4\\
+    3 & 2 & 1 & 2 & 3 & 2 & 1\end{bmatrix}.
+
+After calculating the Laplacian of the image, the focus score is given by the variance of the filtered image :math:`x` is
+defined by 
+
+.. math::
+    
+    \text{variance}(x) = \frac{\sum_i(x_i - \text{mean}(x))}{(\text{length}(x) - 1)}
+
+Local Focus Score
+-----------
+
+Local Focus Score is calculated the same way as Focus Score. However, the focus score is calculated for each tile 
+of the image, where the number of tiles are determined by the input parameter `scale`. The image is divided into
+:math:`\text{scale}^2` non-overlapping tiles and the mean and median values of the tiles are returned.
+
+GLCM Correlation and Dissimilarity
+----------------------------------
+
+See the documentation for these features in :doc:`GLCM math documentation <f_glcm>`. 
+
+Power Spectrum Slope
+--------------------
+
+This feature calculates the slope of the image's log-log power spectrum. The power spectrum of an :math:`m\times n` image is obtained 
+through Fourier transform
+
+.. math::
+
+    F(u,v) = \sum(x,y)l(x,y)-\mu\mu w(x,y)e^{2\pi i(ux/m+vy/n)},
+
+where :math:`u` and :math:`v` are spatial frequency coordinate. The power spectrum :math:`S`can then be found by the 
+square amplitude of the Fourier transform
+
+.. math::
+
+    S(u,v) = L|F(u,v)|^2\Gamma(u,v),
+
+where :math:`\Gamma` is the correction factor and :math:`L` is the number of pixels in the image. The resulting power spectrum
+of the image is then converted to log-log scale and the slope of log-log power spectrum is obtained through the least squares solution.
+
+Saturation Metrics
+------------------
+
+The min saturation and max saturation metrics are obtained by counting the number of pixels at the maximum and minimum intensities of the image and returning 
+the percent of pixels found at both intensity levels.
+
+Sharpness
+---------
+
+The Sharpness calculation aims to determine the blurriness of edges in an image. To accomplish this estimation,
+a difference of differences in gray scale value of the median-filtered image, :math:`\Delta \text{DoM}`, is used[1]. 
+To begin, median-filtering is applied to the image to reduce noise while preserving edges. The median-filtered image :math:`I_m` is calculated
+using a 3x3 averaging filter, given by,
+
+.. math::
+
+    I_m = \frac{1}{9}\sum_{j=-1}^1\sum_{i=-1}^1 \text{I}(x+i, x+j).
+
+The :math:`\Delta`DoM is calculated separately for horizontal and vertical directions. The horizontal direction is calculated as,
+
+.. math::
+
+    \Delta\text{DoM}_x(i,j) = [I_m(i+2, j) - I_m(i,j)] - [I_m(i, j) - I_m(i-2, j)].
+
+
+The sharpness at a pixel :math:`S_x(i,j)` is computed as
+
+.. math::
+
+    S_x(i,j) = \frac{\sum_{i-w\leq k \leq i+w} |\Delta \text{DoM}(k,j)|}{\sum_{i-w\leq k \leq i+w} |I(k,j)-I(k-l,j)|}.
+
+
+This ratio is high at sharp locations and low at blurred locations. The sharpness in the vertical or :math:`y` direction is calculated in a similar way.
+The sharpness for the entire image is calculated as,
+
+
+.. math::
+
+    R_x = \frac{\text{#}SharpPixels_x}{\text{#}EdgePixels_x} \\ \\
+    R_y= \frac{\text{#}SharpPixels_y}{\text{#}EdgePixels_y}.
+
+
+The sharpness of the image is then given by the Frobenius norm,
+
+.. math::
+
+    S_I = \sqrt{R_x^2+R_y^2}.
+
+
+[1] J. Kumar, F. Chen and D. Doermann, "Sharpness estimation for document and scene images," Proceedings of the 21st International Conference on Pattern Recognition (ICPR2012), Tsukuba, Japan, 2012, pp. 3292-3295.
+

docs/source/imagequality.rst

@@ -0,0 +1,27 @@
+Image Quality features
+======================
+
+Image quality features are available in Nyxus to determine how blurry an image is. 
+These features are available in the `IMAGE_QUALITY` feature group. This group is not included
+in the `ALL` group and must be enabled separately. The following features are in the image quality
+group.
+
+.. list-table::
+   :header-rows: 1
+
+   * - Image Quality feature code
+     - Description
+   * - FOCUS_SCORE
+     - Uses edge detection to highlight regions where intesnity changes rapidly. Higher focus score means lower blurriness
+   * - LOCAL_FOCUS_SCORE
+     - Tiles image into non-overlapping regions and calculates the FOCUS_SCORE for each region. Higher local focus score means lower blurriness
+   * - GLCM_DIS
+     - Blurry images low dissimilarity
+   * - GLCM_CORRELATION
+     - Blurry images have a high correlation
+   * - POWER_SPECTRUM
+     - The slope of the image log-log power spectrum. A low score means a blurry image
+   * - SATURATION
+     - Percent of pixels at minimum and maximum pixel values
+   * - SHARPNESS
+     - Uses median-filtered image as indicator of edge scharpness. Values range from 0 to sqrt(2). Low scores indicate blurriness.

docs/source/index.rst

@@ -22,4 +22,5 @@ to operate at any scale, its highly validated algorithms, and its modular nature
    devguide
    Math
    supported_formats
+   imagequality
 

src/nyx/env_features.cpp

@@ -166,6 +166,8 @@ bool Environment::spellcheck_raw_featurelist(const std::string& comma_separated_
 	return success;
 }
 
+
+
 // Returns:
 //		true if s is recognized as a group name
 //		false if not recognized (s may be an individual feature name)
@@ -400,6 +402,22 @@ bool Environment::expand_3D_featuregroup (const std::string& s)
 	return false;
 }
 
+// Returns:
+//		true if s is recognized as a group name
+//		false if not recognized (s may be an individual feature name)
+//
+bool Environment::expand_IMQ_featuregroup (const std::string & s)
+{
+	if (s == Nyxus::theFeatureSet.findGroupNameByCode(FgroupIMQ::ALL_IMQ))
+	{
+		Nyxus::theFeatureSet.enableAllIMQ();
+		return true; 
+	}
+
+	return false;
+}
+
+
 void Environment::expand_featuregroups()
 {
 	theFeatureSet.enableAll(false); 
@@ -408,6 +426,18 @@ void Environment::expand_featuregroups()
 		// Enforce the feature names to be in uppercase
 		s = Nyxus::toupper(s);
 
+		if (is_imq()) {
+			if (expand_IMQ_featuregroup (s)) 
+				return;
+
+			FeatureIMQ a;
+			if (!theFeatureSet.find_IMQ_FeatureByString(s, a))
+				throw std::invalid_argument("Error: '" + s + "' is not a valid Image Quality feature name \n");
+
+			theFeatureSet.enableFeature (int(a));
+			continue;
+		}
+
 		if (dim() == 2)
 		{
 			if (expand_2D_featuregroup (s))

src/nyx/environment.h

@@ -107,6 +107,8 @@ class Environment: public BasicEnvironment
 	// Returns the expected dataset dimensionality based on the command line options
 	int dim() { return dim_; }
 	void set_dim(int d) { dim_ = d; }
+	bool is_imq() {return is_imq_;};
+	void set_imq(bool is_imq) {is_imq_ = is_imq;}
 
 	bool singleROI = false; // Applies to dim()==2: singleROI is set to 'true' parse_cmdline() if labels_dir==intensity_dir
 
@@ -149,6 +151,9 @@ class Environment: public BasicEnvironment
 	size_t get_ram_limit();
 	void expand_featuregroups();
 
+	void expand_IMQ_featuregroups();
+
+
 	static bool gpu_is_available();
 
 	static bool ibsi_compliance;
@@ -238,6 +243,9 @@ class Environment: public BasicEnvironment
 	std::string raw_dim = "";
 	bool expand_2D_featuregroup (const std::string& name);
 	bool expand_3D_featuregroup (const std::string& name);
+	bool expand_IMQ_featuregroup (const std::string & s);
+
+	bool is_imq_ = false;
 
 	// data members implementing exclusive-inclusive timing switch
 	#ifdef CHECKTIMING

src/nyx/feature_method.cpp

@@ -13,6 +13,12 @@ void FeatureMethod::provide_features (const std::initializer_list<Nyxus::Feature
 		provided_features.push_back((int)f);
 }
 
+void FeatureMethod::provide_features (const std::initializer_list<Nyxus::FeatureIMQ> & F)
+{
+	for (auto f : F)
+		provided_features.push_back((int)f);
+}
+
 void FeatureMethod::provide_features (const std::initializer_list<Nyxus::Feature3D> & F)
 {
 	for (auto f : F)
@@ -25,6 +31,12 @@ void FeatureMethod::add_dependencies (const std::initializer_list<Nyxus::Feature
 		dependencies.push_back((int)f);
 }
 
+void FeatureMethod::add_dependencies (const std::initializer_list<Nyxus::FeatureIMQ>& F)
+{
+	for (auto f : F)
+		dependencies.push_back((int)f);
+}
+
 void FeatureMethod::add_dependencies(const std::initializer_list<Nyxus::Feature3D>& F)
 {
 	for (auto f : F)
@@ -47,6 +59,11 @@ bool FeatureMethod::depends (Nyxus::Feature2D fcode)
 	return std::find (dependencies.begin(), dependencies.end(), (int)fcode) != dependencies.end();
 }
 
+bool FeatureMethod::depends (Nyxus::FeatureIMQ fcode)
+{
+	return std::find (dependencies.begin(), dependencies.end(), (int)fcode) != dependencies.end();
+}
+
 bool FeatureMethod::depends (Nyxus::Feature3D fcode)
 {
 	return std::find(dependencies.begin(), dependencies.end(), (int)fcode) != dependencies.end();

src/nyx/feature_method.h

@@ -39,6 +39,10 @@ class FeatureMethod
 	/// @param F 
 	void provide_features (const std::initializer_list<Nyxus::Feature2D> & F);	// Queried by provides()
 
+	/// @brief 
+	/// @param F 
+	void provide_features (const std::initializer_list<Nyxus::FeatureIMQ> & F);	// Queried by provides()
+
 	/// @brief 
 	/// @param F 
 	void provide_features (const std::initializer_list<Nyxus::Feature3D> & F);	// Queried by provides()
@@ -52,12 +56,14 @@ class FeatureMethod
 	/// @param F 
 	void add_dependencies (const std::initializer_list<Nyxus::Feature2D>& F);	// Queried with depends()
 	void add_dependencies (const std::initializer_list<Nyxus::Feature3D>& F);	// Queried with depends()
+	void add_dependencies (const std::initializer_list<Nyxus::FeatureIMQ>& F);	// Queried with depends()
 
 	/// @brief Checks if a feature code is among dependencies
 	/// @param  
 	/// @return 
 	bool depends (Nyxus::Feature2D);	// Looks up in 'dependencies'
 	bool depends (Nyxus::Feature3D);	// Looks up in 'dependencies'
+	bool depends (Nyxus::FeatureIMQ);	// Looks up in 'dependencies'
 
 private:
 	// Dependency manager support

src/nyx/feature_mgr_init.cpp

@@ -28,6 +28,10 @@
 #include "features/caliper.h"
 #include "features/roi_radius.h"
 #include "features/zernike.h"
+#include "features/focus_score.h"
+#include "features/power_spectrum.h"
+#include "features/saturation.h"
+#include "features/sharpness.h"
 
 FeatureManager::FeatureManager()
 {
@@ -60,10 +64,12 @@ FeatureManager::FeatureManager()
 	register_feature (new GaborFeature());
 	register_feature (new ZernikeFeature());
 	register_feature (new RadialDistributionFeature());
-	register_feature(new PixelIntensityFeatures_3D());
-}
-
-bool FeatureManager::init_feature_classes()
+	register_feature (new PixelIntensityFeatures_3D());
+	register_feature (new FocusScoreFeature());
+	register_feature (new PowerSpectrumFeature());
+	register_feature (new SaturationFeature());
+	register_feature (new SharpnessFeature());
+}bool FeatureManager::init_feature_classes()
 {
 	return GaborFeature::init_class();
 }