tea_cv/utils.h

#pragma once

#include <opencv2\imgproc\imgproc.hpp>
//#include <boost\math\distributions\binomial.hpp>
#include <vector>
#include <time.h>
#include <string>
#include <numeric>

#include "data_def_api.h"
#include "data_def.h"

//using namespace cv;
using namespace std;

namespace graft_cv{

	// convert ImgInfo to cv mat 
	cv::Mat imginfo2mat(ImgInfo*);
	ImgInfo* mat2imginfo(const cv::Mat&);
	void imginfo_release(ImgInfo**);
	bool isvalid(const ImgInfo*);
	string currTime();
	string getImgId(int im_type);
	string getImgIdOa(string iid, int im_idx);
	inline void throw_msg(string& msg);

	// set bottom image rows_cnt rows to "value", for xiaoxu camera
	void image_set_bottom(cv::Mat& img,unsigned char value,int rows_cnt);
	void image_set_top(cv::Mat& img,	unsigned char value,	int rows_cnt);
	void image_draw_line(cv::Mat& img,int x0,int y0,int x1,int y1);
	// histogram
	void mat_histogram(
		cv::Mat&,
		vector<int>&hist,
		int axis=0, 
		int r0=-1, 
		int r1=-1, 
		int c0=-1, 
		int c1=-1
		);
	// weighted histogram
	void mat_histogram_w(
		cv::Mat& img,
		vector<int>&hist,
		int axis=0,
		int r0=-1, 
		int r1=-1, 
		int c0=-1, 
		int c1=-1,
		bool asc_w=true //true---ascending weight with x/y coordinate, [0-1.0], false --- descending weight [1.0--0.0]
	);
	void mat_histogram_yfork(
		cv::Mat& img,
		vector<int>&hist,	
		int r0, 
		int r1	
	);

	// histogram to image	
	void hist2image(
		vector<int>& hist,
		cv::Mat&img,
		int aix=1,
		int grid_col=50,
		int grid_row=50
		);
	
	template<class T>
	void hist2image_line(
	vector<T>& hist, 
		cv::Mat&img,
	int grid_col/*=50*/,
	int grid_row/*=50*/
	)
{
	if(hist.size()<2){return;}
	vector<T>::iterator maxit = max_element(begin(hist),end(hist)); 
	if( grid_col<10){grid_col=10;}
	if( grid_row<10){grid_row=10;}	
	
	img = Mat::zeros(*maxit,hist.size(), CV_8UC1);
	drawgrid(img,grid_col, grid_row,1);
	for(size_t i=1; i<hist.size();++i){
		line(img,Point(i-1,(int)*maxit - (int)hist[i-1]),Point(i,(int)*maxit - (int)hist[i]),Scalar(255));
	}
	
};

	template<class T>
	void hist2image_scale_line(
	vector<T>& hist, 
		cv::Mat&img,
	T ymin,
	T ymax,
	T yratio,//数值像素比，一个像素代表的数值量
	T th,
	int grid_col/*=50*/,
	int grid_row/*=50*/
	)
{
	if(hist.size()<2){return;}	
	if( grid_col<10){grid_col=10;}
	if( grid_row<10){grid_row=10;}	
	int height =(int)((ymax-ymin)/yratio);
	
	img = cv::Mat::zeros(height,hist.size(), CV_8UC1);
	drawgrid(img,grid_col, grid_row,1);
	for(size_t i=1; i<hist.size();++i){
		int y_1 = height-(int)((hist[i-1]-ymin)/yratio);
		int y_0 = height-(int)((hist[i]-ymin)/yratio);
		cv::line(img, cv::Point(i-1,y_1), cv::Point(i,y_0), cv::Scalar(255));
	}
	int h1 = height-(int)((th-ymin)/yratio);
	int h2 = height-(int)((1.0/th-ymin)/yratio);
	cv::line(img, cv::Point(0,h1), cv::Point(hist.size()-1,h1), cv::Scalar(200));
	cv::line(img, cv::Point(0,h2), cv::Point(hist.size()-1,h2), cv::Scalar(200));
	
};

	void drawgrid(		
		cv::Mat&img,
		int grid_col=50,
		int grid_row=50,
		int line_thick=1
		);

	//r2 index calculation for line regression
	//         y = beta0 + beta1 * x
	template<class T>
	double r2index(
		const vector<T>& x, 
		const vector<T>& y, 
		size_t i0, 
		size_t i1,
		double& beta0,
		double& beta1)
	{
		double r2 = 0.0;
		assert(x.size() == y.size());
		assert(i0<i1);
		assert(i1-i0>1);
		assert(i0>=0 && i1<y.size());
    
		// y = beta0 + beta1 * x
		//double beta0=0, beta1=0;
		double t1=0, t2=0,t3=0,t4=0;
		double n = (double)(i1-i0+1);
		for(size_t i=i0; i<=i1; ++i)
		{
			t1 += (double)x[i] * (double)x[i];
			t2 += (double)x[i];
			t3 += (double)x[i] * (double)y[i];
			t4 += (double)y[i];
		}
		beta0 = (t1*t4 - t2*t3) / (t1*n - t2*t2);
		beta1 = (t3*n - t2*t4) / (t1*n -t2*t2);
    
		double sst=0, ssr=0, sse=0;
		double y_mu = t4/n;
		double y_hat;
		for(size_t i=i0; i<=i1; ++i)
		{
			y_hat = beta0 + beta1*(double)x[i];
			ssr += (y_hat - y_mu)*(y_hat - y_mu);
			sse += ((double)y[i] - y_hat)* ((double)y[i] - y_hat);
		}
		sst = ssr + sse;
		r2 = 1.0 - sse/sst;   
    
		return r2;
	};

	void get_stem_x_range_oa(
		const 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
		);
	void get_stem_x_range_rscut(
		const 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
		);
	
	void get_stem_x_range_scion(
		const vector<int>& hist_h, 
		double h_ratio, 
		int padding, 
		int&x0, 
		int&x1, 
		int&stem_x0,
		int&stem_x1
		);

	////////////////////////////////////////////////////////////////////////
	// y-fork detect
	void get_stem_y_fork(
		const vector<int>& hist, 
		double ratio,
		int stem_dia_min, 
		int stem_fork_y_min,
		double stem_dia_mp,
		int& fork_y,     //output, 茎分叉y像素位置
		int& end_y,      //output，茎根部y像素位置
		int& stem_dia    //output， end_y和fork_y间茎粗
	);

	void get_stem_y_fork_rs(
	const 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
	);
	void get_stem_y_fork_rs_update(
	const cv::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,
		cv::Point& fork_cent,
	double& max_radius,
	double& stem_angle
	);
	double yfork_validity_position(
		int max_y, 
		int end_y, 
		int fork_y);
	double yfork_validity_stemdiaratio(
		int stem_dia, 
		int dia_y,
		double opt_dia_ratio);
	void yfork_validity_stemdiachange(
	    const vector<int>& hist,
		int stem_dia_min,
		vector<double>& var_ratio		
		);

	//通过计算茎均值，运用统计过程控制方法检测显著点
	void get_stem_y_fork_3sigma(
		const vector<int>& hist, 
		double ratio,
		int stem_dia_min, 
		int stem_fork_y_min,
		double stem_dia_mp,
		int& fork_y,     //output, 茎分叉y像素位置
		int& end_y,      //output，茎根部y像素位置
		int& stem_dia    //output， end_y和fork_y间茎粗
	);

	double stem_y_fork_validity(
		const vector<int>& hist,
		int fork_y,
		int end_y);

	void get_hist_segment(
		const vector<int>& hist, 
		int threshold,
		vector<gcv_point<int>>& segments
	);
	void get_hist_mean_var_local(
		const vector<int>& hist, 
		int x0,
		int x1,
		double& mean,
		double& var
	);
	//void get_stem_y_index(
	//	Mat& bin_img,
	//	int x0,
	//	int x1,
	//	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
	//);
	void imshow_wait(
		const char* winname, 
		cv::Mat&,
		int waittype=-1
		);

	int get_stem_fork_right(
		cv::Mat&stem_img,
		int stem_fork_y, 
		int stem_x0, 
		int stem_x1,
		int detect_window);
	int get_stem_fork_left(
		cv::Mat&stem_img,
		int stem_fork_y, 
		int stem_x0, 
		int stem_x1,
		int detect_window);
	void get_stem_fork_xs(
		cv::Mat&stem_img,
		int stem_fork_y, 
		int stem_x0, 
		int stem_x1,
		int x0, 
		int x1,		
		int & fork_x0,
		int & fork_x1
	);

	void get_next_pt(
		cv::Mat& edge_img,
	    gcv_point<int>&start_pt, 
	    gcv_point<int>&pre_pt, 
		gcv_point<int>&nxt_pt, //output
	    bool is_up=true);

	// 二次曲线拟合，返回二次曲线对称轴x的坐标
	double qua_fitting(vector<double>&xx, vector<double>&yy);

	template <typename T>
    vector<size_t> sort_indexes_e(vector<T> &v, bool ascending=true)
    {
        vector<size_t> idx(v.size());
        iota(idx.begin(), idx.end(), 0);
		if(ascending){
        sort(idx.begin(), idx.end(),
            [&v](size_t i1, size_t i2) {return v[i1] < v[i2]; });
		}
		else{
			sort(idx.begin(), idx.end(),
            [&v](size_t i1, size_t i2) {return v[i1] > v[i2]; });		
		}
        return idx;
    }
	template<typename T>
	void hist_filter(vector<T>& hist, int method=0, int radius=2)
	{
		//mothod: 0---mid; 1--mean
		assert(radius>=1);
		if(hist.size()<3){return;}
		vector<T>tmp;
		
		if(method==0){
			T* buff = new T[2+2*radius];
			for(int i=0; i<hist.size();++i){
				int cnt=0;
				for(int r = -radius;r<=radius;++r){
					int idx=i+r;
					if(idx<0 || idx>=hist.size()){continue;}
					buff[cnt++]=hist[idx];
				}
				sort(buff,buff+cnt);
				int cent = cnt/2;
				tmp.push_back(buff[cent]);
			}
			delete []buff;
		}
		if(method==1){
			for(int i=0; i<hist.size();++i){
				int cnt=0;
				T mean = T(0);
				for(int r = -radius;r<=radius;++r){
					int idx=i+r;
					if(idx<0 || idx>=hist.size()){continue;}
					mean+=hist[idx];
				}
				mean =mean/cnt;
				tmp.push_back(mean);
			}
		}
		
		for(size_t j=0;j<hist.size();++j){hist[j]=tmp[j];}
	}
	template<typename T>
	double otsu(vector<T>&data)
	{
		if (data.size() == 0) {
			return 0.0;
		}
		std::vector<T>data_sort(data);
		std::sort(data_sort.begin(), data_sort.end());
		double t0 = data_sort[int(data_sort.size() / 2)];
		double t1 = data_sort.back();
		double t = t0;
		double max_g = 0;
		double max_t = t0;
		std::vector<T>d0;
		std::vector<T>d1;
		double n = data.size();
		while (t < t1) {
			d0.clear();
			d1.clear();
			for (auto& v : data_sort) {
				if (v < t) { d0.push_back(v); }
				else { d1.push_back(v); }
			}
			double g = get_hist_mean(d0) - get_hist_mean(d1);
			g *= g;
			g = g* d0.size() * d1.size() / n / n;
			if (g > max_g) {
				max_g = g;
				max_t = t;
			}
			t += 1;
		}
		return max_t;

	}
	template<typename T>
	double get_hist_mean(vector<T>&data) {
		double mu = 0;
		for (auto &v : data) { mu += v; }
		if (data.size() > 0) {
			mu /= data.size();
		}
		return mu;
	}

	template<typename T>
	double get_hist_std(vector<T>&data, T mu) {
		double std_v = 0;
		for (auto &v : data) { std_v += (v-mu)*(v - mu); }
		if (data.size() > 0) {
			std_v /= data.size();
			std_v = sqrt(std_v);
		}
		return std_v;
	}

	template<typename T>
	int trend_detect_pos(const vector<T>&data, int data_th, int step=15) {
		int pos = -1;
		if (data.size() < 2 * step) {
			return pos;
		}		
		int radius = int(step/2);
		boost::math::binomial_distribution<> binomal(radius, 0.5);
		pos = radius;
		int max_v = 0;
		int pos_n, neg_n;
		for (size_t i = radius; i < data.size() - radius; ++i) {
			pos_n = neg_n = 0;
			for (int j = 1; j <= radius; ++j) {
				int diff = data[i + j] - data[i + j -radius];
				if (diff > 0) {
					pos_n++;
				}
				else {
					neg_n++;
				}
			}
			int k = min(pos_n, neg_n);			
			double p = 2.0 * boost::math::cdf(binomal, k);
			std::cout <<"idx="<<i<< ", posn=" << pos_n << ", negn=" << neg_n << ", pvalue=" << p << endl;
			if (p > 0.05) {
				if(fabs(data.at(i) - data_th)<3)
				{
					pos = i;
					break;
				}
			}
		}
		return pos;

	}

	//nms
	int nms_bbox(
		vector<Bbox>& candidate_droplets,	//input
		float th,							//input
		vector<int>& keep					//output
	);
	//iou
	float iou_bbox(
		const Bbox & det_a,
		const Bbox & det_b);
	


	
};