|
|
@@ -0,0 +1,1910 @@
|
|
|
+#include "utils.h"
|
|
|
+#include <opencv.hpp>
|
|
|
+#include <time.h>
|
|
|
+#include <algorithm>
|
|
|
+#include <math.h>
|
|
|
+
|
|
|
+namespace graft_cv{
|
|
|
+
|
|
|
+ cv::Mat imginfo2mat(ImgInfo* imginfo)
|
|
|
+ {
|
|
|
+ assert(imginfo->data);
|
|
|
+ assert(imginfo->height>0 && imginfo->width>0);
|
|
|
+ assert(imginfo->channel == 1 || imginfo->channel == 3);
|
|
|
+ cv::Mat m;
|
|
|
+ if (imginfo->channel == 1) {
|
|
|
+ m = cv::Mat(imginfo->height, imginfo->width, CV_8UC1);
|
|
|
+ }
|
|
|
+ if (imginfo->channel == 3) {
|
|
|
+ m = cv::Mat(imginfo->height, imginfo->width, CV_8UC3);
|
|
|
+ }
|
|
|
+
|
|
|
+ for(int h=0; h<imginfo->height; ++h)
|
|
|
+ {
|
|
|
+ memcpy((void*)(m.ptr(h)),
|
|
|
+ (void*)(imginfo->data+h*imginfo->width * imginfo->channel),
|
|
|
+ imginfo->width * imginfo->channel);
|
|
|
+ };
|
|
|
+ return m;
|
|
|
+ };
|
|
|
+
|
|
|
+ ImgInfo* mat2imginfo(const cv::Mat&img){
|
|
|
+ assert(img.channels()==1 || img.channels() == 3);
|
|
|
+ ImgInfo* imginfo = new ImgInfo();
|
|
|
+ imginfo->channel= img.channels();
|
|
|
+ imginfo->height=img.rows;
|
|
|
+ imginfo->width = img.cols;
|
|
|
+ imginfo->data = 0;
|
|
|
+ imginfo->data = new byte[imginfo->height * imginfo->width * imginfo->channel];
|
|
|
+ memset(imginfo->data, 0,imginfo->height * imginfo->width * imginfo->channel);
|
|
|
+ //IplImage ipl_img = cvIplImage(img);
|
|
|
+ for(int i =0; i<imginfo->height;++i){
|
|
|
+ memcpy(imginfo->data+i*imginfo->width * imginfo->channel, img.ptr(i), imginfo->width * imginfo->channel);
|
|
|
+ }
|
|
|
+
|
|
|
+ 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(
|
|
|
+ cv::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(
|
|
|
+ cv::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(
|
|
|
+ cv::Mat& img,
|
|
|
+ int x0,int y0,
|
|
|
+ int x1,int y1)
|
|
|
+ {
|
|
|
+ cv::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);
|
|
|
+ }
|
|
|
+ cv::line(img,p0,p1, cv::Scalar(255,255,255));
|
|
|
+
|
|
|
+ }
|
|
|
+ string currTime(){
|
|
|
+ char tmp[64];
|
|
|
+ struct tm ptime;
|
|
|
+ time_t time_seconds = time(0);
|
|
|
+ localtime_s(&ptime, &time_seconds);
|
|
|
+ strftime(
|
|
|
+ tmp,
|
|
|
+ sizeof(tmp),
|
|
|
+ "%Y-%m-%d %H:%M:%S",
|
|
|
+ &ptime
|
|
|
+ );
|
|
|
+ return tmp;
|
|
|
+ }
|
|
|
+ // function getImgId
|
|
|
+ // input im_type, img_type
|
|
|
+ static unsigned int ID_COUNTER = 0;
|
|
|
+ string getImgId(int im_type)
|
|
|
+ {
|
|
|
+ char tmp[64];
|
|
|
+ struct tm ptime;
|
|
|
+ time_t time_seconds = time(0);
|
|
|
+ localtime_s(&ptime, &time_seconds);
|
|
|
+ strftime(
|
|
|
+ tmp,
|
|
|
+ sizeof(tmp),
|
|
|
+ "%Y%m%d%H%M%S",
|
|
|
+ &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(
|
|
|
+ cv::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){
|
|
|
+ cv::line(img, cv::Point(0,row), cv::Point(img.cols-1,row), cv::Scalar(100),line_thick);
|
|
|
+ }
|
|
|
+ //veritcal grid
|
|
|
+ for(int col=0; col<img.cols; col+=grid_col){
|
|
|
+ cv::line(img, cv::Point(col,0), cv::Point(col,img.rows-1), cv::Scalar(100),line_thick);
|
|
|
+ }
|
|
|
+
|
|
|
+ }
|
|
|
+ void hist2image(
|
|
|
+ vector<int>& hist,
|
|
|
+ cv::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 = cv::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 = cv::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(
|
|
|
+ cv::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(
|
|
|
+ cv::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(
|
|
|
+ cv::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);
|
|
|
+ }
|
|
|
+ cv::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 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
|
|
|
+ )
|
|
|
+ {
|
|
|
+ //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;
|
|
|
+ cv::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){
|
|
|
+ cv::Mat tmp = roi_img.clone();
|
|
|
+ cv::Point p0((int)(beta0),0);
|
|
|
+ cv::Point p1((int)((double)tmp.rows *beta1 +beta0),tmp.rows);
|
|
|
+ cv::line(tmp,p0,p1, cv::Scalar(128,0,0));
|
|
|
+
|
|
|
+ cv::line(tmp, cv::Point(cent_xs[0],cent_ys[0]),
|
|
|
+ cv::Point(cent_xs.back(),cent_ys.back()), cv::Scalar(128,0,0));
|
|
|
+
|
|
|
+
|
|
|
+ imshow_wait("rs_bin_roi", tmp);
|
|
|
+ }
|
|
|
+ ///////////////////////////////////////////////////////////
|
|
|
+ //兴趣区域图片findContours
|
|
|
+ vector<vector<cv::Point>> contours;
|
|
|
+ vector<cv::Vec4i> hierarchy;
|
|
|
+ findContours(roi_img,contours,hierarchy, cv::RETR_EXTERNAL, cv::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 = cv::pointPolygonTest( contours[max_length_idx], cv::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){
|
|
|
+ cv::Mat edge;
|
|
|
+ edge = cv::Mat::zeros(roi_img.rows, roi_img.cols, CV_8UC1);
|
|
|
+ cv::drawContours(edge, contours, -1, cv::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]);
|
|
|
+ cv::line(edge, cv::Point(c-rad,r), cv::Point(c+rad,r), cv::Scalar(80,0,0));
|
|
|
+ }
|
|
|
+ cv::Point fork_cent_origin;
|
|
|
+ fork_cent_origin.x = xs[max_idx];
|
|
|
+ fork_cent_origin.y =max_idx;
|
|
|
+ cv::circle(edge,fork_cent_origin,3, cv::Scalar(200,0,0));
|
|
|
+ cv::circle(edge,fork_cent_origin,(int)(max_radius), cv::Scalar(200,0,0));
|
|
|
+
|
|
|
+ //circle(edge,Point(xs[opt_max_idx],opt_max_idx),5,Scalar(200,0,0));
|
|
|
+
|
|
|
+ cv::circle(edge,fork_cent,3, cv::Scalar(128,0,0));
|
|
|
+ cv::circle(edge,fork_cent,(int)(max_radius), cv::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,
|
|
|
+ cv::Mat&img,
|
|
|
+ int waittype/*=-1*/
|
|
|
+ )
|
|
|
+ {
|
|
|
+ cv::namedWindow(winname, cv::WINDOW_NORMAL);
|
|
|
+ cv::imshow(winname, img);
|
|
|
+ cv::waitKey(waittype);
|
|
|
+ };
|
|
|
+
|
|
|
+ int get_stem_fork_right(
|
|
|
+ cv::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(
|
|
|
+ cv::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(
|
|
|
+ cv::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(
|
|
|
+ cv::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
|
|
|
+ cv::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];
|
|
|
+ }
|
|
|
+
|
|
|
+ cv::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];
|
|
|
+ }
|
|
|
+ cv::Mat c, d;
|
|
|
+ c = A.t() * A;
|
|
|
+ d = A.t() * b;
|
|
|
+ cv::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;
|
|
|
+ }
|
|
|
+
|
|
|
+ int nms_bbox(
|
|
|
+ vector<Bbox>& candidate_droplets,
|
|
|
+ float th,
|
|
|
+ vector<int>& keep
|
|
|
+ )
|
|
|
+ {
|
|
|
+ keep.clear();
|
|
|
+ if (candidate_droplets.size() == 1) {
|
|
|
+ keep.push_back(0);
|
|
|
+ return 0;
|
|
|
+ }
|
|
|
+
|
|
|
+ sort(candidate_droplets.begin(), candidate_droplets.end(),
|
|
|
+ [=](const Bbox& left, const Bbox& right) {
|
|
|
+ return left.score > right.score;
|
|
|
+ });
|
|
|
+
|
|
|
+ vector<int> order(candidate_droplets.size());
|
|
|
+ for (size_t i = 0; i < candidate_droplets.size(); ++i) {
|
|
|
+ order[i] = (int)i;
|
|
|
+ }
|
|
|
+
|
|
|
+ while (order.size()) {
|
|
|
+ int i = order[0];
|
|
|
+ keep.push_back(i);
|
|
|
+ vector<int> del_idx;
|
|
|
+ for (size_t j = 1; j < order.size(); ++j) {
|
|
|
+ float iou = iou_bbox(
|
|
|
+ candidate_droplets[i],
|
|
|
+ candidate_droplets[order[j]]);
|
|
|
+ if (iou > th) {
|
|
|
+ del_idx.push_back((int)j);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ vector<int> order_update;
|
|
|
+ for (size_t j = 1; j < order.size(); ++j) {
|
|
|
+ vector<int>::iterator it = find(del_idx.begin(), del_idx.end(), j);
|
|
|
+ if (it == del_idx.end()) {
|
|
|
+ order_update.push_back(order[j]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ order.clear();
|
|
|
+ order.assign(order_update.begin(), order_update.end());
|
|
|
+ }
|
|
|
+ return 0;
|
|
|
+ }
|
|
|
+
|
|
|
+ float iou_bbox(
|
|
|
+ const Bbox & det_a,
|
|
|
+ const Bbox & det_b)
|
|
|
+ {
|
|
|
+ float aa = (float)(det_a.x2 - det_a.x1 + 1) * (det_a.y2 - det_a.y1 + 1);
|
|
|
+ float ab = (float)(det_b.x2 - det_b.x1 + 1) * (det_b.y2 - det_b.y1 + 1);
|
|
|
+
|
|
|
+ float xx1 = (float)max(det_a.x1, det_b.x1);
|
|
|
+ float yy1 = (float)max(det_a.y1, det_b.y1);
|
|
|
+ float xx2 = (float)min(det_a.x2, det_b.x2);
|
|
|
+ float yy2 = (float)min(det_a.y2, det_b.y2);
|
|
|
+
|
|
|
+ float w = (float)max(0.0, xx2 - xx1 + 1.0);
|
|
|
+ float h = (float)max(0.0, yy2 - yy1 + 1.0);
|
|
|
+ float inter = w * h;
|
|
|
+ float ovr = inter / (aa + ab - inter);
|
|
|
+ return ovr;
|
|
|
+ }
|
|
|
+
|
|
|
+};
|
chenhongjiang
1 年間 前