POL5826 un-hounsfelding 2D montage (PolusAI#297)

Author: friskluft

README.md

@@ -15,7 +15,7 @@ Nyxus is a feature-rich, optimized, Python/C++ application capable of analyzing
 
 Nyxus can be used via Python or command line and is available in containerized form for reproducible execution. Nyxus computes over 450 combined intensity, texture, and morphological features at the ROI or whole image level with more in development. Key features that make Nyxus unique among other image feature extraction applications is its ability to operate at any scale, its highly validated algorithms, and its modular nature that makes the addition of new features straightforward.
 
-Currently, Nyxus can read image data from OME-TIFF, OME-Zarr and DICOM 2D Grayscale images. It also has a Python API to support in-memory image data via Numpy array. 
+Currently, Nyxus can read 2D image data from OME-TIFF, OME-Zarr, and DICOM 2D Grayscale images. Nyxus also reads compressed and uncompressed NIFTI 3D files. Nyxus Python API supports featurizing in-memory 2D image data represented by NumPy arrays. 
 
 The docs can be found at [Read the Docs](https://nyxus.readthedocs.io/en/latest/).
 
@@ -40,6 +40,8 @@ The library provides class `Nyxus` for 2-dimensional TIFF, OME.TIFF, OME.ZARR, a
 
 Given `intensities` and `labels` folders, Nyxus pairs up intensity-segmentation mask images and extracts features from all of them. A summary of the available feature are [listed below](#available-features).
 
+#### featurizing data in file system directories
+
 ```python 
 from nyxus import Nyxus
 nyx = Nyxus (["*ALL*"])
@@ -48,6 +50,7 @@ maskDir = "/path/to/images/labels/"
 features = nyx.featurize_directory (intensityDir, maskDir) # selecting all the .ome.tif slides (default)
 ```
 
+#### featurizing explicitly defined lists of files
 Alternatively, Nyxus can process explicitly defined pairs of intensity-mask images thus specifying custom 1:N and M:N mapping between segmentation mask and intensity image files. 
 The following example extracts all the features (note parameter "*ALL*") from intensity images 'i1', 'i2', and 'i3' related with mask images 'm1' and 'm2' via a custom mapping:
 
@@ -68,7 +71,7 @@ features = nyx.featurize_files(
 	False) # pass True to featurize intensity files as whole segments
 ```
 
-The result `features` variable is a Pandas dataframe similar to what is shown below. Note that if multiple segments are stored in a segmentation mask file, each segment's features in the resultcan be identified by the mask file name and segment mask label.
+The result variable `features` is a Pandas dataframe similar to what is shown below. Note that if multiple segments are stored in a segmentation mask file, each segment's features in the resultcan be identified by the mask file name and segment mask label.
 
 |     | mask_image           | intensity_image      |   label |    MEAN |   MEDIAN |...|    GABOR_6 |
 |----:|:---------------------|:---------------------|--------:|--------:|---------:|--:|-----------:|
@@ -80,7 +83,9 @@ The result `features` variable is a Pandas dataframe similar to what is shown be
 | ... | ...                  | ...                  |     ... | ...     |  ...     |...|   ...      |
 | 734 | p5_y0_r51_c0.ome.tif | p5_y0_r51_c0.ome.tif |     223 | 54573.3 |  54573.3 |...|   0.980769 |
 
-Nyxus can also process intensity-mask pairs that are loaded as NumPy arrays using the `featurize` method. This method takes in either a single pair of 2D intensity-mask pairs
+#### featurizing in-memory 2D images; featurizing a montage
+
+Nyxus can also featurize in-memory intensity-mask pairs that are loaded as NumPy arrays using the `featurize` method. This method takes in either a single pair of 2D intensity-mask pairs
 or a pair of 3D arrays containing 2D intensity and mask images. There is also two optional parameters to supply names to the resulting dataframe, . 
 
 ```python 
@@ -106,6 +111,34 @@ seg = np.array([
 features = nyx.featurize(intens, seg)
 ```
 
+<u>Note:</u> if array `intens` contains negative values similarly to intensities in Hounsfeld units observed in CT-scan datasets, method `featurize()` automatically adjusts values of of array `intens` while passing them to Nyxus backend so as to make them zero-based, like in the following example:
+
+```
+import numpy as np
+import nyxus
+I = np.array([
+  [-1024.74,      -1019.67,      -1005.70,       -998.60,       -998.66,      -1005.82,      -1019.65,      -1024.72],
+  [-1019.44,      -1001.22,      -1023.82,      -1034.34,      -1035.81,      -1027.00,      -1001.89,      -1019.62],
+  [-1011.86,      -1002.17,       -724.06,       -521.43,       -471.04,       -671.30,      -1006.98,      -1010.62],
+  [-1008.78,       -703.58,         21.66,         44.32,        130.35,        113.37,       -608.11,      -1056.33],
+  [-415.46,       -106.08,         69.80,         59.70,         97.64,        120.62,        -49.77,       -480.57],
+  [-464.06,       -176.81,         76.79,         93.34,        131.99,         73.16,       -106.70,       -348.06],
+  [-1012.75,       -740.21,       -502.72,       -370.36,       -377.42,       -497.65,       -719.82,      -1000.82],
+  [-1032.57,       -979.63,       -867.71,       -815.30,       -830.90,       -875.08,       -983.76,      -1033.78]], np.float64)
+M = np.array([
+  [0,      0,      0,       0,       0,      0,      0,      0],
+  [0,      0,      1,       1,       1,      1,      0,      0],
+  [0,      1,      1,       1,       1,      1,      1,      0],
+  [0,      1,      1,       1,       1,      1,      1,      0],
+  [0,      1,      1,       1,       1,      1,      1,      0],
+  [0,      1,      1,       1,       1,      1,      1,      0],
+  [0,      0,      1,       1,       1,      1,      0,      0],
+  [0,      0,      0,       0,       0,      0,      0,      0]], np.uint16)
+nyx = nyxus.Nyxus (["*ALL_INTENSITY*"])
+f = nyx.featurize (I, M, intensity_names=['I'], label_names=['M'])
+```
+
+
 The `features` variable is a Pandas dataframe similar to what is shown below.
 
 |     | mask_image    | intensity_image | label | MEAN    |   MEDIAN |...|    GABOR_6 |

src/nyx/features/contour.cpp

@@ -374,13 +374,64 @@ void ContourFeature::buildRegularContour(LR& r)
 	r.contour.clear();
 
 	// gather contour pixels undecorating their intensities back to original values
+	Pixel2 lastNonzeroPx (0, 0, 0);
 	for (int y = 0; y < height + 2; y++)
 		for (int x = 0; x < width + 2; x++)
 		{
 			size_t idx = x + y * (width + 2);
 			auto inte = borderImage.at(idx);
 			if (inte)
 			{
+				// this pixel may happen to be isolated (a speckle), nonetheless, remember it 
+				// as we'll need to report it as a degenerate contour if no properly neighbored 
+				// pixel group is found
+				lastNonzeroPx = { x, y, inte - 1 };
+				
+				// register a pixel only if it has any immediate neighbor
+				bool hasNeig = false;
+				if (x > 0)	// left neighbor
+				{
+					size_t idxNeig = (x-1) + y * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (x < width-1)	// right neighbor
+				{
+					size_t idxNeig = (x+1) + y * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (y > 0)	// upper neighbor
+				{
+					size_t idxNeig = x + (y-1) * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (y < height-1)	// lower neighbor
+				{
+					size_t idxNeig = x + (y+1) * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (x>0 && y > 0)	// upper left neighbor
+				{
+					size_t idxNeig = (x-1) + (y-1) * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (x < width-1 && y > 0)	// upper right neighbor
+				{
+					size_t idxNeig = (x+1) + (y-1) * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (x>0 && y < height-1)	// lower left neighbor
+				{
+					size_t idxNeig = (x-1) + (y+1) * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (x < width-1 && y < height-1)	// lower right neighbor
+				{
+					size_t idxNeig = (x+1) + (y+1) * (width+2);
+					hasNeig = hasNeig || borderImage.at(idxNeig) != 0;
+				}
+				if (!hasNeig)
+					continue;
+				// pixel is good, save it
 				Pixel2 p(x, y, inte - 1);
 				r.contour.push_back(p);
 			}
@@ -390,65 +441,71 @@ void ContourFeature::buildRegularContour(LR& r)
 
 	//==== Reorder the contour cloud
 
-	//	--containers for unordered (temp) and ordered (result) pixels
-	std::list<Pixel2> unordered(r.contour.begin(), r.contour.end());
-	std::vector<Pixel2> ordered;
-	ordered.reserve(unordered.size());
-	std::vector<Pixel2> pants;
-
-	//	--initialize vector 'ordered' with 1st pixel of 'unordered'
-	auto itBeg = unordered.begin();
-	Pixel2 pxTip = *itBeg;
-	ordered.push_back(pxTip);
-	unordered.remove(pxTip);
-
-	//	-- tip of the ordered contour
-	pxTip = ordered.at(0);
-
-	//	-- harvest items of 'unordered' 
-	while (unordered.size())
+	// are there any good candidate pixels ?
+	if (r.contour.size())
 	{
-		//	--find tip's neighbors 
-		std::vector<Pixel2> cands = find_cands (unordered, pxTip);
-		if (cands.empty())
+		//	--containers for unordered (temp) and ordered (result) pixels
+		std::list<Pixel2> unordered (r.contour.begin(), r.contour.end());
+		std::vector<Pixel2> ordered;
+		ordered.reserve (unordered.size());
+		std::vector<Pixel2> pants;
+
+		//	--initialize vector 'ordered' with 1st pixel of 'unordered'
+		auto itBeg = unordered.begin();
+		Pixel2 pxTip = *itBeg;
+		ordered.push_back(pxTip);
+		unordered.remove(pxTip);
+
+		//	-- tip of the ordered contour
+		pxTip = ordered.at(0);
+
+		//	-- harvest items of 'unordered' 
+		while (unordered.size())
 		{
-			// -- we have a gap and need to fix it
-			VERBOSLVL4(dump_2d_image_with_halfcontour(borderImage, unordered, ordered, pxTip, width + 2, height + 2, "\nhalfcontour:\n", ""));
-				
-			// -- no 'break;' ,instead, jump the tip to the closest U-pixel
-			Pixel2 pxPants;
-			pxPants = pants.back();
-			pxTip = pxPants;
-			Pixel2 closest = find_closest (unordered, pxTip);
-
-			// -- discharge
-			ordered.push_back (closest);
-			unordered.remove (closest);
-			pxTip = ordered.at(ordered.size() - 1);
+			//	--find tip's neighbors 
+			std::vector<Pixel2> cands = find_cands(unordered, pxTip);
+			if (cands.empty())
+			{
+				// -- we have a gap and need to fix it
+				VERBOSLVL4(dump_2d_image_with_halfcontour(borderImage, unordered, ordered, pxTip, width + 2, height + 2, "\nhalfcontour:\n", ""));
+
+				// -- no 'break;' ,instead, jump the tip to the closest U-pixel
+				Pixel2 pxPants;
+				pxPants = pants.back();
+				pxTip = pxPants;
+				Pixel2 closest = find_closest(unordered, pxTip);
+
+				// -- discharge
+				ordered.push_back(closest);
+				unordered.remove(closest);
+				pxTip = ordered.at(ordered.size() - 1);
+			}
+			else
+			{
+				// -- register pants
+				if (cands.size() >= 2)
+					pants.push_back(pxTip);
+
+				// -- score thems
+				std::vector<int> candScores = score_cands(cands, pxTip);
+
+				// -- choose the best
+				auto itBest = std::min_element(candScores.begin(), candScores.end());
+				int idxBest = (int)std::distance(candScores.begin(), itBest);
+
+				// -- discharge the found pixel from set 'unordered' and update the tip
+				Pixel2& px = cands.at(idxBest);
+				ordered.push_back(px);
+				unordered.remove(px);
+				pxTip = ordered.at(ordered.size() - 1);
+			}
 		}
-		else
-		{
-			// -- register pants
-			if (cands.size() >= 2)
-				pants.push_back(pxTip);
 
-			// -- score thems
-			std::vector<int> candScores = score_cands(cands, pxTip);
-
-			// -- choose the best
-			auto itBest = std::min_element (candScores.begin(), candScores.end());
-			int idxBest = (int)std::distance (candScores.begin(), itBest);
-
-			// -- discharge the found pixel from set 'unordered' and update the tip
-			Pixel2& px = cands.at(idxBest);
-			ordered.push_back(px);
-			unordered.remove(px);
-			pxTip = ordered.at(ordered.size() - 1);
-		}
+		// done sorting. Now set the ordered contour in the ROI
+		r.contour = ordered;
 	}
-
-	// done sorting. Now set the ordered contour in the ROI
-	r.contour = ordered;
+	else
+		r.contour.push_back(lastNonzeroPx);	// just use the last speckle as a contour because we have no legit contour
 
 	VERBOSLVL4(dump_2d_image_with_vertex_chain(borderImage, r.contour, width + 2, height + 2, "\n\n-- ContourFeature / buildRegularContour / Padded contour image + sorted contour--\n", "\n\n"));
 

src/nyx/python/nyxus/nyxus.py

@@ -291,9 +291,9 @@ def featurize_directory(
                 axis=1,
             )
 
-            # Labels should always be uint.
+            # Labels should always be uint
             if "ROI_label" in df.columns:
-                df["ROI_label"] = df.ROI_label.astype(np.uint32)
+                df.ROI_label = df.ROI_label.astype(np.uint32)
 
             return df
 
@@ -399,13 +399,24 @@ def featurize(
         if (label_images.shape[0] != len(label_names)):
             raise ValueError("Number of segmentation names must be the same as the number of images.")
 
+        # (1) check if the intensities are represented in Hounsfeld type scale, and (2) adjust them into the unsigned integer scale
+        I = intensity_images
+        min_raw_I = np.min (intensity_images)
+        if (min_raw_I < 0):
+            I -= min_raw_I
+        if (not isinstance(I.flat[0], np.uint32)):
+            I = I.astype (np.uint32)
+
+        # cast mask data to unsigned integer, too
+        M = label_images.astype (np.uint32)
+
+        # featurize
         if (output_type == 'pandas'):
 
-            header, string_data, numeric_data, error_message = featurize_montage_imp (intensity_images, label_images, intensity_names, label_names, output_type, "")
-            
+            header, string_data, numeric_data, error_message = featurize_montage_imp (I, M, intensity_names, label_names, output_type, "")
             self.error_message = error_message
-            if(error_message != ''):
-                print(error_message)
+            if error_message != '':
+                print (error_message)
 
             df = pd.concat(
                 [
@@ -415,18 +426,16 @@ def featurize(
                 axis=1,
             )
 
-            # Labels should always be uint.
-            if "label" in df.columns:
-                df["label"] = df.label.astype(np.uint32)
+            # labels should always be uint
+            if "ROI_label" in df.columns:
+                df.ROI_label = df.ROI_label.astype(np.uint32)
 
             return df
 
         else:
 
-            ret = featurize_montage_imp (intensity_images, label_images, intensity_names, label_names, output_type, output_path)
-            
+            ret = featurize_montage_imp (I, M, intensity_names, label_names, output_type, output_path)
             self.error_message = ret[0]
-            
             if(self.error_message != ''):
                 raise RuntimeError('Error calculating features: ' + error_message[0])
 
@@ -494,9 +503,9 @@ def featurize_files (
                 axis=1,
             )
 
-            # Labels should always be uint.
-            if "label" in df.columns:
-                df["label"] = df.label.astype(np.uint32)
+            # Labels should always be uint
+            if "ROI_label" in df.columns:
+                df.ROI_label = df.ROI_label.astype(np.uint32)
 
             return df
 
@@ -1057,9 +1066,9 @@ def featurize_directory(
                 ],
                 axis=1,
             )
-            # Labels should always be uint.
-            if "label" in df.columns:
-                df["label"] = df.label.astype(np.uint32)
+            # Labels should always be uint
+            if "ROI_label" in df.columns:
+                df.ROI_label = df.ROI_label.astype(np.uint32)
             return df
         else:
             featurize_directory_3D_imp (intensity_dir, label_dir, file_pattern, output_type, output_path)
@@ -1124,9 +1133,9 @@ def featurize_files (
                 axis=1,
             )
 
-            # Labels should always be uint.
-            if "label" in df.columns:
-                df["label"] = df.label.astype(np.uint32)
+            # Labels should always be uint
+            if "ROI_label" in df.columns:
+                df.ROI_label = df.ROI_label.astype(np.uint32)
 
             return df