graft_cv/optimal_angle.cpp

//#include "StdAfx.h"
#include <opencv2\imgproc\imgproc.hpp>
#include "optimal_angle.h"
#include "utils.h"
#include <string.h>

namespace graft_cv{

//COptimalAngle::COptimalAngle(ConfigParam&cp, CGcvLogger* pLog)
//	:m_cparam(cp),
//	m_pLogger(pLog),
//	m_pImginfo(0),
//	m_pImginfoBase(0),
//	m_height(0),
//	m_width(0),
//	m_patch_id(""),
//	m_imgIndex(0),
//	m_imgId(""),
//	m_ppImgSaver(0)
//{
//}
//COptimalAngle::~COptimalAngle(void)
//{
//	this->clear_imginfo();
//}
//void COptimalAngle::clear_imginfo(){
//	if (m_pImginfo){
//		imginfo_release(&m_pImginfo);
//		m_pImginfo=0;
//	}
//	if (m_pImginfoBase){
//		imginfo_release(&m_pImginfoBase);
//		m_pImginfoBase=0;
//	}
//}
//int COptimalAngle::reset()
//{
//	if(m_pLogger){
//		m_pLogger->DEBUG("optimal angle reset begin.");
//	}
//	m_fork_ys.clear();
//	m_end_ys.clear();
//	m_widths.clear();
//	m_an2width.clear();
//	m_an2widthBase.clear();
//	m_width=0;
//	m_height=0;
//	m_patch_id="";
//	m_imgIndex=0;
//	if(m_pLogger){
//		m_pLogger->DEBUG("optimal angle reset finished.");
//	}
//	return 0;
//}
//
//int COptimalAngle::append(
//	ImgInfo* imginfo,
//	PositionInfo& posinfo	
//	)
//{	
//	//clear opencv windows
//	if(m_cparam.image_show){destroyAllWindows();}
//
//	clock_t t;
//	clock_t t0 = clock();
//	if(m_imgIndex==0){
//		//new patch
//		m_patch_id = getImgId(img_type::oa);	
//	}
//	m_imgId = getImgIdOa(m_patch_id,m_imgIndex);	
//	m_imgIndex+=1;
//
//	if(m_pLogger){
//		m_pLogger->INFO(m_imgId + " append image begin.");
//	}	
//	if(!isvalid(imginfo)){
//		throw_msg(m_imgId+" invalid input image");
//	}
//	if(m_width==0 || m_height==0){
//		m_width = imginfo->width;
//		m_height = imginfo->height;
//	}
//	int angle = imginfo->angle;
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_imgId<<" append image, angle="<<angle
//			<<"\twidth="<<imginfo->width
//			<<"\theight="<<imginfo->height;
//		m_pLogger->INFO(buff.str());
//	}
//	Mat img = imginfo2mat(imginfo);	
//	
//	//imshow_wait("src", img);
//	if(m_cparam.self_camera){
//		image_set_bottom(img,20,8);
//		if(m_pLogger){				
//			m_pLogger->DEBUG(m_imgId+" oa_append image set bottom with pixel value 20.");				
//		}
//	}
//	if(m_cparam.oa_y_flip){
//		flip(img,img,0);
//		if(m_pLogger){				
//			m_pLogger->DEBUG(m_imgId+" oa_append image y fliped.");				
//		}
//	}
//	if(m_ppImgSaver && *m_ppImgSaver){
//		(*m_ppImgSaver)->saveImage(img, m_imgId);
//		if(m_pLogger){				
//			m_pLogger->DEBUG(m_imgId+" oa_append after image save.");				
//		}
//	}
//	if(m_cparam.image_show){
//		Mat tmp_b = img.clone();
//		cv::line(tmp_b,Point(0,m_cparam.oa_clip_y_min), Point(tmp_b.cols,m_cparam.oa_clip_y_min),Scalar(200),3);
//		cv::line(tmp_b,Point(0,m_cparam.oa_clip_y_max), Point(tmp_b.cols,m_cparam.oa_clip_y_max),Scalar(200),3);			
//		imshow_wait("oa_sub_image_range", tmp_b);
//	}
//
//	Rect rect_rs(Point(0,0),Point(img.cols,m_cparam.oa_clip_y_min));
//	//Rect rect_base(Point(0,m_cparam.oa_clip_y_max),Point(img.cols,img.rows));
//
//	Mat rs_img = img(rect_rs);
//	//Mat bs_img = img(rect_base);
//
//	
//	// 1 rootstock image processing 
//	// 1.1 calculate the image, get binary image and get object width
//	int min_idx, max_idx;
//	vector<int> hist_col;
//	int width = imgproc(rs_img, hist_col, min_idx, max_idx);
//	m_an2width.insert(make_pair<int,int>(angle, width));
//	posinfo.rs_oa_width=width;
//
//	if(m_cparam.image_show){			
//		Mat hist_col_img;
//		hist2image(hist_col,hist_col_img,1,m_cparam.image_col_grid,m_cparam.image_row_grid);
//		imshow_wait("oa_hist_col", hist_col_img);			
//	}
//
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_imgId<<" append image, rootstock width="<<width;
//		m_pLogger->INFO(buff.str());
//	}
//
//	if(!m_cparam.image_show){
//		t = clock();		
//		if(1000.0*((float)(t-t0))/CLOCKS_PER_SEC>(float)m_cparam.timeout_append){
//			throw_msg(m_imgId+" time out");
//		}
//	}
//
//	//1.2 茎粗计算
//	//1.2.1 茎x方向位置：x方向位置范围	
//	int r0,r1;
//	r0=r1=-1;	
//	std::vector<int>hist_row_leaf;
//	mat_histogram(m_binImg,hist_row_leaf,1);
//	for(size_t i=0;i<hist_row_leaf.size();++i){
//		if(hist_row_leaf[i]>=m_cparam.oa_min_hist_value && r0<0){r0=i;}
//		if(hist_row_leaf[i]>=m_cparam.oa_min_hist_value){r1=i;}
//	}
//	if(r0<0 || r1<0){throw_msg(string(m_imgId+" invalid image, with no object in binary image"));}
//		
//	vector<int> hist_col_yfork;
//	mat_histogram_yfork(m_binImg,hist_col_yfork,r0,r1);
//	if(m_cparam.image_show){			
//		Mat tmp;
//		hist2image(hist_col_yfork,tmp,1,m_cparam.image_col_grid,m_cparam.image_row_grid);
//		imshow_wait("oa_hist_col_yfork", tmp);			
//	}	
//
//	int x0,x1,stem_x0, stem_x1;
//	int centx = (int)((min_idx+ max_idx)/2.0);
//	get_stem_x_range_oa(
//		hist_col_yfork,
//		m_cparam.oa_col_th_ratio,
//		m_cparam.oa_stem_x_padding,
//		centx,
//		width,
//		x0,
//		x1,
//		stem_x0,
//		stem_x1);
//
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_imgId<<" append image, x0="<<x0
//			<<"\tx1="<<x1
//			<<"\tstem_x0="<<stem_x0
//			<<"\tstem_x1="<<stem_x1;
//		m_pLogger->INFO(buff.str());
//	}
//
//	if(m_cparam.image_show){
//		Mat tmp_b = m_binImg.clone();
//		cv::line(tmp_b,Point(x0,0), Point(x0,tmp_b.rows),Scalar(200),3);
//		cv::line(tmp_b,Point(x1,0), Point(x1,tmp_b.rows),Scalar(200),3);			
//		imshow_wait("oa_x_range", tmp_b);
//	}
//
//	if(!m_cparam.image_show){
//		t = clock();
//		if(1000.0*((float)(t-t0))/CLOCKS_PER_SEC>(float)m_cparam.timeout_append){
//			throw_msg(m_imgId+" time out");
//		}		
//	}
//
//	vector<int> hist_row;
//	mat_histogram(m_binImg,hist_row,1,-1,-1,x0,x1+1);
//
//	if(m_cparam.image_show){		
//		Mat hist_row_img;
//		hist2image(hist_row,hist_row_img, 0,m_cparam.image_col_grid,m_cparam.image_row_grid);		
//		imshow_wait("oa_hist_row", hist_row_img);			
//	}
//	// 1.2.2 茎分叉点检测，y方向检测
//	int stem_fork_y=-1,stem_end_y=-1, stem_dia=-1;
//	get_stem_y_fork(
//		hist_row,
//		m_cparam.oa_row_th_ratio,
//		m_cparam.oa_stem_dia_min,
//		m_cparam.oa_stem_fork_y_min,
//		m_cparam.oa_stem_dia_mp,
//		stem_fork_y,
//		stem_end_y,
//		stem_dia);
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_imgId<<" append image, stem_fork_y="<<stem_fork_y
//			<<"\tstem_end_y="<<stem_end_y
//			<<"\tstem_dia="<<stem_dia;
//		m_pLogger->INFO(buff.str());
//	}
//	
//
//	if(m_cparam.image_show){
//		Mat tmp_bin = m_binImg.clone();
//		cv::line(tmp_bin,Point(min_idx,0), Point(min_idx,tmp_bin.rows),Scalar(200),3);
//		cv::line(tmp_bin,Point(max_idx,0), Point(max_idx,tmp_bin.rows),Scalar(200),3);		
//		cv::line(tmp_bin,Point(x0,stem_fork_y), Point(x1,stem_fork_y),Scalar(200),3);
//		cv::line(tmp_bin,Point(x0,stem_end_y), Point(x1,stem_end_y),Scalar(200),3);
//		imshow_wait("oa_result", tmp_bin);
//	}
//	
//	m_fork_ys.push_back(stem_fork_y);
//	m_end_ys.push_back(stem_end_y);
//	m_widths.push_back(width);
//
//
//	//2 基质图像处理
//	//int min_idx_base, max_idx_base;
//	//vector<int> hist_col_base;
//	//int width_base = imgprocBase(bs_img, hist_col_base, min_idx_base, max_idx_base);
//	//m_an2widthBase.insert(make_pair<int,int>(angle, width_base));
//	//posinfo.rs_oa_width_base = width_base;
//
//	//if(m_pLogger){
//	//	stringstream buff;
//	//	buff<<m_imgId<<" append image, base width="<<width_base;
//	//	m_pLogger->INFO(buff.str());
//	//
//
//	// 3 返回结果
//	double fork_y = (double)stem_fork_y;
//	double end_y = (double)stem_end_y;
//
//	fork_y = ((double)m_height/2.0 - fork_y)*m_cparam.rs_oa_pixel_ratio;
//	end_y =  ((double)m_height/2.0 - end_y)*m_cparam.rs_oa_pixel_ratio;
//	posinfo.rs_oa_stem_y_fork=fork_y;
//    posinfo.rs_oa_clamp_y_end=end_y;
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_imgId<<" append image, stem_fork_y(mm)="<<fork_y
//			<<"\tstem_end_y(mm)="<<end_y;
//		m_pLogger->INFO(buff.str());
//	}
//	//4 返回图像
//	if(m_cparam.image_return){
//		this->clear_imginfo();
//		cv::line(m_binImg,Point(min_idx,0), Point(min_idx,rs_img.rows),Scalar(200),3);
//		cv::line(m_binImg,Point(max_idx,0), Point(max_idx,rs_img.rows),Scalar(200),3);
//
//		cv::line(m_binImg,Point(x0,stem_fork_y), Point(x1,stem_fork_y),Scalar(200),3);
//		cv::line(m_binImg,Point(x0,stem_end_y), Point(x1,stem_end_y),Scalar(200),3);
//
//		clear_imginfo();
//		m_pImginfo = mat2imginfo(m_binImg);
//		//m_pImginfoBase = mat2imginfo(m_binImgBase);
//		posinfo.pp_images[0]=m_pImginfo;
//		//posinfo.pp_images[1]=m_pImginfoBase;
//
//		if(m_ppImgSaver && *m_ppImgSaver){
//			(*m_ppImgSaver)->saveImage(m_binImg, m_imgId+"_rst_0");			
//		}
//	}
//	if(m_pLogger){
//		m_pLogger->INFO(m_imgId + " append image finished.");
//	}	
//	return 0;
//}
//
//double COptimalAngle::infer(PositionInfo& posinfo){
//	
//	if(m_pLogger){
//		m_pLogger->INFO(m_patch_id + " optimal angle infer begin.");
//	}	
//	double oa = this->angle_fit(this->m_an2width);
//	//double oa_base = 0.0;
//	//try{
//	//	oa_base = this->angle_fit_base(this->m_an2widthBase);
//	//}
//	//catch(...){
//	//	if(m_pLogger){		
//	//		m_pLogger->WARNING("angle_fit_base() error");
//	//	}
//	//}
//	posinfo.rs_oa = oa;
//	//posinfo.rs_oa_base=oa_base;
//
//	if(m_fork_ys.size()==0 ||m_end_ys.size()==0){
//		throw_msg(m_patch_id+ " empty fork_ys or end_ys, NEED append image");
//	}
//	vector<size_t>des_idx = sort_indexes_e(m_widths,false);	
//	int e_idx = (int)((float)m_end_ys.size() * 0.5);	
//	sort(m_end_ys.begin(), m_end_ys.end());
//	double fork_y = ((m_fork_ys[des_idx[0]] +  m_fork_ys[des_idx[1]])/2.0);
//	double end_y = m_end_ys[e_idx];
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_patch_id<<" angle fit result, stem_fork_y(pixel)="<<(int)(fork_y)
//			<<"\tstem_end_y(pixel)="<<(int)(end_y);
//		m_pLogger->INFO(buff.str());
//	}
//
//	fork_y = ((double)m_height/2.0 - fork_y)*m_cparam.rs_oa_pixel_ratio;
//	end_y =  ((double)m_height/2.0 - end_y)*m_cparam.rs_oa_pixel_ratio;
//	posinfo.rs_oa_stem_y_fork=fork_y;
//    posinfo.rs_oa_clamp_y_end=end_y;
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_patch_id<<" angle fit result, stem_fork_y(mm)="<<fork_y
//			<<"\tstem_end_y(mm)="<<end_y;
//		m_pLogger->INFO(buff.str());
//	
//		m_pLogger->INFO(m_patch_id + " optimal angle infer finished.");
//	}	
//	return oa;
//}
//
//int COptimalAngle::imgproc(
//	Mat& img, 
//	vector<int>& hist,
//	int& min_idx, 
//	int& max_idx
//	)
//{
//	//int morph_size = 1;
//	//int min_hist_value = 10;
//	hist.clear();
//
//	Mat b_img;
//	double th = threshold(
//		img, 
//		b_img, 
//		255,
//		255,
//		THRESH_OTSU);//+THRESH_BINARY_INV
//
//	Mat kernel = getStructuringElement(
//		MORPH_RECT, 
//		Size( 2*m_cparam.oa_morph_radius+1, 2*m_cparam.oa_morph_radius+1), 
//		Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
//		);
//    
//    morphologyEx(
//		b_img, 
//		m_binImg, 
//		MORPH_CLOSE,
//		kernel,
//		Point(-1,-1),
//		m_cparam.oa_morph_iteration);
//
//	
//	
//	mat_histogram(m_binImg, hist, 0);		
//
//	/*int min_idx, max_idx;*/
//	min_idx = hist.size();
//	max_idx = 0;
//	vector<int>::const_iterator it0 = hist.begin();
//	for(vector<int>::const_iterator it = hist.begin();it!=hist.end();++it){
//		if(*it>=m_cparam.oa_min_hist_value){
//			int idx = it - it0;
//			if(idx<min_idx){min_idx = idx;}
//			if(idx > max_idx) {max_idx = idx;}
//		}
//	}
//	if(max_idx<=min_idx){
//		throw_msg(m_imgId+" invalid binary image, not exist valid objects");		
//	}
//	if(m_cparam.image_show){		
//		Mat tmp = m_binImg.clone();
//		cv::line(tmp,Point(min_idx,0), Point(min_idx,tmp.rows),Scalar(156),3);
//		cv::line(tmp,Point(max_idx,0), Point(max_idx,tmp.rows),Scalar(156),3);
//		imshow_wait("oa_bin", tmp);		
//	}
//    int width = max_idx-min_idx+1;	
//	return width;
//}
//
//int COptimalAngle::imgprocBase(
//	Mat& img, 
//	vector<int>& hist,
//	int& min_idx, 
//	int& max_idx
//	)
//{
//	//int morph_size = 1;
//	//int min_hist_value = 10;
//	hist.clear();
//
//	Mat b_img;
//	double th = threshold(
//		img, 
//		b_img, 
//		255,
//		255,
//		THRESH_OTSU+THRESH_BINARY_INV);
//
//	Mat kernel = getStructuringElement(
//		MORPH_RECT, 
//		Size( 2*m_cparam.oa_morph_radius_base+1, 2*m_cparam.oa_morph_radius_base+1), 
//		Point( m_cparam.oa_morph_radius_base, m_cparam.oa_morph_radius_base)
//		);
//    
//    morphologyEx(
//		b_img, 
//		m_binImgBase, 
//		MORPH_OPEN,
//		kernel,
//		Point(-1,-1),
//		m_cparam.oa_morph_iteration_base);
//
//	if(m_cparam.image_show){	
//		destroyAllWindows();
//		imshow_wait("oa_bin_base", m_binImgBase);		
//	}
//	
//	mat_histogram(m_binImgBase, hist, 0);	
//
//	/*int min_idx, max_idx;*/
//	min_idx = hist.size();
//	max_idx = 0;
//	vector<int>::const_iterator it0 = hist.begin();
//	for(vector<int>::const_iterator it = hist.begin()+5;it!=hist.end();++it){
//		if(*it>=m_cparam.oa_min_hist_value_base){
//			int idx = it - it0;
//			if(idx<min_idx){min_idx = idx;}
//			if(idx > max_idx) {max_idx = idx;}
//		}
//	}
//	if(max_idx<=min_idx){
//		throw_msg(m_imgId+" invalid binary image, not exist valid objects");		
//	}
//    int width = max_idx-min_idx+1;	
//	return width;
//}
////angle_fit（）方法
////append()加载的图片，满足下述条件：0-180度采样，其中0和180度都有数值
//double COptimalAngle::angle_fit(map<int,int>& an2width,bool is_base)
//{
//	if(an2width.size()<3){
//		throw_msg(m_patch_id+" not enough valid images, NEED append image");
//	}
//	double oa = 0.0;
//	std::vector<double> x;//angle
//	std::vector<double> y;//width
//
//	map<int,int>::iterator it = an2width.begin();
//	for(it=an2width.begin(); it!=an2width.end(); ++it)
//	{
//		x.push_back((double)(it->first));
//		y.push_back((double)(it->second));
//	}
//	//mean value of first and latest y
//	double y_mu = 0.5*(y[0]+ y[y.size()-1]);
//	y[0] = y[y.size()-1] = y_mu;
//	size_t ap_times = x.size()-1;
//	for (size_t i=0; i<ap_times;++i)
//	{
//		x.push_back(x[i+1]+180.0);
//		y.push_back(y[i+1]);
//	}
//	vector<double>::iterator smallest = min_element(begin(y), end(y));
//	size_t min_idx, max_idx;
//	min_idx = distance(y.begin(), smallest);
//
//	for(size_t i = min_idx+1;i<y.size();++i){
//		if(fabs(y[i]-*smallest)<0.1){
//			max_idx = i;
//			break;
//		}
//	}
//
//	// slice x and y the convex segment
//	vector<double>xx;
//	vector<double>yy;
//	for(size_t i=min_idx;i<=max_idx;++i){
//		xx.push_back(x[i]);
//		yy.push_back(y[i]);
//	}
//
//	// 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(m_patch_id+" a2 = 0 in solver");
//	}
//	oa = -a1 / a2 / 2;
//	if(oa >180.0){oa -= 180.0;}
//	return oa;
//}
//double COptimalAngle::angle_fit_base(map<int, int>&){
//	return 0.0;
//}
//
/////////////////////////////////////////////////////////////////////////////////
//COptimalAnglePart::COptimalAnglePart(ConfigParam&cp, CGcvLogger* pLog)
//	:COptimalAngle(cp,pLog)
//{
//}
//COptimalAnglePart::~COptimalAnglePart()
//{
//	COptimalAngle::clear_imginfo();
//}
//// append()加载的图片，满足不是0-180度采样，但保证部分采样中一定存在最大值，
//// 直接利用这些数值进行二次函数的拟合
//double COptimalAnglePart::angle_fit(map<int,int>& an2width,bool is_base)
//{ 
//	
//	double oa = 0.0;
//	std::vector<double> xx;//angle
//	std::vector<double> yy;//width
//
//	map<int,int>::iterator it = an2width.begin();
//	for(it=an2width.begin(); it!=an2width.end(); ++it)
//	{
//		xx.push_back((double)(it->first));
//		yy.push_back((double)(it->second));
//	}
//	//sort 
//	for(size_t i =0;i<xx.size()-1;++i){
//		for(size_t j = i+1; j<xx.size();++j){
//			if(xx[j]<xx[i]){
//				double tmp = xx[i];
//				xx[i]=xx[j];
//				xx[j]=tmp;
//				tmp = yy[i];
//				yy[i]=yy[j];
//				yy[j]=tmp;
//			}
//		}
//	}
//
//	if(m_pLogger){			
//		stringstream buff;
//		buff<<m_patch_id;
//		if(is_base){
//			buff<<" rootstock_base ";
//		}
//		else{
//			buff<<" rootstock_leaf ";
//		}
//		buff<<"opt-angle, (xi,yi): ";
//		for(size_t i=0;i<xx.size();++i){
//			buff<<"("<<xx[i]<<","<<yy[i]<<"), ";
//		}		
//		m_pLogger->INFO(buff.str());
//	}
//
//	if(an2width.size()<3){
//		throw_msg(m_patch_id+" not enough valid images, NEED append image");
//	}
//	// angle field check
//	if((xx.back()-xx.front())<90){
//		if(m_pLogger){			
//			stringstream buff;
//			buff<<m_patch_id<<" x_end="<<xx.back()<<" x_front="<<xx.front();				
//			buff<<" not enough angle field, small than 90 degree.";
//			m_pLogger->ERRORINFO(buff.str());
//		}
//		throw_msg(m_patch_id+" not enough angle field, small than 90 degree");
//	}
//	oa = opt_angle_part_impl(xx,yy,is_base);
//	if(m_pLogger){
//		stringstream buff;
//		buff<<m_patch_id;
//		if(is_base){
//			buff<<" rootstock_base ";
//		}
//		else{
//			buff<<" rootstock_leaf ";
//		}
//        buff<<" opt-angle: "<<oa;
//		m_pLogger->INFO(buff.str());
//	}
//	return oa;	
//}
//double COptimalAnglePart::angle_fit_base(map<int, int>& an2width){
//	
//	double oa = 0.0;
//	//find the maximun angle
//	oa = angle_fit(an2width,true);
//
//	//get the minimum angle
//	if(oa>45.0){oa -= 45.0;}
//	else{oa +=45.0;}
//	return oa;	
//}
//void COptimalAnglePart::zero_cross_detect(
//	vector<int>&d_sign, //input
//	bool& is_local_max, // whether exists local max 
//	int& local_limit_idx // local limit index of original vector
//)
//{
//	int min_zc_idx = -1;
//    int max_zc_idx = -1;
//	bool b_local_min =false;
//	bool b_local_max =false;
//    bool is_all_down=true;
//    bool is_all_up = true;
//    for (size_t i=0; i<d_sign.size();++i){
//		if (d_sign[i] >0){is_all_down=false;}
//		if (d_sign[i] <0){is_all_up=false;}
//		if (i == d_sign.size()-1){continue;}
//
//        if (d_sign[i] >=0 && d_sign[i+1] <=0){
//			//is_local_max=true;
//			max_zc_idx = i+1;
//			b_local_max=true;
//            
//		}
//        if (d_sign[i] <=0 && d_sign[i+1] >=0){
//			//is_local_max=false;
//            min_zc_idx=i+1;
//			b_local_min=true;
//		}
//	}
//    if(b_local_min){
//		is_local_max=false;
//        local_limit_idx = min_zc_idx;
//		return;
//	}
//	if(b_local_max){
//		is_local_max=true;
//        local_limit_idx = max_zc_idx;
//		return;
//	}
//    if (is_all_down){
//		is_local_max=true;
//        local_limit_idx=0;
//        return;
//	}
//    if( is_all_up){
//		is_local_max=true;
//        local_limit_idx=d_sign.size();
//        return;        
//	}
//	is_local_max=true;
//    local_limit_idx=-1;
//    return;
//}
//double COptimalAnglePart::limit_min_angle(
//		vector<double>&x, //sorted, ascending
//		vector<double>&y,
//		size_t min_idx
//		)
//{
//	size_t i = min_idx;
//	if(fabs(x[i+1]-x[i])<1.0e-6 || fabs(x[i]-x[i-1])<1.0e-6){
//	 throw_msg(m_patch_id+" limit_min_angle, diff_x < 1.0e-6");
//	}
//    double k_r = (y[i+1]-y[i])/(x[i+1]-x[i]);
//    double k_l = (y[i]-y[i-1])/(x[i]-x[i-1]);
//    double oa = 0.0;
//	double b_r = 0.0;
//	double b_l = 0.0;
//	double cross_x = 0;
//    if (fabs(k_r) > fabs(k_l)){
//        // right line
//        b_r = y[i]-k_r*x[i];
//        b_l = y[i-1] + k_r*x[i-1];
//        cross_x = line_cross(k_r,b_r,-k_r,b_l);
//        oa = cross_x;
//	}
//    else{
//        //# left line
//        b_l = y[i]-k_l*x[i];
//        b_r =  y[i+1] + k_l*x[i+1];
//        cross_x = line_cross(k_l,b_l,-k_l,b_r);
//        oa = cross_x;
//	}
//    return oa;
//}
//
//double COptimalAnglePart::line_cross(double k1,double b1,double k2,double b2)
//{
//	double x = (b2-b1) / (k1-k2);
//    return x;
//}
//
//double COptimalAnglePart::opt_angle_part_impl(
//		vector<double>&x, //sorted, ascending
//		vector<double>&y,
//		bool is_base)
//{
//	// sign of difference of y
//	vector<int> dy_sign;
//	for(size_t i=1; i<y.size();++i){
//		double dd = y[i] - y[i-1];
//		if(dd>0){
//			dy_sign.push_back(1);
//		}
//		else{
//			if(dd<0){dy_sign.push_back(-1);}
//			else{dy_sign.push_back(0);}
//		}
//	}
//	//local limit search
//	bool is_local_max = false; 
//	int local_limit_idx = -1;
//	zero_cross_detect(dy_sign,is_local_max,local_limit_idx);
//	if(local_limit_idx<0){
//	    throw_msg(m_patch_id+" not found local limit");
//	}
//	// optimal angle calculation
//	vector<double> xx;
//	vector<double> yy;
//	double oa=0.0;
//	if(is_local_max){
//		//存在局部最大
//        if( local_limit_idx == 0 || local_limit_idx==y.size()-1){
//            //#单调  
//			oa = x[local_limit_idx];
//			//if (local_limit_idx == 0){				
//			//	for(size_t i=0;i<3;++i){
//			///		xx.push_back(x[i]);
//			//		yy.push_back(y[i]);
//			//	}
//			//}
//			//else	{
//			//	for(size_t i=y.size()-3;i<y.size();++i){
//			//		xx.push_back(x[i]);
//			//		yy.push_back(y[i]);
//			//	}
//			//}             
//			//oa = qua_fitting(xx,yy); 
//			if(m_pLogger){				
//				m_pLogger->INFO(m_patch_id+ " angle fit -- signle --- max_angle");
//			}
//			
//		}
//        else{
//            // 极大值        
//			for(size_t i=local_limit_idx-1;i<local_limit_idx+2;++i){
//					xx.push_back(x[i]);
//					yy.push_back(-y[i]);
//				}
//			oa = limit_min_angle(xx,yy,1);  
//			if(m_pLogger){				
//				m_pLogger->INFO(m_patch_id+" angle fit -- local max --- 3p insert");
//			}
//		}
//	}
//    else{
//		// 极小值   
//        oa = limit_min_angle(x,y,local_limit_idx);
//		if(is_base){
//			//for base, offset 45 degree
//			if(oa >45){oa = oa-45.0;}
//			else{oa = oa+45.0;}  
//		}
//		else{  
//			// for rootstock, offset 90 degree
//			if(oa >90){oa = oa-90.0;}
//			else{oa = oa+90.0;}   
//		}
//		if(m_pLogger){				
//			m_pLogger->INFO(m_patch_id+" angle fit -- local min --- 3p insert");
//		}
//	}
//    return oa;
//}

CCotyledonAngle::CCotyledonAngle(ConfigParam&cp, CGcvLogger* pLog)
	:m_cparam(cp),
	m_pLogger(pLog),	
	m_imgId(""),
	m_ppImgSaver(0),
	m_height(0),
	m_width(0)
{
}

CCotyledonAngle::~CCotyledonAngle(void)
{	
}

void CCotyledonAngle::clear_imginfo(){
	if (m_pImginfoOpen){
		imginfo_release(&m_pImginfoOpen);
		m_pImginfoOpen=0;
	}
	if (m_pImginfoAngle){
		imginfo_release(&m_pImginfoAngle);
		m_pImginfoAngle=0;
	}
}
int CCotyledonAngle::angle_recognize(
	ImgInfo* imginfo,
	PositionInfo& posinfo	
	)
{
	//clear opencv windows
	if(m_cparam.image_show){destroyAllWindows();}

	clock_t t;
	clock_t t0 = clock();	
	m_imgId = getImgId(img_type::oa);	

	if(m_pLogger){
		m_pLogger->INFO(m_imgId + " cotyledon angle recognize begin.");
	}	
	if(!isvalid(imginfo)){
		throw_msg(m_imgId+" invalid input image");
	}
	if(m_width==0 || m_height==0){
		m_width = imginfo->width;
		m_height = imginfo->height;
	}	
	if(m_pLogger){
		stringstream buff;
		buff<<m_imgId<<" load image\twidth="<<imginfo->width
			<<"\theight="<<imginfo->height;
		m_pLogger->INFO(buff.str());
	}
	Mat img = imginfo2mat(imginfo);



	//for DEBUG
	//Mat img = *(Mat*)(imginfo);
	//m_width = img.cols/2;
	//m_height = img.rows/2;
	//cv::resize(img, img, cv::Size(m_width,m_height ));
	//end DEBUG


	//imshow_wait("src", img);
	if(m_cparam.self_camera){
		image_set_bottom(img,20,8);
		if(m_pLogger){				
			m_pLogger->DEBUG(m_imgId+" oa_recognize image set bottom with pixel value 20.");				
		}
	}
	if(m_cparam.oa_y_flip){
		flip(img,img,0);
		if(m_pLogger){				
			m_pLogger->DEBUG(m_imgId+" oa_recognize image y fliped.");				
		}
	}
	if(m_ppImgSaver && *m_ppImgSaver){
		(*m_ppImgSaver)->saveImage(img, m_imgId);
		if(m_pLogger){				
			m_pLogger->DEBUG(m_imgId+" oa_recognize after image save.");				
		}
	}
	if(m_cparam.image_show){
		Mat tmp_b = img.clone();			
		imshow_wait("oa_image", tmp_b);
	}

	//1图像分割
	this->img_segment(img);
	if(m_cparam.image_show){					
		imshow_wait("oa_image_bin", m_binImg);
	}
	Mat img_median, img_base, img_open, img_xor;
	cv::medianBlur(m_binImg,img_median, 5);	

	Mat kernel = getStructuringElement(
			MORPH_ELLIPSE,
			Size( 2*m_cparam.oa_morph_radius + 1, 2*m_cparam.oa_morph_radius+1),
			Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
			);  
	morphologyEx(
		img_median,
		img_base,
		MORPH_CLOSE,
		kernel,
		Point(-1,-1),
		1);

	if(m_cparam.image_show){					
		imshow_wait("img_base", img_base);
	}

	//2 找到2片叶子，计算各自重心，主方向，面积，横纵比，椭圆度，选择最优，计算需要旋转角度
	int iter_step = 3;
	int iter_open = 10;
	
	double leaf_area_th = (double)m_cparam.oa_min_leaf_area;
	vector<vector<Point>> contours;
	vector<Vec4i> hierarchy;
	vector<t_leaf> leafs;

	while( iter_open >0 && iter_open < 50){

		morphologyEx(
			img_base,
			img_open,
			MORPH_OPEN,
			kernel,
			Point(-1,-1),
			iter_open);
		if(m_cparam.image_show){
			imshow_wait("img_base", img_base);
			imshow_wait("img_open", img_open);
		}
		m_openImg = img_open.clone();
	
		contours.clear();
		hierarchy.clear();
		leafs.clear();
		findContours(img_open,contours,hierarchy,RETR_EXTERNAL,CHAIN_APPROX_NONE);
		for(int i=0;i<contours.size();++i){
			double area = contourArea(contours[i]);
			if (area < leaf_area_th){continue;}			
			RotatedRect minrect_ellips = fitEllipse(contours[i]);
			Moments moment = moments(contours[i]);
			double centx = moment.m10 / moment.m00;
			double centy = moment.m01 / moment.m00;

			t_leaf leaf;
			leaf.area = area;
			leaf.centx = centx;
			leaf.centy = centy;			
			leaf.minrect_ellips = minrect_ellips;
			leaf.contours = contours[i];

			leafs.push_back(leaf);
		}	

		if(leafs.size()<2){
			iter_open += iter_step;
		}
		if(leafs.size()>2){
			iter_open -= iter_step;
		}
		if(leafs.size()==2){
			break;
		}

	}
	if(leafs.size()!=2){
		throw_msg(m_imgId+" leafs detect failed");
	}

	//select one leaf for angle detection
	double min_err = 1.0e16;
	int min_err_idx = -1;
	for(size_t i =0; i<leafs.size();++i){		
		double err = fit_error_ellipse(leafs[i].minrect_ellips, leafs[i].contours);
		if(err<min_err){
			min_err = err;
			min_err_idx = i;
		}
	}
	if(min_err_idx<0){
		throw_msg(m_imgId+" no normal leaf");
	}

	//计算旋转角度
	double all_area = leafs[0].area+leafs[1].area;
    double r0 = leafs[0].area/all_area;
    double all_centx = leafs[0].centx * r0 + leafs[1].centx * (1-r0);
    double all_centy = leafs[0].centy * r0 + leafs[1].centy * (1 - r0);

	float tar_angle1 = 90 + leafs[min_err_idx].minrect_ellips.angle;//长轴与x轴的夹角（图像坐标系下），90-270间
	float tar_angle2 =  leafs[min_err_idx].minrect_ellips.angle - 90;//-90 ~ 90
	if(tar_angle2<0){
		tar_angle2+=360;
	}
	//叶尖方向（相对于两个叶片中心）
	float leaf_apex_angle = atan2(leafs[min_err_idx].centy - all_centy, leafs[min_err_idx].centx - all_centx) * 180 / CV_PI;
	if(leaf_apex_angle<0){
		leaf_apex_angle += 360;
	}

	float diff1 = fabs(leaf_apex_angle - tar_angle1);
	if(diff1>180){
		diff1 = 360 - diff1;
	}
	float diff2 = fabs(leaf_apex_angle - tar_angle2);
	if(diff2>180){
		diff2 = 360 - diff2;
	}
	float rot_angle = tar_angle1;
	if(diff2 < diff1){
		rot_angle = tar_angle2;
	}

	rot_angle -= 180;//需要旋转的角度

	//draw image
	Mat ellipse_img = m_openImg.clone();
	for(size_t j=0; j<leafs.size(); ++j){
		stringstream buff;
		buff<<"wa="<<(float)((int)(leafs[j].minrect_ellips.angle*10.0)/10.0);	
		Scalar color = Scalar(100);
		if(j == min_err_idx){
			color[0] = 180;
			buff<<",rot="<<(float)((int)(rot_angle*10.0)/10.0);
		}
			
		Point2f vertices[4];
		leafs[j].minrect_ellips.points(vertices);
		for (int i = 0; i < 4; i++)
			line(ellipse_img, vertices[i], vertices[(i+1)%4], color);
					
		putText(ellipse_img,buff.str(),leafs[j].minrect_ellips.center, 0,0.5,color);
			
	}
	if(m_cparam.image_show){		
		imshow_wait("ellipse_img", ellipse_img);			
	}

	// 3 返回结果	
	posinfo.rs_oa=(double)rot_angle;
    
	if(m_pLogger){
		stringstream buff;
		buff<<m_imgId<<" detect rotat angle "<<posinfo.rs_oa;
		m_pLogger->INFO(buff.str());
	}
	//4 返回图像
	if(m_cparam.image_return){
		this->clear_imginfo();
		m_pImginfoOpen = mat2imginfo(m_openImg);
		m_pImginfoAngle = mat2imginfo(ellipse_img);
		posinfo.pp_images[0]=m_pImginfoOpen;
		posinfo.pp_images[1]=m_pImginfoAngle;

		if(m_ppImgSaver && *m_ppImgSaver){
			(*m_ppImgSaver)->saveImage(m_binImg, m_imgId+"_rst_0");	
			(*m_ppImgSaver)->saveImage(ellipse_img, m_imgId+"_rst_1");
			if(!m_grayImg.empty()){
				(*m_ppImgSaver)->saveImage(m_grayImg, m_imgId+"_rst_2");
			}
		}
	}
	
	if(m_pLogger){
		m_pLogger->INFO(m_imgId + " oa_recognize image finished.");
	}	
	return 0;
}
double CCotyledonAngle::fit_error_ellipse(
	RotatedRect& ellipse_rect,
	vector<Point>& pts)
{
	/*
	  ref:https://blog.csdn.net/fangyan90617/article/details/89486331
	  ellipse common function: A * x^2 + B * x*y + C * y^2 + D * x + E * y + F = 0
	  reformed function: 
	              A * x'^2 + B * x'*y' + C * y'^2 + F = 0
				  x' = x - x0
				  y' = y - y0
	  from ellipse_rect, calculate A,B,C,F

	*/
	if(pts.size()==0){
		return 1.0e6;
	}
	float A, B,C, F;
	get_ellipse_param(ellipse_rect, A, B,C, F);
	double mean_err = 0.0;
	double x,y,err;
	for(size_t i=0; i<pts.size(); ++i){
		x = (double)pts[i].x - ellipse_rect.center.x;
		y = (double)pts[i].y - ellipse_rect.center.y;
		err = A * x * x + B * x * y + C * y * y + F;
		mean_err += err*err;
	}
	mean_err /= (double)pts.size();
	mean_err = sqrt(mean_err);
	return mean_err;	
}
void CCotyledonAngle::get_ellipse_param(
	RotatedRect&ellipse_rect,
	float&A, 
	float&B, 
	float&C, 
	float&F)
{
	float theta = ellipse_rect.angle * CV_PI / 180.0;
	float a = ellipse_rect.size.width / 2.0;
	float b = ellipse_rect.size.height / 2.0;
	A = a * a * sin(theta) * sin(theta) + b * b * cos(theta) * cos(theta);
	B = (-2.0) * (a * a - b * b) * sin(theta) * cos(theta);
	C = a * a * cos(theta) * cos(theta) + b * b * sin(theta) * sin(theta);
    F = (-1.0) * a * a * b * b;
}
void CCotyledonAngle::img_segment(Mat&img)
{
	if(img.channels()!=1){
		//color image ,bgr, for testing
		/*Mat img_gray,b_img;
		cvtColor(img,img_gray,COLOR_BGR2GRAY);
		double th = threshold(
			img_gray,
			b_img,
			255,
			255,
			THRESH_OTSU
			);
	
		b_img/= 255;*/

		Mat img_hsv, b_img_hsv;		
		cvtColor(
			img, 
			img_hsv,
			COLOR_BGR2HSV
			);

		inRange(
			img_hsv,
			Scalar(30 , 30, 30),
			Scalar(75, 255, 255),
			b_img_hsv
			);

		/*imshow_wait("oa_image_range", b_img_hsv);
		b_img_hsv/=255;
		b_img = b_img.mul(b_img_hsv);	*/

		//int morph_size = 1;
		Mat kernel = getStructuringElement(
			MORPH_RECT,
			Size( 2*m_cparam.oa_morph_radius + 1, 2*m_cparam.oa_morph_radius+1),
			Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
			);  
		morphologyEx(
			b_img_hsv,
			m_binImg,
			MORPH_CLOSE,
			kernel,
			Point(-1,-1),
			m_cparam.oa_morph_iteration);
		
	}
	else{
		//from imginfo image, gray image 
		//int morph_size = 1;	

		Mat b_img;
		m_grayImg = img.clone();
		double th = threshold(img, b_img, 255, 255,THRESH_OTSU);
		Mat kernel = getStructuringElement( 
			MORPH_RECT, 
			Size( 2*m_cparam.oa_morph_radius + 1, 2*m_cparam.oa_morph_radius+1 ), 
			Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
			);    
		morphologyEx(
			b_img, 
			m_binImg, 
			MORPH_OPEN,
			kernel,
			Point(-1,-1),
			m_cparam.oa_morph_iteration
			);
	}	
}


};