12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205 |
- //#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/*=0*/)
- :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/*=false*/)
- {
- 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/*=0*/)
- :COptimalAngle(cp,pLog)
- {
- }
- COptimalAnglePart::~COptimalAnglePart()
- {
- COptimalAngle::clear_imginfo();
- }
- // append()加载的图片,满足不是0-180度采样,但保证部分采样中一定存在最大值,
- // 直接利用这些数值进行二次函数的拟合
- double COptimalAnglePart::angle_fit(map<int,int>& an2width,bool is_base/*=false*/)
- {
-
- 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/*=false*/)
- {
- // 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);
- 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
- );
- }
- }
- };
|