Use RLE to calculate morphological features

CMakeLists.txt

@@ -69,7 +69,7 @@ if(NOT CMAKE_BUILD_TYPE)
   set(CMAKE_BUILD_TYPE Release)
 endif()
 
-set(CMAKE_CXX_FLAGS_RELEASE "-O2")
+set(CMAKE_CXX_FLAGS_RELEASE "-O2 -march=native")
 
 #== GTest 
 option(RUN_GTEST "Downloads google unit test API and runs google test scripts to test Nyxus" OFF)
@@ -137,11 +137,13 @@ set(SOURCE
 	src/nyx/output_2_csv.cpp
 	src/nyx/parallel.cpp
 	src/nyx/phase1.cpp
+	src/nyx/phase1_fastloop.cpp
 	src/nyx/phase2.cpp
 	src/nyx/phase3.cpp
 	src/nyx/pixel_feed.cpp
 	src/nyx/reduce_by_feature.cpp
 	src/nyx/reduce_trivial_rois.cpp
+	src/nyx/rle.cpp
 	src/nyx/roi_cache.cpp
 	src/nyx/roi_cache_basic.cpp
 	src/nyx/scan_fastloader_way.cpp

src/nyx/features/aabb.h

@@ -26,6 +26,14 @@ class AABB
 		xmin = std::min(xmin, x);
 		xmax = std::max(xmax, x);
 	}
+	void update_xmax(StatsInt x)
+	{
+		xmax = std::max(xmax, x);
+	}
+	void update_xmin(StatsInt x)
+	{
+		xmin = std::min(xmin, x);
+	}
 	void update_y(StatsInt y)
 	{
 		ymin = std::min(ymin, y);
@@ -52,4 +60,4 @@ class AABB
 
 private:
 	StatsInt xmin = INT32_MAX, xmax = INT32_MIN, ymin = INT32_MAX, ymax = INT32_MIN;
-};
+};

src/nyx/features/basic_morphology.cpp

@@ -36,10 +36,19 @@ void BasicMorphologyFeatures::calculate(LR& r)
 	// --CENTROID_XY
 	double cen_x = 0.0,
 		cen_y = 0.0;
-	for (auto& px : r.raw_pixels)
-	{
-		cen_x += px.x;
-		cen_y += px.y;
+	if (r.rle_pixels.empty()) {
+		for (auto& px : r.raw_pixels)
+		{
+			cen_x += px.x;
+			cen_y += px.y;
+		}
+	} else {
+		for (int i = 0; i < r.rle_pixels.size(); i+=2) {
+			double n = r.rle_pixels[i+1].x - r.rle_pixels[i].x + 1;
+			double s = r.rle_pixels[i+1].x + r.rle_pixels[i].x;
+			cen_x += n/2 * s;
+			cen_y += r.rle_pixels[i].y * n;
+		}
 	}
 
 	val_CENTROID_X = cen_x;
@@ -62,10 +71,19 @@ void BasicMorphologyFeatures::calculate(LR& r)
 
 	//==== Basic morphology :: Centroids
 	val_CENTROID_X = val_CENTROID_Y = 0;
-	for (auto& px : r.raw_pixels)
-	{
-		val_CENTROID_X += px.x;
-		val_CENTROID_Y += px.y;
+	if (r.rle_pixels.empty()) {
+		for (auto& px : r.raw_pixels)
+		{
+			val_CENTROID_X += px.x;
+			val_CENTROID_Y += px.y;
+		}
+	} else {
+		for (int i = 0; i < r.rle_pixels.size(); i+=2) {
+			double n = r.rle_pixels[i+1].x - r.rle_pixels[i].x + 1;
+			double s = r.rle_pixels[i+1].x + r.rle_pixels[i].x;
+			val_CENTROID_X += n/2 * s;
+			val_CENTROID_Y += r.rle_pixels[i].y * n;
+		}
 	}
 	val_CENTROID_X /= n;
 	val_CENTROID_Y /= n;

src/nyx/features/contour.cpp

@@ -1,3 +1,6 @@
+#include "pixel.h"
+#include <cstdint>
+#include <vector>
 #define _USE_MATH_DEFINES	// For M_PI, etc.
 #include <cmath>
 #include <memory>
@@ -30,8 +33,169 @@ ContourFeature::ContourFeature() : FeatureMethod("ContourFeature")
 		});
 }
 
+void ContourFeature::buildRegularContourFast(LR& r)
+{
+
+	r.contour.clear();
+
+	// Order the RLE
+	std::vector<Pixel2> rle = r.rle_pixels;
+	std::sort (rle.begin(), rle.end(), compare_locations);
+
+	// Get starting points of each row
+	// Also merge overlapping RLE
+	std::vector<uint32_t> rows;
+	rows.push_back(0);
+	for (uint32_t i = 2; i < rle.size(); i+=2) {
+
+		// Finds row boundary
+		if (rle[i].y != rle[i-1].y) {
+			rows.push_back(i);
+			continue;
+		}
+
+		// Find overlapping RLE and merge them
+		if (rle[i].x-1 == rle[i-1].x) {
+			rle[i-1] = rle[i+1];
+			rle[i].x = rle[i+1].x+1;
+		}
+	}
+	rows.push_back(rle.size());
+
+	// The first row must be entirely on the edge, so store it in the contour
+	uint32_t pix_ind = 0;
+	for (uint32_t i = 0; i < rows[1]; i+=2) {
+		uint32_t n = rle[i+1].x - rle[i].x + 1;
+		for (uint32_t ind = pix_ind; ind < pix_ind + n; ++ind) {
+			Pixel2 p(r.raw_pixels[ind].x,r.raw_pixels[ind].y,r.raw_pixels[ind].inten);
+			r.contour.push_back(p);
+		}
+		pix_ind += n;
+	}
+
+	// Loop over rows from the 2nd row until the 2nd to last row
+	for (uint32_t p = 1; p < rows.size()-2; ++p) {
+		uint32_t prev = rows[p-1];
+		uint32_t prev_end = rows[p];
+		uint32_t next = rows[p+1];
+		uint32_t next_end = rows[p+2];
+
+		// Find the intersection of previous and next rows
+		std::vector<uint32_t> overlap;
+		while (prev < prev_end && next < next_end) {
+			if (rle[prev].x > rle[next+1].x) {
+				next += 2;
+				continue;
+			}
+			if (rle[next].x > rle[prev+1].x) {
+				prev += 2;
+				continue;
+			}
+
+			// If we get to this point, we have found overlap
+			overlap.push_back(std::max(rle[prev].x,rle[next].x));
+			overlap.push_back(std::min(rle[prev+1].x,rle[next+1].x));
+
+			// Determine whether to advance next or previous
+			if (rle[prev+1].x > rle[next+1].x) {
+				next += 2;
+			} else {
+				prev += 2;
+			}
+		}
+
+		// If there is no overlap, all pixels are edge pixels
+		if (overlap.empty()) {
+			for (uint32_t current = rows[p]; current < rows[p+1]; current+=2) {
+				uint32_t n = rle[current+1].x - rle[current].x + 1;
+				for (uint32_t ind = pix_ind; ind < pix_ind + n; ++ind) {
+					Pixel2 p(r.raw_pixels[ind].x,r.raw_pixels[ind].y,r.raw_pixels[ind].inten);
+					r.contour.push_back(p);
+				}
+				pix_ind += n;
+			}
+		} else {
+			// Add pixels to non-overlapping regions
+			uint32_t oind = 0;
+			uint32_t offset = 0;
+			std::vector<uint32_t> nonoverlap;
+			for (uint32_t current = rows[p]; current < rows[p+1]; current+=2) {
+
+				// skip over invalid rle
+				if (rle[current+1].x < rle[current].x) {
+					continue;
+				}
+
+				// The first pixel is always on the edge
+				nonoverlap.push_back(offset);
+
+				uint32_t x_start = rle[current].x;
+
+				while (oind < overlap.size()) {
+
+					if (overlap[oind+1] < rle[current].x) {
+						oind += 2;
+						continue;
+					}
+
+					if (overlap[oind] < rle[current+1].x) {
+						long end = std::max(rle[current].x + 1, (long) overlap[oind]);
+						long start = std::min(rle[current + 1].x, (long) overlap[oind + 1]+1);
+						nonoverlap.push_back(end - x_start + offset);
+						nonoverlap.push_back(start - x_start + offset);
+					}
+
+					if (overlap[oind+1] < rle[current+1].x) {
+						oind += 2;
+					} else {
+						break;
+					}
+				}
+
+				// The last pixel is always on the edge
+				// In the case of a single pixel gap, it might already be there
+				uint32_t last = rle[current+1].x - x_start + 1 + offset;
+				if (nonoverlap[nonoverlap.size()-2] == last) {
+					nonoverlap.pop_back();
+				} else {
+					nonoverlap.push_back(last);
+				}
+
+				offset += rle[current+1].x - rle[current].x + 1;
+			}
+
+			// Add pixels to the contour
+			for (uint32_t ind = 0; ind < nonoverlap.size(); ind += 2) {
+				for (uint32_t pind = nonoverlap[ind]; pind < nonoverlap[ind+1]; ++pind) {
+					Pixel2 p(r.raw_pixels[pix_ind+pind].x, r.raw_pixels[pix_ind+pind].y, r.raw_pixels[pix_ind+pind].inten);
+					r.contour.push_back(p);
+				}
+			}
+			pix_ind += offset;
+		}
+	}
+
+	// The last row must be entirely on the edge, so store it in the contour
+	pix_ind = r.raw_pixels.size();
+	for (uint32_t i = rows[rows.size()-2]; i < rows[rows.size()-1]; i+=2) {
+		uint32_t n = rle[i+1].x - rle[i].x + 1;
+		pix_ind -= n;
+		for (uint32_t ind = pix_ind; ind < pix_ind + n; ++ind) {
+			Pixel2 p(r.raw_pixels[ind].x,r.raw_pixels[ind].y,r.raw_pixels[ind].inten);
+			r.contour.push_back(p);
+		}
+	}
+}
+
 void ContourFeature::buildRegularContour(LR& r)
 {
+
+	// If RLE is present, use it
+	if (!r.rle_pixels.empty()) {
+		buildRegularContourFast(r);
+		return;
+	}
+
 	//==== Pad the image
 
 	int width = r.aux_image_matrix.width,
@@ -179,6 +343,11 @@ void ContourFeature::buildRegularContour(LR& r)
 	}
 }
 
+bool ContourFeature::compare_locations (const Pixel2& lhs, const Pixel2& rhs)
+{
+	return (lhs.y < rhs.y) || (lhs.y == rhs.y && lhs.x < rhs.x);
+}
+
 void ContourFeature::buildWholeSlideContour(LR& r)
 {
 	// Push the 4 slide vertices of dummy intensity 999
@@ -265,8 +434,10 @@ namespace Nyxus
 		if (r.roi_disabled)
 			return;
 
-		//==== Calculate ROI's image matrix
-		r.aux_image_matrix.calculate_from_pixelcloud (r.raw_pixels, r.aabb);
+		//==== Calculate ROI's image matrix if no RLE is present
+		if (r.rle_pixels.empty()) {
+			r.aux_image_matrix.calculate_from_pixelcloud (r.raw_pixels, r.aabb);
+		}
 
 		//==== Contour, ROI perimeter, equivalent circle diameter
 		ContourFeature f;

src/nyx/features/contour.h

@@ -113,6 +113,8 @@ class ContourFeature: public FeatureMethod
 
 private:
 	void buildRegularContour(LR& r);
+	void buildRegularContourFast(LR& r);
+	static bool compare_locations (const Pixel2& lhs, const Pixel2& rhs);
 	void buildWholeSlideContour(LR& r);
 	std::tuple<double, double, double, double> calc_min_max_mean_stddev_intensity (const std::vector<Pixel2> & contour_pixels);
 	double

src/nyx/features/convex_hull_nontriv.cpp

@@ -24,7 +24,11 @@ void ConvexHullFeature::cleanup_instance()
 void ConvexHullFeature::calculate (LR& r)
 {
 	// Build the convex hull
-	build_convex_hull(r.raw_pixels, r.convHull_CH);
+	if (r.rle_pixels.empty()) {
+		build_convex_hull(r.raw_pixels, r.convHull_CH);
+	} else {
+		build_convex_hull(r.rle_pixels, r.convHull_CH);
+	}
 
 	// Calculate related features
 	double s_hull = polygon_area(r.convHull_CH),

src/nyx/features_calc_workflow.cpp

@@ -134,6 +134,27 @@ namespace Nyxus
 		r.intFname = intFile;
 	}
 
+	void init_label_record_2_fast (LR& r, const std::string& segFile, const std::string& intFile, int x1, int x2, int y, int label, PixIntens maxInt, PixIntens minInt, unsigned int tile_index)
+	{
+		// Cache the host tile's index
+		r.host_tiles.insert (tile_index);
+
+		// Initialize basic counters
+		r.aux_area = x2-x1;
+		r.aux_min = minInt;
+		r.aux_max = maxInt;
+		r.aabb.init_y (y);
+		r.aabb.init_x(x1);
+		r.aabb.update_xmax(x2-1);
+
+		// Cache the ROI label
+		r.label = label;
+
+		// File names
+		r.segFname = segFile;
+		r.intFname = intFile;
+	}
+
 	// This function 'digests' the 2nd and the following pixel of a label and updates the label's feature calculation state - the instance of structure 'LR'
 	void update_label_record(LR& lr, int x, int y, int label, PixIntens intensity)
 	{
@@ -163,4 +184,18 @@ namespace Nyxus
 		lr.update_aabb (x,y);
 	}
 
-}
+	void update_label_record_2_fast (LR& lr, int x1, int x2, int y, int label, PixIntens maxInt, PixIntens minInt, unsigned int tile_index)
+	{
+		// Cache the host tile's index
+		lr.host_tiles.insert(tile_index);
+
+		// Initialize basic counters
+		lr.aux_area += x2-x1;
+		lr.aux_min = std::min(lr.aux_min, minInt);
+		lr.aux_max = std::max(lr.aux_max, maxInt);
+		lr.aabb.update_xmin (x1);
+		lr.aabb.update_xmax (x2-1);
+		lr.aabb.update_y (y);
+	}
+
+}