نحوه عملکرد الگوریتم MSER

Question

سلام.آیا کسی با این الگوریتم کار کرده و از جزئیات پیاده سازی آن اطلاع داره یکم مطالعه کردم خیلی سخت بود.

Answer 1

#pragma once

#include "MserType.h"
#include <opencv2/core/core.hpp>

enum class ErState { none, max_var, max_area, min_area, min_div };
class EcvMser {
public:
EcvMser();
~EcvMser();
void newObject(int ndims, int const* dims);
//EcvMserFilt*  vl_mser_new(int ndims, int const* dims);

void  process(uchar  *im);

void oneOrder(int d);
void twoOrder(int d);
void integrate();
void integrateMoments();
void computeCentralMmoments();
void  ellFit();


uchar delta();
void setDelta(uchar x);
double minDiversity();
void setMinDiversity(double x);
EcvMserStats const* stats();
double maxArea();
void setMaxArea(double x);
double minArea();
void setMinArea(double x);
double max_variation();
void setMaxVariation(double x);
uint const * regions();
uint regionsNum();
float const * ell();
uint ellDof();
uint ellNum();
void getHistogram(uchar* img,int pix_count,uint* hist);
void getCDF(uint* hist);
void setPixOrder( uint* hist);
void initForest();
void addAsRoot(int idx);
void computeRegion();

bool isExtermal(int idx);
void reg2ExtermalTree(int idx, int er_count);
void createExtermalTree();
void linkExtermalRegions2Tree();
void computeVaribility();
void SelectMaximallyStableBranches();
void furtherFiltering();
void saveResult();

//void vl_mser_ell_fit();
//void computeRegion2();

private:

void drawImages();
const cv::Point n_coords[9] = { cv::Point(-1,-1),cv::Point(0,-1) ,cv::Point(1,-1),
cv::Point(-1,0),cv::Point(0,0) ,cv::Point(1,0) ,
cv::Point(-1,1),cv::Point(0,1) ,cv::Point(1,1) };

cv::Mat parent_img_, shortcut_img_, area_img_,  height_img_,join_img_,order_img_,extrmal_img_,temp_img_;
uchar* img_;
cv::Size image_size_;
cv::Rect image_rect_;

/** @name Image data and meta data @internal */
/*@{*/
int                ndims_;   /**< number of dimensions                    */
int               *dims_;    /**< dimensions                              */
int                pixel_count_;     /**< number of image elements (pixels)       */
int               *subs_;    /**< N-dimensional subscript                 */
int               *dsubs_;   /**< another subscript                       */
int               *strides_; /**< strides to move in image data           */
/*@}*/

uint           *pix_order_;    /**< pixel ordering                          */
uint           *joins_;   /**< sequence of join ops                    */
int                join_count_;  /**< number of join ops                      */

/** @name Regions */
/*@{*/
EcvMserReg         *r_;       /**< basic regions                           */
EcvMserExtrReg     *e_tree_;      /**< extremal tree                           */
uint           *mer_;     /**< maximally stable extremal regions       */
int                er_count_;     /**< number of extremal regions              */
int                mer_count_;    /**< number of maximally stable extr. reg.   */
int                rer_;     /**< size of er buffer                       */
int                rmer_;    /**< size of mer buffer                      */
/*@}*/

/** @name Ellipsoids fitting */
/*@{*/
float             *acc_;     /**< moment accumulator.                    */
float             *ell_;     /**< ellipsoids list.                       */
int                rell_;    /**< size of ell buffer                     */
int                nell_;    /**< number of ellipsoids extracted         */
int                dof_;     /**< number of dof of ellipsoids.           */

/*@}*/

/** @name Configuration */
/*@{*/
bool   verbose_;          /**< be verbose                             */
int       delta_;            /**< delta filter parameter                 */
double    max_area_;         /**< badness test parameter                 */
double    min_area_;         /**< badness test parameter                 */
double    max_variation_;    /**< badness test parameter                 */
double    min_diversity_;    /**< minimum diversity                      */
/*@}*/

EcvMserStats stats_;          /** run statistic                           */
cv::Mat view_img_;
uint upPass(uint& idx);
void downPass(uint& idx);
uint climb( uint idx);
void adv(int ndims, int const *dims, int *subs);
cv::Point index2Coord(int idx);
bool getNeighbor(int pos, const cv::Point& pnt, int index,  int& n_index);
void processNeighbor(bool good, int n_idx, int idx, int r_idx);
void rootIdxChild(int r_idx, int nr_idx, int hgt, int n_hgt);
void rootIdxParent(int r_idx, int nr_idx, int hgt, int n_hgt);
void removeDuplicate( int i, int min_div, int& dup_count, ErState& state);

};

 

 

Answer 2

من قبلا یک کد ansi c از الگوریتم MSER را به c++ تبدیل کرده بودم .کلیت عملکردش به صورت زیر هستش :

 

محاسبه هیستوگرام

محاسبه توزیع تجمعی هیستوگرام

مرتب سازی پیکسل ها با توزیع تجمعی :

برای مرتب سازی پیکسل ها استفاده کنیم که به آن bucketsort هم اطلاق میشه. تصویر را برحسب شدت نور به صورت صعودی مرتب می کنیم نقاطی که دارای شدت نور  یکسان هستند به ترتیب از بالا به پایین و چپ به راست مرتب می شوند که از طریق آن آرایه اولویت  متناظر با تصویر ساخته می شود. در این روش پیکسل ها به صورت یک بعدی از انتهای تصویر پردازش می شوند در هر مرحله با توجه به مقدار هر پیکسل به ستون تجمعی مورد نظر رفته و از مقدار آن ستون یک واحد کم کرده و مقدار بدست آمده  را در آرایه اولویت در نقطه متناظر قرار می دهیم.

ایجاد و اتصال نواحی :

در هر مرحله به ازای هر پیکسل از idx متناظر در آرایه اولویت  استفاده می کنیم.

کلیه پیکسل های تصویر را پردازش می کنیم:

در هر مرحله :

• افزودن ناحیه:

به ازای هر idx یک ناحیه اضافه می کنیم مساحت و عمق آن را یک قرار می دهیم . همچنین parent و shortcut آن را idx قرار میدیم.

• بررسی تمامی نقاط همسایه :

هر پیکسل 8 همسایه داره .

شرایط همسایگی :

1. مختصات آنها داخل تصویر باشد.
2. از قبل ناحیه ای برای آنها ایجاد شده باشد

 شرایط  ترکیب نواحی :

تابع climb جهت محاسبه idx ریشه هر ناحیه مورد استفاده قرار می گیره که در ادامه به توضیح آن خواهیم پرداخت.

ریشه ناحیه جاری و ریشه ناحیه همسایه را با climb محاسبه می کنیم .

اگر idx  ریشه ها نابرابر باشند آنگاه می توانیم عملیات ترکیب ریشه ها را انجام بدیم.

اگر idx ریشه همسایه با idx ناحیه جاری یکسان باشه و مقدار جاری از مقدار ریشه همسایه کمتر باشه آنگاه ریشه ناحیه جاری، فرزند ریشه همسایه میشه اگر هم شرط بالا برقرار نباشه ریشه همسایه ،فرزند ریشه ناحیه جاری میشه.

در عملیات ترکیب ، ناحیه ای ریشه ی که فرزند میشه، idx ناحیه ریشه مقابل را به عنوان parent و shortcut خود درنظر می گیرد.

مساحت فرزند به والد اضافه میشه

محاسبه عمق جدید والد به صورت زیر است :

ماکزیمم  مقدار بین عمق فرزند یک واحد بیشتر و عمق والد

 

شمارنده ای به نام njoins داریم که از صفر شروع میشه. که به ازای هر عمیلات ترکیب یک واحد به آن اضافه میشه .

آرایه ای به نام joins داریم  که idx فرزندان را در خود نگهداری می کنه .idx مورد نظر را در موقعیت njoins قرار می دهیم.

 

تابع climb :

برای بدست آوردن idx ریشه یک ناحیه از این تابع استفاده می کنیم .

این تابع از 2 بخش تشکیل شده است:

1.  از طریق shortcut  ناحیه به صورت متوالی به سمت بالا می رویم ه تا به ریشه برسیم.

در هر مرحله پس از استفاده از shortcut مقدار جدید آن را با ناحیه قبلی(پایینی) مقداردهی می کنیم.

1.  shortcut را به حالت اولیه آن بر می گردانیم.

 

ایجاد نواحی extermal :

تا اینجا نواحی بدست آمدند در این مرحله  قصد داریم نواحی extermal  را ایجاد کنیم.

کلیه نواحی را پردازش می کنیم اگر ناحیه شرایط extermal داشته باشد یک ناحیه جدید به آرایه نواحی extermal اضافه می کنیم .

 

 اتصال ناحیه به ناحیه extermal  :

 برای اتصال شاخص ناحیه extermal را در  shortcut ناحیه قرار می دهیم.

اگر هم ناحیه ای دارای ناحیه extermal نباشد shortcut آن را void قرار می دهیم.

 

مقدار دهی ناحیه extermal:

 index آن را idx ناحیه  قرار می دهیم.

shortcut آن شاخص مربوط به ناحیه extermal قرار می دهیم.

از مقدار پیکسل و مساحت ناحیه ، برای ناحیه extermal استفاده می کنیم.

 

اتصال نواحی extermal به صورت درخت :

 

 

کلیه نواحی extermal موجود را پردازش می کنیم در هر مرحله :

 

Answer 3

یافتن والد غیر void یک ناحیه :

با توجه به فیلد index در نواحی extermal که اشاره دارد به ناحیه idx ، به صورت متوالی والد ها ناحیه مورد نظر را پردازش می کنیم تا زمانیکه  shortcut نواحی والد برابر void باشد به جستجو ادامه می دهیم.

پس از پیدا شدن ناحیه والد ، shortcut آن را در parent ناحیه extermal قرار می دهیم.

 shortcut ناحیه extermal را با شاخص فعلی قرار میدیم.

 

محاسبه variability :

کلیه نواحی extermal موجود را پردازش می کنیم در هر مرحله :

برای هر ناحیه extermal مقدار variation آن را محاسبه می کنیم برای اینکار نیاز به یافتن ناحیه top هر ناحیه extermal داریم .

Variation = (area_top –area)/area

چطور به ناحیه top برسیم؟

برای هر ناحیه extermal مقدار top_value را به صورت زیر محاسبه می کنیم.

Top_value=value+delta

   delta  :پارامتر ورودی الگوریتم.

سپس از طریق مقدار parent  هر ناحیه extermal به صورت متوالی عملیات را تکرار می کنیم تا parent را پیدا کنیم که دارای بزرگترین مقدار باشه این مقدار باید کوچکتر از top_value باشه یا اینکه parent ی باشه که خودش top باشه(اصطلاحاً خودش parent خودشه)

 

 

تعیین max_stable برای نواحی extermal  :

کلیه نواحی extermal را max_stable در نظر می گیریم که در ادامه با شرایطی آن را اصلاح می کنیم.

اگر مقدار والد یک ناحیه extermal از خودش بیشتر باشه این ناحیه extermal را max_stable در نظر می گیریم در غیر اینصورت باید به variation ها نگاه کنیم هر کدام که variation بیشتری داشته باشد دیگر max_stable نیست.

 

فیلتر کردن نواحی extermalی که  max_stable :

نواحی extermalی که شرایط زیر را داشته باشند max_stable باقی می مانند.

1. نواحی extermal باید مساحت آنها در بازه مشخصی باشد تا max_stable باشند دو  پارامتر max_area و min_area داریم که این بازه را مشخص می کنند.
2. چنانچه variation ناحیه extermal ی بزرگتر از max_vari باشد آن ناحیه  max_stable نیست.
3. اگر diversity یک ناحیه extermal کوچکتر از پارامتر min_div باشد آن ناحیه دیگر max_stable نیست.

محاسبه diversity  به استفاده از parent به شرح زیر است :

Diversity= (p_area – area)/p_area

یافتن parent  :

برای هر ناحیه extermal  به صورت متوالی با فیلد parent به دنبال والدی می گردیم که max_stable باشد .

 

 

 

Answer 4

#include "mser2.h"
/** -------------------------------------------------------------------
** @brief Create a new MSER filter
**
** Initializes a new MSER filter for images of the specified
** dimensions. Images are @a ndims -dimensional arrays of dimensions
** @a dims.
**
** @param ndims number of dimensions.
** @param dims  dimensions.
**/
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

EcvMser::EcvMser()
{


}


EcvMser::~EcvMser()
{
	free(acc_);
	free(ell_);

	free(e_tree_);
	free(r_);
	free(joins_);
	free(pix_order_);

	free(strides_);
	free(dsubs_);
	free(subs_);
	free(dims_);

	free(mer_);


}

void EcvMser::newObject(int ndims, int const* dims)
{
	image_size_ = cv::Size(dims[0], dims[1]);
	image_rect_ = cv::Rect(cv::Point(0, 0), image_size_);

	int *strides, k;


	ndims_ = ndims;
	dims_ = (int*)malloc(sizeof(int) * ndims);
	subs_ = (int*)malloc(sizeof(int) * ndims);
	dsubs_ = (int*)malloc(sizeof(int) * ndims);
	strides_ = (int*)malloc(sizeof(int) * ndims);

	for (k = 0; k < ndims_; ++k) {
		dims_[k] = dims[k];
	}

	/* compute strides to move into the N-dimensional image array */
	strides_[0] = 1;
	for (k = 1; k < ndims_; ++k) {
		strides_[k] = strides_[k - 1] * dims_[k - 1];
	}


	pixel_count_ = strides_[ndims_ - 1] * dims_[ndims_ - 1];

	/* dof of ellipsoids */
	dof_ = ndims_ * (ndims_ + 1) / 2 + ndims_;

	/* more buffers */
	pix_order_ = (uint*)malloc(sizeof(uint)   * pixel_count_);
	joins_ = (uint*)malloc(sizeof(uint)   * pixel_count_);
	r_ = (EcvMserReg*)malloc(sizeof(EcvMserReg) * pixel_count_);

	e_tree_ = 0;
	rer_ = 0;
	mer_ = 0;
	rmer_ = 0;
	ell_ = 0;
	rell_ = 0;


	delta_ = 5;
	max_area_ = 0.75;
	min_area_ = 3.0 / pixel_count_;
	max_variation_ = 0.25;
	min_diversity_ = 0.2;




}

void EcvMser::process(uchar *img)
{
	img_ = img;
	nell_ = 0;

	uint buckets[ECV_MSER_PIX_MAXVAL];
	getHistogram(img_, pixel_count_, buckets);
	getCDF(buckets);
	setPixOrder(buckets);

	initForest();
	computeRegion();
	createExtermalTree();
	linkExtermalRegions2Tree();
	computeVaribility();
	SelectMaximallyStableBranches();
	furtherFiltering();
	saveResult();



}


void EcvMser::oneOrder(int d)
{
	for (int index = 0; index < pixel_count_; ++index) {
		acc_[index] = subs_[d];
		adv(ndims_, dims_, subs_);
	}

}


void EcvMser::twoOrder(int d)
{

	// map the dof d to a second order moment E[x_i x_j] 
	int i = d - ndims_;
	int j = 0;
	while (i > j) {
		i -= j + 1;
		j++;
	}
	// initialize acc with  x_i * x_j 
	for (int index = 0; index < pixel_count_; ++index) {
		acc_[index] = subs_[i] * subs_[j];
		adv(ndims_, dims_, subs_);
	}
}

void EcvMser::integrate()
{
	for (int i = 0; i < join_count_; ++i) {
		uint index = joins_[i];
		uint parent = r_[index].parent;
		acc_[parent] += acc_[index];
	}

}

void EcvMser::integrateMoments()
{

	for (int d = 0; d < dof_; ++d) {

		// start from the upper-left pixel (0,0,...,0) 
		memset(subs_, 0, sizeof(int) * ndims_);

		// step 1: fill acc pretending that each region has only one pixel 
		if (d < ndims_)
			oneOrder(d);
		else
			twoOrder(d);

		integrate();

		// step 3: save back to ellpises 
		for (int i = 0; i < mer_count_; ++i) {
			uint idx = mer_[i];
			ell_[d + dof_*i] = acc_[idx];
		}

	}

}

void EcvMser::computeCentralMmoments()
{

	for (int index = 0; index < mer_count_; ++index) {
		float  *pt = ell_ + index * dof_;
		uint    idx = mer_[index];
		float  area = r_[idx].area;

		for (int d = 0; d < dof_; ++d) {

			pt[d] /= area;

			if (d >= ndims_) {
				// remove squared mean from moment to get variance 
				int i = d - ndims_;
				int j = 0;
				while (i > j) {
					i -= j + 1;
					j++;
				}
				pt[d] -= pt[i] * pt[j];
			}

		}
	}
}


void EcvMser::drawImages()
{
	

	for (int i = 0; i < pixel_count_; ++i)
	{
		auto idx = pix_order_[i];
		auto cur_point = index2Coord(idx);


			parent_img_.at<float>(cur_point) = r_[idx].parent;
			shortcut_img_.at<float>(cur_point) = r_[idx].shortcut;
			area_img_.at<float>(cur_point) = r_[idx].area;
			height_img_.at<float>(cur_point) = r_[idx].height;
			join_img_.at<float>(cur_point) = joins_[idx];
			order_img_.at<float>(cur_point) = idx;
		

	}
}

uint EcvMser::upPass(uint& idx)
{
	auto prev_idx = idx;
	int next_idx;
	// move towards root to find it 
	while (true) {

		// next jump to the root 
		next_idx = r_[idx].shortcut;

		// recycle shortcut to remember how we came here 
		r_[idx].shortcut = prev_idx;

		// stop if the root is found
		if (next_idx == idx)
			break;

		// next guy 
		prev_idx = idx;
		idx = next_idx;
	}
	return idx;

}


void EcvMser::downPass(uint& idx)
{
	uint prev_idx;
	auto root_idx = idx;
	// move backward to update shortcuts 
	while (true) {

		// get previously visited one 
		prev_idx = r_[idx].shortcut;

		// update shortcut to point to the new root 
		r_[idx].shortcut = root_idx;

		// stop if the first visited node is reached 
		if (prev_idx == idx) break;

		// next guy 
		idx = prev_idx;
	}

}

uint EcvMser::climb(uint idx)
{
	auto root_idx = upPass(idx);
	downPass(idx);
	return root_idx;
}



uchar EcvMser::delta()
{
	return delta_;
}


void EcvMser::setDelta(uchar x)
{
	delta_ = x;
}


double EcvMser::minDiversity()
{
	return min_diversity_;
}


void EcvMser::setMinDiversity(double x)
{
	min_diversity_ = x;
}

EcvMserStats const* EcvMser::stats()
{
	return &stats_;
}

double EcvMser::maxArea()
{
	return max_area_;
}

void EcvMser::setMaxArea(double x)
{
	max_area_ = x;
}


double EcvMser::minArea()
{
	return min_area_;
}


void EcvMser::setMinArea(double x)
{
	min_area_ = x;
}


double EcvMser::max_variation()
{
	return max_variation_;
}


void EcvMser::setMaxVariation(double x)
{
	max_variation_ = x;
}


uint const * EcvMser::regions()
{
	return mer_;
}


uint EcvMser::regionsNum()
{
	return mer_count_;
}


float const * EcvMser::ell()
{
	return ell_;
}


uint EcvMser::ellDof()
{
	return dof_;
}


uint EcvMser::ellNum()
{
	return nell_;
}


void EcvMser::ellFit()
{

	int d, index, i, j;

	// already fit ? 
	if (nell_ == mer_count_)
		return;

	// make room 
	if (rell_ < mer_count_) {
		if (ell_)
			free(ell_);
		ell_ = (float*)malloc(sizeof(float) * mer_count_ * dof_);
		rell_ = mer_count_;
	}

	if (acc_ == 0)
		acc_ = (float*)malloc(sizeof(float) * pixel_count_);


	view_img_ = cv::Mat(dims_[1], dims_[0], CV_32FC1, acc_, dims_[0] * sizeof(float));

	integrateMoments();
	computeCentralMmoments();


	nell_ = mer_count_;
}

void EcvMser::adv(int ndims, int const *dims, int *subs)
{
	int d = 0;
	while (d < ndims) {
		if (++subs[d] < dims[d]) 
			return;
		subs[d++] = 0;
	}
}


cv::Point EcvMser::index2Coord(int idx)
{
	/* convert the index IDX into the subscript SUBS; also initialize
	DSUBS to (-1,-1,...,-1) */
	return cv::Point(idx % image_size_.width, int(idx / image_size_.width));
	/*{
		uint temp = idx;
		for (int k = ndims_ - 1; k >= 0; --k) {

			subs_[k] = temp / strides_[k];
			temp = temp % strides_[k];
		}
	}*/

}


bool EcvMser::getNeighbor(int pos, const cv::Point& pnt, int index, int& n_index)
{
	cv::Point n_pnt = pnt + n_coords[pos];
	n_index = index + n_coords[pos].y * image_size_.width + n_coords[pos].x;
	return image_rect_.contains(n_pnt);
}


void EcvMser::processNeighbor(bool good, int n_idx, int idx, int r_idx)
{
	auto val = img_[idx];
	if (good &&
		n_idx != idx &&
		r_[n_idx].parent != ECV_MSER_VOID_NODE) {

		uchar nr_val = 0;
		uint nr_idx = 0;

		auto hgt = r_[r_idx].height;
		auto n_hgt = r_[nr_idx].height;

		r_idx = climb(idx);
		nr_idx = climb(n_idx);


		if (r_idx != nr_idx) {

			nr_val = img_[nr_idx];

			if (nr_val == val && hgt < n_hgt) {

				rootIdxChild(r_idx, nr_idx, hgt, n_hgt);
				
				joins_[join_count_++] = r_idx;

			}
			else {
				rootIdxParent(r_idx, nr_idx, hgt, n_hgt);
				
				joins_[join_count_++] = nr_idx;

				// count if extremal 
				if (nr_val != val)
					++er_count_;

			}
		}
	}

}


//(B)I(ROOT(IDX)) == I(ROOT(NR_IDX)).In this case the pixel
//IDX is extending an extremal region with the same
//intensity value.Since ROOT(NR_IDX) will NOT be an
//extremal region of the full image, ROOT(IDX) can be
//safely added as children of ROOT(NR_IDX) if this
//reduces the height according to the union rank
//heuristic.

void EcvMser::rootIdxChild(int r_idx, int nr_idx, int hgt, int n_hgt)
{
	/* ROOT(IDX) becomes the child */
	r_[r_idx].parent = nr_idx;
	r_[r_idx].shortcut = nr_idx;
	r_[nr_idx].area += r_[r_idx].area;
	r_[nr_idx].height = std::max(n_hgt, hgt + 1);

}

//(C)I(ROOT(IDX)) > I(ROOT(NR_IDX)).In this case the pixel
//IDX is starting a new extremal region.Thus ROOT(NR_IDX)
//WILL be an extremal region of the final image and the
//only possibility is to add ROOT(NR_IDX) as children of
//ROOT(IDX), which becomes parent.


void EcvMser::rootIdxParent(int r_idx, int nr_idx, int hgt, int n_hgt)
{
	/* cases ROOT(IDX) becomes the parent */
	r_[nr_idx].parent = r_idx;
	r_[nr_idx].shortcut = r_idx;
	r_[r_idx].area += r_[nr_idx].area;
	r_[r_idx].height = std::max(hgt, n_hgt + 1);

}


void EcvMser::removeDuplicate(int i, int min_div, int& dup_count, ErState& state)
{
	uint   parent = e_tree_[i].parent;
	int       area, p_area;
	float div;

	// check all but the root mser 
	if ((int)parent != i) {

		// search for the maximally stable parent region 
		while (!e_tree_[parent].max_stable) {
			uint next = e_tree_[parent].parent;
			if (next == parent) 
				break;
			parent = next;
		}

		// Compare with the parent region; if the current and parent
		// regions are too similar, keep only the parent. 
		area = e_tree_[i].area;
		p_area = e_tree_[parent].area;
		div = (float)(p_area - area) / (float)p_area;

		if (div < min_div) {
			++dup_count;
			state = ErState::min_div;
		}
	}
}


void EcvMser::getHistogram(uchar* img, int pix_count, uint* hist)
{

	memset(hist, 0, sizeof(uint) * ECV_MSER_PIX_MAXVAL);

	/* compute bucket size (how many pixels for each intensity value) */
	for (int i = 0; i < (int)pixel_count_; ++i) {
		auto v = img[i];
		++hist[v];
	}


}

void EcvMser::getCDF(uint* hist)
{
	/* cumulatively add bucket sizes */
	for (int i = 1; i < ECV_MSER_PIX_MAXVAL; ++i) {
		hist[i] += hist[i - 1];
	}

}

void EcvMser::setPixOrder(uint* hist)
{
	/* empty buckets computing pixel ordering */
	for (int i = pixel_count_; i >= 1; ) {
		auto v = img_[--i];
		auto j = --hist[v];
		pix_order_[j] = i;
	}

}

void EcvMser::initForest()
{
	// initialize the forest with all void nodes 
	for (int i = 0; i < pixel_count_; ++i)
		r_[i].parent = ECV_MSER_VOID_NODE;


}


void EcvMser::addAsRoot(int idx)
{
	// add the pixel to the forest as a root for now 
	r_[idx].parent = idx;
	r_[idx].shortcut = idx;
	r_[idx].area = 1;
	r_[idx].height = 1;

}

void EcvMser::computeRegion()
{
	parent_img_ = cv::Mat(image_size_, CV_32FC1, cv::Scalar::all(0));
	shortcut_img_ = parent_img_.clone();
	area_img_ = parent_img_.clone();
	height_img_ = parent_img_.clone();
	join_img_ = parent_img_.clone();
	order_img_ = parent_img_.clone();
	extrmal_img_ = parent_img_.clone();
	temp_img_ = parent_img_.clone();

	er_count_ = 0;
	join_count_ = 0;

	for (int i = 0; i < pixel_count_; ++i)
	{
		auto idx = pix_order_[i];
		addAsRoot(idx);
		auto r_idx = idx;
		auto cur_point = index2Coord(idx);
		
		// examine the neighbors of the current pixel 
		for (int kernel_index = 0; kernel_index < 9; kernel_index++) {
			int n_idx;

			bool good = getNeighbor(kernel_index, cur_point, idx, n_idx);

			processNeighbor(good, n_idx, idx, r_idx);

			

		}
	}
	
	
	drawImages();

	cv::Mat a1 = cv::Mat(image_size_, CV_32FC1,cv::Scalar::all(0));
	int counter = 1;
	int cur_val = 0;
	while(1){
		auto coord = index2Coord(cur_val);

		auto parent_val = shortcut_img_.at<float>(coord);
		if (parent_val == cur_val)
			break;
		a1.at<float>(coord) = counter;
		cur_val = parent_val;
		
		imshow("a", a1);
		cv::waitKey(1);

	}



	// the last root is extremal too 
	++er_count_;



	stats_.num_extremal = er_count_;

	/* -----------------------------------------------------------------
	*                                          Extract extremal regions
	* -------------------------------------------------------------- */

	
}

Answer 5

bool EcvMser::isExtermal(int idx)
{
	auto val = img_[idx];
	auto p_idx = r_[idx].parent;
	auto p_val = img_[p_idx];

	bool is_extremal = (p_val > val) || (idx == p_idx);
	return is_extremal;
}


void EcvMser::reg2ExtermalTree(int idx, int er_count)
{
	e_tree_[er_count].index = idx;
	e_tree_[er_count].parent = er_count;
	e_tree_[er_count].value = img_[idx];
	e_tree_[er_count].area = r_[idx].area;

	


}

void EcvMser::createExtermalTree()
{
	/*
	Extremal regions are extracted and stored into the array ER.  The
	structure R is also updated so that .SHORTCUT indexes the
	corresponding extremal region if any (otherwise it is set to
	VOID).
	*/


	if (rer_ < er_count_) {
		if (e_tree_)
			free(e_tree_);

		e_tree_ = (EcvMserExtrReg*)malloc(sizeof(EcvMserExtrReg) * er_count_);
		rer_ = er_count_;
	};


	mer_count_ = er_count_;
	er_count_ = 0;

	for (int i = 0; i < pixel_count_; ++i) {
		

		auto idx = pix_order_[i];
		auto cur_point = index2Coord(idx);

		if (isExtermal(idx)) {

			extrmal_img_.at<float>(cur_point) = 255;

			reg2ExtermalTree(idx, er_count_);

			r_[idx].shortcut = er_count_;	// link this region to this extremal region 
			
			++er_count_;// increase count 
		}
		else {
			r_[idx].shortcut = ECV_MSER_VOID_NODE;// link this region to void 
			extrmal_img_.at<float>(cur_point) = 128;
		}

	}

}

void EcvMser::linkExtermalRegions2Tree()
{

	for (int i = 0; i < er_count_; ++i) {

		uint idx = e_tree_[i].index;

		do {
			idx = r_[idx].parent;
		} while (r_[idx].shortcut == ECV_MSER_VOID_NODE);

		e_tree_[i].parent = r_[idx].shortcut;
		e_tree_[i].shortcut = i;
	}

}

void EcvMser::computeVaribility()
{
	/* -----------------------------------------------------------------
	*                            Compute variability of +DELTA branches
	* -------------------------------------------------------------- */
	/* For each extremal region Xi of value VAL we look for the biggest
	* parent that has value not greater than VAL+DELTA. This is dubbed
	* `top parent'. */

	for (int i = 0; i < er_count_; ++i) {

		/* Xj is the current region the region and Xj are the parents */
		int     top_val = e_tree_[i].value + delta_;
		int     top = e_tree_[i].shortcut;

		/* examine all parents */
		while (1) {
			int next = e_tree_[top].parent;
			int next_val = e_tree_[next].value;

			/* Break if:
			* - there is no node above the top or
			* - the next node is above the top value.
			*/
			if (next == top || next_val > top_val) break;

			/* so next could be the top */
			top = next;
		}

		/* calculate branch variation */
		{
			int area = e_tree_[i].area;
			int area_top = e_tree_[top].area;
			e_tree_[i].variation = (float)(area_top - area) / area;
			e_tree_[i].max_stable = 1;
		}

		/* Optimization: since extremal regions are processed by
		* increasing intensity, all next extremal regions being processed
		* have value at least equal to the one of Xi. If any of them has
		* parent the parent of Xi (this comprises the parent itself), we
		* can safely skip most intermediate node along the branch and
		* skip directly to the top to start our search. */
		{
			int parent = e_tree_[i].parent;
			int curr = e_tree_[parent].shortcut;
			e_tree_[parent].shortcut = std::max(top, curr);
		}
	}
}

void EcvMser::SelectMaximallyStableBranches()
{

	mer_count_ = er_count_;
	for (int i = 0; i < er_count_; ++i) {
		auto parent = e_tree_[i].parent;
		auto val = e_tree_[i].value;
		float var = e_tree_[i].variation;
		uchar p_val = e_tree_[parent].value;
		float   p_var = e_tree_[parent].variation;



		//Notice that R_parent = R_{l+1} only if p_val = val + 1. If not,
		//this and the parent region coincide and there is nothing to do.

		if (p_val > val + 1)
			continue;

		// decide which one to keep and put that in loser
		uint loser = var < p_var ? parent : i;

		// make loser NON maximally stable 
		if (e_tree_[loser].max_stable) {
			--mer_count_;
			e_tree_[loser].max_stable = 0;
		}
	}

	stats_.num_unstable = er_count_ - mer_count_;
}

void EcvMser::furtherFiltering()
{
	// It is critical for correct duplicate detection to remove regions
	// from the bottom (smallest one first).                          

	int big_count = 0;
	int small_count = 0;
	int bad_count = 0;
	int dup_count = 0;


	float max_area = (float)max_area_ * pixel_count_;
	float min_area = (float)min_area_ * pixel_count_;
	float max_var = (float)max_variation_;
	float min_div = (float)min_diversity_;

	// scan all extremal regions (intensity value order) 
	for (int i = er_count_ - 1; i >= 0; --i) {

		// process only maximally stable extremal regions 
		if (!e_tree_[i].max_stable)
			continue;
		ErState er_state = ErState::none;

		if (e_tree_[i].variation >= max_var) {
			++bad_count;
			er_state = ErState::max_var;
		}
		else if (e_tree_[i].area > max_area) {
			++big_count;
			er_state = ErState::max_area;
		}
		else if (e_tree_[i].area < min_area) {
			++small_count;
			er_state = ErState::min_area;
		}
		else if (min_div < 1.0) {
			removeDuplicate(i, min_div, dup_count, er_state);
			

		}
		if (er_state != ErState::none) {
			e_tree_[i].max_stable = 0;
			--mer_count_;
		}
	}

	stats_.num_abs_unstable = bad_count;
	stats_.num_too_big = big_count;
	stats_.num_too_small = small_count;
	stats_.num_duplicates = dup_count;

}

void EcvMser::saveResult()
{

	// make room 
	if (rmer_ < mer_count_) {
		if (mer_)
			free(mer_);
		mer_ = (uint*)malloc(sizeof(uint) * mer_count_);
		rmer_ = mer_count_;
	}

	// save back 
	mer_count_ = mer_count_;

	int j = 0;
	for (int i = 0; i < er_count_; ++i)
		if (e_tree_[i].max_stable)
			mer_[j++] = e_tree_[i].index;

}