Added RLE fast pass to phase1/2

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_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);
+	}
+
+}

src/nyx/globals.h

@@ -21,9 +21,11 @@ namespace Nyxus
 
 	bool scanFilePairParallel(const std::string& intens_fpath, const std::string& label_fpath, int num_fastloader_threads, int num_sensemaker_threads, int filepair_index, int tot_num_filepairs);
 	std::string getPureFname(const std::string& fpath);
+	int processDataset(const std::vector<std::string>& intensFiles, const std::vector<std::string>& labelFiles, int numFastloaderThreads, int numSensemakerThreads, int numReduceThreads, int min_online_roi_size, bool save2csv, const std::string& csvOutputDir, bool use_fastloop);
 	int processDataset(const std::vector<std::string>& intensFiles, const std::vector<std::string>& labelFiles, int numFastloaderThreads, int numSensemakerThreads, int numReduceThreads, int min_online_roi_size, bool save2csv, const std::string& csvOutputDir);
 	bool gatherRoisMetrics(const std::string& intens_fpath, const std::string& label_fpath, int num_FL_threads);
-	bool processTrivialRois (const std::vector<int>& trivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, int num_FL_threads, size_t memory_limit);
+	bool gatherRoisMetricsFast(const std::string& intens_fpath, const std::string& label_fpath, int num_FL_threads);
+	bool processTrivialRois (const std::vector<int>& trivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, int num_FL_threads, size_t memory_limit, bool use_fastloop);
 	bool processNontrivialRois (const std::vector<int>& nontrivRoiLabels, const std::string& intens_fpath, const std::string& label_fpath, int num_FL_threads);
 	void dump_roi_metrics(const std::string & label_fpath);
 
@@ -49,8 +51,11 @@ namespace Nyxus
 
 	void init_label_record(LR& lr, const std::string& segFile, const std::string& intFile, int x, int y, int label, PixIntens intensity);
 	void init_label_record_2(LR& lr, const std::string& segFile, const std::string& intFile, int x, int y, int label, PixIntens intensity, unsigned int tile_index);
+	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);
 	void update_label_record(LR& lr, int x, int y, int label, PixIntens intensity);
+	void update_label_record_fast(LR& lr, int x1, int x2, int y, int label, PixIntens maxInt, PixIntens minInt);
 	void update_label_record_2(LR& lr, int x, int y, int label, PixIntens intensity, unsigned int tile_index);
+	void update_label_record_2_fast(LR& lr, int x1, int x2, int y, int label, PixIntens maxInt, PixIntens minInt, unsigned int tile_index);
 	void reduce_neighbors(int labels_collision_radius);
 
 	void allocateTrivialRoisBuffers(const std::vector<int>& Pending);
@@ -69,6 +74,7 @@ namespace Nyxus
 	/// @param intensity -- pixel's intensity
 	/// @param tile_index -- index of pixel's tile in the image
 	void feed_pixel_2_metrics(int x, int y, PixIntens intensity, int label, unsigned int tile_index);
+	void feed_pixel_2_metrics_fast(int x1, int x2, int y, PixIntens maxInt, PixIntens minInt, int label, unsigned int tile_index);
 
 	/// @brief Copies a pixel to the ROI's cache. 
 	/// @param x -- x-coordinate of the pixel in the image

src/nyx/phase1_fastloop.cpp

@@ -0,0 +1,140 @@
+#include <cstdint>
+#include <popcntintrin.h>
+#include <ratio>
+#include <smmintrin.h>
+#include <string>
+#include <sstream>
+#include <vector>
+#include <map>
+#include <array>
+#ifdef WITH_PYTHON_H
+#include <pybind11/pybind11.h>
+#endif
+#include "environment.h"
+#include "globals.h"
+#include "helpers/timing.h"
+#include "rle.hpp"
+
+namespace Nyxus
+{
+	bool gatherRoisMetricsFast (const std::string& intens_fpath, const std::string& label_fpath, int num_FL_threads)
+	{
+
+		int lvl = 0, // Pyramid level
+			lyr = 0; //	Layer
+
+		// Read the tiff. The image loader is put in the open state in processDataset()
+		size_t nth = theImLoader.get_num_tiles_hor(),
+			ntv = theImLoader.get_num_tiles_vert(),
+			fw = theImLoader.get_tile_width(),
+			th = theImLoader.get_tile_height(),
+			tw = theImLoader.get_tile_width(),
+			tileSize = theImLoader.get_tile_size(),
+			fullwidth = theImLoader.get_full_width(),
+			fullheight = theImLoader.get_full_height();
+
+		int cnt = 1;
+		for (unsigned int row = 0; row < nth; row++)
+			for (unsigned int col = 0; col < ntv; col++)
+			{
+				// Fetch the tile 
+				bool ok = theImLoader.load_tile(row, col);
+				if (!ok)
+				{
+					std::stringstream ss;
+					ss << "Error fetching tile row=" << row << " col=" << col;
+					#ifdef WITH_PYTHON_H
+						throw ss.str();
+					#endif	
+					std::cerr << ss.str() << "\n";
+					return false;
+				}
+
+				// Get ahold of tile's pixel buffer
+				auto tileIdx = row * nth + col;
+				auto dataI = theImLoader.get_int_tile_buffer(),
+					dataL = theImLoader.get_seg_tile_buffer();
+
+				// Iterate pixels
+                int y = 0;
+				int y_max = (fullheight < y+th) ? fullheight - row*th : th;
+				while (y < y_max)
+				{
+
+                    size_t i = y*tw;
+
+                    int x = 0;
+					int x_max = (fullwidth < x+tw) ? fullwidth - col*tw : tw;
+
+
+					int x_stream = 32 * (x_max / 32);
+
+					// Compress row to an RLE stream
+					RLEStream<uint32_t, uint16_t> stream = rle_encode_long_stream_32(dataL.data()+i,x_stream);
+
+					// Loop over row objects in RLE stream
+					unsigned int ind = 0;
+					while (ind+1 < stream.offsets.size()) {
+						uint16_t x1 = stream.offsets[ind++];
+						auto label = dataL[i+x1];
+
+						// If not background, store the object
+						if (label > 0) {
+							// Collapse all the labels to one if single-ROI mde is requested
+							if (theEnvironment.singleROI)
+								label = 1;
+
+							uint16_t x2 = stream.offsets[ind];
+							auto minInt = dataI[i+x1];
+							auto maxInt = dataI[i+x1];
+
+							// Find the min and max intensity
+							for (int xi=i+x1; xi<i+x2; xi++) {
+								if (maxInt < dataI[xi]) {
+									maxInt = dataI[xi];
+								} else if (minInt > dataI[xi]) {
+									minInt = dataI[xi];
+								}
+
+							}
+
+							feed_pixel_2_metrics_fast (x1, x2, y, maxInt, minInt, label, tileIdx); // Updates 'uniqueLabels' and 'roiData'
+						}
+					}
+
+					i += x_stream;
+					x = x_stream;
+					while (x < x_max) {
+
+                        // Skip non-mask pixels
+                        auto label = dataL[i];
+                        if (label == 0) {
+                            ++x, ++i;
+                            continue;
+                        }
+
+                        // Collapse all the labels to one if single-ROI mde is requested
+                        if (theEnvironment.singleROI)
+                            label = 1;
+
+                        // Update pixel's ROI metrics
+                        feed_pixel_2_metrics (x, y, dataI[i], label, tileIdx); // Updates 'uniqueLabels' and 'roiData'
+                        ++x, ++i;
+                    }
+                    ++y;
+				}
+
+#ifdef WITH_PYTHON_H
+				if (PyErr_CheckSignals() != 0)
+					throw pybind11::error_already_set();
+#endif
+
+				// Show stayalive progress info
+				if (cnt++ % 4 == 0)
+					std::cout << "\t" << int((row * nth + col) * 100 / float(nth * ntv) * 100) / 100. << "%\t" << uniqueLabels.size() << " ROIs" << "\n";
+			}
+
+		return true;
+	}
+
+}