graft_cv/utils.cpp

#include "utils.h"
#include <opencv.hpp>
#include <time.h>
#include <algorithm>
#include <math.h>
//#include "fork_rs.h"

namespace graft_cv{

	Mat imginfo2mat(ImgInfo* imginfo)
	{
		assert(imginfo->data);
		assert(imginfo->height>0 && imginfo->width>0);
		Mat m = Mat(imginfo->height,imginfo->width,CV_8UC1);
		for(int h=0; h<imginfo->height; ++h)
		{
			memcpy((void*)(m.ptr(h)),
				(void*)(imginfo->data+h*imginfo->width),
				imginfo->width);
		};
		return m;
	};

	ImgInfo* mat2imginfo(const Mat&img){
		assert(img.channels()==1);
		ImgInfo* imginfo = new ImgInfo();
		imginfo->angle=0;
		imginfo->height=img.rows;
		imginfo->width = img.cols;
		imginfo->data = 0;
		imginfo->data = new byte[imginfo->height * imginfo->width];
		memset(imginfo->data, 0,imginfo->height * imginfo->width);	
		//IplImage ipl_img = cvIplImage(img);
		for(int i =0; i<imginfo->height;++i){
			memcpy(imginfo->data+i*imginfo->width, img.ptr(i), imginfo->width);
		}

		return imginfo;
	};
	void imginfo_release(ImgInfo**ppImgInfo){
		ImgInfo* pImgInfo = *ppImgInfo;
		if(pImgInfo->data){
			delete [] pImgInfo->data;
			pImgInfo->data=0;
		}
		delete pImgInfo;
		pImgInfo = 0;
	};
	bool isvalid(const ImgInfo*imginfo)
	{
		if(!imginfo){return false;}
		if(!imginfo->data){return false;}
		if(imginfo->height*imginfo->width<=0){return false;}
		return true;
	}

	void image_set_bottom(
		Mat& img,
		unsigned char value,
		int rows_cnt)
	{
		for(size_t r=(img.rows-rows_cnt);r<img.rows;++r){
			for(size_t c=0;c<img.cols;++c){
				img.at<unsigned char>(r,c) = value;			
			}
		}
	}
	void image_set_top(
		Mat& img,
		unsigned char value,
		int rows_cnt)
	{
		for(size_t r=0;r<rows_cnt;++r){
			for(size_t c=0;c<img.cols;++c){
				img.at<unsigned char>(r,c) = value;			
			}
		}
	}

	void image_draw_line(
		Mat& img,
		int x0,int y0,
		int x1,int y1)
	{
		Point p0,p1;
		if(x0==x1){
			p0.x=x0;
			p0.y=0;
			p1.x=x0;
			p1.y=img.rows-1;
		}
		else{
			double k = (double)(y1-y0)/(double)(x1-x0);
			double b = (double)y0 - k * x0;

			p0.x=0;
			p0.y=b;

			p1.x=img.cols;
			p1.y = (int)(b + k * p1.x);
		}
		line(img,p0,p1,Scalar(255,255,255));
	
	}
	string currTime(){
			char tmp[64];
			time_t ptime;
			time(&ptime);
			strftime(
				tmp, 
				sizeof(tmp), 
				"%Y-%m-%d %H:%M:%S",
				localtime(&ptime)
				);
			return tmp;
		}
	// function getImgId
	// input im_type, img_type
	static unsigned int ID_COUNTER = 0;
	string getImgId(int im_type)
	{
		char tmp[64];
		time_t ptime;
		time(&ptime);
		strftime(
			tmp, 
			sizeof(tmp), 
			"%Y%m%d%H%M%S",
			localtime(&ptime)
			);

		unsigned int id_serial = ID_COUNTER % 100;	
		stringstream buff;
		buff<<tmp;
		buff.width(2);
		buff.fill('0');	
		buff<<id_serial;
		buff<<"_"<<im_type;
		ID_COUNTER+=1;	
		return buff.str();
	}
	string getImgIdOa(string iid, int im_idx)
	{	
		stringstream buff;
		buff<<im_idx;
		return iid+"_"+buff.str();		
	}

	inline void throw_msg(string& msg){	
			stringstream buff;
			buff<<"Error:"<<msg<<"\tfile:"<<__FILE__<<"\tline:"<<__LINE__;			
			throw(buff.str());		
		}

	void drawgrid(		
			Mat&img,		
			int grid_col/*=50*/,
			int grid_row/*=50*/,
			int line_thick/*=2*/
			)
	{
		if(img.empty()){return;}
		//horizontal grid
		for(int row=img.rows-1; row>=0; row-=grid_row){
			line(img,Point(0,row),Point(img.cols-1,row),Scalar(100),line_thick);
		}
		//veritcal grid
		for(int col=0; col<img.cols; col+=grid_col){
			line(img,Point(col,0),Point(col,img.rows-1),Scalar(100),line_thick);
		}

	}
	void hist2image(
		vector<int>& hist, 
		Mat&img, 
		int aix, /*=1*/
		int grid_col/*=50*/,
		int grid_row/*=50*/
		)
	{
		vector<int>::iterator maxit = max_element(begin(hist),end(hist)); 
		if( grid_col<10){grid_col=10;}
		if( grid_row<10){grid_row=10;}
		// aix = 0, column hist, aix=1, row histogtam
		if (aix==0){
			img = Mat::zeros(hist.size(),*maxit, CV_8UC1);
			drawgrid(img,grid_col, grid_row,1);
			for(size_t i=0; i<hist.size();++i){
				if(i>=img.rows){break;}
				for(size_t j=0;j<=hist[i];++j){
					if(j>=img.cols){break;}
					img.at<unsigned char>(i,j)=255;
				}
			}
		}
		else{
			img = Mat::zeros(*maxit,hist.size(), CV_8UC1);
			drawgrid(img,grid_col, grid_row,1);
			for(size_t i=0; i<hist.size();++i){
				int y = *maxit-hist[i]-1;
				if (y<0){y=0;}
				for(size_t j=y;j<img.rows;++j){
					img.at<unsigned char>(j,i)=255;
				}
			}
		}
	};


	//matHistogram
	// axis == 0, statistic histogram for columns; 
	// others, statistic histogram for rows
	void mat_histogram(
		Mat& img,	
		std::vector<int>&hist,
		int axis,
		int r0, 
		int r1, 
		int c0, 
		int c1
		)
	{
		hist.clear();
		if(r0<0){r0 = 0;}
		if(r1<0){r1 = img.rows;}
		if(c0<0){c0 = 0;}
		if(c1<0){c1 = img.cols;}

		if (axis==0){
			for(int c=c0;c<c1;++c){
				int cnt = 0;
				for(int r=r0;r<r1; ++r){
					if (img.at<unsigned char>(r,c)>0){cnt+=1;}
				}
				hist.push_back(cnt);
			}
		}
		else{
			for(int r=r0;r<r1; ++r){
				int cnt = 0;
				for(int c=c0;c<c1;++c){
					if (img.at<unsigned char>(r,c)>0){cnt+=1;}
				}
				hist.push_back(cnt);
			}
		}
	};

	//matHistogram
	// axis == 0, statistic histogram for columns; 
	// others, statistic histogram for rows
	// with weight
	void mat_histogram_w(
		Mat& img,	
		std::vector<int>&hist,
		int axis,
		int r0, 
		int r1, 
		int c0, 
		int c1,
		bool asc_w //true---ascending weight, [0-1.0], false --- descending weight [1.0--0.0]
		)
	{
		hist.clear();
		if(r0<0){r0 = 0;}
		if(r1<0){r1 = img.rows;}
		if(c0<0){c0 = 0;}
		if(c1<0){c1 = img.cols;}
	
		if (axis==0){
			double step = 1.0/(double)(r1-1);
			for(int c=c0;c<c1;++c){
				double cnt = 0.0;
				double hi=0.0;
				for(int r=r0;r<r1; ++r){
					if (img.at<unsigned char>(r,c)>0){
						if(asc_w){
							hi = (r-r0)*step;
							if(hi>=0.66666){
								cnt+=hi;
							}
						}
						else{
							hi = (r1-r+r0)*step;
							if(hi>=0.66666){
								cnt+=hi;
							}
						}
					}
				}
				hist.push_back((int)cnt);
			}
		}
		else{
			double step = 1.0/(double)(c1-1);
			for(int r=r0;r<r1; ++r){
				double cnt = 0.0;
				double hi=0.0;
				for(int c=c0;c<c1;++c){
					if (img.at<unsigned char>(r,c)>0){
						if(asc_w){
							hi = (c-c0)*step;
							if(hi>=0.66666){
								cnt+=hi;
							}
						}
						else{
							hi = (c1-c+c0)*step;
							if(hi>=0.66666){
								cnt+=hi;
							}
						}
					}
				}
				hist.push_back((int)cnt);
			}
		}
	};

	void mat_histogram_yfork(
		Mat& img,	
		std::vector<int>&hist,	
		int r0, 
		int r1	
		)
	{
		hist.clear();
		if(r0<0){r0 = 0;}
		if(r1<0){r1 = img.rows;}	
	
		double step = 1.0/(double)(r1-r0);
		for(int c=0;c<img.cols;++c){
			double cnt = 0.0;
			double hi=0.0;
			for(int r=r0;r<r1; ++r){
				if (img.at<unsigned char>(r,c)>0){				
					hi = (r-r0)*step;
					//if(hi>=0.66666){
						cnt+=hi;
					//}
				
				
				}
			}
			hist.push_back((int)cnt);
		}
	
	};

	//get_stem_x_range() 确定茎的位置，x方向
	// 2022/1/4		chenhj	修改，大叶子下垂，将叶子识别成茎，增加histogram满足阈值情况下，两侧数据方差检测 
	void  get_stem_x_range_oa(
		const std::vector<int>&hist_h,
		double h_ratio, 
		int padding, 
		int cent_x,
		int width_x,
		int&x0, 
		int&x1, 
		int&stem_x0,
		int&stem_x1
		)
	{
		//hist_h: histogram in column
		//h_ratio: ratio for generating threshold
		// x0,x1: sub-image x-range
		// stem_x0,1: stem x-range
		// padding: pad for stem x-range

		// calculate the max-value of hist_h
		int min_segment_len = 0;//最小线段长度，如果小于此数值，不计为有效数据，跳过
		int min_segment_dist = 20;//多峰情况下，间隔像素数量
		int var_padding = 20;
		auto max_it = max_element(hist_h.begin(), hist_h.end());
		size_t max_idx = max_it - hist_h.begin();
		int th = int(h_ratio * (double)(*max_it));

		if(th<=0){
			throw_msg(string("error: threshold is 0 in get_stem_x_range()"));
		}

		// 1 计算小线段的起始位置
		vector<gcv_point<int>>segments;	
		get_hist_segment(hist_h,th,segments);

		// 2 获取多峰位置
		vector<gcv_point<int>>peaks;	
		
		if(segments.size()==1){
			peaks = segments; //单峰
		}
		else{// 小线段删除
			vector<gcv_point<int>>big_segments;
			for(size_t i=0;i<segments.size();++i){
				if(segments[i].y - segments[i].x>=min_segment_len){
					big_segments.push_back(segments[i]);
				}
			}

			if(big_segments.size()==1){
				peaks = big_segments;//单峰
			}
			else{// 合并
				peaks.push_back(big_segments[0]);
				for(size_t j=1;j<big_segments.size();++j){
					int dist = big_segments[j].x - peaks.back().y;
					if(dist>min_segment_dist){
						peaks.push_back(big_segments[j]);
					}
					else{
						peaks.back().y = big_segments[j].y;
					}
				}

			}
		}

		// 3 多峰情况，先计算center_ratio, 统计峰段两侧数据的方差，方差大的为茎
		int opt_index=0;
		if(peaks.size()>1){
			vector<double>center_r;
			for(size_t i=0;i<peaks.size();++i){
				double r = fabs(0.5*(double)(peaks[i].x+ peaks[i].y) - cent_x)/(double)(width_x/2.0);
				center_r.push_back(1.0-r);
			}
			vector<size_t>sort_idx = sort_indexes_e(center_r,false);

			double cr_1 = center_r[sort_idx[0]];
			double cr_2 = center_r[sort_idx[1]];		

			if(cr_1-cr_2>0.2 || cr_1>=0.65){//如果存在居中显著的，按居中优先
				opt_index = sort_idx[0];
			}
			else{
				vector<int> candidate_idx;
				candidate_idx.push_back(sort_idx[0]);
				candidate_idx.push_back(sort_idx[1]);
				for(size_t id = 2;id<sort_idx.size();++id){
					if(cr_1 - center_r[sort_idx[id]]<=0.3){
						candidate_idx.push_back(sort_idx[id]);
					}
				}
				vector<double>vars;
				int vx0,vx1;
				for(size_t i=0;i<peaks.size();++i){
					vector<int>::iterator it = find(candidate_idx.begin(),candidate_idx.end(),i);
					if(it==candidate_idx.end()){
						vars.push_back(0.0);
						continue;
					}
					vx0 = peaks[i].x - var_padding;
					vx1 = peaks[i].y + var_padding;
					if(vx0<0){vx0=0;}
					if(vx1>hist_h.size()-1){vx1=hist_h.size()-1;}
					double mean,var;
					mean = 0.0;
					var=-1;
					get_hist_mean_var_local(hist_h,vx0,vx1,mean,var);
					if(var<0){var=0.0;}
					vars.push_back(var);
				}
				auto var_max_it = max_element(vars.begin(), vars.end());
				size_t var_max_idx = var_max_it - vars.begin();
				opt_index = var_max_idx;
			}		
		}

		stem_x0 = peaks[opt_index].x;
		stem_x1 = peaks[opt_index].y;
		x0 = stem_x0 - padding;
		x1 = stem_x1 + padding;
		if(x0<0){x0=0;};
		if(x1>hist_h.size()-1){x1 = hist_h.size()-1;};    


		//////////////////////////////////////

		//int global_start_idx=-1;
		//int global_end_idx=-1;
		//int start_idx=-1;
		//int end_idx=-1;

		//for(size_t i=0;i<hist_h.size();++i){
		//	if(i==0){
		//		if(hist_h[i]>=th){
		//			start_idx=i;
		//		}
		//		continue;
		//	}
		//	if(hist_h[i]>=th && hist_h[i-1]<th){
		//		start_idx=i;
		//		continue;
		//	}
		//	if(i==hist_h.size()-1){
		//		if(hist_h[i]>=th && hist_h[i-1]>=th){
		//			if((i-start_idx+1)>=min_segment_len){
		//				global_end_idx=i;
		//			}
		//		}
		//	}
		//	if(hist_h[i]<th && hist_h[i-1]>=th){
		//		if((i-start_idx+1)>=min_segment_len){
		//			if( global_start_idx<0){
		//				global_start_idx = start_idx;
		//			}
		//			global_end_idx = i;
		//		}
		//	}

		//}
		///*stem_x0 = 0;
		//stem_x1 = hist_h.size()-1;

		//for(size_t i=max_idx; i < hist_h.size(); ++i){
		//	if(hist_h[i]>=th){
		//		stem_x1 = i;
		//	}
		//	else{
		//		break;
		//	}
		//}
		//for(int i=max_idx; i >= 0; i--){
		//	if(hist_h[i]>=th){
		//		stem_x0 = i;
		//	}
		//	else{
		//		break;
		//	}
		//}*/
		//stem_x0 = global_start_idx;
		//stem_x1 = global_end_idx;
		//x0 = stem_x0 - padding;
		//x1 = stem_x1 + padding;
		//if(x0<0){x0=0;};
		//if(x1>hist_h.size()-1){x1 = hist_h.size()-1;};    
	};

	void  get_stem_x_range_rscut(
		const std::vector<int>&hist_h,
		double h_ratio, 
		int padding, 
		int cent_x,
		int width_x,
		int&x0, 
		int&x1, 
		int&stem_x0,
		int&stem_x1
		)
	{
		//hist_h: histogram in column
		//h_ratio: ratio for generating threshold
		// x0,x1: sub-image x-range
		// stem_x0,1: stem x-range
		// padding: pad for stem x-range

		// calculate the max-value of hist_h
		int min_segment_len = 0;//最小线段长度，如果小于此数值，不计为有效数据，跳过
		int min_segment_dist = 20;//多峰情况下，间隔像素数量
		int var_padding = 20;
		auto max_it = max_element(hist_h.begin(), hist_h.end());
		size_t max_idx = max_it - hist_h.begin();
		int th = int(h_ratio * (double)(*max_it));

		if(th<=0){
			throw_msg(string("error: threshold is 0 in get_stem_x_range()"));
		}

		// 1 计算小线段的起始位置
		vector<gcv_point<int>>segments;	
		get_hist_segment(hist_h,th,segments);

		// 2 获取多峰位置
		vector<gcv_point<int>>peaks;	
		
		if(segments.size()==1){
			peaks = segments; //单峰
		}
		else{// 小线段删除
			vector<gcv_point<int>>big_segments;
			for(size_t i=0;i<segments.size();++i){
				if(segments[i].y - segments[i].x>=min_segment_len){
					big_segments.push_back(segments[i]);
				}
			}

			if(big_segments.size()==1){
				peaks = big_segments;//单峰
			}
			else{// 合并
				peaks.push_back(big_segments[0]);
				for(size_t j=1;j<big_segments.size();++j){
					int dist = big_segments[j].x - peaks.back().y;
					if(dist>min_segment_dist){
						peaks.push_back(big_segments[j]);
					}
					else{
						peaks.back().y = big_segments[j].y;
					}
				}

			}
		}

		// 3 多峰情况，先计算center_ratio, 统计峰段两侧数据的方差，方差大的为茎
		int opt_index=0;
		if(peaks.size()>1){
			vector<double>center_r;
			for(size_t i=0;i<peaks.size();++i){
				double r = fabs(0.5*(double)(peaks[i].x+ peaks[i].y) - cent_x)/(double)(width_x/2.0);
				center_r.push_back(1.0-r);
			}
			vector<size_t>sort_idx = sort_indexes_e(center_r,false);

			double cr_1 = center_r[sort_idx[0]];
			double cr_2 = center_r[sort_idx[1]];		

			if(cr_1-cr_2>0.2 || cr_1>=0.65){//如果存在居中显著的，按居中优先
				opt_index = sort_idx[0];
			}
			else{
				vector<int> candidate_idx;
				candidate_idx.push_back(sort_idx[0]);
				candidate_idx.push_back(sort_idx[1]);
				for(size_t id = 2;id<sort_idx.size();++id){
					if(cr_1 - center_r[sort_idx[id]]<=0.3){
						candidate_idx.push_back(sort_idx[id]);
					}
				}
				vector<double>vars;
				int vx0,vx1;
				for(size_t i=0;i<peaks.size();++i){
					vector<int>::iterator it = find(candidate_idx.begin(),candidate_idx.end(),i);
					if(it==candidate_idx.end()){
						vars.push_back(0.0);
						continue;
					}
					vx0 = peaks[i].x - var_padding;
					vx1 = peaks[i].y + var_padding;
					if(vx0<0){vx0=0;}
					if(vx1>hist_h.size()-1){vx1=hist_h.size()-1;}
					double mean,var;
					mean = 0.0;
					var=-1;
					get_hist_mean_var_local(hist_h,vx0,vx1,mean,var);
					if(var<0){var=0.0;}
					vars.push_back(var);
				}
				auto var_max_it = max_element(vars.begin(), vars.end());
				size_t var_max_idx = var_max_it - vars.begin();
				opt_index = var_max_idx;
			}		
		}

		stem_x0 = peaks[opt_index].x;
		stem_x1 = peaks[opt_index].y;
		x0 = stem_x0 - padding;
		x1 = stem_x1 + padding;
		if(x0<0){x0=0;};
		if(x1>hist_h.size()-1){x1 = hist_h.size()-1;};    


		//////////////////////////////////////

		//int global_start_idx=-1;
		//int global_end_idx=-1;
		//int start_idx=-1;
		//int end_idx=-1;

		//for(size_t i=0;i<hist_h.size();++i){
		//	if(i==0){
		//		if(hist_h[i]>=th){
		//			start_idx=i;
		//		}
		//		continue;
		//	}
		//	if(hist_h[i]>=th && hist_h[i-1]<th){
		//		start_idx=i;
		//		continue;
		//	}
		//	if(i==hist_h.size()-1){
		//		if(hist_h[i]>=th && hist_h[i-1]>=th){
		//			if((i-start_idx+1)>=min_segment_len){
		//				global_end_idx=i;
		//			}
		//		}
		//	}
		//	if(hist_h[i]<th && hist_h[i-1]>=th){
		//		if((i-start_idx+1)>=min_segment_len){
		//			if( global_start_idx<0){
		//				global_start_idx = start_idx;
		//			}
		//			global_end_idx = i;
		//		}
		//	}

		//}
		///*stem_x0 = 0;
		//stem_x1 = hist_h.size()-1;

		//for(size_t i=max_idx; i < hist_h.size(); ++i){
		//	if(hist_h[i]>=th){
		//		stem_x1 = i;
		//	}
		//	else{
		//		break;
		//	}
		//}
		//for(int i=max_idx; i >= 0; i--){
		//	if(hist_h[i]>=th){
		//		stem_x0 = i;
		//	}
		//	else{
		//		break;
		//	}
		//}*/
		//stem_x0 = global_start_idx;
		//stem_x1 = global_end_idx;
		//x0 = stem_x0 - padding;
		//x1 = stem_x1 + padding;
		//if(x0<0){x0=0;};
		//if(x1>hist_h.size()-1){x1 = hist_h.size()-1;};    
	};


	void  get_stem_x_range_scion(
		const std::vector<int>&hist_h,
		double h_ratio, 
		int padding, 
		int&x0, 
		int&x1, 
		int&stem_x0,
		int&stem_x1
		)
	{
		//hist_h: histogram in column
		//h_ratio: ratio for generating threshold
		// x0,x1: sub-image x-range
		// stem_x0,1: stem x-range
		// padding: pad for stem x-range

		// calculate the max-value of hist_h
		auto max_it = max_element(hist_h.begin(), hist_h.end());
		size_t max_idx = max_it - hist_h.begin();
		int th = int(h_ratio * (double)(*max_it));

		stem_x0 = 0;
		stem_x1 = hist_h.size()-1;

		for(size_t i=max_idx; i < hist_h.size(); ++i){
			if(hist_h[i]>=th){
				stem_x1 = i;
			}
			else{
				break;
			}
		}
		for(int i=max_idx; i >= 0; i--){
			if(hist_h[i]>=th){
				stem_x0 = i;
			}
			else{
				break;
			}
		}

		x0 = stem_x0 - padding;
		x1 = stem_x1 + padding;
		if(x0<0){x0=0;};
		if(x1>hist_h.size()-1){x1 = hist_h.size()-1;};    
	};

	void get_hist_segment(
			const std::vector<int>& hist, 
			int th,//threshold
			std::vector<gcv_point<int>>& segments
		)
	{
		segments.clear();	
		int start_idx=-1;
		int end_idx=-1;

		for(size_t i=0;i<hist.size();++i){
			if(i==0){
				if(hist[i]>=th){
					start_idx=i;
				}
				continue;
			}
			if(hist[i]>=th && hist[i-1]<th){
				start_idx=i;
				continue;
			}
			if(i==hist.size()-1){
				if(hist[i]>=th && hist[i-1]>=th){
					segments.push_back(gcv_point<int>(start_idx,i));
				}
			}
			if(hist[i]<th && hist[i-1]>=th){
				segments.push_back(gcv_point<int>(start_idx,i-1));
			}
		}
	}

	void get_hist_mean_var_local(
			const std::vector<int>& hist, 
			int x0,
			int x1,
			double& mean,
			double& var
		)
	{
		mean=0.0;
		var = 0.0;
		if(hist.size()==0){
			return;
		}
		if(x0<0){x0=0;}
		if(x1<0 || x1>hist.size()-1){hist.size()-1;}
	
		for(int i=x0;i<=x1;++i){
			mean+=hist[i];
		}
		mean /=(double)(x1-x0+1);
	
		for(int i=x0;i<=x1;++i){
			var+=((double)(hist[i])-mean) * ((double)(hist[i])-mean);
		}
		var /= (double)(x1-x0+1);
		
	}
	void get_stem_y_fork(
		const std::vector<int>& hist, 
		double ratio,
		int stem_dia_min, 
		int stem_fork_y_min,
		double stem_dia_mp,
		int& fork_y, //output
		int& end_y,  //output
		int& stem_dia    //output
		)
	{
		//# 由y的底部向上检测是否有茎,有的话计算移动均值
		//# 通过当前值和移动均值的比值识别分叉点

		//# stem_dia_min = 10 # 茎粗最小值
		//# stem_fork_y_min = 10 # 茎分叉位置与茎根的最小像素高度
		//# stem_dia_mp = 0.8 # moving power

		fork_y = end_y = stem_dia = -1;

		std::vector<int>reversed_hist;
		for(int ii=hist.size()-1; ii>=0;ii--){
			reversed_hist.push_back(hist[ii]);
		}

		int stem_root_y = -1;
		double stem_dia_ma = -1;//# stem diameter moving-average

		size_t i=0;
		for(; i<hist.size(); ++i){
			if (i==0){continue;}
			if ((reversed_hist[i-1]<stem_dia_min && reversed_hist[i]>=stem_dia_min && stem_root_y<0)
				||(reversed_hist[i-1]>=stem_dia_min && reversed_hist[i]>=stem_dia_min  && stem_root_y<0))
			{
				//# stem root begin
				stem_root_y = i;
				stem_dia_ma = (double)reversed_hist[i];
			}
			else {
				if (reversed_hist[i-1]>=stem_dia_min &&
					reversed_hist[i]>=stem_dia_min)
				{
					double fork_ratio = reversed_hist[i] / stem_dia_ma;
					//if (i<150){
					//	cout<<fork_ratio<<"\t"<<reversed_hist[i]<<"\t"<<stem_dia_ma<<endl;
					//}
					if( i - stem_root_y > stem_fork_y_min && 
						fork_ratio >= ratio)
					{
						//# found the fork y level
						fork_y = hist.size() - i;
						stem_dia = reversed_hist[i-1];
						break;
					}
					//# update stem_dia_ma
					stem_dia_ma = stem_dia_mp *stem_dia_ma + (1.0-stem_dia_mp) * reversed_hist[i];
				}
			}

		}; 
		if(stem_root_y<0){
		  throw_msg(string("not exists diameter bigger than stem_dia_min"));
		}
		if(fork_y<0){
			throw_msg(string("not exists diameter fork_ratio bigger than threshold"));
		}

		//end_y
		end_y =  hist.size() - stem_root_y-1;

		//统计stem_root_y 到 i 间reversed_hist的数值，得到stem_dia。方法可用平均值，75%分位数
		vector<int> tmp;
		for(size_t j=stem_root_y; j<=i;++j){
			tmp.push_back(reversed_hist[j]);
		}
		sort(tmp.begin(), tmp.end());
		int tar_idx = (int)((float)tmp.size()*0.75);
		stem_dia = tmp[tar_idx];
	};

	void get_stem_y_fork_3sigma(
		const std::vector<int>& hist, 
		double ratio,
		int stem_dia_min, 
		int stem_fork_y_min,
		double stem_dia_mp,
		int& fork_y, //output
		int& end_y,  //output
		int& stem_dia    //output
		)
	{
	
    

		fork_y = end_y = stem_dia = -1;

		std::vector<int>reversed_hist;
		for(int ii=hist.size()-1; ii>=0;ii--){
			reversed_hist.push_back(hist[ii]);
		}

		int mean_radius =2;
		double mean_dia = 2.0*mean_radius +1.0;
		double mean = 0.0;
		double min_m, max_m;
		min_m = max_m = 0.0;
		std::vector<double>reversed_mean;
		for(int ii=0; ii<reversed_hist.size();ii++){
			if(ii-mean_radius<0 || ii + mean_radius>=reversed_hist.size()){
				reversed_mean.push_back(0.0);
				continue;
			}
			mean = 0.0;
			for(int k = -mean_radius;k<=mean_radius;++k){
				mean+=(double)reversed_hist[ii+k];
			}
			mean /= mean_dia;
			if(mean>max_m){max_m = mean;}
			reversed_mean.push_back(mean);
		}
		Mat tmp;
		hist2image_scale_line(reversed_mean,tmp,min_m,max_m,1.0,min_m,50,50);
		imshow_wait("mean5",tmp);
		return;

		int stem_root_y = -1;
		double stem_dia_ma = -1;//# stem diameter moving-average

		size_t i=0;
		for(; i<reversed_mean.size(); ++i){
			if (i==0){continue;}
			if ((reversed_hist[i-1]<stem_dia_min && reversed_hist[i]>=stem_dia_min && stem_root_y<0)
				||(reversed_hist[i-1]>=stem_dia_min && reversed_hist[i]>=stem_dia_min  && stem_root_y<0))
			{
				//# stem root begin
				stem_root_y = i;
				stem_dia_ma = (double)reversed_hist[i];
			}
			else {
				if (reversed_hist[i-1]>=stem_dia_min &&
					reversed_hist[i]>=stem_dia_min)
				{
					double fork_ratio = reversed_hist[i] / stem_dia_ma;
					//if (i<150){
					//	cout<<fork_ratio<<"\t"<<reversed_hist[i]<<"\t"<<stem_dia_ma<<endl;
					//}
					if( i - stem_root_y > stem_fork_y_min && 
						fork_ratio >= ratio)
					{
						//# found the fork y level
						fork_y = hist.size() - i;
						stem_dia = reversed_hist[i-1];
						break;
					}
					//# update stem_dia_ma
					stem_dia_ma = stem_dia_mp *stem_dia_ma + (1.0-stem_dia_mp) * reversed_hist[i];
				}
			}

		}; 

		//r2index(

		/*
		int stem_root_y = -1;
		double stem_dia_ma = -1;//# stem diameter moving-average

		size_t i=0;
		for(; i<hist.size(); ++i){
			if (i==0){continue;}
			if ((reversed_hist[i-1]<stem_dia_min && reversed_hist[i]>=stem_dia_min && stem_root_y<0)
				||(reversed_hist[i-1]>=stem_dia_min && reversed_hist[i]>=stem_dia_min  && stem_root_y<0))
			{
				//# stem root begin
				stem_root_y = i;
				stem_dia_ma = (double)reversed_hist[i];
			}
			else {
				if (reversed_hist[i-1]>=stem_dia_min &&
					reversed_hist[i]>=stem_dia_min)
				{
					double fork_ratio = reversed_hist[i] / stem_dia_ma;
					//if (i<150){
					//	cout<<fork_ratio<<"\t"<<reversed_hist[i]<<"\t"<<stem_dia_ma<<endl;
					//}
					if( i - stem_root_y > stem_fork_y_min && 
						fork_ratio >= ratio)
					{
						//# found the fork y level
						fork_y = hist.size() - i;
						stem_dia = reversed_hist[i-1];
						break;
					}
					//# update stem_dia_ma
					stem_dia_ma = stem_dia_mp *stem_dia_ma + (1.0-stem_dia_mp) * reversed_hist[i];
				}
			}

		}; 
		if(stem_root_y<0){
		  throw_msg(string("not exists diameter bigger than stem_dia_min"));
		}
		if(fork_y<0){
			throw_msg(string("not exists diameter fork_ratio bigger than threshold"));
		}

		//end_y
		end_y =  hist.size() - stem_root_y-1;

		//统计stem_root_y 到 i 间reversed_hist的数值，得到stem_dia。方法可用平均值，75%分位数
		vector<int> tmp;
		for(size_t j=stem_root_y; j<=i;++j){
			tmp.push_back(reversed_hist[j]);
		}
		sort(tmp.begin(), tmp.end());
		int tar_idx = (int)((float)tmp.size()*0.75);
		stem_dia = tmp[tar_idx];
		*/
	};

	void get_stem_y_fork_rs(
		const std::vector<int>& hist, 
		double ratio,
		int stem_dia_min, 
		int stem_fork_y_min,
		double stem_dia_mp,
		int& fork_y, //output
		int& end_y,  //output
		int& stem_dia,    //output
		int& roi_max_y  //output
		)
	{
		//# 由y的底部向上检测是否有茎,有的话计算移动均值
		//# 通过当前值和移动均值的比值识别分叉点

		//# stem_dia_min = 10 # 茎粗最小值
		//# stem_fork_y_min = 10 # 茎分叉位置与茎根的最小像素高度
		//# stem_dia_mp = 0.8 # moving power

		fork_y = end_y = stem_dia = -1;

		std::vector<int>reversed_hist;
		for(int ii=hist.size()-1; ii>=0;ii--){
			reversed_hist.push_back(hist[ii]);
		}

		int stem_root_y = -1;
		double stem_dia_ma = -1;//# stem diameter moving-average
		int max_y = 0;
		int max_val = reversed_hist[0];

		size_t i=0;
		vector<double> index_of_diachange;	
		vector<double> index_of_diameter;
		for(; i<reversed_hist.size(); ++i){
			//find max value index
			if(reversed_hist[i]>max_val){
				max_val = reversed_hist[i];
				max_y = i;
			}

			if (i==0){
				index_of_diachange.push_back(0.0);
				continue;
			}
			if ((reversed_hist[i-1]<stem_dia_min && reversed_hist[i]>=stem_dia_min && stem_root_y<0)
				||(reversed_hist[i-1]>=stem_dia_min && reversed_hist[i]>=stem_dia_min  && stem_root_y<0))
			{
				//# stem root begin
				stem_root_y = i;
				stem_dia_ma = (double)reversed_hist[i];
				index_of_diachange.push_back(0.0);
			}
			else {
				if (reversed_hist[i-1]>=stem_dia_min &&
					reversed_hist[i]>=stem_dia_min)
				{
					double fork_ratio = reversed_hist[i] / stem_dia_ma;
					if(fork_ratio>1.5){fork_ratio=1.5;}
					index_of_diachange.push_back(fork_ratio);				
					stem_dia_ma = stem_dia_mp *stem_dia_ma + (1.0-stem_dia_mp) * reversed_hist[i];
				}
				else{
					index_of_diachange.push_back(0.0);
				}
			}

		}; 
		if(stem_root_y<0){
		  throw_msg(string("not exists diameter bigger than stem_dia_min"));
		}
	
		//end_y
		end_y =  hist.size() - stem_root_y-1;

		//统计stem_root_y 到 max_y 间reversed_hist的数值，得到stem_dia。方法可用平均值，40%分位数
		vector<int> tmp;
		for(size_t j=stem_root_y; j<max_y;++j){
			if(j<0 || j>=reversed_hist.size()){continue;}
			tmp.push_back(reversed_hist[j]);
		}	
		sort(tmp.begin(), tmp.end());
		int tar_idx = (int)((float)tmp.size()*0.40);
		if(tmp.size()==0 || tar_idx>=tmp.size()){
			throw_msg(string("not exists y fork point"));
		}
		stem_dia = tmp[tar_idx];

		////////////////////////////////////////
		//update 重新检验max_y，因为max_y是全局最大，由于苗的倾斜，可能出现局部最大为分叉点位置
		// 2022-2-21 将比例参数1.5调整至2.0
		int local_max_y = max_y;
		for(size_t j=stem_root_y; j<max_y;++j){
			if((double)(reversed_hist[j]/(double)stem_dia) >=2.0){
				local_max_y = j;
				break;
			}
		}
		max_y = local_max_y ;	
		roi_max_y = hist.size() - max_y-1;
	
		//////////////////////////////////////////////////////////////////
		// 计算vector<double>index_of_diameter 
		for(int idx = 0;idx<reversed_hist.size();++idx){
		
			if(reversed_hist[idx]==0){
				index_of_diameter.push_back(0.0);
			}
			else {
				if(idx>=max_y){
					index_of_diameter.push_back(0.0);
				}
				else{
					double r = yfork_validity_stemdiaratio(stem_dia,reversed_hist[idx],ratio);
					index_of_diameter.push_back(r);
				}
			}
		}
		//////////////////////////////////////////////////////////////////
		// 计算fork_y
		vector<double> index_yfork;
		for(int idx = 0;idx<reversed_hist.size();++idx){
			index_yfork.push_back(index_of_diachange[idx] *	index_of_diameter[idx]);
		}
		vector<double>::iterator max_it = max_element(begin(index_yfork),end(index_yfork));
		int max_idx_y = max_it -  index_yfork.begin();
		fork_y = hist.size() - max_idx_y;
	};
	void get_stem_y_fork_rs_update(
		const Mat& bin_img,
		int stem_x_padding,
		int x0,
		int x1,
		int max_y,
		int stem_root_y,
		int stem_dia,
		bool image_show,
		double cut_point_offset_ratio,
		Point& fork_cent,
		double& max_radius,
		double& stem_angle
		)
	{
		//1 通过bin_img的二值图像，找到茎中心线
		//2 通过中心线，更新二值图像采集区域（有可能倾斜）
		//3 在max_y附近找最大内接圆中心（在茎中心线上）

		fork_cent.x = fork_cent.y = -1;
		max_radius = 0.0;
		stem_angle = 90.0;

		// 茎中心点提取， 在图像max_y 和 stem_root_y之间
		int max_len,start_x,cent_x, pre_val, cur_val;
		vector<int>cent_xs;
		vector<int>cent_ys;
		for(int i=0;i<bin_img.rows;++i){
			if(i<=max_y || i>=stem_root_y){continue;}
			max_len = start_x = cent_x = -1;
			for(int j=x0;j<=x1;++j){
				if(j==x0 && bin_img.at<unsigned char>(i,j)==0){continue;}
				if(j==x0 && bin_img.at<unsigned char>(i,j)>0){start_x=j;continue;}

				pre_val = bin_img.at<unsigned char>(i,j-1);
				cur_val = bin_img.at<unsigned char>(i,j);
				if(pre_val==0 && cur_val>0){
					start_x = j;
				}

				if((pre_val>0 && cur_val==0) ||(j==x1 && cur_val>0)){
					int len = j-start_x +1;
					if(len>max_len){
						max_len = len;
						cent_x = (start_x+j)/2;
						//cout<<i<<"\t"<<"x0="<<start_x<<", x1="<<j<<", centx="<<cent_x<<endl;
					}
				}

			}
			if(cent_x>0){
				cent_xs.push_back(cent_x);
				cent_ys.push_back(i);
			}
		}
		if(cent_xs.size()!=cent_ys.size() || cent_xs.size()<3){
			throw_msg(string("stem length < 3, in function get_stem_y_fork_rs_update()"));
		}

		//茎中心线拟合
		double beta0, beta1;
		beta0 = beta1 = 0.0;
		// x = beta0 + beta1 * y
		double r2 = r2index(cent_ys,cent_xs,0,cent_xs.size()-1, beta0, beta1);
		double th = (double)stem_dia/2.0 + stem_x_padding;
		Mat roi_img = bin_img.clone();

		//提取兴趣区域图片
		stem_angle = atan(beta1)*57.2957795130 + 90.0;
		if(fabs(stem_angle-90.0)<2.0){
			for(int r=0;r<roi_img.rows;++r){		
				for(int c=0;c<roi_img.cols;++c){
					if(c<x0 || c>x1){
						roi_img.at<unsigned char>(r,c) = 0;
					}
				}
			}
		}
		else{
			for(int r=0;r<roi_img.rows;++r){	
				double cx = beta0 + beta1 * r;
				int cx0 = (int)(cx-th);
				int cx1 = (int)(cx+th);				
				for(int c=0;c<roi_img.cols;++c){
					if(c<cx0 || c>cx1){
						roi_img.at<unsigned char>(r,c) = 0;
					}
				}
			}
		}

		////////////////////////////////////////////////////////////
		//
		if(image_show){
			Mat tmp = roi_img.clone();
			Point p0((int)(beta0),0);
			Point p1((int)((double)tmp.rows *beta1 +beta0),tmp.rows);
			line(tmp,p0,p1,Scalar(128,0,0));

			line(tmp,Point(cent_xs[0],cent_ys[0]),
				Point(cent_xs.back(),cent_ys.back()),Scalar(128,0,0));
	

			imshow_wait("rs_bin_roi", tmp);	
		}
		///////////////////////////////////////////////////////////
		//兴趣区域图片findContours
		vector<vector<Point>> contours;
		vector<Vec4i> hierarchy;
		findContours(roi_img,contours,hierarchy,RETR_EXTERNAL,CHAIN_APPROX_NONE);
		//找到周长最长的区域作为目标区域
		int max_length_idx = 0;
		int max_length_val = contours[0].size();
		for(int i=0;i<contours.size();++i){
			if(i==0){continue;}
			if(contours[i].size()>max_length_val){
				max_length_val = contours[i].size();
				max_length_idx = i;
			}
		}

		vector<float>inner_max_radius;
		vector<int> nnindex;
		vector<int> xs;

		/*
		//////////////////////////////////////////
		//基于flann的方法
		vector<Point> contours_serial;
		for(size_t i=0;i<contours.size();++i){
			if(contours_serial.size()==0){
				contours_serial.assign(contours[i].begin(),contours[i].end());
			}
			else{
				for(size_t j=0;j<contours[i].size();++j){
				contours_serial.push_back(contours[i][j]);
				}
			}
		}
		//find inner max radius
		Mat source = Mat(contours_serial).reshape(1);
		source.convertTo(source,CV_32F);
		flann::KDTreeIndexParams indexParams(2);
		flann::Index kdtree(source, indexParams);

		unsigned queryNum=1;
		vector<float> vecQuery(2);
		vector<int> vecIndex(queryNum);
		vector<float> vecDist(queryNum);
		flann::SearchParams params(32);		

		//在茎中心线上计算最优点指标：
		//   1) inner_max_radius  --- 中心点对应的最小内接圆半径
		for(size_t r =0; r<roi_img.rows;++r){	
			if(r < max_y - 50){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(-1);
				continue;
			}
			// x = beta0 + beta1 * y
			float x = (float) (beta0 + beta1 * r);
			int xi = (int)(x+0.5);
			if(xi<0 || xi>=roi_img.cols){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(-1);
				continue;
			}
			if(roi_img.at<unsigned char>(r,xi)==0){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(xi);
				continue;
			}
			vecQuery[0] = x;//x
			vecQuery[1] = (float)r;//y
			kdtree.knnSearch(vecQuery, vecIndex, vecDist,queryNum, params);
			if(vecDist.size()<1 || vecIndex.size()<1){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(xi);
				continue;
			}
			inner_max_radius.push_back(vecDist[0]);
			nnindex.push_back(vecIndex[0]);
			xs.push_back(xi);			
		}
		*/

		////////////////////////////////////////////////////////////////
		//基于pointPolygonTest的方法
		//在茎中心线上计算最优点指标：
		//   1) inner_max_radius  --- 中心点对应的最小内接圆半径
		for(size_t r =0; r<roi_img.rows;++r){	
			if(r < max_y - 50){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(-1);
				continue;
			}
			// x = beta0 + beta1 * y
			float x = (float) (beta0 + beta1 * r);
			int xi = (int)(x+0.5);
			if(xi<0 || xi>=roi_img.cols){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(-1);
				continue;
			}
			if(roi_img.at<unsigned char>(r,xi)==0){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(xi);
				continue;
			}
			float d = pointPolygonTest( contours[max_length_idx], Point2f(x,r), true );			
			if(d<0){
				inner_max_radius.push_back(0.0);
				nnindex.push_back(-1);
				xs.push_back(xi);
				continue;
			}
			inner_max_radius.push_back(d);
			nnindex.push_back(-1);
			xs.push_back(xi);			
		}	

		//find max radius point
		vector<float>::iterator max_it = max_element(inner_max_radius.begin(),inner_max_radius.end());
		int max_idx = max_it - inner_max_radius.begin();
		//基于flann的方法
		//max_radius = sqrt(inner_max_radius[max_idx]);	
		max_radius = inner_max_radius[max_idx];

		//offset  
		double r_offset = cut_point_offset_ratio * max_radius;
		double dy = r_offset * (1.0 / sqrt(1.0 + beta1 * beta1));
		double x_offset =  beta0 + beta1 * (max_idx + dy);
		fork_cent.x = (int)(x_offset+0.5);
		fork_cent.y =(int)(max_idx + dy + 0.5);

		///////////////////////////////////////////////
		//test CForkDetectOptimal 2022-2-17
		/*CForkDetectOptimal fdo(100,20);
		vector<float>opt_max;
		int opt_max_idx = 0;
		double optm=0;
		for(size_t r =0; r<xs.size();++r){	
			if(abs((int)r-max_idx)>50){
				opt_max.push_back(0.0);
				continue;
			}
			else{
				if(xs[r]<0){
					opt_max.push_back(0.0);
					continue;
				}
				Point cent_pt = Point2f(xs[r],r);
				double rt = fdo.center_pt_index(bin_img,cent_pt,-60.0);
				opt_max.push_back((float)rt);
				if(rt>optm){
					optm=rt;
					opt_max_idx=r;
				}
			}
		}*/

		///////////////////////////////////////////////

		///////////////////////////////////////////////
		//test
		if(image_show){
			Mat edge;
			edge = Mat::zeros(roi_img.rows, roi_img.cols, CV_8UC1);
			drawContours(edge, contours, -1, Scalar(255, 255, 0), 1);			
			for(int r=0;r<xs.size();++r){
				int c = xs[r];
				if(c<0){continue;}
				if(inner_max_radius[r]<=0){continue;}
				/*
				//基于flann的方法
				if(nnindex[r]<0){continue;}
				line(edge,Point(c,r),contours_serial[nnindex[r]],Scalar(128,0,0));
				*/
				int rad = int(inner_max_radius[r]);
				line(edge,Point(c-rad,r),Point(c+rad,r),Scalar(80,0,0));		
			}
			Point fork_cent_origin;
			fork_cent_origin.x = xs[max_idx];
			fork_cent_origin.y =max_idx;
			circle(edge,fork_cent_origin,3,Scalar(200,0,0));
			circle(edge,fork_cent_origin,(int)(max_radius),Scalar(200,0,0));

			//circle(edge,Point(xs[opt_max_idx],opt_max_idx),5,Scalar(200,0,0));

			circle(edge,fork_cent,3,Scalar(128,0,0));
			circle(edge,fork_cent,(int)(max_radius),Scalar(128,0,0));
			imshow_wait("rs_bin_roi_radius_circle", edge);	
		}
		///////////////////////////////////////////////   

	};

	//分叉位置指标 yfork_validity_position()
	// max_y---hist中最大点的y-index,
	// end_y ----hist中茎最底端点的y-index
	// fork_y --- 预期分叉点y-index
	// max_y > fork_y > end_y
	double yfork_validity_position(int max_y, int end_y, int fork_y){
		//确定hist最大点位置向下到end_y点最优分叉点位置系数
		double validity = 0.0;
		double opt_pos = 0.85;
		if(fork_y<end_y || fork_y>max_y){return validity;}
		double opt_y = (max_y-end_y+1)*opt_pos + end_y;
		if(fork_y<=end_y || fork_y>=max_y){return validity;}
		if(fork_y>end_y && fork_y<opt_y){
			validity =(double)( (fork_y-end_y+1) / (opt_y-end_y+1));
		}
		else{

			validity =(double)( (max_y -fork_y+1) / (max_y -opt_y+1));
		}
		return validity;
	}

	//茎粗比指标  yfork_validity_stemdiaratio()
	// stem_dia---茎粗,
	// dia_y ----hist中y-index对应的茎粗
	double yfork_validity_stemdiaratio(int stem_dia, int dia_y,double opt_dia_ratio){
		//确定hist最大点位置向下到end_y点最优分叉点位置系数
		double validity = 0.0;
		//double opt_pos = 1.2;
		double opt_dia = stem_dia*opt_dia_ratio;
		if(dia_y<opt_dia){return dia_y/opt_dia;}
		else{
			return opt_dia/dia_y;
		}
	}
	//计算上下窗口茎粗变化率
	void yfork_validity_stemdiachange(
			const std::vector<int>& hist,
			int stem_dia_min,
			std::vector<double>& var_ratio		
			)
	{
		var_ratio.clear();
		int smooth_width = 5;
		double mean_pre=0.0;
		double mean_nxt=0.0;

		int low_y,up_y;
		low_y=up_y=-1;
		for(int i=0;i<hist.size();++i){
			if(hist[i]>=stem_dia_min){
				if( low_y<0){low_y=i;}
				if(i>up_y){up_y = i;}
			}
		}
		for(int i=0;i<hist.size();++i){
			if(i-smooth_width<0 || i+smooth_width>=hist.size()){
				var_ratio.push_back(0.0);
				continue;
			}
			if(i<low_y+smooth_width || i>up_y-smooth_width){
				var_ratio.push_back(0.0);
				continue;
			}

			mean_pre = mean_nxt=0.0;
			//low sum
			for(int di = -smooth_width;di<0;++di){
				mean_pre += hist[i+di];
			}
			//up sum
			for(int di = 0;di<smooth_width;++di){
				mean_nxt += hist[i+di];
			}
			if(mean_pre<1.0e-5){
				var_ratio.push_back(0.0);
				continue;
			}
			else{
				double ratio = mean_nxt / mean_pre;
				if(ratio>3.0){ratio=3.0;}
				var_ratio.push_back(ratio);
			}
		}

	}

	void imshow_wait(
		const char* winname,
		Mat&img, 
		int waittype/*=-1*/
		)
	{	
		namedWindow(winname, CV_WINDOW_NORMAL);
		imshow(winname, img);
		waitKey(waittype);
	};

	int get_stem_fork_right(
		Mat&stem_img, 
		int stem_fork_y, 
		int stem_x0, 
		int stem_x1,
		int detect_window
		)
	{
		/*
		通过茎分叉坐标找到右侧茎分叉点坐标
		stem_img, 二值图像，茎部分的局部图像
		stem_fork_y, 检测到的茎分叉的y坐标
		stem_x0, 茎左侧边缘（通过hist检测，针对stem_fork_y不一定精确）
		stem_x1， 茎右侧边缘（通过hist检测，针对stem_fork_y不一定精确）
		*/

    
		// 注释掉 2021-11-3 开始
		/*double min_dist = 1.0;
		int min_x1 = -1;
		for(size_t i=stem_x0; i<stem_x1+2*detect_window; ++i){
			double left_cnt = 0.0;
			double right_cnt = 0.0;
			for(size_t j=i-detect_window+1; j<i+1; ++j){
				if (stem_img.at<unsigned char>(stem_fork_y,j)>0){
					left_cnt+=1;
				}
			}
			for(size_t j=i+1; j<i+detect_window+1; ++j){
				if( stem_img.at<unsigned char>(stem_fork_y,j)==0){
					right_cnt+=1;
				}
			}

			if( left_cnt==0 || right_cnt ==0){continue;}
			double dist = 1.0 - min(left_cnt/right_cnt, right_cnt/left_cnt);
			if (dist < min_dist){
				min_dist = dist;
				min_x1 = i;
			}
		}
		return min_x1;
		*/
		// 注释掉 2021-11-3 结束
		size_t left_min = stem_x0-4*detect_window;
		size_t right_max = stem_x1+4*detect_window;
		if(left_min<0){left_min = 0;}
		if(right_max>stem_img.cols){right_max = stem_img.cols;}

		int max_len = -1;
		int max_idx = -1;
		int start_x=left_min;
		for(size_t i=left_min; i<right_max; ++i){
			if(i==left_min){ 
				if(stem_img.at<unsigned char>(stem_fork_y,i)>0){
					start_x =i;					
				}
				continue;
			}			
			if (stem_img.at<unsigned char>(stem_fork_y,i-1)==0 && 
				stem_img.at<unsigned char>(stem_fork_y,i)>0){
					start_x =i;
					continue;
			}
		
			if ((stem_img.at<unsigned char>(stem_fork_y,i-1)>0 && stem_img.at<unsigned char>(stem_fork_y,i)==0) ||
				(stem_img.at<unsigned char>(stem_fork_y,i)>0 && i==right_max-1)){
				if(((int)i-start_x) > max_len){
					max_len = i-start_x;
					max_idx = i;
				}
			}		
		}
		if(max_idx<0){max_idx = stem_x1;}
		return max_idx;
	};

	int get_stem_fork_left(
		Mat&stem_img, 
		int stem_fork_y, 
		int stem_x0, 
		int stem_x1,
		int detect_window
		)
	{
		/*
		通过茎分叉坐标找到右侧茎分叉点坐标
		stem_img, 二值图像，茎部分的局部图像
		stem_fork_y, 检测到的茎分叉的y坐标
		stem_x0, 茎左侧边缘（通过hist检测，针对stem_fork_y不一定精确）
		stem_x1， 茎右侧边缘（通过hist检测，针对stem_fork_y不一定精确）
		*/    
	
		size_t left_min = stem_x0-8*detect_window;
		size_t right_max = stem_x1+8*detect_window;
		if(left_min<0){left_min = 0;}
		if(right_max>stem_img.cols-1){right_max=stem_img.cols;}

		int max_len = -1;
		int max_idx = -1;
		int start_x=left_min;
		for(size_t i=left_min; i<right_max; ++i){
			if(i==left_min){ 
				if(stem_img.at<unsigned char>(stem_fork_y,i)>0){
					start_x =i;					
				}
				continue;
			}			
			if (stem_img.at<unsigned char>(stem_fork_y,i-1)==0 && 
				stem_img.at<unsigned char>(stem_fork_y,i)>0)
			{
					start_x =i;
					continue;
			}		
			if ((stem_img.at<unsigned char>(stem_fork_y,i-1)>0 && stem_img.at<unsigned char>(stem_fork_y,i)==0) ||
				(stem_img.at<unsigned char>(stem_fork_y,i)>0 && i==right_max-1)){
				if(((int)i-start_x) > max_len){
					max_len = i-start_x;
					max_idx = start_x;
				}
			}		
		}
		if(max_idx<0){max_idx = stem_x0;}
		return max_idx;
	}

	void get_stem_fork_xs(
		Mat&stem_img, 
		int stem_fork_y, 
		int stem_x0, 
		int stem_x1,	
		int x0, 
		int x1,	
		int & fork_x0,
		int & fork_x1
		)
	{
		size_t left_min = x0;
		size_t right_max = x1;
		if(left_min<0){left_min = 0;}
		if(right_max>=stem_img.cols){right_max = stem_img.cols-1;}

		int max_len = -1;
		int max_idx = -1;
		int max_idx_right = -1;

		int start_x=left_min;
		for(size_t i=left_min; i<right_max; ++i){
			if ( i==0 && 
				stem_img.at<unsigned char>(stem_fork_y,i)>0)
			{
					start_x =i;
					continue;
			}
			if (stem_img.at<unsigned char>(stem_fork_y,i)==0 &&
				stem_img.at<unsigned char>(stem_fork_y,i+1)>0)
			{
					start_x =i;
			}
		
			if ((stem_img.at<unsigned char>(stem_fork_y,i)>0 && stem_img.at<unsigned char>(stem_fork_y,i+1)==0) ||
				(stem_img.at<unsigned char>(stem_fork_y,i)>0 && i==right_max-1)){
				if(((int)i-start_x) > max_len)
				{
					max_len = i-start_x;
					max_idx = start_x;
					max_idx_right = i;
				}
			}		
		}
		if(max_idx<0){max_idx = stem_x0;}
		fork_x0 = max_idx;	
		if(max_idx_right<0){max_idx_right = stem_x1;}
		fork_x1 = max_idx_right;
	}

	// 根据边缘上的一点，爬取下一点
	void get_next_pt(
		Mat& edge_img,
		gcv_point<int>&start_pt, 
		gcv_point<int>&pre_pt, 	
		gcv_point<int>&nxt_pt,//output
		bool is_up //true 向上爬取， false 向下爬取
		)
	{
		nxt_pt.x=-1;
		nxt_pt.y=-1;
		int dy_range[3] = {-1,0,1};
		int dx_range[3] = {1,0,-1};
		if (!is_up){
			dy_range[0] = 1;
			dy_range[2] = -1;
		}
    
		for (int dyi =0; dyi<3;++dyi){
			int dy = dy_range[dyi];
			bool break_x = false;
			for( int dxi=0;dxi<3;++dxi){
				int dx = dx_range[dxi];
				int x = start_pt.x+dx;
				int y = start_pt.y+dy;
				if (dx==0 && dy == 0){continue;}
				if (x==pre_pt.x && y==pre_pt.y){continue;}
				if (edge_img.at<unsigned char>(y,x)>0){
					nxt_pt.x = x;
					nxt_pt.y = y;
					break_x = true;
					break;
				}
			}
			if (break_x){
				break;
			} 
		}
	}

	// qua_fitting
	// input: xx  ---- sorted x values, ascending
	//        yy  ----- y values
	//
	// return: oa --- optimal angle, which is the x coordinates of center line of quadratic fitting line
	double qua_fitting(vector<double>&xx, vector<double>&yy)
	{
		double oa = 0.0;
		// fit quadratic function with opencv
		Mat A = cv::Mat::zeros(cv::Size(3,xx.size()), CV_64FC1);
		for(int i=0; i<xx.size();++i){
			A.at<double>(i,0) = 1;
			A.at<double>(i,1) = xx[i];
			A.at<double>(i,2) = xx[i] * xx[i];
		}

		Mat b = cv::Mat::zeros(cv::Size(1,yy.size()), CV_64FC1);
		for(int i = 0; i<yy.size(); ++i){
			b.at<double>(i,0) = yy[i];
		}
		Mat c, d;
		c = A.t() * A;
		d = A.t() * b;
		Mat result = cv::Mat::zeros(cv::Size(1,3),CV_64FC1);
		cv::solve(c,d,result);
	

		double a0 = result.at<double>(0, 0);
		double a1 = result.at<double>(1, 0);
		double a2 = result.at<double>(2, 0);
		if(fabs(a2)<1.0e-6){
			throw_msg(string("Error:  a2 = 0 in solver"));
		}
		oa = -a1 / a2 / 2;
		if(oa >180.0){oa -= 180.0;}
		return oa;
	}

};