text: /*
  * File:   ROI2D.cpp
  * Author: justin
  *
  * Created on February 11, 2015, 12:02 AM
  */
  
  #include "ROI2D.h"
  
  namespace ncorr {
  
  	// Additional Constructors ---------------------------------------------------//    
  	ROI2D::ROI2D(Array2D<bool> mask, difference_type cutoff) : mask_ptr(std::make_shared<Array2D<bool>>(std::move(mask))), regions_ptr(std::make_shared<std::vector<region>>()) {
  		if (cutoff < 0) {
  			throw std::invalid_argument("Attempted to form ROI2D with cutoff of: " + std::to_string(cutoff) +
  				". cutoff must be an integer of value 0 or greater.");
  		}
  
  		// Given an input of a mask, form ROI2D ----------------------------------//    
  		// Form the nlinfos first
  		auto nlinfos_pair = ROI2D::region_nlinfo::form_nlinfos(*this->mask_ptr, cutoff);
  
  		// Cycle over each nlinfo and create a corresponding region_boundary
  		Array2D<bool> mask_buf(this->mask_ptr->height(), this->mask_ptr->width());
  		for (const auto &nlinfo : nlinfos_pair.first) {
  			auto boundary = nlinfo.to_boundary(mask_buf);
  			this->regions_ptr->emplace_back(std::move(nlinfo), std::move(boundary));
  		}
  
  		if (nlinfos_pair.second) {
  			// If regions were removed, update the mask
  			draw_mask();
  		}
  
  		// Set points
  		set_points();
  	}
  
  	ROI2D::ROI2D(region_nlinfo nlinfo, difference_type h, difference_type w) : mask_ptr(std::make_shared<Array2D<bool>>(h, w)), regions_ptr(std::make_shared<std::vector<region>>()) {
  		// Given an input of an nlinfo, form ROI2D. .
  
  		// Set regions - can use mask_ptr as buffer since it's empty at this point
  		auto boundary = nlinfo.to_boundary(*this->mask_ptr);
  		this->regions_ptr->emplace_back(std::move(nlinfo), std::move(boundary));
  
  		// Draw mask
  		draw_mask();
  
  		// Set points
  		set_points();
  	}
  
  	ROI2D::ROI2D(region_boundary boundary, difference_type h, difference_type w) : mask_ptr(std::make_shared<Array2D<bool>>(h, w)), regions_ptr(std::make_shared<std::vector<region>>()) {
  		// Given an input of a boundary, form ROI2D.
  
  		// Set region - can use mask_ptr as buffer since it's empty at this point
  		auto nlinfo = boundary.to_nlinfo(*this->mask_ptr);
  		this->regions_ptr->emplace_back(std::move(nlinfo), std::move(boundary));
  
  		// Draw mask
  		draw_mask();
  
  		// Set points
  		set_points();
  	}
  
  	ROI2D::ROI2D(std::vector<region_nlinfo> nlinfos, difference_type h, difference_type w) : mask_ptr(std::make_shared<Array2D<bool>>(h, w)), regions_ptr(std::make_shared<std::vector<region>>()) {
  		// Given input of nlinfos, form ROI2D.
  
  		// Set region(s) - can use mask_ptr as buffer since it's empty at this point
  		for (const auto &nlinfo : nlinfos) {
  			auto boundary = nlinfo.to_boundary(*this->mask_ptr);
  			this->regions_ptr->emplace_back(std::move(nlinfo), std::move(boundary));
  		}
  
  		// Draw mask
  		draw_mask();
  
  		// Set points
  		set_points();
  	}
  
  	ROI2D::ROI2D(std::vector<region_boundary> boundaries, difference_type h, difference_type w) : mask_ptr(std::make_shared<Array2D<bool>>(h, w)), regions_ptr(std::make_shared<std::vector<region>>()) {
  		// Given an input of boundaries, form ROI2D.
  
  		// Set region(s) - can use mask_ptr as buffer since it's empty at this point
  		for (const auto &boundary : boundaries) {
  			auto nlinfo = boundary.to_nlinfo(*this->mask_ptr);
  			this->regions_ptr->emplace_back(std::move(nlinfo), std::move(boundary));
  		}
  
  		// Draw mask
  		draw_mask();
  
  		// Set points
  		set_points();
  	}
  
  	// Static factory methods ----------------------------------------------------//
  
  	// Note that 'simple' ROIs require:
  	//  1) Contains 1 region/nlinfo
  	//  2) nlinfo must have two nodes per column
  	//  3) Mask must have odd width
  	//  4) Must contain its center point. 
  	ROI2D ROI2D::simple_circle(difference_type r) {
  		difference_type length = 2 * r + 1;
  		double r_squared = std::pow(r, 2);
  
  		Array2D<bool> mask(length, length);
  		for (difference_type p2 = 0; p2 < length; ++p2) {
  			double h_squared = std::pow(p2 - r, 2);
  			difference_type np_top = std::ceil(-std::sqrt(r_squared - h_squared)) + r;
  			difference_type np_bottom = std::floor(std::sqrt(r_squared - h_squared)) + r;
  			for (difference_type p1 = np_top; p1 <= np_bottom; ++p1) {
  				mask(p1, p2) = true;
  			}
  		}
  
  		return ROI2D(std::move(mask));
  	}
  
  	ROI2D ROI2D::simple_square(difference_type r) {
  		Array2D<bool> mask(2 * r + 1, 2 * r + 1, true);
  
  		return ROI2D(std::move(mask));
  	}
  
  	ROI2D ROI2D::load(std::ifstream &is) {
  		// Form empty ROI2D then fill in values in accordance to how they are saved
  		ROI2D roi;
  
  		// Load mask
  		roi.mask_ptr = std::make_shared<Array2D<bool>>(Array2D<bool>::load(is));
  
  		// Load regions
  		difference_type num_regions = 0;
  		is.read(reinterpret_cast<char*>(&num_regions), std::streamsize(sizeof(difference_type)));
  		roi.regions_ptr = std::make_shared<std::vector<ROI2D::region>>(num_regions);
  		for (auto &reg : *roi.regions_ptr) {
  			reg = ROI2D::region::load(is);
  		}
  
  		// Load points
  		is.read(reinterpret_cast<char*>(&roi.points), std::streamsize(sizeof(difference_type)));
  
  		return roi;
  	}
  
  	// Operators Interface -------------------------------------------------------//  
  	std::ostream& operator<<(std::ostream &os, const ROI2D &roi) {
  		os << "Mask :" << '\n';
  		os << *roi.mask_ptr;
  
  		for (const auto &region : *roi.regions_ptr) {
  			os << region << '\n';
  		}
  
  		os << "Total points: " << roi.points;
  
  		return os;
  	}
  
  	bool isequal(const ROI2D &roi1, const ROI2D &roi2) {
  		typedef ROI2D::difference_type                              difference_type;
  
  		if (roi1.points == roi2.points && isequal(*roi1.mask_ptr, *roi2.mask_ptr) && roi1.regions_ptr->size() == roi2.regions_ptr->size()) {
  			for (difference_type region_idx = 0; region_idx < difference_type(roi1.regions_ptr->size()); ++region_idx) {
  				if (!isequal((*roi1.regions_ptr)[region_idx], (*roi2.regions_ptr)[region_idx])) {
  					return false;
  				}
  			}
  
  			// At this point, the points, masks, and all the regions in roi1 and roi2
  			// are the same
  			return true;
  		}
  
  		return false;
  	}
  
  	void save(const ROI2D &roi, std::ofstream &os) {
  		typedef ROI2D::difference_type                              difference_type;
  
  		// Save mask -> regions -> points
  		save(*roi.mask_ptr, os);
  
  		difference_type num_regions = roi.regions_ptr->size();
  		os.write(reinterpret_cast<const char*>(&num_regions), std::streamsize(sizeof(difference_type)));
  		for (const auto &reg : *roi.regions_ptr) {
  			save(reg, os);
  		}
  
  		os.write(reinterpret_cast<const char*>(&roi.points), std::streamsize(sizeof(difference_type)));
  	}
  
  	// Access --------------------------------------------------------------------//
  	const ROI2D::region_nlinfo& ROI2D::get_nlinfo(difference_type region_idx) const {
  		chk_region_idx_in_bounds(region_idx, "get_nlinfo()");
  
  		return (*regions_ptr)[region_idx].nlinfo;
  	}
  
  	const ROI2D::region_boundary& ROI2D::get_boundary(difference_type region_idx) const {
  		chk_region_idx_in_bounds(region_idx, "get_boundary()");
  
  		return (*regions_ptr)[region_idx].boundary;
  	}
  
  	std::pair<ROI2D::difference_type, ROI2D::difference_type> ROI2D::get_region_idx(difference_type p1, difference_type p2) const {
  		// .first refers to idx of region - will be -1 if it does not exist
  		// .second refers to idx of nodepair which contains p2 - will be -1 if it does not exists
  		if (!mask_ptr->in_bounds(p1, p2)) {
  			return{ -1,-1 };
  		}
  
  		for (difference_type region_idx = 0; region_idx < difference_type(regions_ptr->size()); ++region_idx) {
  			const auto &nlinfo = (*regions_ptr)[region_idx].nlinfo;
  			difference_type nl_idx = p2 - nlinfo.left_nl;
  			if (nlinfo.noderange.in_bounds(nl_idx)) {
  				for (difference_type np_idx = 0; np_idx < nlinfo.noderange(nl_idx); np_idx += 2) {
  					difference_type np_top = nlinfo.nodelist(np_idx, nl_idx);
  					difference_type np_bottom = nlinfo.nodelist(np_idx + 1, nl_idx);
  					if (p1 < np_top) {
  						break; // p1 comes before top of nodepair
  					}
  					if (p1 <= np_bottom) {
  						return{ region_idx, np_idx }; // p1 is contained in this nodepair
  					}
  				}
  			}
  		}
  		// No nlinfo contained p1 and p2
  		return{ -1,-1 };
  	}
  
  	// Arithmetic operations ------------------------------------------------------//
  	ROI2D ROI2D::reduce(difference_type scalefactor) const {
  		if (scalefactor < 1) {
  			throw std::invalid_argument("Attempted to reduce ROI2D with scalefactor of: " + std::to_string(scalefactor) +
  				". scalefactor must be an integer of value 1 or greater.");
  		}
  
  		// Reduction is done by reducing nlinfos. It is possible for nlinfo to 
  		// disappear after reduction, but an empty nlinfo is kept in order to 
  		// maintain a 1 to 1 correspondence of regions. If nlinfo becomes 
  		// non-contiguous after reduction, reduce() returns the largest contiguous 
  		// region.
  		std::vector<region_nlinfo> nlinfos_reduced;
  		Array2D<bool> mask_buf(height(), width());
  		for (difference_type region_idx = 0; region_idx < difference_type(regions_ptr->size()); ++region_idx) {
  			nlinfos_reduced.push_back((*regions_ptr)[region_idx].nlinfo.reduce(scalefactor, mask_buf));
  		}
  
  		return ROI2D(std::move(nlinfos_reduced), std::ceil(double(height()) / scalefactor), std::ceil(double(width()) / scalefactor));
  	}
  
  	ROI2D ROI2D::form_union(const Array2D<bool> &mask) const {
  		if (mask.height() != height() || mask.width() != width()) {
  			throw std::invalid_argument("Attempted to form union of ROI2D of size: " + size_2D_string() +
  				" with mask of size : " + mask.size_2D_string() + ". Mask and ROI2D must have the same size.");
  		}
  
  		// Union is done by finding new nodepair values depending on the values of 
  		// the mask between the original nodepair values. It is possible for nlinfo to 
  		// disappear after union, but an empty nlinfo is kept in order to maintain a 
  		// 1 to 1 correspondence of regions. If nlinfo becomes non-contiguous after 
  		// union, form_union() returns the largest contiguous region.
  		std::vector<region_nlinfo> nlinfos_unioned;
  		Array2D<bool> mask_buf(height(), width());
  		for (difference_type region_idx = 0; region_idx < difference_type(regions_ptr->size()); ++region_idx) {
  			nlinfos_unioned.push_back((*regions_ptr)[region_idx].nlinfo.form_union(mask, mask_buf));
  		}
  
  		return ROI2D(std::move(nlinfos_unioned), height(), width());
  	}
  
  	// Incrementor ---------------------------------------------------------------//
  	ROI2D::incrementor ROI2D::begin_inc() const {
  		return (regions_ptr->empty() ? incrementor(*this, 0, {}) : incrementor(*this, 0, regions_ptr->front().nlinfo.begin_inc()));
  	}
  
  	ROI2D::incrementor ROI2D::end_inc() const {
  		return (regions_ptr->empty() ? incrementor(*this, 0, {}) : incrementor(*this, regions_ptr->size(), regions_ptr->back().nlinfo.end_inc()));
  	}
  
  	// contig_subregion_generator ------------------------------------------------//
  	ROI2D::contig_subregion_generator ROI2D::get_contig_subregion_generator(SUBREGION subregion_type, difference_type r) const {
  		return contig_subregion_generator(*this, subregion_type, r);
  	}
  
  	// Utility methods -----------------------------------------------------------//
  	void ROI2D::set_points() {
  		points = 0;
  		for (const auto &region : *regions_ptr) {
  			points += region.nlinfo.points;
  		}
  	}
  
  	void ROI2D::draw_mask() {
  		(*this->mask_ptr)() = false;
  		for (const auto &region : *regions_ptr) {
  			fill(*this->mask_ptr, region.nlinfo, true);
  		}
  	}
  
  	// Checks --------------------------------------------------------------------//
  	void ROI2D::chk_region_idx_in_bounds(difference_type idx, const std::string &op) const {
  		if (idx < 0 || idx >= difference_type(regions_ptr->size())) {
  			throw std::invalid_argument("Attempted to access region " + std::to_string(idx) + " for " + op + " operation, but there are only " + std::to_string(regions_ptr->size()) + " regions.");
  		}
  	}
  
  	// Static factory methods ----------------------------------------------------//
  	ROI2D::region_nlinfo ROI2D::region_nlinfo::load(std::ifstream &is) {
  		// Form empty nlinfo then fill in values in accordance to how they are saved
  		region_nlinfo nlinfo;
  
  		// Load bounds
  		is.read(reinterpret_cast<char*>(&nlinfo.top), std::streamsize(sizeof(difference_type)));
  		is.read(reinterpret_cast<char*>(&nlinfo.bottom), std::streamsize(sizeof(difference_type)));
  		is.read(reinterpret_cast<char*>(&nlinfo.left), std::streamsize(sizeof(difference_type)));
  		is.read(reinterpret_cast<char*>(&nlinfo.right), std::streamsize(sizeof(difference_type)));
  		is.read(reinterpret_cast<char*>(&nlinfo.left_nl), std::streamsize(sizeof(difference_type)));
  		is.read(reinterpret_cast<char*>(&nlinfo.right_nl), std::streamsize(sizeof(difference_type)));
  
  		// Load nodelist and noderange
  		nlinfo.nodelist = Array2D<difference_type>::load(is);
  		nlinfo.noderange = Array2D<difference_type>::load(is);
  
  		// Load points
  		is.read(reinterpret_cast<char*>(&nlinfo.points), std::streamsize(sizeof(difference_type)));
  
  		return nlinfo;
  	}
  
  	// Operations interface ------------------------------------------------------//
  	std::ostream& operator<<(std::ostream &os, const ROI2D::region_nlinfo &nlinfo) {
  		os << "nlinfo: " << '\n';
  		os << "Top bound of nlinfo: " << std::to_string(nlinfo.top) << "." << '\n';
  		os << "Bottom bound of nlinfo: " << std::to_string(nlinfo.bottom) << "." << '\n';
  		os << "Left bound of nlinfo: " << std::to_string(nlinfo.left) << "." << '\n';
  		os << "Right bound of nlinfo: " << std::to_string(nlinfo.right) << "." << '\n';
  		os << "Left position of nodelist: " << std::to_string(nlinfo.left_nl) << "." << '\n';
  		os << "Right position of nodelist: " << std::to_string(nlinfo.right_nl) << "." << '\n';
  
  		os << '\n' << "nodelist: " << '\n';
  		os << nlinfo.nodelist << '\n';
  		os << "noderange: " << '\n';
  		os << nlinfo.noderange << '\n';
  
  		os << '\n' << "nlinfo points: " << std::to_string(nlinfo.points) << "." << '\n';
  
  		return os;
  	}
  
  	bool isequal(const ROI2D::region_nlinfo &nlinfo1, const ROI2D::region_nlinfo &nlinfo2) {
  		return nlinfo1.top == nlinfo2.top &&
  			nlinfo1.bottom == nlinfo2.bottom &&
  			nlinfo1.left == nlinfo2.left &&
  			nlinfo1.right == nlinfo2.right &&
  			nlinfo1.left_nl == nlinfo2.left_nl &&
  			nlinfo1.right_nl == nlinfo2.right_nl &&
  			isequal(nlinfo1.nodelist, nlinfo2.nodelist) &&
  			isequal(nlinfo1.noderange, nlinfo2.noderange) &&
  			nlinfo1.points == nlinfo2.points;
  	}
  
  	void save(const ROI2D::region_nlinfo &nlinfo, std::ofstream &os) {
  		typedef ROI2D::difference_type                              difference_type;
  
  		// Save bounds -> nodelist -> noderange -> points
  		os.write(reinterpret_cast<const char*>(&nlinfo.top), std::streamsize(sizeof(difference_type)));
  		os.write(reinterpret_cast<const char*>(&nlinfo.bottom), std::streamsize(sizeof(difference_type)));
  		os.write(reinterpret_cast<const char*>(&nlinfo.left), std::streamsize(sizeof(difference_type)));
  		os.write(reinterpret_cast<const char*>(&nlinfo.right), std::streamsize(sizeof(difference_type)));
  		os.write(reinterpret_cast<const char*>(&nlinfo.left_nl), std::streamsize(sizeof(difference_type)));
  		os.write(reinterpret_cast<const char*>(&nlinfo.right_nl), std::streamsize(sizeof(difference_type)));
  
  		save(nlinfo.nodelist, os);
  		save(nlinfo.noderange, os);
  
  		os.write(reinterpret_cast<const char*>(&nlinfo.points), std::streamsize(sizeof(difference_type)));
  	}
  
  	// Arithmetic methods --------------------------------------------------------//
  	bool ROI2D::region_nlinfo::in_nlinfo(difference_type p1, difference_type p2) const {
  		difference_type nl_idx = p2 - left_nl;
  		if (noderange.in_bounds(nl_idx)) {
  			for (difference_type np_idx = 0; np_idx < noderange(nl_idx); np_idx += 2) {
  				difference_type np_top = nodelist(np_idx, nl_idx);
  				difference_type np_bottom = nodelist(np_idx + 1, nl_idx);
  				if (p1 < np_top) {
  					return false; // p1 comes before top of nodepair
  				}
  				if (p1 <= np_bottom) {
  					return true;  // p1 is contained in this nodepair
  				}
  			}
  		}
  
  		return false;
  	}
  
  	ROI2D::region_nlinfo& ROI2D::region_nlinfo::shift(difference_type p1_shift, difference_type p2_shift) {
  		// Note that no checking is done for shifting out of bounds so be careful.
  
  		// Shift bounds first    
  		top += p1_shift;
  		bottom += p1_shift;
  		left += p2_shift;
  		right += p2_shift;
  		left_nl += p2_shift;
  		right_nl += p2_shift;
  
  		// Shift nodelist
  		for (difference_type nl_idx = 0; nl_idx < nodelist.width(); ++nl_idx) {
  			for (difference_type np_idx = 0; np_idx < noderange(nl_idx); ++np_idx) {
  				nodelist(np_idx, nl_idx) += p1_shift;
  			}
  		}
  
  		return *this;
  	}
  
  	// Incrementor ---------------------------------------------------------------//
  	ROI2D::region_nlinfo::incrementor ROI2D::region_nlinfo::begin_inc() const {
  		return (empty() ? incrementor(*this, 0, 0, 0) : incrementor(*this, first_pos_idx(), 0, first_pos_p1()));
  	}
  
  	ROI2D::region_nlinfo::incrementor ROI2D::region_nlinfo::end_inc() const {
  		return (empty() ? incrementor(*this, 0, 0, 0) : incrementor(*this, nodelist.width(), noderange(last_pos_idx()), last_pos_p1() + 1));
  	}
  
  	// Static factory methods ----------------------------------------------------//
  	namespace details {
  		void add_interacting_nodes(ROI2D::difference_type np_adj_p2,
  			ROI2D::difference_type np_loaded_top,
  			ROI2D::difference_type np_loaded_bottom,
  			Array2D<std::vector<ROI2D::difference_type>> &overall_nodelist,
  			Array2D<std::vector<bool>> &overall_active_nodepairs,
  			std::stack<ROI2D::difference_type> &queue_np_idx) {
  			typedef ROI2D::difference_type                          difference_type;
  
  			// Make sure adjacent nodepair(s) position is in range.
  			if (overall_nodelist.in_bounds(np_adj_p2)) {
  				// Scans nodes from top to bottom
  				for (difference_type np_adj_idx = 0; np_adj_idx < difference_type(overall_nodelist(np_adj_p2).size()); np_adj_idx += 2) {
  					difference_type np_adj_top = overall_nodelist(np_adj_p2)[np_adj_idx];
  					difference_type np_adj_bottom = overall_nodelist(np_adj_p2)[np_adj_idx + 1];
  					if (np_loaded_bottom < np_adj_top) {
  						return; // top node of adjacent nodepair is below bottom node of loaded nodepair
  					}
  					if (overall_active_nodepairs(np_adj_p2)[np_adj_idx / 2] && np_loaded_top <= np_adj_bottom) {
  						// Inactivate node pair, and then insert into the queue
  						overall_active_nodepairs(np_adj_p2)[np_adj_idx / 2] = false;
  						queue_np_idx.push(np_adj_top);
  						queue_np_idx.push(np_adj_bottom);
  						queue_np_idx.push(np_adj_p2);
  					}
  				}
  			}
  		}
  	}
  
  	std::pair<std::vector<ROI2D::region_nlinfo>, bool> ROI2D::region_nlinfo::form_nlinfos(const Array2D<bool> &mask, difference_type cutoff) {
  		if (cutoff < 0) {
  			throw std::invalid_argument("Attempted to form nlinfos with cutoff of: " + std::to_string(cutoff) +
  				". cutoff must be an integer of value 0 or greater.");
  		}
  
  		// Form overall_nodelist and overall_active_nodepairs --------------------//
  		Array2D<std::vector<difference_type>> overall_nodelist(1, mask.width());
  		Array2D<std::vector<bool>> overall_active_nodepairs(1, mask.width());
  		for (difference_type p2 = 0; p2 < mask.width(); ++p2) {
  			bool in_nodepair = false;
  			for (difference_type p1 = 0; p1 < mask.height(); ++p1) {
  				if (!in_nodepair && mask(p1, p2)) {
  					in_nodepair = true;
  					overall_nodelist(p2).push_back(p1); // sets top node
  				}
  				if (in_nodepair && (!mask(p1, p2) || p1 == mask.height() - 1)) {
  					in_nodepair = false;
  					overall_nodelist(p2).push_back((p1 == mask.height() - 1 && mask(p1, p2)) ? p1 : p1 - 1); // Sets bottom node
  				}
  			}
  
  			// Update overall_active_nodepairs - this keeps track of which node
  			// pairs have been analyzed when doing contiguous separation; set to 
  			// true initially.
  			overall_active_nodepairs(p2).resize(overall_nodelist(p2).size() / 2, true);
  		}
  
  		// Separate regions ------------------------------------------------------//   
  		// Regions are made 4-way contiguous. Scan over columns and separate 
  		// contiguous nodelists in overall_nodelist.
  		std::vector<ROI2D::region_nlinfo> nlinfos; // This will get updated and returned
  		bool removed = false;                      // Keeps track if regions are removed due to "cutoff" parameter
  		for (difference_type p2_sweep = 0; p2_sweep < overall_nodelist.width(); ++p2_sweep) {
  			// Find first active node pair
  			difference_type np_init_idx = -1;
  			for (difference_type np_idx = 0; np_idx < difference_type(overall_nodelist(p2_sweep).size()); np_idx += 2) {
  				if (overall_active_nodepairs(p2_sweep)[np_idx / 2]) {       // Test if nodepair is active
  					overall_active_nodepairs(p2_sweep)[np_idx / 2] = false; // Inactivate node pair
  					np_init_idx = np_idx;                                 // Store nodepair idx
  					break;
  				}
  			}
  
  			// If there are no active node pairs, then continue to next column
  			if (np_init_idx == -1) {
  				continue;
  			}
  
  			// nlinfo_buf to be updated, then inserted into nlinfos
  			ROI2D::region_nlinfo nlinfo_buf(mask.height() - 1,    // Top bound of region                (this gets updated)
  				0,                  // Bottom bound of region             (this gets updated)
  				p2_sweep,           // Left bound of region               (this is correct)
  				0,                  // Right bound of region              (this gets updated)
  				p2_sweep,           // Left position of nodelist          (this is correct)
  				0,                  // Right position of nodelist         (this gets updated)
  				0,                  // h of nodelist                      (this gets updated)  
  				0,                  // w of nodelist                      (this gets updated)     
  				0);                 // Number of points in region         (this gets updated) 
  
  									// Keep track of nodes with separate_nodelist
  			Array2D<std::vector<difference_type>> separate_nodelist(1, mask.width());
  
  			// Initialize queue and enter while loop - exit when queue is empty
  			std::stack<difference_type> queue_np_idx;                       // Holds all nodepairs (along with their index) which need to be processed
  			queue_np_idx.push(overall_nodelist(p2_sweep)[np_init_idx]);     // Top
  			queue_np_idx.push(overall_nodelist(p2_sweep)[np_init_idx + 1]); // Bottom
  			queue_np_idx.push(p2_sweep);                                    // position of nodepair
  			while (!queue_np_idx.empty()) {
  				// Pop nodepair and its position out of queue and compare 
  				// it to adjacent nodepairs (left and right of np_loaded_p2)
  				difference_type np_loaded_p2 = queue_np_idx.top(); queue_np_idx.pop();
  				difference_type np_loaded_bottom = queue_np_idx.top(); queue_np_idx.pop();
  				difference_type np_loaded_top = queue_np_idx.top(); queue_np_idx.pop();
  
  				// Compare to node pairs LEFT. Any node pairs which interact are added to the queue
  				details::add_interacting_nodes(np_loaded_p2 - 1,
  					np_loaded_top,
  					np_loaded_bottom,
  					overall_nodelist,
  					overall_active_nodepairs,
  					queue_np_idx);
  
  				// Compare to node pairs RIGHT. Any node pairs which interact are added to the queue
  				details::add_interacting_nodes(np_loaded_p2 + 1,
  					np_loaded_top,
  					np_loaded_bottom,
  					overall_nodelist,
  					overall_active_nodepairs,
  					queue_np_idx);
  
  				// Update points
  				nlinfo_buf.points += np_loaded_bottom - np_loaded_top + 1;
  
  				// Update bounds - note that "right_nl" and "right" are the same
  				if (np_loaded_top < nlinfo_buf.top) { nlinfo_buf.top = np_loaded_top; }                         // Top
  				if (np_loaded_bottom > nlinfo_buf.bottom) { nlinfo_buf.bottom = np_loaded_bottom; }             // Bottom
  				if (np_loaded_p2 > nlinfo_buf.right) { nlinfo_buf.right_nl = nlinfo_buf.right = np_loaded_p2; } // Right
  
  																												// Insert node pairs and then sort - usually very small so BST isn't
  																												// necessary.
  				separate_nodelist(np_loaded_p2).push_back(np_loaded_top);
  				separate_nodelist(np_loaded_p2).push_back(np_loaded_bottom);
  				std::sort(separate_nodelist(np_loaded_p2).begin(), separate_nodelist(np_loaded_p2).end());
  			}
  
  			// Now finish setting nodelist and noderange for this region.
  			// Find max nodes first so we can use it to set the correct height
  			// for nodelist.
  			difference_type max_nodes = 0;
  			for (const auto &nodes : separate_nodelist) {
  				if (difference_type(nodes.size()) > max_nodes) {
  					max_nodes = nodes.size();
  				}
  			}
  
  			// Set and fill nodelist and noderange
  			nlinfo_buf.nodelist = Array2D<difference_type>(max_nodes, nlinfo_buf.right_nl - nlinfo_buf.left_nl + 1);
  			nlinfo_buf.noderange = Array2D<difference_type>(1, nlinfo_buf.right_nl - nlinfo_buf.left_nl + 1);
  			for (difference_type nl_idx = 0; nl_idx < nlinfo_buf.nodelist.width(); ++nl_idx) {
  				difference_type p2 = nl_idx + nlinfo_buf.left_nl;
  				// noderange:
  				nlinfo_buf.noderange(nl_idx) = separate_nodelist(p2).size();
  				// nodelist:
  				for (difference_type np_idx = 0; np_idx < difference_type(separate_nodelist(p2).size()); ++np_idx) {
  					nlinfo_buf.nodelist(np_idx, nl_idx) = separate_nodelist(p2)[np_idx];
  				}
  			}
  
  			// Subtract one from p2_sweep in order to recheck the column to 
  			// ensure all nodes are deactivated before proceeding
  			--p2_sweep;
  
  			// Make sure number of points in nlinfo is greater than or equal to
  			// the cutoff
  			if (nlinfo_buf.points >= cutoff) {
  				nlinfos.push_back(nlinfo_buf);
  			}
  			else {
  				removed = true; // Parameter lets caller know regions were removed
  			}
  		}
  
  		return{ std::move(nlinfos), removed };
  	}
  
  	// Arithmetic methods --------------------------------------------------------//
  	ROI2D::region_nlinfo ROI2D::region_nlinfo::reduce(difference_type scalefactor, Array2D<bool> &mask_buf) const {
  		if (scalefactor < 1) {
  			throw std::invalid_argument("Attempted to reduce nlinfo with scalefactor of: " + std::to_string(scalefactor) +
  				". scalefactor must be an integer of value 1 or greater.");
  		}
  
  		// This function will return a copy of a reduced nlinfo based on input
  		// scalefactor.
  
  		if (empty()) {
  			// Return empty nlinfo
  			return region_nlinfo();
  		}
  
  		// Only *_nl bounds can be set now. Other bounds must be calculated since 
  		// regions can disappear if a large scalefactor is used. Note that since
  		// nlinfo is nonempty at this point, using the bounds is valid.
  		difference_type top_reduced = std::ceil(double(top) / scalefactor);
  		difference_type bottom_reduced = std::floor(double(bottom) / scalefactor);
  		difference_type left_reduced_nl = std::ceil(double(left_nl) / scalefactor);
  		difference_type right_reduced_nl = std::floor(double(right_nl) / scalefactor);
  		region_nlinfo nlinfo_reduced(bottom_reduced,                              // top               (gets updated)
  			top_reduced,                                 // bottom            (gets updated)
  			right_reduced_nl,                            // left              (gets updated)
  			left_reduced_nl,                             // right             (gets updated)
  			left_reduced_nl,                             // left_nl           (correct)
  			right_reduced_nl,                            // right_nl          (correct)
  			nodelist.height(),                           // nodelist height   (correct)
  			right_reduced_nl - left_reduced_nl + 1,      // nodelist width    (correct)
  			0);                                          // points            (gets updated)
  
  		for (difference_type nl_idx_reduced = 0; nl_idx_reduced < nlinfo_reduced.nodelist.width(); ++nl_idx_reduced) {
  			difference_type nl_idx = (nl_idx_reduced + nlinfo_reduced.left_nl) * scalefactor - left_nl;
  			for (difference_type np_idx = 0; np_idx < noderange(nl_idx); np_idx += 2) {
  				difference_type np_top_reduced = std::ceil(double(nodelist(np_idx, nl_idx)) / scalefactor);
  				difference_type np_bottom_reduced = std::floor(double(nodelist(np_idx + 1, nl_idx)) / scalefactor);
  				if (np_bottom_reduced >= np_top_reduced) {
  					// Update points
  					nlinfo_reduced.points += np_bottom_reduced - np_top_reduced + 1;
  
  					// Update bounds
  					if (np_top_reduced < nlinfo_reduced.top) { nlinfo_reduced.top = np_top_reduced; }
  					if (np_bottom_reduced > nlinfo_reduced.bottom) { nlinfo_reduced.bottom = np_bottom_reduced; }
  					if (nl_idx_reduced + nlinfo_reduced.left_nl < nlinfo_reduced.left) { nlinfo_reduced.left = nl_idx_reduced + nlinfo_reduced.left_nl; }
  					if (nl_idx_reduced + nlinfo_reduced.left_nl > nlinfo_reduced.right) { nlinfo_reduced.right = nl_idx_reduced + nlinfo_reduced.left_nl; }
  
  					// Insert nodepairs
  					nlinfo_reduced.nodelist(nlinfo_reduced.noderange(nl_idx_reduced), nl_idx_reduced) = np_top_reduced;
  					nlinfo_reduced.nodelist(nlinfo_reduced.noderange(nl_idx_reduced) + 1, nl_idx_reduced) = np_bottom_reduced;
  
  					// Update noderange
  					nlinfo_reduced.noderange(nl_idx_reduced) += 2;
  				}
  			}
  		}
  		// Return largest contiguous region since reduced nlinfo may no longer be 
  		// 4-way contiguous
  		return nlinfo_reduced.largest_contig_nlinfo(mask_buf);
  	}
  
  	ROI2D::region_nlinfo ROI2D::region_nlinfo::form_union(const Array2D<bool> &mask, Array2D<bool> &mask_buf) const {
  		if (!mask.same_size(mask_buf)) {
  			throw std::invalid_argument("Attempted to form union with mask of size: " + mask.size_2D_string() +
  				" with mask_buf of size : " + mask_buf.size_2D_string() + ". Mask and mask_buf must have the same size.");
  		}
  
  		// This function will return a copy of a unioned nlinfo based on input mask.
  
  		if (empty()) {
  			// Return empty nlinfo
  			return region_nlinfo();
  		}
  
  		// Only *_nl bounds can be set now. Other bounds must be calculated since 
  		// it is indeterminate how many new nodes will be added.
  		region_nlinfo nlinfo_union(bottom,           // top               (gets updated)
  			top,              // bottom            (gets updated)
  			right_nl,         // left              (gets updated)
  			left_nl,          // right             (gets updated)
  			left_nl,          // left_nl           (correct)
  			right_nl,         // right_nl          (correct)
  			0,                // nodelist height   (gets updated)
  			nodelist.width(), // nodelist width    (correct)
  			0);               // points            (gets updated)
  
  
  							  // Use a vector buffer to hold nodepairs
  		Array2D<std::vector<difference_type>> nodelist_buf_vec(nodelist.width(), 1);
  		difference_type max_height = 0; // will be used to allocate array for nodelist
  
  										// Cycle over nlinfo and find new nodepairs
  		for (difference_type nl_idx = 0; nl_idx < nodelist.width(); ++nl_idx) {
  			difference_type p2 = nl_idx + left_nl;
  			for (difference_type np_idx = 0; np_idx < noderange(nl_idx); np_idx += 2) {
  				difference_type np_top = nodelist(np_idx, nl_idx);
  				difference_type np_bottom = nodelist(np_idx + 1, nl_idx);
  				bool in_nodepair = false;
  				for (difference_type p1 = np_top; p1 <= np_bottom; p1++) {
  					difference_type np_top_new, np_bottom_new;
  					if (!in_nodepair && mask(p1, p2)) {
  						in_nodepair = true;
  						np_top_new = p1;
  						nodelist_buf_vec(nl_idx).push_back(np_top_new);
  					}
  					if (in_nodepair && (!mask(p1, p2) || p1 == np_bottom)) {
  						in_nodepair = false;
  						np_bottom_new = ((p1 == np_bottom && mask(p1, p2)) ? p1 : p1 - 1);
  
  						// update points
  						nlinfo_union.points += np_bottom_new - np_top_new + 1;
  
  						// update bounds
  						if (np_top_new < nlinfo_union.top) { nlinfo_union.top = np_top_new; }
  						if (np_bottom_new > nlinfo_union.bottom) { nlinfo_union.bottom = np_bottom_new; }
  						if (nl_idx + nlinfo_union.left_nl < nlinfo_union.left) { nlinfo_union.left = nl_idx + nlinfo_union.left_nl; }
  						if (nl_idx + nlinfo_union.left_nl > nlinfo_union.right) { nlinfo_union.right = nl_idx + nlinfo_union.left_nl; }
  
  						// insert bottom of nodepair
  						nodelist_buf_vec(nl_idx).push_back(np_bottom_new);
  
  						// Update noderange
  						nlinfo_union.noderange(nl_idx) += 2;
  
  						// Update max height
  						if (nlinfo_union.noderange(nl_idx) > max_height) { max_height = nlinfo_union.noderange(nl_idx); }
  					}
  				}
  			}
  		}
  
  		// Store nodepairs in nodelist_buf_vec in nlinfo_union
  		nlinfo_union.nodelist = Array2D<difference_type>(max_height, nlinfo_union.nodelist.width());
  		for (difference_type nl_idx = 0; nl_idx < nodelist.width(); ++nl_idx) {
  			for (difference_type np_idx = 0; np_idx < difference_type(nodelist_buf_vec(nl_idx).size()); ++np_idx) {
  				nlinfo_union.nodelist(np_idx, nl_idx) = nodelist_buf_vec(nl_idx)[np_idx];
  			}
  		}
  		// Return largest contiguous region since unioned nlinfo may no longer be 
  		// 4-way contiguous
  		return nlinfo_union.largest_contig_nlinfo(mask_buf);
  	}
  
  	namespace details {
  		Array2D<double> calc_boundary(const Array2D<bool> &mask, double init_p1, double init_p2, ROI2D::difference_type init_direc = 0) {
  			typedef ROI2D::difference_type                          difference_type;
  
  			if (init_direc < 0 || init_direc > 3) {
  				throw std::invalid_argument("Initial direction of: " + std::to_string(init_direc) + " was provided to calc_boundary function. Initial direction must be between (inclusive) 0 and 3.");
  			}
  
  			// Note that the returned boundary is "open" (does not contain duplicate 
  			// initial and end points) and is the "outer boundary" of the mask (i.e. 
  			// if the mask contains 1 true pixel, a boundary containing the 4 points 
  			// of the corners of the pixel is returned). This will return a boundary 
  			// for an 8-way connected region. init_p1 and init_p2 need to be the 
  			// location WRT the outer boundary (i.e (1.5, 2.5) for a pixel located 
  			// at (2, 3) assuming its the left-most top pixel). If this function is
  			// called correctly, boundary is guaranteed to contain at least 4 points.
  
  			// direc is based on:
  			//
  			//        1
  			//        |
  			//   2 ------- 0   
  			//        |
  			//        3
  			//
  			// Initial direction of 0 is correct for inputs containing the left-most
  			// top pixel in an 8-way connected region.
  
  			// Initialize p1, p2, and direc, which keep track of the current point 
  			// location and the direction taken to get there, respectively.
  			double p1 = init_p1;
  			double p2 = init_p2;
  			difference_type direc = init_direc;
  
  			// Initialize boundary
  			std::vector<double> boundary_vec = { p1,p2 };
  			bool initpoint_encountered = false;
  			while (!initpoint_encountered) {
  				// Get initial new direction
  				direc = (direc + 3) % 4;
  
  				// Cycle clockwise to find new direction
  				for (difference_type it = 0; it < 4; ++it) {
  					// Increment point based on direc
  					double p1_inc, p2_inc;
  					difference_type p1_inc_mask, p2_inc_mask;
  					switch (direc) {
  					case 0: p1_inc = p1; p2_inc = p2 + 1;
  						p1_inc_mask = std::round(p1_inc - 0.5); p2_inc_mask = std::round(p2_inc - 0.5);
  						break;
  					case 1: p1_inc = p1 - 1; p2_inc = p2;
  						p1_inc_mask = std::round(p1_inc + 0.5); p2_inc_mask = std::round(p2_inc - 0.5);
  						break;
  					case 2: p1_inc = p1; p2_inc = p2 - 1;
  						p1_inc_mask = std::round(p1_inc + 0.5); p2_inc_mask = std::round(p2_inc + 0.5);
  						break;
  					case 3: p1_inc = p1 + 1; p2_inc = p2;
  						p1_inc_mask = std::round(p1_inc - 0.5); p2_inc_mask = std::round(p2_inc + 0.5);
  						break;
  					}
  
  					if (mask.in_bounds(p1_inc_mask, p2_inc_mask) && mask(p1_inc_mask, p2_inc_mask)) {
  						// This is the next point in the boundary
  						if (boundary_vec[0] != p1_inc || boundary_vec[1] != p2_inc) {
  							// Insert point into boundary and then set current point.
  							boundary_vec.push_back(p1_inc); boundary_vec.push_back(p2_inc);
  							p1 = p1_inc; p2 = p2_inc;
  						}
  						else {
  							initpoint_encountered = true;
  						}
  						// break from for loop
  						break;
  					}
  					else {
  						// Increment clockwise one position, then re-check
  						direc = (direc + 1) % 4;
  					}
  
  					if (it == 3) {
  						// This means an empty mask was provided and is a programmer 
  						// error
  						throw std::invalid_argument("Invalid initial point or empty mask was provided to calc_boundary() function.");
  					}
  				}
  			}
  
  			// Convert vec_boundary to Array2D
  			Array2D<double> boundary(boundary_vec.size() / 2, 2);
  			for (difference_type idx = 0; idx < difference_type(boundary_vec.size()); idx += 2) {
  				boundary(idx / 2, 0) = boundary_vec[idx];
  				boundary(idx / 2, 1) = boundary_vec[idx + 1];
  			}
  
  			return boundary;
  		}
  	}
  
  	ROI2D::region_nlinfo ROI2D::region_nlinfo::largest_contig_nlinfo(Array2D<bool> &mask_buf) const {
  		if (empty()) {
  			// nlinfo is empty, so return empty nlinfo.
  			return ROI2D::region_nlinfo();
  		}
  
  		// Clear mask buffer and then fill with nlinfo
  		mask_buf() = false;
  		fill(mask_buf, *this, true);
  
  		// Form new nlinfos and then return the one with the largest number of points
  		auto nlinfos_pair = ROI2D::region_nlinfo::form_nlinfos(mask_buf);
  		// nlinfo is non-empty at this point, so max element is guaranteed to exist.
  		return *std::max_element(nlinfos_pair.first.begin(), nlinfos_pair.first.end(), [](const ROI2D::region_nlinfo &a, const ROI2D::region_nlinfo &b) { return a.points < b.points; });
  	}
  
  	ROI2D::region_boundary ROI2D::region_nlinfo::to_boundary(Array2D<bool> &mask_buf) const {
  		// Converts input nlinfo to a region_boundary. Note that no bounds checking 
  		// nor contiguity checking is done for input nlinfo's, so they must be 
  		// correct.        
  
  		if (empty()) {
  			// If nlinfo is empty, return an empty boundary.
  			return ROI2D::region_boundary();
  		}
  
  		// Form 1 "add" boundary and varying "sub" boundaries per nlinfo.
  		// Create mask_buf and then fill nlinfo.
  		mask_buf() = false;
  		fill(mask_buf, *this, true);
  
  		// Get the "add" boundary ------------------------------------------------//
  		// Use the left-most top pixel in the region; must use the "outer" position, 
  		// so subtract by 0.5.
  		auto add_boundary = details::calc_boundary(mask_buf, first_pos_p1() - 0.5, first_pos_p2() - 0.5);
  
  		// Find the "sub" boundaries ---------------------------------------------//  
  		// Use mask_inv_buf to keep track of which sub boundaries have been analyzed.
  		Array2D<bool> mask_inv_buf(mask_buf.height(), mask_buf.width());
  		fill(mask_inv_buf, add_boundary, true);
  		mask_inv_buf = mask_inv_buf & ~mask_buf; // This shows "holes" as white regions
  
  												 // Cycle over nlinfo and check for holes in columns with more than 1 nodepair
  		std::vector<Array2D<double>> sub_boundaries;
  		for (difference_type nl_idx = 0; nl_idx < nodelist.width(); ++nl_idx) {
  			if (noderange(nl_idx) > 2) {
  				difference_type p2 = nl_idx + left_nl;
  				// Test one pixel below bottom node of every node pair except 
  				// for the last node pair, since this will be the bottom of the ROI
  				for (difference_type np_idx = 0; np_idx < noderange(nl_idx) - 2; np_idx += 2) {
  					difference_type np_bottom = nodelist(np_idx + 1, nl_idx);
  					if (mask_inv_buf(np_bottom + 1, p2)) {
  						// This is a boundary which hasn't been analyzed yet
  						sub_boundaries.push_back(details::calc_boundary(mask_inv_buf, np_bottom + 0.5, p2 - 0.5));
  
  						// Fill in this boundary so it does not get analyzed again
  						// Note that sub_boundary cannot have a hole, so it is
  						// safe to use boundary fill. Also, boundary is guaranteed
  						// to fill region covered by nlinfo.
  						fill(mask_inv_buf, sub_boundaries.back(), false);
  					}
  				}
  			}
  		}
  
  		return ROI2D::region_boundary(std::move(add_boundary), std::move(sub_boundaries));
  	}
  
  	// Utility -------------------------------------------------------------------//
  	void ROI2D::region_nlinfo::chk_nonempty_op(const std::string &op) const {
  		if (this->empty()) {
  			throw std::invalid_argument("Attempted to use " + op + " operator on empty nlinfo" +
  				". nlinfo must be nonempty.");
  		}
  	}
  
  	// Static factory methods ----------------------------------------------------//
  	ROI2D::region_boundary ROI2D::region_boundary::load(std::ifstream &is) {
  		// Form empty boundary then fill in values in accordance to how they are saved
  		region_boundary boundary;
  
  		// Load add boundary
  		boundary.add = Array2D<double>::load(is);
  
  		// Allocate new sub boundaries and then load them
  		difference_type num_sub_boundaries = 0;
  		is.read(reinterpret_cast<char*>(&num_sub_boundaries), std::streamsize(sizeof(difference_type)));
  		boundary.sub.resize(num_sub_boundaries);
  		for (auto &sub_boundary : boundary.sub) {
  			sub_boundary = Array2D<double>::load(is);
  		}
  
  		return boundary;
  	}
  
  	// Operations interface ------------------------------------------------------//
  	std::ostream& operator<<(std::ostream &os, const ROI2D::region_boundary &boundary) {
  		typedef ROI2D::difference_type                              difference_type;
  
  		os << "Add Boundary: " << '\n' << boundary.add;
  		for (difference_type boundary_idx = 0; boundary_idx < difference_type(boundary.sub.size()); ++boundary_idx) {
  			os << '\n' << "Sub Boundary(" << boundary_idx << ") :" << '\n' << boundary.sub[boundary_idx];
  		}
  
  		return os;
  	}
  
  	bool isequal(const ROI2D::region_boundary &boundary1, const ROI2D::region_boundary &boundary2) {
  		typedef ROI2D::difference_type                              difference_type;
  
  		if (isequal(boundary1.add, boundary2.add) && boundary1.sub.size() == boundary2.sub.size()) {
  			for (difference_type boundary_idx = 0; boundary_idx < difference_type(boundary1.sub.size()); ++boundary_idx) {
  				if (!isequal(boundary1.sub[boundary_idx], boundary2.sub[boundary_idx])) {
  					return false;
  				}
  			}
  			// Add boundary and sub boundaries are the same
  			return true;
  		}
  
  		return false;
  	}
  
  	void save(const ROI2D::region_boundary &boundary, std::ofstream &os) {
  		typedef ROI2D::difference_type                              difference_type;
  
  		// Save add boundary -> vec of sub boundaries
  		save(boundary.add, os);
  
  		// Save length of vector
  		difference_type num_sub_boundaries = boundary.sub.size();
  		os.write(reinterpret_cast<const char*>(&num_sub_boundaries), std::streamsize(sizeof(difference_type)));
  		// Save each sub boundary
  		for (const auto &sub_boundary : boundary.sub) {
  			save(sub_boundary, os);
  		}
  	}
  
  	// Arithmetic methods --------------------------------------------------------//
  	ROI2D::region_nlinfo ROI2D::region_boundary::to_nlinfo(Array2D<bool> &mask_buf) const {
  		// Converts boundary to nlinfo. This will form an nlinfo corresponding to 
  		// the largest contiguous enclosing region of the boundary.
  
  		// Clear mask_buf and then draw boundary
  		mask_buf() = false;
  		fill(mask_buf, add, true);
  
  		//crashed here
  		for (const auto &sub_boundary : sub) {
  			fill(mask_buf, sub_boundary, false);
  		}
  
  		// Get nlinfos for this boundary
  		auto nlinfos_pair = ROI2D::region_nlinfo::form_nlinfos(mask_buf);
  
  		if (nlinfos_pair.first.empty()) {
  			// Boundary is empty, so return empty nlinfo
  			return ROI2D::region_nlinfo();
  		}
  
  		// Its possible for boundary to produce multiple nlinfos if it gets "pinched",
  		// so return the nlinfo with the largest number of points. nlinfo is 
  		// non-empty at this point, so max element is guaranteed to exist.
  		return *std::max_element(nlinfos_pair.first.begin(), nlinfos_pair.first.end(), [](const ROI2D::region_nlinfo &a, const ROI2D::region_nlinfo &b) { return a.points < b.points; });
  	}
  
  	// Static factory methods ----------------------------------------------------//
  	ROI2D::region ROI2D::region::load(std::ifstream &is) {
  		// Form empty region then fill in values in accordance to how they are saved
  		region reg;
  
  		// Load nlinfo
  		reg.nlinfo = ROI2D::region_nlinfo::load(is);
  
  		// Load boundary
  		reg.boundary = ROI2D::region_boundary::load(is);
  
  		return reg;
  	}
  
  	// Operations interface ------------------------------------------------------//
  	std::ostream& operator<<(std::ostream &os, const ROI2D::region &reg) {
  		os << reg.nlinfo;
  		os << reg.boundary;
  
  		return os;
  	}
  
  	void save(const ROI2D::region &reg, std::ofstream &os) {
  		// Save nlinfo -> boundary        
  		save(reg.nlinfo, os);
  		save(reg.boundary, os);
  	}
  
  	namespace details {
  		// Arithmetic methods ----------------------------------------------------//
  		nlinfo_incrementor& nlinfo_incrementor::operator++() {
  			if (p1 == nlinfo_ptr->nodelist(np_idx + 1, nl_idx)) {
  				// incremented position will go beyond bottom node
  				if (np_idx == nlinfo_ptr->noderange(nl_idx) - 2) {
  					// Incremented nodepair index will go beyond nodepairs in this column,
  					// so search columns to find next one with nodepairs
  					for (++nl_idx; nl_idx < nlinfo_ptr->nodelist.width() && nlinfo_ptr->noderange(nl_idx) == 0; ++nl_idx) {}
  					if (nl_idx == nlinfo_ptr->nodelist.width()) {
  						// End has been reached, increment nodepair_idx and p1 to one
  						// beyond the end to match the end incrementor.
  						np_idx += 2;
  						++p1;
  					}
  					else {
  						np_idx = 0;
  						p1 = nlinfo_ptr->nodelist(np_idx, nl_idx); // top node
  					}
  				}
  				else {
  					np_idx += 2;
  					p1 = nlinfo_ptr->nodelist(np_idx, nl_idx); // top node
  				}
  			}
  			else {
  				++p1;
  			}
  
  			return *this;
  		}
  
  		// Additional Constructors -----------------------------------------------//
  		ROI2D_incrementor::ROI2D_incrementor(const ROI2D &roi, difference_type region_idx, const ROI2D::region_nlinfo::incrementor &nlinfo_inc) : roi(roi), region_idx(region_idx), nlinfo_inc(nlinfo_inc) {
  			// If the initial nlinfo_inc is the beginning incrementor for the first
  			// region, increment until the first position is found. This is done for 
  			// convenience so the begin_inc() can be set easily by the caller.
  			if (this->roi.size_regions() > 0 && this->nlinfo_inc == this->roi.get_nlinfo(0).begin_inc() && this->roi.get_nlinfo(0).empty()) {
  				// Set to first non-empty region
  				for (++this->region_idx; this->region_idx < this->roi.size_regions() && this->roi.get_nlinfo(this->region_idx).empty(); ++this->region_idx) {}
  				if (this->region_idx == this->roi.size_regions()) {
  					// All regions are empty - set nlinfo incrementor to the end
  					// incrementor of the last region in order to match the end 
  					// incrementor
  					this->nlinfo_inc = this->roi.get_nlinfo(this->region_idx - 1).end_inc();
  				}
  				else {
  					this->nlinfo_inc = this->roi.get_nlinfo(region_idx).begin_inc();
  				}
  			}
  		}
  
  		// Arithmetic methods ----------------------------------------------------//
  		ROI2D_incrementor& ROI2D_incrementor::operator++() {
  			if (nlinfo_inc.pos_2D().first == roi.get_nlinfo(region_idx).last_pos_p1() && nlinfo_inc.pos_2D().second == roi.get_nlinfo(region_idx).last_pos_p2()) {
  				// Reached last position in this region, increment region_idx until 
  				// nonempty region is found - or until all regions have been checked
  				for (++region_idx; region_idx < roi.size_regions() && roi.get_nlinfo(region_idx).empty(); ++region_idx) {}
  				if (region_idx == roi.size_regions()) {
  					// This is the end, set nlinfo_inc to end_incrementor of last 
  					// region to match the end incrementor.
  					nlinfo_inc = roi.get_nlinfo(region_idx - 1).end_inc();
  				}
  				else {
  					nlinfo_inc = roi.get_nlinfo(region_idx).begin_inc();
  				}
  			}
  			else {
  				++nlinfo_inc;
  			}
  
  			return *this;
  		}
  
  		// Additional Constructors -----------------------------------------------//
  		ROI2D_contig_subregion_generator::ROI2D_contig_subregion_generator(const ROI2D &roi, SUBREGION subregion_type, difference_type r) : roi(roi), r(r) {
  			// Find max height of nodelist in roi so nodelist can be set.
  			difference_type max_height = 0;
  			for (difference_type region_idx = 0; region_idx < this->roi.size_regions(); ++region_idx) {
  				if (this->roi.get_nlinfo(region_idx).nodelist.height() > max_height) {
  					max_height = this->roi.get_nlinfo(region_idx).nodelist.height();
  				}
  			}
  
  			// Set buffers. These buffers are the upperbound for the sizes, so they 
  			// never need to be resized when forming contiguous region. 
  			this->nlinfo_subregion.nodelist = Array2D<difference_type>(max_height, 2 * this->r + 1);
  			this->nlinfo_subregion.noderange = Array2D<difference_type>(1, 2 * this->r + 1);
  			this->active_nodepairs = Array2D<bool>(max_height / 2, 2 * this->r + 1);
  
  			// Note that simple_nlinfo must be 'simple' (look at the requirements
  			// for a 'simple' ROI at the simple_* static factory method definitions).
  			switch (subregion_type) {
  			case SUBREGION::CIRCLE:
  				this->nlinfo_simple = ROI2D::simple_circle(r).get_nlinfo(0); // Safe since ROI is simple
  				break;
  			case SUBREGION::SQUARE:
  				this->nlinfo_simple = ROI2D::simple_square(r).get_nlinfo(0); // Safe since ROI is simple
  				break;
  			}
  		}
  
  		// Local function --------------------------------------------------------//
  		void simple_contig_subregion_add_interacting_nodes(ROI2D::difference_type np_adj_p2,
  			ROI2D::difference_type np_loaded_top,
  			ROI2D::difference_type np_loaded_bottom,
  			const ROI2D::ROI2D::region_nlinfo &nlinfo_roi,
  			const ROI2D::ROI2D::region_nlinfo &nlinfo_simple,
  			Array2D<bool> &active_nodepairs,
  			std::stack<ROI2D::difference_type> &queue_np_idx) {
  			typedef ROI2D::difference_type                          difference_type;
  
  			// Make sure idx's are in range of both nlinfo_roi and nlinfo_simple
  			difference_type nl_adj_idx = np_adj_p2 - nlinfo_roi.left_nl;
  			difference_type nl_simple_idx = np_adj_p2 - nlinfo_simple.left_nl;
  			if (nlinfo_roi.noderange.in_bounds(nl_adj_idx) && nlinfo_simple.noderange.in_bounds(nl_simple_idx)) {
  				// Get simple nodepair at adjacent position - safe since nlinfo is 'simple'
  				difference_type np_simple_top = nlinfo_simple.nodelist(0, nl_simple_idx);
  				difference_type np_simple_bottom = nlinfo_simple.nodelist(1, nl_simple_idx);
  				// Scans nodes from top to bottom
  				for (difference_type np_adj_idx = 0; np_adj_idx < nlinfo_roi.noderange(nl_adj_idx); np_adj_idx += 2) {
  					difference_type np_adj_top = nlinfo_roi.nodelist(np_adj_idx, nl_adj_idx);
  					difference_type np_adj_bottom = nlinfo_roi.nodelist(np_adj_idx + 1, nl_adj_idx);
  					if (np_loaded_bottom < np_adj_top || np_loaded_bottom < np_simple_top) {
  						return; // top node of adjacent nodepair or simple nodepair is below bottom node of loaded nodepair
  					}
  					if (active_nodepairs(np_adj_idx / 2, nl_simple_idx) && np_loaded_top <= np_adj_bottom && np_loaded_top <= np_simple_bottom) {
  						// At this point, loaded nodepair interacts with both adjacent
  						// nodepair and simple nodepair. Take the union of the 
  						// intersection. Note that it is possible for adjacent nodepair
  						// and simple nodepair to be disjoint, which results in a 
  						// "flipped" nodepair, so test for it.
  						difference_type np_top = std::max(np_adj_top, np_simple_top);
  						difference_type np_bottom = std::min(np_adj_bottom, np_simple_bottom);
  						if (np_top <= np_bottom) {
  							// Inactivate node pair, and then insert into the queue.
  							// Since simple_nlinfo only contains two nodes per column, 
  							// a nodepair in nlinfo_roi can only interact with it
  							// once, so deactivating it is safe.
  							active_nodepairs(np_adj_idx / 2, nl_simple_idx) = false;
  							queue_np_idx.push(np_top);
  							queue_np_idx.push(np_bottom);
  							queue_np_idx.push(np_adj_p2);
  						}
  					}
  				}
  			}
  		}
  
  		// Arithmetic methods ----------------------------------------------------//
  		const ROI2D::region_nlinfo& ROI2D_contig_subregion_generator::operator()(difference_type p1, difference_type p2) const {
  			// Clear/initialize values in nlinfo_output - must do this here in case 
  			// empty nlinfo_output is returned
  			nlinfo_subregion.top = roi.height() - 1;           // Gets updated
  			nlinfo_subregion.bottom = 0;                     // Gets updated
  			nlinfo_subregion.left = roi.width() - 1;           // Gets updated
  			nlinfo_subregion.right = 0;                      // Gets updated
  			nlinfo_subregion.left_nl = p2 - r;               // Correct
  			nlinfo_subregion.right_nl = p2 + r;              // Correct
  			nlinfo_subregion.noderange() = 0;                // Gets updated
  			nlinfo_subregion.points = 0;                     // Gets updated
  
  															 // Get region idx containing (p1,p2)
  			auto region_idx_pair = roi.get_region_idx(p1, p2);
  			if (region_idx_pair.first == -1) {
  				// ROI does not contain the x,y coordinate - return the empty nlinfo_output
  				return nlinfo_subregion;
  			}
  
  			// Set active_nodepairs to true
  			active_nodepairs() = true;
  
  			// Get nlinfo corresponding to p1 and p2
  			auto &nlinfo_roi = roi.get_nlinfo(region_idx_pair.first);
  
  			// Shift nlinfo_simple's position in-place
  			nlinfo_simple.shift(p1 - r, p2 - r);
  
  			// Get node pairs containing x and y, take their union with simple 
  			// nodepair, and then add to queue.
  			active_nodepairs(region_idx_pair.second / 2, r) = false; // Inactivate nodepair    
  			std::stack<difference_type> queue_np_idx; // Holds all nodepairs (along with their index) which need to be processed
  			queue_np_idx.push(std::max(nlinfo_simple.nodelist(0, r), nlinfo_roi.nodelist(region_idx_pair.second, p2 - nlinfo_roi.left_nl)));     // Top
  			queue_np_idx.push(std::min(nlinfo_simple.nodelist(1, r), nlinfo_roi.nodelist(region_idx_pair.second + 1, p2 - nlinfo_roi.left_nl))); // Bottom
  			queue_np_idx.push(p2);                                                                                                               // idx
  			while (!queue_np_idx.empty()) {
  				// Pop nodepair and its position out of queue and compare 
  				// it to adjacent nodepairs (left and right of np_loaded_p2)
  				difference_type np_loaded_p2 = queue_np_idx.top(); queue_np_idx.pop();
  				difference_type np_loaded_bottom = queue_np_idx.top(); queue_np_idx.pop();
  				difference_type np_loaded_top = queue_np_idx.top(); queue_np_idx.pop();
  
  				// Compare to node pairs LEFT. Any node pairs which interact are added to the queue
  				details::simple_contig_subregion_add_interacting_nodes(np_loaded_p2 - 1,
  					np_loaded_top,
  					np_loaded_bottom,
  					nlinfo_roi,
  					nlinfo_simple,
  					active_nodepairs,
  					queue_np_idx);
  
  				// Compare to node pairs RIGHT. Any node pairs which interact are added to the queue
  				details::simple_contig_subregion_add_interacting_nodes(np_loaded_p2 + 1,
  					np_loaded_top,
  					np_loaded_bottom,
  					nlinfo_roi,
  					nlinfo_simple,
  					active_nodepairs,
  					queue_np_idx);
  
  				// Update points
  				nlinfo_subregion.points += np_loaded_bottom - np_loaded_top + 1;
  
  				// Update bounds
  				if (np_loaded_top < nlinfo_subregion.top) { nlinfo_subregion.top = np_loaded_top; }               // Top
  				if (np_loaded_bottom > nlinfo_subregion.bottom) { nlinfo_subregion.bottom = np_loaded_bottom; }   // Bottom
  				if (np_loaded_p2 < nlinfo_subregion.left) { nlinfo_subregion.left = np_loaded_p2; }               // Left
  				if (np_loaded_p2 > nlinfo_subregion.right) { nlinfo_subregion.right = np_loaded_p2; }             // Right
  
  																												  // Insert node pairs and then sort 
  				difference_type nl_output_idx = np_loaded_p2 - p2 + r;
  				nlinfo_subregion.nodelist(nlinfo_subregion.noderange(nl_output_idx), nl_output_idx) = np_loaded_top;
  				nlinfo_subregion.nodelist(nlinfo_subregion.noderange(nl_output_idx) + 1, nl_output_idx) = np_loaded_bottom;
  				std::sort(&nlinfo_subregion.nodelist(0, nl_output_idx), &nlinfo_subregion.nodelist(0, nl_output_idx) + nlinfo_subregion.noderange(nl_output_idx) + 2);
  
  				// Update noderange
  				nlinfo_subregion.noderange(nl_output_idx) += 2;
  			}
  
  			// Shift nlinfo_simple's position back to original place
  			nlinfo_simple.shift(r - p1, r - p2);
  
  			return nlinfo_subregion;
  		}
  	}
  
  }