optimal_angle.cpp 29 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205
  1. //#include "StdAfx.h"
  2. #include <opencv2\imgproc\imgproc.hpp>
  3. #include "optimal_angle.h"
  4. #include "utils.h"
  5. #include <string.h>
  6. namespace graft_cv{
  7. COptimalAngle::COptimalAngle(ConfigParam&cp, CGcvLogger* pLog/*=0*/)
  8. :m_cparam(cp),
  9. m_pLogger(pLog),
  10. m_pImginfo(0),
  11. m_pImginfoBase(0),
  12. m_height(0),
  13. m_width(0),
  14. m_patch_id(""),
  15. m_imgIndex(0),
  16. m_imgId(""),
  17. m_ppImgSaver(0)
  18. {
  19. }
  20. COptimalAngle::~COptimalAngle(void)
  21. {
  22. this->clear_imginfo();
  23. }
  24. void COptimalAngle::clear_imginfo(){
  25. if (m_pImginfo){
  26. imginfo_release(&m_pImginfo);
  27. m_pImginfo=0;
  28. }
  29. if (m_pImginfoBase){
  30. imginfo_release(&m_pImginfoBase);
  31. m_pImginfoBase=0;
  32. }
  33. }
  34. int COptimalAngle::reset()
  35. {
  36. if(m_pLogger){
  37. m_pLogger->DEBUG("optimal angle reset begin.");
  38. }
  39. m_fork_ys.clear();
  40. m_end_ys.clear();
  41. m_widths.clear();
  42. m_an2width.clear();
  43. m_an2widthBase.clear();
  44. m_width=0;
  45. m_height=0;
  46. m_patch_id="";
  47. m_imgIndex=0;
  48. if(m_pLogger){
  49. m_pLogger->DEBUG("optimal angle reset finished.");
  50. }
  51. return 0;
  52. }
  53. int COptimalAngle::append(
  54. ImgInfo* imginfo,
  55. PositionInfo& posinfo
  56. )
  57. {
  58. //clear opencv windows
  59. if(m_cparam.image_show){destroyAllWindows();}
  60. clock_t t;
  61. clock_t t0 = clock();
  62. if(m_imgIndex==0){
  63. //new patch
  64. m_patch_id = getImgId(img_type::oa);
  65. }
  66. m_imgId = getImgIdOa(m_patch_id,m_imgIndex);
  67. m_imgIndex+=1;
  68. if(m_pLogger){
  69. m_pLogger->INFO(m_imgId + " append image begin.");
  70. }
  71. if(!isvalid(imginfo)){
  72. throw_msg(m_imgId+" invalid input image");
  73. }
  74. if(m_width==0 || m_height==0){
  75. m_width = imginfo->width;
  76. m_height = imginfo->height;
  77. }
  78. int angle = imginfo->angle;
  79. if(m_pLogger){
  80. stringstream buff;
  81. buff<<m_imgId<<" append image, angle="<<angle
  82. <<"\twidth="<<imginfo->width
  83. <<"\theight="<<imginfo->height;
  84. m_pLogger->INFO(buff.str());
  85. }
  86. Mat img = imginfo2mat(imginfo);
  87. //imshow_wait("src", img);
  88. if(m_cparam.self_camera){
  89. image_set_bottom(img,20,8);
  90. if(m_pLogger){
  91. m_pLogger->DEBUG(m_imgId+" oa_append image set bottom with pixel value 20.");
  92. }
  93. }
  94. if(m_cparam.oa_y_flip){
  95. flip(img,img,0);
  96. if(m_pLogger){
  97. m_pLogger->DEBUG(m_imgId+" oa_append image y fliped.");
  98. }
  99. }
  100. if(m_ppImgSaver && *m_ppImgSaver){
  101. (*m_ppImgSaver)->saveImage(img, m_imgId);
  102. if(m_pLogger){
  103. m_pLogger->DEBUG(m_imgId+" oa_append after image save.");
  104. }
  105. }
  106. if(m_cparam.image_show){
  107. Mat tmp_b = img.clone();
  108. 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);
  109. 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);
  110. imshow_wait("oa_sub_image_range", tmp_b);
  111. }
  112. Rect rect_rs(Point(0,0),Point(img.cols,m_cparam.oa_clip_y_min));
  113. //Rect rect_base(Point(0,m_cparam.oa_clip_y_max),Point(img.cols,img.rows));
  114. Mat rs_img = img(rect_rs);
  115. //Mat bs_img = img(rect_base);
  116. // 1 rootstock image processing
  117. // 1.1 calculate the image, get binary image and get object width
  118. int min_idx, max_idx;
  119. vector<int> hist_col;
  120. int width = imgproc(rs_img, hist_col, min_idx, max_idx);
  121. m_an2width.insert(make_pair<int,int>(angle, width));
  122. posinfo.rs_oa_width=width;
  123. if(m_cparam.image_show){
  124. Mat hist_col_img;
  125. hist2image(hist_col,hist_col_img,1,m_cparam.image_col_grid,m_cparam.image_row_grid);
  126. imshow_wait("oa_hist_col", hist_col_img);
  127. }
  128. if(m_pLogger){
  129. stringstream buff;
  130. buff<<m_imgId<<" append image, rootstock width="<<width;
  131. m_pLogger->INFO(buff.str());
  132. }
  133. if(!m_cparam.image_show){
  134. t = clock();
  135. if(1000.0*((float)(t-t0))/CLOCKS_PER_SEC>(float)m_cparam.timeout_append){
  136. throw_msg(m_imgId+" time out");
  137. }
  138. }
  139. //1.2 茎粗计算
  140. //1.2.1 茎x方向位置:x方向位置范围
  141. int r0,r1;
  142. r0=r1=-1;
  143. std::vector<int>hist_row_leaf;
  144. mat_histogram(m_binImg,hist_row_leaf,1);
  145. for(size_t i=0;i<hist_row_leaf.size();++i){
  146. if(hist_row_leaf[i]>=m_cparam.oa_min_hist_value && r0<0){r0=i;}
  147. if(hist_row_leaf[i]>=m_cparam.oa_min_hist_value){r1=i;}
  148. }
  149. if(r0<0 || r1<0){throw_msg(string(m_imgId+" invalid image, with no object in binary image"));}
  150. vector<int> hist_col_yfork;
  151. mat_histogram_yfork(m_binImg,hist_col_yfork,r0,r1);
  152. if(m_cparam.image_show){
  153. Mat tmp;
  154. hist2image(hist_col_yfork,tmp,1,m_cparam.image_col_grid,m_cparam.image_row_grid);
  155. imshow_wait("oa_hist_col_yfork", tmp);
  156. }
  157. int x0,x1,stem_x0, stem_x1;
  158. int centx = (int)((min_idx+ max_idx)/2.0);
  159. get_stem_x_range_oa(
  160. hist_col_yfork,
  161. m_cparam.oa_col_th_ratio,
  162. m_cparam.oa_stem_x_padding,
  163. centx,
  164. width,
  165. x0,
  166. x1,
  167. stem_x0,
  168. stem_x1);
  169. if(m_pLogger){
  170. stringstream buff;
  171. buff<<m_imgId<<" append image, x0="<<x0
  172. <<"\tx1="<<x1
  173. <<"\tstem_x0="<<stem_x0
  174. <<"\tstem_x1="<<stem_x1;
  175. m_pLogger->INFO(buff.str());
  176. }
  177. if(m_cparam.image_show){
  178. Mat tmp_b = m_binImg.clone();
  179. cv::line(tmp_b,Point(x0,0), Point(x0,tmp_b.rows),Scalar(200),3);
  180. cv::line(tmp_b,Point(x1,0), Point(x1,tmp_b.rows),Scalar(200),3);
  181. imshow_wait("oa_x_range", tmp_b);
  182. }
  183. if(!m_cparam.image_show){
  184. t = clock();
  185. if(1000.0*((float)(t-t0))/CLOCKS_PER_SEC>(float)m_cparam.timeout_append){
  186. throw_msg(m_imgId+" time out");
  187. }
  188. }
  189. vector<int> hist_row;
  190. mat_histogram(m_binImg,hist_row,1,-1,-1,x0,x1+1);
  191. if(m_cparam.image_show){
  192. Mat hist_row_img;
  193. hist2image(hist_row,hist_row_img, 0,m_cparam.image_col_grid,m_cparam.image_row_grid);
  194. imshow_wait("oa_hist_row", hist_row_img);
  195. }
  196. // 1.2.2 茎分叉点检测,y方向检测
  197. int stem_fork_y=-1,stem_end_y=-1, stem_dia=-1;
  198. get_stem_y_fork(
  199. hist_row,
  200. m_cparam.oa_row_th_ratio,
  201. m_cparam.oa_stem_dia_min,
  202. m_cparam.oa_stem_fork_y_min,
  203. m_cparam.oa_stem_dia_mp,
  204. stem_fork_y,
  205. stem_end_y,
  206. stem_dia);
  207. if(m_pLogger){
  208. stringstream buff;
  209. buff<<m_imgId<<" append image, stem_fork_y="<<stem_fork_y
  210. <<"\tstem_end_y="<<stem_end_y
  211. <<"\tstem_dia="<<stem_dia;
  212. m_pLogger->INFO(buff.str());
  213. }
  214. if(m_cparam.image_show){
  215. Mat tmp_bin = m_binImg.clone();
  216. cv::line(tmp_bin,Point(min_idx,0), Point(min_idx,tmp_bin.rows),Scalar(200),3);
  217. cv::line(tmp_bin,Point(max_idx,0), Point(max_idx,tmp_bin.rows),Scalar(200),3);
  218. cv::line(tmp_bin,Point(x0,stem_fork_y), Point(x1,stem_fork_y),Scalar(200),3);
  219. cv::line(tmp_bin,Point(x0,stem_end_y), Point(x1,stem_end_y),Scalar(200),3);
  220. imshow_wait("oa_result", tmp_bin);
  221. }
  222. m_fork_ys.push_back(stem_fork_y);
  223. m_end_ys.push_back(stem_end_y);
  224. m_widths.push_back(width);
  225. //2 基质图像处理
  226. /*int min_idx_base, max_idx_base;
  227. vector<int> hist_col_base;
  228. int width_base = imgprocBase(bs_img, hist_col_base, min_idx_base, max_idx_base);
  229. m_an2widthBase.insert(make_pair<int,int>(angle, width_base));
  230. posinfo.rs_oa_width_base = width_base;
  231. if(m_pLogger){
  232. stringstream buff;
  233. buff<<m_imgId<<" append image, base width="<<width_base;
  234. m_pLogger->INFO(buff.str());
  235. }*/
  236. // 3 返回结果
  237. double fork_y = (double)stem_fork_y;
  238. double end_y = (double)stem_end_y;
  239. fork_y = ((double)m_height/2.0 - fork_y)*m_cparam.rs_oa_pixel_ratio;
  240. end_y = ((double)m_height/2.0 - end_y)*m_cparam.rs_oa_pixel_ratio;
  241. posinfo.rs_oa_stem_y_fork=fork_y;
  242. posinfo.rs_oa_clamp_y_end=end_y;
  243. if(m_pLogger){
  244. stringstream buff;
  245. buff<<m_imgId<<" append image, stem_fork_y(mm)="<<fork_y
  246. <<"\tstem_end_y(mm)="<<end_y;
  247. m_pLogger->INFO(buff.str());
  248. }
  249. //4 返回图像
  250. if(m_cparam.image_return){
  251. this->clear_imginfo();
  252. cv::line(m_binImg,Point(min_idx,0), Point(min_idx,rs_img.rows),Scalar(200),3);
  253. cv::line(m_binImg,Point(max_idx,0), Point(max_idx,rs_img.rows),Scalar(200),3);
  254. cv::line(m_binImg,Point(x0,stem_fork_y), Point(x1,stem_fork_y),Scalar(200),3);
  255. cv::line(m_binImg,Point(x0,stem_end_y), Point(x1,stem_end_y),Scalar(200),3);
  256. clear_imginfo();
  257. m_pImginfo = mat2imginfo(m_binImg);
  258. //m_pImginfoBase = mat2imginfo(m_binImgBase);
  259. posinfo.pp_images[0]=m_pImginfo;
  260. //posinfo.pp_images[1]=m_pImginfoBase;
  261. if(m_ppImgSaver && *m_ppImgSaver){
  262. (*m_ppImgSaver)->saveImage(m_binImg, m_imgId+"_rst_0");
  263. }
  264. }
  265. if(m_pLogger){
  266. m_pLogger->INFO(m_imgId + " append image finished.");
  267. }
  268. return 0;
  269. }
  270. double COptimalAngle::infer(PositionInfo& posinfo){
  271. if(m_pLogger){
  272. m_pLogger->INFO(m_patch_id + " optimal angle infer begin.");
  273. }
  274. double oa = this->angle_fit(this->m_an2width);
  275. /*double oa_base = 0.0;
  276. try{
  277. oa_base = this->angle_fit_base(this->m_an2widthBase);
  278. }
  279. catch(...){
  280. if(m_pLogger){
  281. m_pLogger->WARNING("angle_fit_base() error");
  282. }
  283. }*/
  284. posinfo.rs_oa = oa;
  285. //posinfo.rs_oa_base=oa_base;
  286. if(m_fork_ys.size()==0 ||m_end_ys.size()==0){
  287. throw_msg(m_patch_id+ " empty fork_ys or end_ys, NEED append image");
  288. }
  289. vector<size_t>des_idx = sort_indexes_e(m_widths,false);
  290. int e_idx = (int)((float)m_end_ys.size() * 0.5);
  291. sort(m_end_ys.begin(), m_end_ys.end());
  292. double fork_y = ((m_fork_ys[des_idx[0]] + m_fork_ys[des_idx[1]])/2.0);
  293. double end_y = m_end_ys[e_idx];
  294. if(m_pLogger){
  295. stringstream buff;
  296. buff<<m_patch_id<<" angle fit result, stem_fork_y(pixel)="<<(int)(fork_y)
  297. <<"\tstem_end_y(pixel)="<<(int)(end_y);
  298. m_pLogger->INFO(buff.str());
  299. }
  300. fork_y = ((double)m_height/2.0 - fork_y)*m_cparam.rs_oa_pixel_ratio;
  301. end_y = ((double)m_height/2.0 - end_y)*m_cparam.rs_oa_pixel_ratio;
  302. posinfo.rs_oa_stem_y_fork=fork_y;
  303. posinfo.rs_oa_clamp_y_end=end_y;
  304. if(m_pLogger){
  305. stringstream buff;
  306. buff<<m_patch_id<<" angle fit result, stem_fork_y(mm)="<<fork_y
  307. <<"\tstem_end_y(mm)="<<end_y;
  308. m_pLogger->INFO(buff.str());
  309. m_pLogger->INFO(m_patch_id + " optimal angle infer finished.");
  310. }
  311. return oa;
  312. }
  313. int COptimalAngle::imgproc(
  314. Mat& img,
  315. vector<int>& hist,
  316. int& min_idx,
  317. int& max_idx
  318. )
  319. {
  320. //int morph_size = 1;
  321. //int min_hist_value = 10;
  322. hist.clear();
  323. Mat b_img;
  324. double th = threshold(
  325. img,
  326. b_img,
  327. 255,
  328. 255,
  329. THRESH_OTSU);//+THRESH_BINARY_INV
  330. Mat kernel = getStructuringElement(
  331. MORPH_RECT,
  332. Size( 2*m_cparam.oa_morph_radius+1, 2*m_cparam.oa_morph_radius+1),
  333. Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
  334. );
  335. morphologyEx(
  336. b_img,
  337. m_binImg,
  338. MORPH_CLOSE,
  339. kernel,
  340. Point(-1,-1),
  341. m_cparam.oa_morph_iteration);
  342. mat_histogram(m_binImg, hist, 0);
  343. /*int min_idx, max_idx;*/
  344. min_idx = hist.size();
  345. max_idx = 0;
  346. vector<int>::const_iterator it0 = hist.begin();
  347. for(vector<int>::const_iterator it = hist.begin();it!=hist.end();++it){
  348. if(*it>=m_cparam.oa_min_hist_value){
  349. int idx = it - it0;
  350. if(idx<min_idx){min_idx = idx;}
  351. if(idx > max_idx) {max_idx = idx;}
  352. }
  353. }
  354. if(max_idx<=min_idx){
  355. throw_msg(m_imgId+" invalid binary image, not exist valid objects");
  356. }
  357. if(m_cparam.image_show){
  358. Mat tmp = m_binImg.clone();
  359. cv::line(tmp,Point(min_idx,0), Point(min_idx,tmp.rows),Scalar(156),3);
  360. cv::line(tmp,Point(max_idx,0), Point(max_idx,tmp.rows),Scalar(156),3);
  361. imshow_wait("oa_bin", tmp);
  362. }
  363. int width = max_idx-min_idx+1;
  364. return width;
  365. }
  366. int COptimalAngle::imgprocBase(
  367. Mat& img,
  368. vector<int>& hist,
  369. int& min_idx,
  370. int& max_idx
  371. )
  372. {
  373. //int morph_size = 1;
  374. //int min_hist_value = 10;
  375. hist.clear();
  376. Mat b_img;
  377. double th = threshold(
  378. img,
  379. b_img,
  380. 255,
  381. 255,
  382. THRESH_OTSU+THRESH_BINARY_INV);
  383. Mat kernel = getStructuringElement(
  384. MORPH_RECT,
  385. Size( 2*m_cparam.oa_morph_radius_base+1, 2*m_cparam.oa_morph_radius_base+1),
  386. Point( m_cparam.oa_morph_radius_base, m_cparam.oa_morph_radius_base)
  387. );
  388. morphologyEx(
  389. b_img,
  390. m_binImgBase,
  391. MORPH_OPEN,
  392. kernel,
  393. Point(-1,-1),
  394. m_cparam.oa_morph_iteration_base);
  395. if(m_cparam.image_show){
  396. destroyAllWindows();
  397. imshow_wait("oa_bin_base", m_binImgBase);
  398. }
  399. mat_histogram(m_binImgBase, hist, 0);
  400. /*int min_idx, max_idx;*/
  401. min_idx = hist.size();
  402. max_idx = 0;
  403. vector<int>::const_iterator it0 = hist.begin();
  404. for(vector<int>::const_iterator it = hist.begin()+5;it!=hist.end();++it){
  405. if(*it>=m_cparam.oa_min_hist_value_base){
  406. int idx = it - it0;
  407. if(idx<min_idx){min_idx = idx;}
  408. if(idx > max_idx) {max_idx = idx;}
  409. }
  410. }
  411. if(max_idx<=min_idx){
  412. throw_msg(m_imgId+" invalid binary image, not exist valid objects");
  413. }
  414. int width = max_idx-min_idx+1;
  415. return width;
  416. }
  417. //angle_fit()方法
  418. //append()加载的图片,满足下述条件:0-180度采样,其中0和180度都有数值
  419. double COptimalAngle::angle_fit(map<int,int>& an2width,bool is_base/*=false*/)
  420. {
  421. if(an2width.size()<3){
  422. throw_msg(m_patch_id+" not enough valid images, NEED append image");
  423. }
  424. double oa = 0.0;
  425. std::vector<double> x;//angle
  426. std::vector<double> y;//width
  427. map<int,int>::iterator it = an2width.begin();
  428. for(it=an2width.begin(); it!=an2width.end(); ++it)
  429. {
  430. x.push_back((double)(it->first));
  431. y.push_back((double)(it->second));
  432. }
  433. //mean value of first and latest y
  434. double y_mu = 0.5*(y[0]+ y[y.size()-1]);
  435. y[0] = y[y.size()-1] = y_mu;
  436. size_t ap_times = x.size()-1;
  437. for (size_t i=0; i<ap_times;++i)
  438. {
  439. x.push_back(x[i+1]+180.0);
  440. y.push_back(y[i+1]);
  441. }
  442. vector<double>::iterator smallest = min_element(begin(y), end(y));
  443. size_t min_idx, max_idx;
  444. min_idx = distance(y.begin(), smallest);
  445. for(size_t i = min_idx+1;i<y.size();++i){
  446. if(fabs(y[i]-*smallest)<0.1){
  447. max_idx = i;
  448. break;
  449. }
  450. }
  451. // slice x and y the convex segment
  452. vector<double>xx;
  453. vector<double>yy;
  454. for(size_t i=min_idx;i<=max_idx;++i){
  455. xx.push_back(x[i]);
  456. yy.push_back(y[i]);
  457. }
  458. // fit quadratic function with opencv
  459. Mat A = cv::Mat::zeros(cv::Size(3,xx.size()), CV_64FC1);
  460. for(int i=0; i<xx.size();++i){
  461. A.at<double>(i,0) = 1;
  462. A.at<double>(i,1) = xx[i];
  463. A.at<double>(i,2) = xx[i] * xx[i];
  464. }
  465. Mat b = cv::Mat::zeros(cv::Size(1,yy.size()), CV_64FC1);
  466. for(int i = 0; i<yy.size(); ++i){
  467. b.at<double>(i,0) = yy[i];
  468. }
  469. Mat c, d;
  470. c = A.t() * A;
  471. d = A.t() * b;
  472. Mat result = cv::Mat::zeros(cv::Size(1,3),CV_64FC1);
  473. cv::solve(c,d,result);
  474. double a0 = result.at<double>(0, 0);
  475. double a1 = result.at<double>(1, 0);
  476. double a2 = result.at<double>(2, 0);
  477. if(fabs(a2)<1.0e-6){
  478. throw_msg(m_patch_id+" a2 = 0 in solver");
  479. }
  480. oa = -a1 / a2 / 2;
  481. if(oa >180.0){oa -= 180.0;}
  482. return oa;
  483. }
  484. double COptimalAngle::angle_fit_base(map<int, int>&){
  485. return 0.0;
  486. }
  487. ///////////////////////////////////////////////////////////////////////////////
  488. COptimalAnglePart::COptimalAnglePart(ConfigParam&cp, CGcvLogger* pLog/*=0*/)
  489. :COptimalAngle(cp,pLog)
  490. {
  491. }
  492. COptimalAnglePart::~COptimalAnglePart()
  493. {
  494. COptimalAngle::clear_imginfo();
  495. }
  496. // append()加载的图片,满足不是0-180度采样,但保证部分采样中一定存在最大值,
  497. // 直接利用这些数值进行二次函数的拟合
  498. double COptimalAnglePart::angle_fit(map<int,int>& an2width,bool is_base/*=false*/)
  499. {
  500. double oa = 0.0;
  501. std::vector<double> xx;//angle
  502. std::vector<double> yy;//width
  503. map<int,int>::iterator it = an2width.begin();
  504. for(it=an2width.begin(); it!=an2width.end(); ++it)
  505. {
  506. xx.push_back((double)(it->first));
  507. yy.push_back((double)(it->second));
  508. }
  509. //sort
  510. for(size_t i =0;i<xx.size()-1;++i){
  511. for(size_t j = i+1; j<xx.size();++j){
  512. if(xx[j]<xx[i]){
  513. double tmp = xx[i];
  514. xx[i]=xx[j];
  515. xx[j]=tmp;
  516. tmp = yy[i];
  517. yy[i]=yy[j];
  518. yy[j]=tmp;
  519. }
  520. }
  521. }
  522. if(m_pLogger){
  523. stringstream buff;
  524. buff<<m_patch_id;
  525. if(is_base){
  526. buff<<" rootstock_base ";
  527. }
  528. else{
  529. buff<<" rootstock_leaf ";
  530. }
  531. buff<<"opt-angle, (xi,yi): ";
  532. for(size_t i=0;i<xx.size();++i){
  533. buff<<"("<<xx[i]<<","<<yy[i]<<"), ";
  534. }
  535. m_pLogger->INFO(buff.str());
  536. }
  537. if(an2width.size()<3){
  538. throw_msg(m_patch_id+" not enough valid images, NEED append image");
  539. }
  540. // angle field check
  541. if((xx.back()-xx.front())<90){
  542. if(m_pLogger){
  543. stringstream buff;
  544. buff<<m_patch_id<<" x_end="<<xx.back()<<" x_front="<<xx.front();
  545. buff<<" not enough angle field, small than 90 degree.";
  546. m_pLogger->ERRORINFO(buff.str());
  547. }
  548. throw_msg(m_patch_id+" not enough angle field, small than 90 degree");
  549. }
  550. oa = opt_angle_part_impl(xx,yy,is_base);
  551. if(m_pLogger){
  552. stringstream buff;
  553. buff<<m_patch_id;
  554. if(is_base){
  555. buff<<" rootstock_base ";
  556. }
  557. else{
  558. buff<<" rootstock_leaf ";
  559. }
  560. buff<<" opt-angle: "<<oa;
  561. m_pLogger->INFO(buff.str());
  562. }
  563. return oa;
  564. }
  565. double COptimalAnglePart::angle_fit_base(map<int, int>& an2width){
  566. double oa = 0.0;
  567. //find the maximun angle
  568. oa = angle_fit(an2width,true);
  569. //get the minimum angle
  570. if(oa>45.0){oa -= 45.0;}
  571. else{oa +=45.0;}
  572. return oa;
  573. }
  574. void COptimalAnglePart::zero_cross_detect(
  575. vector<int>&d_sign, //input
  576. bool& is_local_max, // whether exists local max
  577. int& local_limit_idx // local limit index of original vector
  578. )
  579. {
  580. int min_zc_idx = -1;
  581. int max_zc_idx = -1;
  582. bool b_local_min =false;
  583. bool b_local_max =false;
  584. bool is_all_down=true;
  585. bool is_all_up = true;
  586. for (size_t i=0; i<d_sign.size();++i){
  587. if (d_sign[i] >0){is_all_down=false;}
  588. if (d_sign[i] <0){is_all_up=false;}
  589. if (i == d_sign.size()-1){continue;}
  590. if (d_sign[i] >=0 && d_sign[i+1] <=0){
  591. //is_local_max=true;
  592. max_zc_idx = i+1;
  593. b_local_max=true;
  594. }
  595. if (d_sign[i] <=0 && d_sign[i+1] >=0){
  596. //is_local_max=false;
  597. min_zc_idx=i+1;
  598. b_local_min=true;
  599. }
  600. }
  601. if(b_local_min){
  602. is_local_max=false;
  603. local_limit_idx = min_zc_idx;
  604. return;
  605. }
  606. if(b_local_max){
  607. is_local_max=true;
  608. local_limit_idx = max_zc_idx;
  609. return;
  610. }
  611. if (is_all_down){
  612. is_local_max=true;
  613. local_limit_idx=0;
  614. return;
  615. }
  616. if( is_all_up){
  617. is_local_max=true;
  618. local_limit_idx=d_sign.size();
  619. return;
  620. }
  621. is_local_max=true;
  622. local_limit_idx=-1;
  623. return;
  624. }
  625. double COptimalAnglePart::limit_min_angle(
  626. vector<double>&x, //sorted, ascending
  627. vector<double>&y,
  628. size_t min_idx
  629. )
  630. {
  631. size_t i = min_idx;
  632. if(fabs(x[i+1]-x[i])<1.0e-6 || fabs(x[i]-x[i-1])<1.0e-6){
  633. throw_msg(m_patch_id+" limit_min_angle, diff_x < 1.0e-6");
  634. }
  635. double k_r = (y[i+1]-y[i])/(x[i+1]-x[i]);
  636. double k_l = (y[i]-y[i-1])/(x[i]-x[i-1]);
  637. double oa = 0.0;
  638. double b_r = 0.0;
  639. double b_l = 0.0;
  640. double cross_x = 0;
  641. if (fabs(k_r) > fabs(k_l)){
  642. // right line
  643. b_r = y[i]-k_r*x[i];
  644. b_l = y[i-1] + k_r*x[i-1];
  645. cross_x = line_cross(k_r,b_r,-k_r,b_l);
  646. oa = cross_x;
  647. }
  648. else{
  649. //# left line
  650. b_l = y[i]-k_l*x[i];
  651. b_r = y[i+1] + k_l*x[i+1];
  652. cross_x = line_cross(k_l,b_l,-k_l,b_r);
  653. oa = cross_x;
  654. }
  655. return oa;
  656. }
  657. double COptimalAnglePart::line_cross(double k1,double b1,double k2,double b2)
  658. {
  659. double x = (b2-b1) / (k1-k2);
  660. return x;
  661. }
  662. double COptimalAnglePart::opt_angle_part_impl(
  663. vector<double>&x, //sorted, ascending
  664. vector<double>&y,
  665. bool is_base/*=false*/)
  666. {
  667. // sign of difference of y
  668. vector<int> dy_sign;
  669. for(size_t i=1; i<y.size();++i){
  670. double dd = y[i] - y[i-1];
  671. if(dd>0){
  672. dy_sign.push_back(1);
  673. }
  674. else{
  675. if(dd<0){dy_sign.push_back(-1);}
  676. else{dy_sign.push_back(0);}
  677. }
  678. }
  679. //local limit search
  680. bool is_local_max = false;
  681. int local_limit_idx = -1;
  682. zero_cross_detect(dy_sign,is_local_max,local_limit_idx);
  683. if(local_limit_idx<0){
  684. throw_msg(m_patch_id+" not found local limit");
  685. }
  686. // optimal angle calculation
  687. vector<double> xx;
  688. vector<double> yy;
  689. double oa=0.0;
  690. if(is_local_max){
  691. //存在局部最大
  692. if( local_limit_idx == 0 || local_limit_idx==y.size()-1){
  693. //#单调
  694. oa = x[local_limit_idx];
  695. /*if (local_limit_idx == 0){
  696. for(size_t i=0;i<3;++i){
  697. xx.push_back(x[i]);
  698. yy.push_back(y[i]);
  699. }
  700. }
  701. else {
  702. for(size_t i=y.size()-3;i<y.size();++i){
  703. xx.push_back(x[i]);
  704. yy.push_back(y[i]);
  705. }
  706. }
  707. oa = qua_fitting(xx,yy); */
  708. if(m_pLogger){
  709. m_pLogger->INFO(m_patch_id+ " angle fit -- signle --- max_angle");
  710. }
  711. }
  712. else{
  713. // 极大值
  714. for(size_t i=local_limit_idx-1;i<local_limit_idx+2;++i){
  715. xx.push_back(x[i]);
  716. yy.push_back(-y[i]);
  717. }
  718. oa = limit_min_angle(xx,yy,1);
  719. if(m_pLogger){
  720. m_pLogger->INFO(m_patch_id+" angle fit -- local max --- 3p insert");
  721. }
  722. }
  723. }
  724. else{
  725. // 极小值
  726. oa = limit_min_angle(x,y,local_limit_idx);
  727. if(is_base){
  728. //for base, offset 45 degree
  729. if(oa >45){oa = oa-45.0;}
  730. else{oa = oa+45.0;}
  731. }
  732. else{
  733. // for rootstock, offset 90 degree
  734. if(oa >90){oa = oa-90.0;}
  735. else{oa = oa+90.0;}
  736. }
  737. if(m_pLogger){
  738. m_pLogger->INFO(m_patch_id+" angle fit -- local min --- 3p insert");
  739. }
  740. }
  741. return oa;
  742. }
  743. CCotyledonAngle::CCotyledonAngle(ConfigParam&cp, CGcvLogger* pLog)
  744. :m_cparam(cp),
  745. m_pLogger(pLog),
  746. m_imgId(""),
  747. m_ppImgSaver(0),
  748. m_height(0),
  749. m_width(0)
  750. {
  751. }
  752. CCotyledonAngle::~CCotyledonAngle(void)
  753. {
  754. }
  755. void CCotyledonAngle::clear_imginfo(){
  756. if (m_pImginfoOpen){
  757. imginfo_release(&m_pImginfoOpen);
  758. m_pImginfoOpen=0;
  759. }
  760. if (m_pImginfoAngle){
  761. imginfo_release(&m_pImginfoAngle);
  762. m_pImginfoAngle=0;
  763. }
  764. }
  765. int CCotyledonAngle::angle_recognize(
  766. ImgInfo* imginfo,
  767. PositionInfo& posinfo
  768. )
  769. {
  770. //clear opencv windows
  771. if(m_cparam.image_show){destroyAllWindows();}
  772. clock_t t;
  773. clock_t t0 = clock();
  774. m_imgId = getImgId(img_type::oa);
  775. /*if(m_pLogger){
  776. m_pLogger->INFO(m_imgId + " cotyledon angle recognize begin.");
  777. }
  778. if(!isvalid(imginfo)){
  779. throw_msg(m_imgId+" invalid input image");
  780. }
  781. if(m_width==0 || m_height==0){
  782. m_width = imginfo->width;
  783. m_height = imginfo->height;
  784. }
  785. if(m_pLogger){
  786. stringstream buff;
  787. buff<<m_imgId<<" load image\twidth="<<imginfo->width
  788. <<"\theight="<<imginfo->height;
  789. m_pLogger->INFO(buff.str());
  790. }
  791. Mat img = imginfo2mat(imginfo);*/
  792. //for DEBUG
  793. Mat img = *(Mat*)(imginfo);
  794. m_width = img.cols/2;
  795. m_height = img.rows/2;
  796. cv::resize(img, img, cv::Size(m_width,m_height ));
  797. //end DEBUG
  798. //imshow_wait("src", img);
  799. if(m_cparam.self_camera){
  800. image_set_bottom(img,20,8);
  801. if(m_pLogger){
  802. m_pLogger->DEBUG(m_imgId+" oa_recognize image set bottom with pixel value 20.");
  803. }
  804. }
  805. if(m_cparam.oa_y_flip){
  806. flip(img,img,0);
  807. if(m_pLogger){
  808. m_pLogger->DEBUG(m_imgId+" oa_recognize image y fliped.");
  809. }
  810. }
  811. if(m_ppImgSaver && *m_ppImgSaver){
  812. (*m_ppImgSaver)->saveImage(img, m_imgId);
  813. if(m_pLogger){
  814. m_pLogger->DEBUG(m_imgId+" oa_recognize after image save.");
  815. }
  816. }
  817. if(m_cparam.image_show){
  818. Mat tmp_b = img.clone();
  819. imshow_wait("oa_image", tmp_b);
  820. }
  821. //1图像分割
  822. this->img_segment(img);
  823. if(m_cparam.image_show){
  824. imshow_wait("oa_image_bin", m_binImg);
  825. }
  826. Mat img_median, img_base, img_open, img_xor;
  827. cv::medianBlur(m_binImg,img_median, 5);
  828. Mat kernel = getStructuringElement(
  829. MORPH_ELLIPSE,
  830. Size( 2*m_cparam.oa_morph_radius + 1, 2*m_cparam.oa_morph_radius+1),
  831. Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
  832. );
  833. morphologyEx(
  834. img_median,
  835. img_base,
  836. MORPH_CLOSE,
  837. kernel,
  838. Point(-1,-1),
  839. 1);
  840. if(m_cparam.image_show){
  841. imshow_wait("img_base", img_base);
  842. }
  843. //2 找到2片叶子,计算各自重心,主方向,面积,横纵比,椭圆度,选择最优,计算需要旋转角度
  844. int iter_step = 3;
  845. int iter_open = 10;
  846. double leaf_area_th = (double)m_cparam.oa_min_leaf_area;
  847. vector<vector<Point>> contours;
  848. vector<Vec4i> hierarchy;
  849. vector<t_leaf> leafs;
  850. while( iter_open >0 && iter_open < 50){
  851. morphologyEx(
  852. img_base,
  853. img_open,
  854. MORPH_OPEN,
  855. kernel,
  856. Point(-1,-1),
  857. iter_open);
  858. imshow_wait("img_base", img_base);
  859. imshow_wait("img_open", img_open);
  860. m_openImg = img_open.clone();
  861. contours.clear();
  862. hierarchy.clear();
  863. leafs.clear();
  864. findContours(img_open,contours,hierarchy,RETR_EXTERNAL,CHAIN_APPROX_NONE);
  865. for(int i=0;i<contours.size();++i){
  866. double area = contourArea(contours[i]);
  867. if (area < leaf_area_th){continue;}
  868. RotatedRect minrect_ellips = fitEllipse(contours[i]);
  869. Moments moment = moments(contours[i]);
  870. double centx = moment.m10 / moment.m00;
  871. double centy = moment.m01 / moment.m00;
  872. t_leaf leaf;
  873. leaf.area = area;
  874. leaf.centx = centx;
  875. leaf.centy = centy;
  876. leaf.minrect_ellips = minrect_ellips;
  877. leaf.contours = contours[i];
  878. leafs.push_back(leaf);
  879. }
  880. if(leafs.size()<2){
  881. iter_open += iter_step;
  882. }
  883. if(leafs.size()>2){
  884. iter_open -= iter_step;
  885. }
  886. if(leafs.size()==2){
  887. break;
  888. }
  889. }
  890. if(leafs.size()!=2){
  891. throw_msg(m_imgId+" leafs detect failed");
  892. }
  893. //select one leaf for angle detection
  894. double min_err = 1.0e16;
  895. int min_err_idx = -1;
  896. for(size_t i =0; i<leafs.size();++i){
  897. double err = fit_error_ellipse(leafs[i].minrect_ellips, leafs[i].contours);
  898. if(err<min_err){
  899. min_err = err;
  900. min_err_idx = i;
  901. }
  902. }
  903. if(min_err_idx<0){
  904. throw_msg(m_imgId+" no normal leaf");
  905. }
  906. //计算旋转角度
  907. double all_area = leafs[0].area+leafs[1].area;
  908. double r0 = leafs[0].area/all_area;
  909. double all_centx = leafs[0].centx * r0 + leafs[1].centx * (1-r0);
  910. double all_centy = leafs[0].centy * r0 + leafs[1].centy * (1 - r0);
  911. float tar_angle1 = 90 + leafs[min_err_idx].minrect_ellips.angle;//长轴与x轴的夹角(图像坐标系下),90-270间
  912. float tar_angle2 = leafs[min_err_idx].minrect_ellips.angle - 90;//-90 ~ 90
  913. if(tar_angle2<0){
  914. tar_angle2+=360;
  915. }
  916. //叶尖方向(相对于两个叶片中心)
  917. float leaf_apex_angle = atan2(leafs[min_err_idx].centy - all_centy, leafs[min_err_idx].centx - all_centx) * 180 / CV_PI;
  918. if(leaf_apex_angle<0){
  919. leaf_apex_angle += 360;
  920. }
  921. float diff1 = fabs(leaf_apex_angle - tar_angle1);
  922. if(diff1>180){
  923. diff1 = 360 - diff1;
  924. }
  925. float diff2 = fabs(leaf_apex_angle - tar_angle2);
  926. if(diff2>180){
  927. diff2 = 360 - diff2;
  928. }
  929. float rot_angle = tar_angle1;
  930. if(diff2 < diff1){
  931. rot_angle = tar_angle2;
  932. }
  933. rot_angle -= 180;//需要旋转的角度
  934. //draw image
  935. Mat ellipse_img = m_openImg.clone();
  936. for(size_t j=0; j<leafs.size(); ++j){
  937. stringstream buff;
  938. buff<<"wa="<<(float)((int)(leafs[j].minrect_ellips.angle*10.0)/10.0);
  939. Scalar color = Scalar(100);
  940. if(j == min_err_idx){
  941. color[0] = 180;
  942. buff<<",rot="<<(float)((int)(rot_angle*10.0)/10.0);
  943. }
  944. Point2f vertices[4];
  945. leafs[j].minrect_ellips.points(vertices);
  946. for (int i = 0; i < 4; i++)
  947. line(ellipse_img, vertices[i], vertices[(i+1)%4], color);
  948. putText(ellipse_img,buff.str(),leafs[j].minrect_ellips.center, 0,0.5,color);
  949. }
  950. if(m_cparam.image_show){
  951. imshow_wait("ellipse_img", ellipse_img);
  952. }
  953. // 3 返回结果
  954. posinfo.rs_oa=(double)rot_angle;
  955. if(m_pLogger){
  956. stringstream buff;
  957. buff<<m_imgId<<" detect rotat angle "<<posinfo.rs_oa;
  958. m_pLogger->INFO(buff.str());
  959. }
  960. //4 返回图像
  961. if(m_cparam.image_return){
  962. this->clear_imginfo();
  963. m_pImginfoOpen = mat2imginfo(m_openImg);
  964. m_pImginfoAngle = mat2imginfo(ellipse_img);
  965. posinfo.pp_images[0]=m_pImginfoOpen;
  966. posinfo.pp_images[1]=m_pImginfoAngle;
  967. if(m_ppImgSaver && *m_ppImgSaver){
  968. (*m_ppImgSaver)->saveImage(m_binImg, m_imgId+"_rst_0");
  969. (*m_ppImgSaver)->saveImage(ellipse_img, m_imgId+"_rst_1");
  970. if(!m_grayImg.empty()){
  971. (*m_ppImgSaver)->saveImage(m_grayImg, m_imgId+"_rst_2");
  972. }
  973. }
  974. }
  975. if(m_pLogger){
  976. m_pLogger->INFO(m_imgId + " oa_recognize image finished.");
  977. }
  978. return 0;
  979. }
  980. double CCotyledonAngle::fit_error_ellipse(
  981. RotatedRect& ellipse_rect,
  982. vector<Point>& pts)
  983. {
  984. /*
  985. ref:https://blog.csdn.net/fangyan90617/article/details/89486331
  986. ellipse common function: A * x^2 + B * x*y + C * y^2 + D * x + E * y + F = 0
  987. reformed function:
  988. A * x'^2 + B * x'*y' + C * y'^2 + F = 0
  989. x' = x - x0
  990. y' = y - y0
  991. from ellipse_rect, calculate A,B,C,F
  992. */
  993. if(pts.size()==0){
  994. return 1.0e6;
  995. }
  996. float A, B,C, F;
  997. get_ellipse_param(ellipse_rect, A, B,C, F);
  998. double mean_err = 0.0;
  999. double x,y,err;
  1000. for(size_t i=0; i<pts.size(); ++i){
  1001. x = (double)pts[i].x - ellipse_rect.center.x;
  1002. y = (double)pts[i].y - ellipse_rect.center.y;
  1003. err = A * x * x + B * x * y + C * y * y + F;
  1004. mean_err += err*err;
  1005. }
  1006. mean_err /= (double)pts.size();
  1007. mean_err = sqrt(mean_err);
  1008. return mean_err;
  1009. }
  1010. void CCotyledonAngle::get_ellipse_param(
  1011. RotatedRect&ellipse_rect,
  1012. float&A,
  1013. float&B,
  1014. float&C,
  1015. float&F)
  1016. {
  1017. float theta = ellipse_rect.angle * CV_PI / 180.0;
  1018. float a = ellipse_rect.size.width / 2.0;
  1019. float b = ellipse_rect.size.height / 2.0;
  1020. A = a * a * sin(theta) * sin(theta) + b * b * cos(theta) * cos(theta);
  1021. B = (-2.0) * (a * a - b * b) * sin(theta) * cos(theta);
  1022. C = a * a * cos(theta) * cos(theta) + b * b * sin(theta) * sin(theta);
  1023. F = (-1.0) * a * a * b * b;
  1024. }
  1025. void CCotyledonAngle::img_segment(Mat&img)
  1026. {
  1027. if(img.channels()!=1){
  1028. //color image ,bgr, for testing
  1029. /*Mat img_gray,b_img;
  1030. cvtColor(img,img_gray,COLOR_BGR2GRAY);
  1031. double th = threshold(
  1032. img_gray,
  1033. b_img,
  1034. 255,
  1035. 255,
  1036. THRESH_OTSU
  1037. );
  1038. b_img/= 255;*/
  1039. Mat img_hsv, b_img_hsv;
  1040. cvtColor(
  1041. img,
  1042. img_hsv,
  1043. COLOR_BGR2HSV
  1044. );
  1045. inRange(
  1046. img_hsv,
  1047. Scalar(30 , 30, 30),
  1048. Scalar(75, 255, 255),
  1049. b_img_hsv
  1050. );
  1051. /*imshow_wait("oa_image_range", b_img_hsv);
  1052. b_img_hsv/=255;
  1053. b_img = b_img.mul(b_img_hsv); */
  1054. //int morph_size = 1;
  1055. Mat kernel = getStructuringElement(
  1056. MORPH_RECT,
  1057. Size( 2*m_cparam.oa_morph_radius + 1, 2*m_cparam.oa_morph_radius+1),
  1058. Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
  1059. );
  1060. morphologyEx(
  1061. b_img_hsv,
  1062. m_binImg,
  1063. MORPH_CLOSE,
  1064. kernel,
  1065. Point(-1,-1),
  1066. m_cparam.oa_morph_iteration);
  1067. }
  1068. else{
  1069. //from imginfo image, gray image
  1070. //int morph_size = 1;
  1071. Mat b_img;
  1072. m_grayImg = img.clone();
  1073. double th = threshold(img, b_img, 255, 255,THRESH_OTSU);
  1074. Mat kernel = getStructuringElement(
  1075. MORPH_RECT,
  1076. Size( 2*m_cparam.oa_morph_radius + 1, 2*m_cparam.oa_morph_radius+1 ),
  1077. Point( m_cparam.oa_morph_radius, m_cparam.oa_morph_radius)
  1078. );
  1079. morphologyEx(
  1080. b_img,
  1081. m_binImg,
  1082. MORPH_OPEN,
  1083. kernel,
  1084. Point(-1,-1),
  1085. m_cparam.oa_morph_iteration
  1086. );
  1087. }
  1088. }
  1089. };