text: //--------------------------------------------------------------------------------------------------
  // Linear time Maximally Stable Extremal Regions implementation as described in D. Nistér and
  // H. Stewénius. Linear Time Maximally Stable Extremal Regions. Proceedings of the European
  // Conference on Computer Vision (ECCV), 2008.
  // 
  // Copyright (c) 2011 Idiap Research Institute, http://www.idiap.ch/.
  // Written by Charles Dubout <charles.dubout@idiap.ch>/.
  // 
  // MSER is free software: you can redistribute it and/or modify it under the terms of the GNU
  // General Public License version 3 as published by the Free Software Foundation.
  // 
  // MSER is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
  // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
  // Public License for more details.
  // 
  // You should have received a copy of the GNU General Public License along with MSER. If not, see
  // <http://www.gnu.org/licenses/>.
  //--------------------------------------------------------------------------------------------------
  
  #include "mser.h"
  
  #include <algorithm>
  #include <cassert>
  #include <limits>
  
  using namespace std;
  
  MSER::MSER(bool eight, int delta, double minArea, double maxArea, double maxVariation,
  		   double minDiversity) : eight_(eight), delta_(delta), minArea_(minArea),
  maxArea_(maxArea), maxVariation_(maxVariation), minDiversity_(minDiversity), pool_(256),
  poolIndex_(0)
  {
  	// Parameter check
  	assert(delta > 0);
  	assert(minArea >= 0.0);
  	assert(maxArea <= 1.0);
  	assert(minArea < maxArea);
  	assert(maxVariation > 0.0);
  	assert(minDiversity >= 0.0);
  }
  
  void MSER::operator()(const uint8_t * bits, int width, int height, vector<Region> & regions)
  {
  	// 1. Clear the accessible pixel mask, the heap of boundary pixels and the component stack. Push
  	// a dummy-component onto the stack, with grey-level higher than any allowed in the image.
  	vector<bool> accessible(width * height);
  	vector<int> boundaryPixels[256];
  	int priority = 256;
  	vector<Region *> regionStack;
  	
  	regionStack.push_back(new (&pool_[poolIndex_++]) Region);
  	
  	// 2. Make the source pixel (with its first edge) the current pixel, mark it as accessible and
  	// store the grey-level of it in the variable current level.
  	int curPixel = 0;
  	int curEdge = 0;
  	int curLevel = bits[0];
  	accessible[0] = true;
  	
  	// 3. Push an empty component with current level onto the component stack.
  step_3:
  	regionStack.push_back(new (&pool_[poolIndex_++]) Region(curLevel, curPixel));
  	
  	if (poolIndex_ == pool_.size())
  		doublePool(regionStack);
  	
  	// 4. Explore the remaining edges to the neighbors of the current pixel, in order, as follows:
  	// For each neighbor, check if the neighbor is already accessible. If it is not, mark it as
  	// accessible and retrieve its grey-level. If the grey-level is not lower than the current one,
  	// push it onto the heap of boundary pixels. If on the other hand the grey-level is lower than
  	// the current one, enter the current pixel back into the queue of boundary pixels for later
  	// processing (with the next edge number), consider the new pixel and its grey-level and go to 3.
  	for (;;) {
  		const int x = curPixel % width;
  		const int y = curPixel / width;
  		
  		for (; curEdge < (eight_ ? 8 : 4); ++curEdge) {
  			int neighborPixel = curPixel;
  			
  			if (eight_) {
  				switch (curEdge) {
  					case 0: if (x < width - 1) neighborPixel = curPixel + 1; break;
  					case 1: if ((x < width - 1) && (y > 0)) neighborPixel = curPixel - width + 1; break;
  					case 2: if (y > 0) neighborPixel = curPixel - width; break;
  					case 3: if ((x > 0) && (y > 0)) neighborPixel = curPixel - width - 1; break;
  					case 4: if (x > 0) neighborPixel = curPixel - 1; break;
  					case 5: if ((x > 0) && (y < height - 1)) neighborPixel = curPixel + width - 1; break;
  					case 6: if (y < height - 1) neighborPixel = curPixel + width; break;
  					default: if ((x < width - 1) && (y < height - 1)) neighborPixel = curPixel + width + 1; break;
  				}
  			}
  			else {
  				switch (curEdge) {
  					case 0: if (x < width - 1) neighborPixel = curPixel + 1; break;
  					case 1: if (y < height - 1) neighborPixel = curPixel + width; break;
  					case 2: if (x > 0) neighborPixel = curPixel - 1; break;
  					default: if (y > 0) neighborPixel = curPixel - width; break;
  				}
  			}
  			
  			if (neighborPixel != curPixel && !accessible[neighborPixel]) {
  				const int neighborLevel = bits[neighborPixel];
  				accessible[neighborPixel] = true;
  				
  				if (neighborLevel >= curLevel) {
  					boundaryPixels[neighborLevel].push_back(neighborPixel << 4);
  					
  					if (neighborLevel < priority)
  						priority = neighborLevel;
  				}
  				else {
  					boundaryPixels[curLevel].push_back((curPixel << 4) | (curEdge + 1));
  					
  					if (curLevel < priority)
  						priority = curLevel;
  					
  					curPixel = neighborPixel;
  					curEdge = 0;
  					curLevel = neighborLevel;
  					
  					goto step_3;
  				}
  			}
  		}
  		
  		// 5. Accumulate the current pixel to the component at the top of the stack (water
  		// saturates the current pixel).
  		regionStack.back()->accumulate(x, y);
  		
  		// 6. Pop the heap of boundary pixels. If the heap is empty, we are done. If the returned
  		// pixel is at the same grey-level as the previous, go to 4.
  		if (priority == 256) {
  			regionStack.back()->detect(delta_, minArea_ * width * height,
  									   maxArea_ * width * height, maxVariation_, minDiversity_,
  									   regions);
  			poolIndex_ = 0;
  			return;
  		}
  		
  		curPixel = boundaryPixels[priority].back() >> 4;
  		curEdge = boundaryPixels[priority].back() & 15;
  		
  		boundaryPixels[priority].pop_back();
  		
  		while (boundaryPixels[priority].empty() && (priority < 256))
  			++priority;
  		
  		const int newPixelGreyLevel = bits[curPixel];
  		
  		if (newPixelGreyLevel != curLevel) {
  			curLevel = newPixelGreyLevel;
  			
  			// 7. The returned pixel is at a higher grey-level, so we must now process
  			// all components on the component stack until we reach the higher
  			// grey-level. This is done with the processStack sub-routine, see below.
  			// Then go to 4.
  			processStack(newPixelGreyLevel, curPixel, regionStack);
  		}
  	}
  }
  
  void MSER::processStack(int newPixelGreyLevel, int pixel, vector<Region *> & regionStack)
  {
  	// 1. Process component on the top of the stack. The next grey-level is the minimum of
  	// newPixelGreyLevel and the grey-level for the second component on the stack.
  	do {
  		Region * top = regionStack.back();
  		
  		regionStack.pop_back();
  		
  		// 2. If newPixelGreyLevel is smaller than the grey-level on the second component on the
  		// stack, set the top of stack grey-level to newPixelGreyLevel and return from sub-routine
  		// (This occurs when the new pixel is at a grey-level for which there is not yet a component
  		// instantiated, so we let the top of stack be that level by just changing its grey-level.
  		if (newPixelGreyLevel < regionStack.back()->level_) {
  			regionStack.push_back(new (&pool_[poolIndex_++]) Region(newPixelGreyLevel, pixel));
  			
  			if (poolIndex_ == pool_.size())
  				top = reinterpret_cast<Region *>(reinterpret_cast<char *>(top) +
  												 doublePool(regionStack));
  			
  			regionStack.back()->merge(top);
  			
  			return;
  		}
  		
  		// 3. Remove the top of stack and merge it into the second component on stack as follows:
  		// Add the first and second moment accumulators together and/or join the pixel lists.
  		// Either merge the histories of the components, or take the history from the winner. Note
  		// here that the top of stack should be considered one ’time-step’ back, so its current
  		// size is part of the history. Therefore the top of stack would be the winner if its
  		// current size is larger than the previous size of second on stack.
  		regionStack.back()->merge(top);
  	}
  	// 4. If(newPixelGreyLevel>top of stack grey-level) go to 1.
  	while (newPixelGreyLevel > regionStack.back()->level_);
  }
  
  ptrdiff_t MSER::doublePool(vector<Region *> & regionStack)
  {
  	assert(!pool_.empty()); // Cannot double the size of an empty pool
  	
  	vector<Region> newPool(pool_.size() * 2);
  	copy(pool_.begin(), pool_.end(), newPool.begin());
  	
  	// Cast to char in case the two pointers do not share the same alignment
  	const ptrdiff_t offset = reinterpret_cast<char *>(&newPool[0]) -
  								  reinterpret_cast<char *>(&pool_[0]);
  	
  	for (size_t i = 0; i < pool_.size(); ++i) {
  		if (newPool[i].parent_)
  			newPool[i].parent_ =
  				reinterpret_cast<Region *>(reinterpret_cast<char *>(newPool[i].parent_) + offset);
  		
  		if (newPool[i].child_)
  			newPool[i].child_ =
  				reinterpret_cast<Region *>(reinterpret_cast<char *>(newPool[i].child_) + offset);
  		
  		if (newPool[i].next_)
  			newPool[i].next_ =
  				reinterpret_cast<Region *>(reinterpret_cast<char *>(newPool[i].next_) + offset);
  	}
  	
  	for (size_t i = 0; i < regionStack.size(); ++i)
  		regionStack[i] =
  			reinterpret_cast<Region *>(reinterpret_cast<char *>(regionStack[i]) + offset);
  	
  	pool_.swap(newPool);
  	
  	return offset;
  }
  