grab_point_rs.cpp 41 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280
  1. /*
  2. 通过点云数据识别抓取位置信息
  3. 1)获取点云
  4. 2)剔除离群点
  5. 3)降采样
  6. 4)植株检测
  7. 5)选出最前,最右侧植株
  8. 6)植株抓取位置检测
  9. */
  10. #include <pcl\io\ply_io.h>
  11. #include <pcl\visualization\cloud_viewer.h>
  12. #include <pcl\filters\crop_box.h>
  13. #include <pcl\filters\radius_outlier_removal.h>
  14. #include <pcl\filters\statistical_outlier_removal.h>
  15. #include <pcl\filters\voxel_grid.h>
  16. #include <pcl\common\common.h>
  17. #include <pcl\features\moment_of_inertia_estimation.h>
  18. #include <pcl\segmentation\sac_segmentation.h>
  19. #include <math.h>
  20. #include "grab_point_rs.h"
  21. #include "utils.h"
  22. #define PI std::acos(-1)
  23. namespace graft_cv {
  24. CRootStockGrabPoint::CRootStockGrabPoint(ConfigParam&cp, CGcvLogger*pLog/*=0*/)
  25. :m_cparam(cp),
  26. m_pLogger(pLog),
  27. m_dtype(0),
  28. m_pcdId(""),
  29. m_ppImgSaver(0)
  30. {
  31. }
  32. CRootStockGrabPoint::~CRootStockGrabPoint() {}
  33. float* CRootStockGrabPoint::get_raw_point_cloud(int &data_size)
  34. {
  35. data_size = m_raw_cloud->width * m_raw_cloud->height;
  36. if (data_size == 0) {
  37. return 0;
  38. }
  39. else {
  40. pcl::PointXYZ* pp = m_raw_cloud->points.data();
  41. return (float*)pp->data;
  42. }
  43. }
  44. int CRootStockGrabPoint::load_data(
  45. float*pPoint,
  46. int pixel_size,
  47. int pt_size,
  48. const char* fn/* = 0*/)
  49. {
  50. int rst = 0;
  51. //1 get point cloud data
  52. if (pPoint != 0 && pt_size>0) {
  53. //read point
  54. m_raw_cloud.reset(new pcl::PointCloud<pcl::PointXYZ>);
  55. int step = pixel_size;
  56. for (int i = 0; i < pt_size; ++i) {
  57. pcl::PointXYZ pt = { pPoint[i*step], pPoint[i*step + 1] , pPoint[i*step + 2] };
  58. m_raw_cloud->push_back(pt);
  59. }
  60. rst = m_raw_cloud->width * m_raw_cloud->height;
  61. if (m_pLogger) {
  62. stringstream buff;
  63. buff << "load raw point cloud " << rst << " data points";
  64. m_pLogger->INFO(buff.str());
  65. }
  66. }
  67. else if (fn != 0) {
  68. //read file
  69. rst = this->read_ply_file(fn);
  70. }
  71. else {//error
  72. if (m_pLogger) {
  73. m_pLogger->ERRORINFO("no valid input");
  74. return (-1);
  75. }
  76. }
  77. if (m_dtype == 0) {
  78. m_pcdId = getImgId(img_type::sola_sc_pcd);
  79. }
  80. else {
  81. m_pcdId = getImgId(img_type::sola_rs_pcd);
  82. }
  83. if (m_ppImgSaver && *m_ppImgSaver) {
  84. (*m_ppImgSaver)->saveBinPly(m_raw_cloud, m_pcdId);
  85. }
  86. if (m_cparam.image_show) {
  87. viewer_cloud(m_raw_cloud, std::string("raw point cloud"));
  88. }
  89. return rst;
  90. }
  91. int CRootStockGrabPoint::grab_point_detect(
  92. int dtype,
  93. PositionInfo& posinfo
  94. )
  95. {
  96. // dtype == 0, scion
  97. // dtype != 0, rootstock
  98. // 输入,穗苗--0, 砧木--1
  99. if (m_raw_cloud->width * m_raw_cloud->height < 1) {
  100. if (m_pLogger) {
  101. stringstream buff;
  102. buff << "m_raw_cloud point cloud " << m_raw_cloud->width * m_raw_cloud->height << " data points";
  103. m_pLogger->ERRORINFO(buff.str());
  104. }
  105. return 1;
  106. }
  107. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  108. //1 crop filter
  109. if (m_pLogger) {
  110. m_pLogger->INFO("begin crop");
  111. }
  112. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  113. pcl::CropBox<pcl::PointXYZ> box_filter;
  114. m_dtype = dtype;
  115. if (dtype != 0) {//rootstock
  116. box_filter.setMin(Eigen::Vector4f(m_cparam.rs_grab_xmin, m_cparam.rs_grab_ymin, m_cparam.rs_grab_zmin, 1));
  117. box_filter.setMax(Eigen::Vector4f(m_cparam.rs_grab_xmax, m_cparam.rs_grab_ymax, m_cparam.rs_grab_zmax, 1));
  118. }
  119. else {//scion
  120. box_filter.setMin(Eigen::Vector4f(m_cparam.sc_grab_xmin, m_cparam.sc_grab_ymin, m_cparam.sc_grab_zmin, 1));
  121. box_filter.setMax(Eigen::Vector4f(m_cparam.sc_grab_xmax, m_cparam.sc_grab_ymax, m_cparam.sc_grab_zmax, 1));
  122. }
  123. box_filter.setNegative(false);
  124. box_filter.setInputCloud(m_raw_cloud);
  125. box_filter.filter(*cloud_inbox);
  126. if (m_pLogger) {
  127. stringstream buff;
  128. buff << "CropBox croped point cloud " << cloud_inbox->width * cloud_inbox->height << " data points";
  129. m_pLogger->INFO(buff.str());
  130. }
  131. if (cloud_inbox->width * cloud_inbox->height < 1) {
  132. return 1;
  133. }
  134. if (m_cparam.image_show) {
  135. viewer_cloud(cloud_inbox, std::string("cloud_inbox"));
  136. }
  137. if (m_pLogger) {
  138. m_pLogger->INFO("end crop");
  139. }
  140. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  141. //2 filtler with radius remove
  142. if (m_pLogger) {
  143. m_pLogger->INFO("begin ror");
  144. }
  145. int nb_points = 20;
  146. double stem_radius = m_cparam.rs_grab_stem_diameter / 2.0;
  147. if(dtype == 0){
  148. stem_radius = m_cparam.sc_grab_stem_diameter / 2.0;
  149. }
  150. double remove_radius = stem_radius * 2.0;
  151. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_ror(new pcl::PointCloud<pcl::PointXYZ>);
  152. pcl::RadiusOutlierRemoval<pcl::PointXYZ> ror(false);
  153. ror.setInputCloud(cloud_inbox);
  154. ror.setRadiusSearch(remove_radius);
  155. ror.setMinNeighborsInRadius(nb_points);
  156. ror.filter(*cloud_ror);
  157. if (m_pLogger) {
  158. stringstream buff;
  159. buff << "RadiusOutlierRemoval filtered point cloud " << cloud_ror->width * cloud_ror->height << " data points. param stem_radius="<<
  160. stem_radius<<", nb_points="<< nb_points;
  161. m_pLogger->INFO(buff.str());
  162. }
  163. if (cloud_ror->width * cloud_ror->height < 1) {
  164. return 1;
  165. }
  166. /*if (m_cparam.image_show) {
  167. viewer_cloud(cloud_ror, std::string("cloud_ror"));
  168. }*/
  169. if (m_pLogger) {
  170. m_pLogger->INFO("end ror");
  171. }
  172. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  173. //3 voxel grid down sampling
  174. if (m_pLogger) {
  175. m_pLogger->INFO("begin down");
  176. }
  177. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_dowm_sampled(new pcl::PointCloud<pcl::PointXYZ>);
  178. pcl::VoxelGrid<pcl::PointXYZ> outrem;
  179. outrem.setInputCloud(cloud_ror);
  180. outrem.setLeafSize(stem_radius, stem_radius, stem_radius);
  181. outrem.filter(*cloud_dowm_sampled);
  182. if (m_pLogger) {
  183. stringstream buff;
  184. buff << "VoxelGrid dowm_sampled point cloud " << cloud_dowm_sampled->width * cloud_dowm_sampled->height << " data points (if < 50, return)";
  185. m_pLogger->INFO(buff.str());
  186. }
  187. if (cloud_dowm_sampled->width * cloud_dowm_sampled->height < 50) {
  188. return 1;
  189. }
  190. /*if (m_cparam.image_show) {
  191. viewer_cloud(cloud_dowm_sampled, std::string("cloud_dowm_sampled"));
  192. }*/
  193. if (m_pLogger) {
  194. m_pLogger->INFO("end down");
  195. }
  196. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  197. //4 seedling position,找到第一个位置上的植株
  198. std::vector<int> first_seedling_cloud_idx; //存储找到的第一个位置上植株茎上直线点的index
  199. pcl::PointXYZ xz_center; //存储找到的第一个位置上植株根部坐标
  200. if (m_pLogger) {
  201. m_pLogger->INFO("begin find seedling position");
  202. }
  203. //find_seedling_position_line_based(cloud_ror, first_seedling_cloud_idx, xz_center);
  204. find_seedling_position(cloud_dowm_sampled, first_seedling_cloud_idx, xz_center);
  205. if (m_pLogger) {
  206. stringstream buff;
  207. buff << "after find_seedling_position(), foud first seedling seeds points size " << first_seedling_cloud_idx.size();
  208. m_pLogger->INFO(buff.str());
  209. }
  210. if (m_pLogger) {
  211. m_pLogger->INFO("end find seedling position");
  212. }
  213. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  214. //5 nearest seedling grab point selection,找到最接近指定高度的最优抓取点坐标
  215. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_seedling_seed(new pcl::PointCloud<pcl::PointXYZ>);
  216. pcl::copyPointCloud(*cloud_dowm_sampled, first_seedling_cloud_idx, *cloud_seedling_seed);
  217. //改进思路:将茎直线上cloud_seedling_seed的点,提取周围邻域xz的宽度,在相同邻域cloud_dowm_sampled点云内提取xz的宽度
  218. //对比两个宽度,变化很大的,说明周边有异物(叶柄或叶子),不适合作为抓取位置;否则就是抓取位置
  219. //在众多抓取位置上,选择离指定高度最近的点作为抓取位置
  220. //
  221. pcl::PointXYZ selected_pt;
  222. int selected_pt_idx;
  223. std::vector<int> optimal_seeds_idx;
  224. get_optimal_seed_compare(cloud_dowm_sampled, cloud_seedling_seed, selected_pt, selected_pt_idx, optimal_seeds_idx);
  225. if (selected_pt_idx < 0) {
  226. if (m_pLogger) {
  227. m_pLogger->ERRORINFO("Not found valid grab point");
  228. }
  229. return 1;
  230. }
  231. if (m_pLogger) {
  232. m_pLogger->INFO("end get_optimal_seed");
  233. }
  234. posinfo.rs_grab_x = selected_pt.x;
  235. posinfo.rs_grab_y = selected_pt.y;
  236. posinfo.rs_grab_z = selected_pt.z;
  237. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  238. //6 debug 显示结果
  239. if (m_cparam.image_show) {
  240. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_cand_demo(new pcl::PointCloud<pcl::PointXYZRGB>);
  241. pcl::copyPointCloud(*cloud_dowm_sampled, *cloud_cand_demo);
  242. for (auto& pt : *cloud_cand_demo) {
  243. pt.r = 255;
  244. pt.g = 255;
  245. pt.b = 255;
  246. }
  247. int cnt = 0;
  248. for (auto& idx : optimal_seeds_idx) {
  249. int p_idx = first_seedling_cloud_idx.at(idx);
  250. /*if (p_idx == optimal_seeds_idx[selected_pt_idx]) {
  251. cloud_cand_demo->points[p_idx].r = 0;
  252. cloud_cand_demo->points[p_idx].g = 255;
  253. cloud_cand_demo->points[p_idx].b = 0;
  254. }
  255. else {*/
  256. cloud_cand_demo->points.at(p_idx).r = 255;
  257. cloud_cand_demo->points.at(p_idx).g = 0;
  258. cloud_cand_demo->points.at(p_idx).b = 0;
  259. /*} */
  260. if (cnt > optimal_seeds_idx.size()) { break; }
  261. cnt++;
  262. }
  263. pcl::PointXYZRGB pt_grab = {0,255,0};
  264. pt_grab.x = selected_pt.x;
  265. pt_grab.y = selected_pt.y;
  266. pt_grab.z = selected_pt.z;
  267. pcl::PointXYZRGB pt_grab_ = { 0,255,120 };
  268. pt_grab_.x = selected_pt.x;
  269. pt_grab_.y = selected_pt.y+0.2;
  270. pt_grab_.z = selected_pt.z;
  271. cloud_cand_demo->push_back(pt_grab);
  272. //viewer_cloud(cloud_cand_demo, std::string("cloud_cand_demo"));
  273. viewer_cloud_debug(cloud_cand_demo, selected_pt, xz_center, std::string("cloud_cand_demo"));
  274. }
  275. return 0;
  276. }
  277. int CRootStockGrabPoint::read_ply_file(const char* fn)
  278. {
  279. m_raw_cloud.reset( new pcl::PointCloud<pcl::PointXYZ>);
  280. if (pcl::io::loadPLYFile<pcl::PointXYZ>(fn, *m_raw_cloud) == -1) {
  281. if (m_pLogger) {
  282. m_pLogger->ERRORINFO("could not load file: "+std::string(fn));
  283. }
  284. return (-1);
  285. }
  286. if (m_pLogger) {
  287. stringstream buff;
  288. buff << "load raw point cloud " << m_raw_cloud->width * m_raw_cloud->height << " data points";
  289. m_pLogger->INFO(buff.str());
  290. }
  291. return m_raw_cloud->width * m_raw_cloud->height;
  292. }
  293. double CRootStockGrabPoint::compute_nearest_neighbor_distance(pcl::PointCloud<pcl::PointXYZ>::Ptr pcd)
  294. {
  295. pcl::KdTreeFLANN<pcl::PointXYZ> tree;
  296. tree.setInputCloud(pcd);
  297. int k = 2;
  298. double res = 0.0;
  299. int n_points = 0;
  300. for (size_t i = 0; i < pcd->size(); ++i) {
  301. std::vector<int> idx(k);
  302. std::vector<float> sqr_distances(k);
  303. if (tree.nearestKSearch(i, k, idx, sqr_distances) == k) {
  304. for (int id = 1; id < k; ++id) {
  305. res += sqrt(sqr_distances.at(id));
  306. ++n_points;
  307. }
  308. }
  309. }
  310. if (n_points > 0) {
  311. res /= (double)n_points;
  312. }
  313. return res;
  314. }
  315. //void CRootStockGrabPoint::find_tray_top_edge(
  316. //pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  317. //float & tray_top_edge_y
  318. //)
  319. //{
  320. // tray_top_edge_y = -1000.0;
  321. // //以毫米为单位构建vector
  322. // pcl::PointXYZ min_point_aabb;
  323. // pcl::PointXYZ max_point_aabb;
  324. // pcl::MomentOfInertiaEstimation<pcl::PointXYZ> feature_extractor;
  325. // feature_extractor.setInputCloud(in_cloud);
  326. // feature_extractor.compute();
  327. // feature_extractor.getAABB(min_point_aabb, max_point_aabb);
  328. // float miny = min_point_aabb.y;
  329. // float maxy = max_point_aabb.y;
  330. // if(maxy - miny <5.0){
  331. // tray_top_edge_y = maxy;
  332. // return;
  333. // }
  334. // std::vector<int> y_cnt_hist(int(maxy - miny), 0);
  335. // for(auto& pt : in_cloud->points){
  336. // int idx = (int)(pt.y - miny);
  337. // if(idx<0){continue;}
  338. // if (idx >= y_cnt_hist.size()) { continue; }
  339. // y_cnt_hist.at(idx) += 1;
  340. // }
  341. // //从上半部分找到最大值,作为平面上沿
  342. // int idx_part = (int)(y_cnt_hist.size() / 2);
  343. // int idx_edge = std::max_element(y_cnt_hist.begin(), y_cnt_hist.begin() + idx_part) - y_cnt_hist.begin();
  344. // tray_top_edge_y = miny + (float)(idx_edge + 2.0);
  345. //}
  346. //void CRootStockGrabPoint::find_seedling_position_line_based(
  347. // pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  348. // std::vector<int> &first_seedling_cloud_idx,
  349. // pcl::PointXYZ&xz_center
  350. //)
  351. //{
  352. // //找到穴盘上沿z,最优抓取的z,再在最优抓取的z基础上加上3作为有效范围
  353. // float tray_y = -1000.0;
  354. // float top_y = 0.0;
  355. // find_tray_top_edge(in_cloud, tray_y);
  356. // if (tray_y < -500.0) {
  357. // if (m_dtype == 0) {
  358. // //scion
  359. // tray_y = m_cparam.sc_grab_y_opt-2.0;
  360. // }
  361. // else {
  362. // tray_y = m_cparam.rs_grab_y_opt-2.0;
  363. // }
  364. // }
  365. // //up limit
  366. // if (m_dtype == 0) {
  367. // top_y = m_cparam.sc_grab_zmax;
  368. // if (top_y > m_cparam.sc_grab_y_opt + 3.0) {
  369. // top_y = m_cparam.sc_grab_y_opt + 3.0;
  370. // }
  371. // }
  372. // else {
  373. // top_y = m_cparam.rs_grab_zmax;
  374. // if (top_y > m_cparam.rs_grab_y_opt + 3.0) {
  375. // top_y = m_cparam.rs_grab_y_opt + 3.0;
  376. // }
  377. // }
  378. // //sub cloud
  379. // pcl::PointCloud<pcl::PointXYZ>::Ptr sub_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  380. // pcl::CropBox<pcl::PointXYZ> box_filter;
  381. // if (m_dtype != 0) {//rootstock
  382. // box_filter.setMin(Eigen::Vector4f(m_cparam.rs_grab_xmin, tray_y, m_cparam.rs_grab_zmin, 1));
  383. // box_filter.setMax(Eigen::Vector4f(m_cparam.rs_grab_xmax, top_y, m_cparam.rs_grab_zmax, 1));
  384. // }
  385. // else {//scion
  386. // box_filter.setMin(Eigen::Vector4f(m_cparam.sc_grab_xmin, tray_y, m_cparam.sc_grab_zmin, 1));
  387. // box_filter.setMax(Eigen::Vector4f(m_cparam.sc_grab_xmax, top_y, m_cparam.sc_grab_zmax, 1));
  388. // }
  389. // box_filter.setNegative(false);
  390. // box_filter.setInputCloud(in_cloud);
  391. // box_filter.filter(*sub_cloud);
  392. // if (m_cparam.image_show) {
  393. // viewer_cloud(sub_cloud, std::string("sub inbox"));
  394. // }
  395. // //在sub_cloud中进行直线检测
  396. // segement_line(sub_cloud);
  397. //}
  398. //void CRootStockGrabPoint::segement_line(
  399. // pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud
  400. //)
  401. //{
  402. // pcl::ModelCoefficients::Ptr coefficients(new pcl::ModelCoefficients);
  403. // pcl::PointIndices::Ptr inliers(new pcl::PointIndices);
  404. // pcl::SACSegmentation<pcl::PointXYZ> seg;
  405. // seg.setOptimizeCoefficients(true);
  406. // seg.setModelType(pcl::SACMODEL_LINE);
  407. // seg.setMethodType(pcl::SAC_RANSAC);
  408. // seg.setDistanceThreshold(0.5);
  409. // seg.setInputCloud(in_cloud);
  410. // seg.segment(*inliers, *coefficients);
  411. // if (m_cparam.image_show) {
  412. // pcl::PointCloud<pcl::PointXYZ>::Ptr line_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  413. // pcl::copyPointCloud(*in_cloud, *inliers, *line_cloud);
  414. // viewer_cloud(line_cloud, std::string("cloud_line"));
  415. // }
  416. //}
  417. ////////////////////////////////////////////////////////////
  418. void CRootStockGrabPoint::find_seedling_position(
  419. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  420. std::vector<int> &first_seedling_cloud_idx,
  421. pcl::PointXYZ&xz_center
  422. )
  423. {
  424. float hole_step = 40.0; //穴盘中穴间距
  425. float hole_step_radius = hole_step / 2.0;
  426. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2d(new pcl::PointCloud < pcl::PointXYZ>);
  427. pcl::copyPointCloud(*in_cloud, *cloud2d);
  428. for (auto&pt : *cloud2d) {
  429. pt.y = 0.0;
  430. }
  431. /*if(m_cparam.image_show){
  432. viewer_cloud(cloud2d, std::string("cloud2d"));
  433. } */
  434. double radius = m_cparam.rs_grab_stem_diameter;
  435. if (m_dtype == 0) {
  436. radius = m_cparam.sc_grab_stem_diameter;
  437. }
  438. std::vector<int> counter;
  439. pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;
  440. kdtree.setInputCloud(cloud2d);
  441. std::vector<int>idx;
  442. std::vector<float>dist_sqr;
  443. for (size_t i = 0; i < cloud2d->points.size(); ++i) {
  444. int k = kdtree.radiusSearch(cloud2d->points.at(i), radius, idx, dist_sqr);
  445. counter.push_back(k);
  446. }
  447. int th = (int)(otsu(counter));
  448. std::vector<int> index;
  449. for (size_t i = 0; i < counter.size(); ++i) {
  450. if (counter.at(i) >= th) {
  451. index.push_back(i);
  452. }
  453. }
  454. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2d_pos(new pcl::PointCloud < pcl::PointXYZ>);
  455. pcl::copyPointCloud(*cloud2d, index, *cloud2d_pos);
  456. if (m_pLogger) {
  457. stringstream buff;
  458. buff << "get 2d seedling seed point cloud " << index.size()<< " data points with thrreshold "<<th;
  459. m_pLogger->INFO(buff.str());
  460. }
  461. if (m_cparam.image_show) {
  462. viewer_cloud(cloud2d_pos, std::string("cloud2d_pos"));
  463. }
  464. //clustering
  465. double d1 = m_cparam.rs_grab_stem_diameter;
  466. double d2 = m_cparam.rs_grab_stem_diameter * 3.0;
  467. if (m_dtype == 0) {
  468. d1 = m_cparam.sc_grab_stem_diameter;
  469. d2 = m_cparam.sc_grab_stem_diameter * 3.0;
  470. }
  471. std::vector<pcl::PointXYZ>cluster_center;
  472. std::vector<std::vector<int>> cluster_member;
  473. euclidean_clustering_ttsas(cloud2d_pos, d1, d2, cluster_center, cluster_member);
  474. if (m_pLogger) {
  475. stringstream buff;
  476. buff << "euclidean_clustering_ttsas: " << cluster_center.size() << " t1=" << d1
  477. << " t2=" << d2;
  478. m_pLogger->INFO(buff.str());
  479. }
  480. if (m_cparam.image_show) {
  481. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_cluster(new pcl::PointCloud<pcl::PointXYZRGB>);
  482. pcl::copyPointCloud(*cloud2d_pos, *cloud_cluster);
  483. for (auto& pt : *cloud_cluster) {
  484. pt.r = 255;
  485. pt.g = 255;
  486. pt.b = 255;
  487. }
  488. viewer_cloud_cluster(cloud_cluster, cluster_center, std::string("cloud2d_cluster_pos"));
  489. }
  490. //对每个类别进行检测,剔除不是茎的类别
  491. std::vector<pcl::PointXYZ> clt_root;
  492. std::vector<pcl::PointXYZ> clt_box;
  493. std::vector<std::vector<int> > clt_line_idx;
  494. for (size_t i = 0; i < cluster_center.size(); ++i) {
  495. pcl::PointIndices cluster_idxs;
  496. for (auto idx_clt : cluster_member.at(i)) {
  497. int idx_raw = index.at(idx_clt);
  498. cluster_idxs.indices.push_back(idx_raw);
  499. }
  500. pcl::PointCloud<pcl::PointXYZ>::Ptr clt_cloud(new pcl::PointCloud < pcl::PointXYZ>);
  501. pcl::copyPointCloud(*in_cloud, cluster_idxs, *clt_cloud);
  502. //计算每一个茎的外接矩形
  503. pcl::PointXYZ min_point_aabb;
  504. pcl::PointXYZ max_point_aabb;
  505. pcl::MomentOfInertiaEstimation<pcl::PointXYZ> feature_extractor;
  506. feature_extractor.setInputCloud(clt_cloud);
  507. feature_extractor.compute();
  508. feature_extractor.getAABB(min_point_aabb, max_point_aabb);
  509. //扩展矩形范围
  510. pcl::PointXYZ min_point_aabb_ex;
  511. pcl::PointXYZ max_point_aabb_ex;
  512. min_point_aabb_ex.x = 0.5*(min_point_aabb.x + max_point_aabb.x) - hole_step_radius;
  513. min_point_aabb_ex.y = m_cparam.sc_grab_ymin;
  514. if (m_dtype != 0) { min_point_aabb_ex.y = m_cparam.rs_grab_ymin; }
  515. min_point_aabb_ex.z = 0.5*(min_point_aabb.z + max_point_aabb.z) - hole_step_radius;
  516. max_point_aabb_ex.x = 0.5*(min_point_aabb.x + max_point_aabb.x) + hole_step_radius;
  517. max_point_aabb_ex.y = m_cparam.sc_grab_ymax;
  518. if (m_dtype != 0) { max_point_aabb_ex.y = m_cparam.rs_grab_ymax; }
  519. max_point_aabb_ex.z = 0.5*(min_point_aabb.z + max_point_aabb.z) + hole_step_radius;
  520. /////////////////////////////////
  521. //扩展矩形内直线检测
  522. std::vector<int>line_idx;
  523. find_line_in_cube(in_cloud, min_point_aabb_ex, max_point_aabb_ex, line_idx);
  524. clt_line_idx.push_back(line_idx);
  525. //找到line_idx中y最小的那个点的idx
  526. int line_root_idx = -1;
  527. float line_root_y_min = 10000.0;
  528. for (size_t lidx = 0; lidx < line_idx.size();++lidx) {
  529. int raw_idx = line_idx.at(lidx);
  530. if (in_cloud->at(raw_idx).y < line_root_y_min) {
  531. line_root_y_min = in_cloud->at(raw_idx).y;
  532. line_root_idx = raw_idx;
  533. }
  534. }
  535. //找到底部中心点,然后找到根部坐标
  536. pcl::PointXYZ pt_root;
  537. pt_root = in_cloud->at(line_root_idx);
  538. clt_root.push_back(pt_root);
  539. clt_box.push_back(min_point_aabb_ex);
  540. clt_box.push_back(max_point_aabb_ex);
  541. //viewer_cloud(clt_cloud, std::string("cluster"));
  542. }
  543. if (m_cparam.image_show) {
  544. viewer_cloud_cluster_box(in_cloud, clt_root, clt_box, clt_line_idx, std::string("cloud2d_cluster_box"));
  545. }
  546. //sort cluster center, get the frontest leftest one
  547. std::vector<float> cluster_index;
  548. for (auto&pt : clt_root) {
  549. float idx = pt.x - pt.z;
  550. cluster_index.push_back(idx);
  551. }
  552. int first_idx = std::min_element(cluster_index.begin(), cluster_index.end()) - cluster_index.begin();
  553. if (m_dtype == 0) {
  554. first_idx = std::max_element(cluster_index.begin(), cluster_index.end()) - cluster_index.begin();
  555. }
  556. first_seedling_cloud_idx.clear();
  557. for (auto&v : clt_line_idx.at(first_idx)) {
  558. first_seedling_cloud_idx.push_back(v);
  559. }
  560. xz_center.x = clt_root.at(first_idx).x;
  561. xz_center.y = clt_root.at(first_idx).y;
  562. xz_center.z = clt_root.at(first_idx).z;
  563. if (m_pLogger) {
  564. stringstream buff;
  565. buff << "euclidean_clustering_ttsas, find cluster center(" << xz_center.x
  566. <<", "<< xz_center.y<<", "<< xz_center.z<<")";
  567. m_pLogger->INFO(buff.str());
  568. }
  569. }
  570. void CRootStockGrabPoint::find_line_in_cube(
  571. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //输入,整体点云
  572. pcl::PointXYZ& min_pt_ex, //输入,
  573. pcl::PointXYZ& max_pt_ex, //输入,
  574. std::vector<int>& out_idx //输出,直线上点的index, 基于输入整体点云
  575. )
  576. {
  577. out_idx.clear();
  578. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  579. pcl::CropBox<pcl::PointXYZ> box_filter;
  580. std::vector<int>inbox_idx;
  581. box_filter.setMin(Eigen::Vector4f(min_pt_ex.x, min_pt_ex.y, min_pt_ex.z, 1));
  582. box_filter.setMax(Eigen::Vector4f(max_pt_ex.x, max_pt_ex.y, max_pt_ex.z, 1));
  583. box_filter.setNegative(false);
  584. box_filter.setInputCloud(in_cloud);
  585. box_filter.filter(inbox_idx);
  586. pcl::copyPointCloud(*in_cloud, inbox_idx, *cloud_inbox);
  587. //找到inbox点云中的直线
  588. float stem_radius = m_cparam.rs_grab_stem_diameter / 2.0;
  589. if (m_dtype == 0) {
  590. stem_radius = m_cparam.sc_grab_stem_diameter / 2.0;
  591. }
  592. pcl::ModelCoefficients::Ptr coefficients(new pcl::ModelCoefficients);
  593. pcl::PointIndices::Ptr inliers(new pcl::PointIndices);
  594. pcl::SACSegmentation<pcl::PointXYZ> seg;
  595. seg.setOptimizeCoefficients(true);
  596. seg.setModelType(pcl::SACMODEL_LINE);
  597. seg.setDistanceThreshold(stem_radius);
  598. seg.setInputCloud(cloud_inbox);
  599. seg.segment(*inliers, *coefficients);
  600. int idx_raw = 0;
  601. for (auto& idx : inliers->indices) {
  602. idx_raw = inbox_idx.at(idx);
  603. out_idx.push_back(idx_raw);
  604. }
  605. if (m_cparam.image_show) {
  606. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_line(new pcl::PointCloud<pcl::PointXYZ>);
  607. pcl::copyPointCloud(*cloud_inbox, *inliers, *cloud_line);
  608. pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("line in cluster"));
  609. viewer->addPointCloud<pcl::PointXYZ>(cloud_inbox, "cluster");
  610. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "cluster");
  611. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "cluster");
  612. viewer->addPointCloud(cloud_line, "fit_line");
  613. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0, 1, 0, "fit_line");
  614. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, "fit_line");
  615. while (!viewer->wasStopped()) {
  616. viewer->spinOnce(100);
  617. }
  618. }
  619. }
  620. void CRootStockGrabPoint::crop_nn_analysis(
  621. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  622. pcl::PointCloud<pcl::PointXYZ>::Ptr seed_cloud,
  623. double dist_mean,
  624. std::vector<double>&mass_indices,
  625. std::vector<int>& candidate_idx
  626. )
  627. {
  628. candidate_idx.clear();
  629. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  630. pcl::CropBox<pcl::PointXYZ> box_filter;
  631. box_filter.setNegative(false);
  632. box_filter.setInputCloud(in_cloud);
  633. double radius = 5;
  634. double rx = 0.8;
  635. double ry = 1.0;
  636. double rz = 0.8;
  637. double cx, cy, cz;
  638. double dz;
  639. for (size_t i = 0; i < seed_cloud->points.size(); ++i) {
  640. cx = seed_cloud->points.at(i).x;
  641. cy = seed_cloud->points.at(i).y;
  642. cz = seed_cloud->points.at(i).z;
  643. box_filter.setMin(Eigen::Vector4f(cx - rx*radius, cy - ry*radius, cz - rz*radius, 1));
  644. box_filter.setMax(Eigen::Vector4f(cx + rx*radius, cy + ry*radius, cz + rz*radius, 1));
  645. box_filter.filter(*cloud_inbox);
  646. //dz
  647. Eigen::Vector4f min_point;
  648. Eigen::Vector4f max_point;
  649. pcl::getMinMax3D(*cloud_inbox, min_point, max_point);
  650. dz = max_point(2) - min_point(2);
  651. //project to xy-plane
  652. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox_xy(new pcl::PointCloud<pcl::PointXYZ>);
  653. pcl::copyPointCloud(*cloud_inbox, *cloud_inbox_xy);
  654. for (auto&pt : *cloud_inbox_xy) { pt.z = 0.0; }
  655. //pca
  656. double dx_obb;
  657. double dy_obb;
  658. double angle_obb;
  659. cal_obb_2d(cloud_inbox_xy, 0, dx_obb, dy_obb, angle_obb);
  660. try {
  661. double dx_grid = dx_obb / dist_mean;
  662. double dy_grid = dy_obb / dist_mean;
  663. double dz_grid = dz / dist_mean;
  664. double fill_ratio = cloud_inbox->points.size() / dx_grid / dy_grid / dz_grid;
  665. double y_ratio = dy_obb / dx_obb/2;
  666. y_ratio = pow(y_ratio, 2);
  667. double a_ratio = cos((angle_obb - 90)*PI / 180.0);
  668. a_ratio = pow(a_ratio, 3);
  669. double mass_index = fill_ratio*y_ratio*a_ratio;
  670. mass_indices.push_back(mass_index);
  671. if (m_pLogger) {
  672. stringstream buff;
  673. buff << i << "/" << seed_cloud->points.size() << " dx=" << dx_obb << ", dy=" << dy_obb << ", fill_ratio=" << fill_ratio
  674. << ", y_ratio=" << y_ratio << ", a_ratio=" << a_ratio << ", mass_index=" << mass_index;
  675. m_pLogger->INFO(buff.str());
  676. }
  677. }
  678. catch (...) {
  679. mass_indices.push_back(0);
  680. }
  681. }
  682. // sort by mass_indices
  683. std::vector<size_t> sort_idx = sort_indexes_e(mass_indices, false);
  684. for (auto&v : sort_idx) { candidate_idx.push_back((int)v); }
  685. }
  686. void CRootStockGrabPoint::euclidean_clustering_ttsas(
  687. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  688. double d1, double d2,
  689. std::vector<pcl::PointXYZ>&cluster_center,
  690. std::vector<std::vector<int>> &clustr_member
  691. )
  692. {
  693. if (m_pLogger) {
  694. stringstream buff;
  695. buff << "euclidean_clustering_ttsas() begin...";
  696. m_pLogger->INFO(buff.str());
  697. }
  698. std::vector<int> cluster_weight;
  699. std::vector<int> data_stat;
  700. std::vector<pcl::PointXYZ>cluster_center_raw;
  701. std::vector<std::vector<int>> clustr_member_raw;
  702. for (size_t i = 0; i < in_cloud->points.size(); ++i) { data_stat.push_back(0); }
  703. size_t data_len = in_cloud->points.size();
  704. int exists_change = 0;
  705. int prev_change = 0;
  706. int cur_change = 0;
  707. int m = 0;
  708. while (std::find(data_stat.begin(), data_stat.end(), 0) != data_stat.end()) {
  709. bool new_while_first = true;
  710. for (size_t i = 0; i < data_len; ++i) {
  711. if (data_stat.at(i) == 0 && new_while_first && exists_change == 0) {
  712. new_while_first = false;
  713. std::vector<int> idx;
  714. idx.push_back(i);
  715. clustr_member_raw.push_back(idx);
  716. pcl::PointXYZ center;
  717. center.x = in_cloud->points.at(i).x;
  718. center.y = in_cloud->points.at(i).y;
  719. center.z = in_cloud->points.at(i).z;
  720. cluster_center_raw.push_back(center);
  721. data_stat.at(i) = 1;
  722. cur_change += 1;
  723. cluster_weight.push_back(1);
  724. m += 1;
  725. }
  726. else if (data_stat[i] == 0) {
  727. std::vector<float> distances;
  728. for (size_t j = 0; j < clustr_member_raw.size(); ++j) {
  729. std::vector<float> distances_sub;
  730. for (size_t jj = 0; jj < clustr_member_raw.at(j).size(); ++jj) {
  731. size_t ele_idx = clustr_member_raw.at(j).at(jj);
  732. double d = sqrt(
  733. (in_cloud->points.at(i).x - in_cloud->points.at(ele_idx).x) * (in_cloud->points.at(i).x - in_cloud->points.at(ele_idx).x) +
  734. (in_cloud->points.at(i).y - in_cloud->points.at(ele_idx).y) * (in_cloud->points.at(i).y - in_cloud->points.at(ele_idx).y) +
  735. (in_cloud->points.at(i).z - in_cloud->points.at(ele_idx).z) * (in_cloud->points.at(i).z - in_cloud->points.at(ele_idx).z));
  736. distances_sub.push_back(d);
  737. }
  738. double min_dist = *std::min_element(distances_sub.begin(), distances_sub.end());
  739. distances.push_back(min_dist);
  740. }
  741. int min_idx = std::min_element(distances.begin(), distances.end()) - distances.begin();
  742. if (distances.at(min_idx) < d1) {
  743. data_stat.at(i) = 1;
  744. double w = cluster_weight.at(min_idx);
  745. cluster_weight.at(min_idx) += 1;
  746. clustr_member_raw.at(min_idx).push_back(i);
  747. cluster_center_raw.at(min_idx).x = (cluster_center_raw.at(min_idx).x * w + in_cloud->points.at(i).x) / (w + 1);
  748. cluster_center_raw.at(min_idx).y = (cluster_center_raw.at(min_idx).y * w + in_cloud->points.at(i).y) / (w + 1);
  749. cluster_center_raw.at(min_idx).z = (cluster_center_raw.at(min_idx).z * w + in_cloud->points.at(i).z) / (w + 1);
  750. cur_change += 1;
  751. }
  752. else if (distances.at(min_idx) > d2) {
  753. std::vector<int> idx;
  754. idx.push_back(i);
  755. clustr_member_raw.push_back(idx);
  756. pcl::PointXYZ center;
  757. center.x = in_cloud->points.at(i).x;
  758. center.y = in_cloud->points.at(i).y;
  759. center.z = in_cloud->points.at(i).z;
  760. cluster_center_raw.push_back(center);
  761. data_stat.at(i) = 1;
  762. cur_change += 1;
  763. cluster_weight.push_back(1);
  764. m += 1;
  765. }
  766. }
  767. else if (data_stat.at(i)== 1) {
  768. cur_change += 1;
  769. }
  770. }
  771. exists_change = fabs(cur_change - prev_change);
  772. prev_change = cur_change;
  773. cur_change = 0;
  774. }
  775. // copy result
  776. for (size_t i = 0; i < clustr_member_raw.size(); ++i) {
  777. if (clustr_member_raw.at(i).size() < 20) { continue; }
  778. cluster_center.push_back(cluster_center_raw.at(i));
  779. clustr_member.push_back(clustr_member_raw.at(i));
  780. }
  781. if (m_pLogger) {
  782. stringstream buff;
  783. buff << "euclidean_clustering_ttsas() end";
  784. m_pLogger->INFO(buff.str());
  785. }
  786. }
  787. void CRootStockGrabPoint::cal_obb_2d(
  788. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  789. int axis, //0--xy, 1--zy
  790. double &dx_obb,
  791. double &dy_obb,
  792. double &angle_obb
  793. )
  794. {
  795. // asign value
  796. Eigen::MatrixXd pts(in_cloud->points.size(), 2);
  797. for (size_t i = 0; i < in_cloud->points.size(); ++i) {
  798. if (axis == 0) {
  799. pts(i, 0) = in_cloud->points.at(i).x;
  800. }
  801. else {
  802. pts(i, 0) = in_cloud->points.at(i).z;
  803. }
  804. pts(i, 1) = in_cloud->points.at(i).y;
  805. }
  806. // centerlize
  807. Eigen::MatrixXd mu = pts.colwise().mean();
  808. Eigen::RowVectorXd mu_row = mu;
  809. pts.rowwise() -= mu_row;
  810. //calculate covariance
  811. Eigen::MatrixXd C = pts.adjoint()*pts;
  812. C = C.array() / (pts.cols() - 1);
  813. //compute eig
  814. Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(C);
  815. Eigen::MatrixXd d = eig.eigenvalues();
  816. Eigen::MatrixXd v = eig.eigenvectors();
  817. //projection
  818. Eigen::MatrixXd p = pts * v;
  819. Eigen::VectorXd minv = p.colwise().minCoeff();
  820. Eigen::VectorXd maxv = p.colwise().maxCoeff();
  821. Eigen::VectorXd range = maxv - minv;
  822. dy_obb = range(1);
  823. dx_obb = range(0);
  824. angle_obb = std::atan2(v(1, 1), v(0, 1));
  825. if (angle_obb < 0) { angle_obb = PI + angle_obb; }
  826. angle_obb = angle_obb * 180 / PI;
  827. }
  828. //void CRootStockGrabPoint::get_optimal_seed(
  829. // pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  830. // pcl::PointXYZ&pt,
  831. // int &pt_idx)
  832. //{
  833. // double d1 = m_cparam.rs_grab_stem_diameter;
  834. // double d2 = m_cparam.rs_grab_stem_diameter * 4.0;
  835. // if (m_dtype != 0) {
  836. // d1 = m_cparam.sc_grab_stem_diameter;
  837. // d2 = m_cparam.sc_grab_stem_diameter * 4.0;
  838. // }
  839. // std::vector<pcl::PointXYZ>cluster_center;
  840. // std::vector<std::vector<int>> cluster_member;
  841. // euclidean_clustering_ttsas(in_cloud, d1, d2, cluster_center, cluster_member);
  842. // float min_y_dist = 1.0e6;
  843. // int opt_idx = -1;
  844. // int member_size = 5;
  845. // float y_opt = m_cparam.rs_grab_y_opt;
  846. // if (m_dtype != 0) {
  847. // y_opt = m_cparam.sc_grab_y_opt;
  848. // }
  849. // for (int i = 0; i < cluster_member.size(); ++i) {
  850. // if (cluster_member.at(i).size() < member_size) {
  851. // continue;
  852. // }
  853. // if (fabs(cluster_center.at(i).y -y_opt) < min_y_dist) {
  854. // min_y_dist = fabs(cluster_center.at(i).y - y_opt);
  855. // opt_idx = i;
  856. // }
  857. // }
  858. // if (opt_idx < 0) {
  859. // if (m_pLogger) {
  860. // stringstream buff;
  861. // buff << "get_optimal_seed failed, get invalid optimal cluster id";
  862. // m_pLogger->ERRORINFO(buff.str());
  863. // }
  864. // return;
  865. // }
  866. // //find nearest pt
  867. // float nearest_dist = 1.0e6;
  868. // int nearest_idx = -1;
  869. // for (int i = 0; i < cluster_member.at(opt_idx).size(); ++i) {
  870. // int idx = cluster_member.at(opt_idx).at(i);
  871. // float dist = fabs(in_cloud->points.at(idx).x - cluster_center.at(opt_idx).x) +
  872. // fabs(in_cloud->points.at(idx).y - cluster_center.at(opt_idx).y) +
  873. // fabs(in_cloud->points.at(idx).z - cluster_center.at(opt_idx).z);
  874. // if (dist < nearest_dist) {
  875. // nearest_dist = dist;
  876. // nearest_idx = idx;
  877. // }
  878. // }
  879. // pt.x = in_cloud->points.at(nearest_idx).x;
  880. // pt.y = in_cloud->points.at(nearest_idx).y;
  881. // pt.z = in_cloud->points.at(nearest_idx).z;
  882. // pt_idx = nearest_idx;
  883. //}
  884. void CRootStockGrabPoint::get_optimal_seed(
  885. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  886. pcl::PointXYZ&pt,
  887. int &pt_idx)
  888. {
  889. /*double d1 = m_cparam.rs_grab_stem_diameter;
  890. double d2 = m_cparam.rs_grab_stem_diameter * 4.0;
  891. if (m_dtype != 0) {
  892. d1 = m_cparam.sc_grab_stem_diameter;
  893. d2 = m_cparam.sc_grab_stem_diameter * 4.0;
  894. }
  895. std::vector<pcl::PointXYZ>cluster_center;
  896. std::vector<std::vector<int>> cluster_member;
  897. euclidean_clustering_ttsas(in_cloud, d1, d2, cluster_center, cluster_member);*/
  898. float min_y_dist = 1.0e6;
  899. int opt_idx = -1;
  900. int member_size = 5;
  901. float y_opt = m_cparam.rs_grab_y_opt;
  902. if (m_dtype == 0) {
  903. y_opt = m_cparam.sc_grab_y_opt;
  904. }
  905. for (int i = 0; i < in_cloud->points.size(); ++i) {
  906. /*if (cluster_member.at(i).size() < member_size) {
  907. continue;
  908. }*/
  909. if (fabs(in_cloud->points.at(i).y - y_opt) < min_y_dist) {
  910. min_y_dist = fabs(in_cloud->points.at(i).y - y_opt);
  911. opt_idx = i;
  912. }
  913. }
  914. if (opt_idx < 0) {
  915. if (m_pLogger) {
  916. stringstream buff;
  917. buff << "get_optimal_seed failed, get invalid optimal cluster id";
  918. m_pLogger->ERRORINFO(buff.str());
  919. }
  920. return;
  921. }
  922. //find nearest pt
  923. /*float nearest_dist = 1.0e6;
  924. int nearest_idx = -1;
  925. for (int i = 0; i < cluster_member.at(opt_idx).size(); ++i) {
  926. int idx = cluster_member.at(opt_idx).at(i);
  927. float dist = fabs(in_cloud->points.at(idx).x - cluster_center.at(opt_idx).x) +
  928. fabs(in_cloud->points.at(idx).y - cluster_center.at(opt_idx).y) +
  929. fabs(in_cloud->points.at(idx).z - cluster_center.at(opt_idx).z);
  930. if (dist < nearest_dist) {
  931. nearest_dist = dist;
  932. nearest_idx = idx;
  933. }
  934. }*/
  935. pt.x = in_cloud->points.at(opt_idx).x;
  936. pt.y = in_cloud->points.at(opt_idx).y;
  937. pt.z = in_cloud->points.at(opt_idx).z;
  938. pt_idx = opt_idx;
  939. }
  940. //通过比较直线点云和原始点云相同位置邻域内xz的范围,确定此点是否是无干扰点的茎
  941. void CRootStockGrabPoint::get_optimal_seed_compare(
  942. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, // 全部有效点云
  943. pcl::PointCloud<pcl::PointXYZ>::Ptr in_line_cloud, // 茎上直线点云
  944. pcl::PointXYZ&pt, // 返回抓取的点坐标
  945. int &pt_idx, // 返回点index
  946. std::vector<int>& valid_line_index // 返回有效直线点index
  947. )
  948. {
  949. valid_line_index.clear();
  950. float th = 0.75; //原始点云和直线点云邻域尺寸增加原来的20%,不能当做抓取点
  951. pt_idx = -1;
  952. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  953. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox_line(new pcl::PointCloud<pcl::PointXYZ>);
  954. pcl::CropBox<pcl::PointXYZ> box_filter;
  955. box_filter.setNegative(false);
  956. box_filter.setInputCloud(in_cloud);
  957. pcl::CropBox<pcl::PointXYZ> box_filter_line;
  958. box_filter_line.setNegative(false);
  959. box_filter_line.setInputCloud(in_line_cloud);
  960. float radius = m_cparam.rs_grab_stem_diameter;
  961. float opt_y = m_cparam.rs_grab_y_opt;
  962. if (m_dtype == 0) {
  963. radius = m_cparam.sc_grab_stem_diameter;
  964. opt_y = m_cparam.sc_grab_y_opt;
  965. }
  966. float rx = 1.5;
  967. float ry = 1.5;
  968. float rz = 1.5;
  969. float cx, cy, cz;
  970. float dz,dx, dz_line, dx_line;
  971. float dist_min = 10000.0;
  972. for (size_t i = 0; i < in_line_cloud->points.size(); ++i) {
  973. cx = in_line_cloud->points.at(i).x;
  974. cy = in_line_cloud->points.at(i).y;
  975. cz = in_line_cloud->points.at(i).z;
  976. //////////////////////////////////////////////////////////////////
  977. //原始点云,获取指定区域的dx,dz
  978. box_filter.setMin(Eigen::Vector4f(cx - rx*radius, cy - ry*radius, cz - rz*radius, 1));
  979. box_filter.setMax(Eigen::Vector4f(cx + rx*radius, cy + ry*radius, cz + rz*radius, 1));
  980. box_filter.filter(*cloud_inbox);
  981. Eigen::Vector4f min_point;
  982. Eigen::Vector4f max_point;
  983. pcl::getMinMax3D(*cloud_inbox, min_point, max_point);
  984. dz = max_point(2) - min_point(2);
  985. dx = max_point(0) - min_point(0);
  986. /////////////////////////////////////////////////////////////////////
  987. //直线点云,获取指定区域的dx_line,dz_line
  988. box_filter_line.setMin(Eigen::Vector4f(cx - rx*radius, cy - ry*radius, cz - rz*radius, 1));
  989. box_filter_line.setMax(Eigen::Vector4f(cx + rx*radius, cy + ry*radius, cz + rz*radius, 1));
  990. box_filter_line.filter(*cloud_inbox_line);
  991. pcl::getMinMax3D(*cloud_inbox_line, min_point, max_point);
  992. dz_line = max_point(2) - min_point(2);
  993. dx_line = max_point(0) - min_point(0);
  994. float ratio_ex = (dx + dz - dz_line - dx_line) / (dz_line + dx_line);
  995. if (ratio_ex < th) {
  996. valid_line_index.push_back(i);
  997. if (fabs(cy - opt_y) < dist_min) {
  998. dist_min = fabs(cy - opt_y);
  999. pt_idx = i;
  1000. }
  1001. }
  1002. }
  1003. if (pt_idx >= 0) {
  1004. pt = in_line_cloud->points.at(pt_idx);
  1005. }
  1006. }
  1007. void CRootStockGrabPoint::viewer_cloud(pcl::PointCloud<pcl::PointXYZ>::Ptr cloud, std::string&winname)
  1008. {
  1009. pcl::visualization::CloudViewer viewer(winname);
  1010. //viewer.runOnVisualizationThreadOnce(viewerOneOff);
  1011. viewer.showCloud(cloud);
  1012. while (!viewer.wasStopped()) {
  1013. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  1014. }
  1015. }
  1016. void CRootStockGrabPoint::viewer_cloud(pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud, std::string&winname)
  1017. {
  1018. pcl::visualization::CloudViewer viewer(winname);
  1019. //viewer.runOnVisualizationThreadOnce(viewerOneOff);
  1020. viewer.showCloud(cloud);
  1021. while (!viewer.wasStopped()) {
  1022. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  1023. }
  1024. }
  1025. void CRootStockGrabPoint::viewer_cloud_debug(
  1026. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud,
  1027. pcl::PointXYZ &p,
  1028. pcl::PointXYZ &root_pt,
  1029. std::string&winname)
  1030. {
  1031. pcl::visualization::PCLVisualizer viewer(winname);
  1032. //viewer.runOnVisualizationThreadOnce(viewerOneOff);
  1033. viewer.addPointCloud(cloud);
  1034. viewer.addCoordinateSystem();
  1035. pcl::PointXYZ p0, x1, y1,p00,x0,y0, root_px0, root_px1, root_py0,root_py1;
  1036. p0.x = p.x;
  1037. p0.y = p.y;
  1038. p0.z = p.z;
  1039. x1.x = p.x + 100.0;
  1040. x1.y = p.y;
  1041. x1.z = p.z;
  1042. y1.x = p.x;
  1043. y1.y = p.y + 20.0;
  1044. y1.z = p.z;
  1045. p00.x = 0.0;
  1046. p00.y = 0.0;
  1047. p00.z = p.z;
  1048. x0.x = 600.0;
  1049. x0.y = 0;
  1050. x0.z = p.z;
  1051. y0.x = 0.0;
  1052. y0.y = 300.0;
  1053. y0.z = p.z;
  1054. root_px0.x = root_pt.x - 5.0;
  1055. root_px0.y = root_pt.y;
  1056. root_px0.z = root_pt.z;
  1057. root_px1.x = root_pt.x + 5.0;
  1058. root_px1.y = root_pt.y;
  1059. root_px1.z = root_pt.z;
  1060. root_py0.x = root_pt.x;
  1061. root_py0.y = root_pt.y - 5.0;
  1062. root_py0.z = root_pt.z;
  1063. root_py1.x = root_pt.x;
  1064. root_py1.y = root_pt.y + 5.0;
  1065. root_py1.z = root_pt.z;
  1066. viewer.addLine(p0, x1, 255, 0, 0, "x");
  1067. viewer.addLine(p0, y1, 0, 255, 0, "y");
  1068. viewer.addLine(p00, x0, 255, 0, 0, "x0");
  1069. viewer.addLine(p00, y0, 0, 255, 0, "y0");
  1070. viewer.addLine(root_px0, root_px1, 255, 0, 0, "rootx");
  1071. viewer.addLine(root_py0, root_py1, 0, 255, 0, "rooty");
  1072. while (!viewer.wasStopped()) {
  1073. viewer.spinOnce(100);
  1074. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  1075. }
  1076. }
  1077. void CRootStockGrabPoint::viewer_cloud_cluster(
  1078. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud,
  1079. std::vector<pcl::PointXYZ>cluster_center,
  1080. std::string&winname)
  1081. {
  1082. pcl::visualization::PCLVisualizer viewer(winname);
  1083. viewer.addPointCloud(cloud);
  1084. viewer.addCoordinateSystem();
  1085. int cnt = 0;
  1086. char buf[8];
  1087. for (auto& p : cluster_center) {
  1088. pcl::PointXYZ px0, px1, py1, py0;
  1089. px0.x = p.x-5.0;
  1090. px0.y = p.y;
  1091. px0.z = p.z;
  1092. px1.x = p.x + 5.0;
  1093. px1.y = p.y;
  1094. px1.z = p.z;
  1095. py0.x = p.x;
  1096. py0.y = p.y - 5.0;
  1097. py0.z = p.z;
  1098. py1.x = p.x;
  1099. py1.y = p.y + 5.0;
  1100. py1.z = p.z;
  1101. viewer.addLine(px0, px1, 255, 0, 0, "x" + string(_itoa(cnt, buf,10)));
  1102. viewer.addLine(py0, py1, 0, 255, 0, "y" + string(_itoa(cnt, buf, 10)));
  1103. cnt += 1;
  1104. }
  1105. while (!viewer.wasStopped()) {
  1106. viewer.spinOnce(100);
  1107. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  1108. }
  1109. }
  1110. void CRootStockGrabPoint::viewer_cloud_cluster_box(
  1111. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud,
  1112. std::vector<pcl::PointXYZ>&cluster_center,
  1113. std::vector<pcl::PointXYZ>&cluster_box,
  1114. std::vector<std::vector<int> >& clt_line_idx,
  1115. std::string&winname)
  1116. {
  1117. pcl::visualization::PCLVisualizer viewer(winname);
  1118. viewer.addCoordinateSystem();
  1119. viewer.addPointCloud<pcl::PointXYZ>(cloud, "all_cloud");
  1120. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "all_cloud");
  1121. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "all_cloud");
  1122. int cnt = 0;
  1123. char buf[8];
  1124. for (size_t i = 0; i < cluster_center.size();++i) {
  1125. pcl::PointXYZ &p = cluster_center.at(i);
  1126. pcl::PointXYZ px0, px1, py1, py0;
  1127. px0.x = p.x - 5.0;
  1128. px0.y = p.y;
  1129. px0.z = p.z;
  1130. px1.x = p.x + 5.0;
  1131. px1.y = p.y;
  1132. px1.z = p.z;
  1133. py0.x = p.x;
  1134. py0.y = p.y - 5.0;
  1135. py0.z = p.z;
  1136. py1.x = p.x;
  1137. py1.y = p.y + 5.0;
  1138. py1.z = p.z;
  1139. viewer.addLine(px0, px1, 255, 0, 0, "x" + string(_itoa(cnt, buf, 10)));
  1140. viewer.addLine(py0, py1, 0, 255, 0, "y" + string(_itoa(cnt, buf, 10)));
  1141. //box
  1142. pcl::PointXYZ & min_pt = cluster_box.at(2 * i);
  1143. pcl::PointXYZ & max_pt = cluster_box.at(2 * i + 1);
  1144. viewer.addCube(min_pt.x, max_pt.x, min_pt.y, max_pt.y, min_pt.z, max_pt.z, 0.5,0.5,0.0,"AABB_"+string(_itoa(cnt, buf, 10)));
  1145. viewer.setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_REPRESENTATION,
  1146. pcl::visualization::PCL_VISUALIZER_REPRESENTATION_WIREFRAME, "AABB_" + string(_itoa(cnt, buf, 10)));
  1147. pcl::PointCloud<pcl::PointXYZ>::Ptr line_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  1148. pcl::copyPointCloud(*cloud, clt_line_idx.at(i), *line_cloud);
  1149. viewer.addPointCloud(line_cloud, "fit_line" + string(_itoa(cnt, buf, 10)));
  1150. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0, 1, 0, "fit_line" + string(_itoa(cnt, buf, 10)));
  1151. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, "fit_line" + string(_itoa(cnt, buf, 10)));
  1152. cnt += 1;
  1153. }
  1154. while (!viewer.wasStopped()) {
  1155. viewer.spinOnce(100);
  1156. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  1157. }
  1158. }
  1159. };