grab_point_rs.cpp 93 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836
  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 <pcl\segmentation\extract_clusters.h>
  20. #include <math.h>
  21. #include <pcl\features\normal_3d.h>
  22. #include <pcl\features\boundary.h>
  23. #include "grab_point_rs.h"
  24. #include "utils.h"
  25. #include "peak_finder.h"
  26. #define PI std::acos(-1)
  27. namespace graft_cv {
  28. CRootStockGrabPoint::CRootStockGrabPoint(ConfigParam&cp, CGcvLogger*pLog/*=0*/)
  29. :m_cparam(cp),
  30. m_pLogger(pLog),
  31. m_dtype(0),
  32. m_pcdId(""),
  33. m_ppImgSaver(0),
  34. m_1st_row_zmean_rs(-1.0),
  35. m_1st_row_zmean_sc(-1.0),
  36. m_cloud_mean_dist(0.0),
  37. m_pImginfoResult(0),
  38. m_pStemInfos(0)
  39. {
  40. }
  41. void CRootStockGrabPoint::clear_imginfo() {
  42. if (m_pImginfoResult) {
  43. imginfo_release(&m_pImginfoResult);
  44. m_pImginfoResult = 0;
  45. }
  46. }
  47. CRootStockGrabPoint::~CRootStockGrabPoint() {
  48. this->clear_imginfo();
  49. if (m_pStemInfos) {
  50. delete m_pStemInfos;
  51. m_pStemInfos = 0;
  52. }
  53. }
  54. float* CRootStockGrabPoint::get_raw_point_cloud(int &data_size)
  55. {
  56. data_size = m_raw_cloud->width * m_raw_cloud->height;
  57. if (data_size == 0) {
  58. return 0;
  59. }
  60. else {
  61. pcl::PointXYZ* pp = m_raw_cloud->points.data();
  62. return (float*)pp->data;
  63. }
  64. }
  65. int CRootStockGrabPoint::load_data(
  66. float*pPoint,
  67. int pixel_size,
  68. int pt_size,
  69. int dtype,
  70. const char* fn/* = 0*/)
  71. {
  72. //数据加载功能实现,并生成imageid,保存原始数据到文件
  73. int rst = 0;
  74. m_dtype = dtype;
  75. //generate image id
  76. if (m_dtype == 0) {
  77. m_pcdId = getImgId(img_type::sola_sc_pcd);
  78. }
  79. else {
  80. m_pcdId = getImgId(img_type::sola_rs_pcd);
  81. }
  82. //1 get point cloud data
  83. if (pPoint != 0 && pt_size>0) {
  84. //read point
  85. m_raw_cloud.reset(new pcl::PointCloud<pcl::PointXYZ>);
  86. int step = pixel_size;
  87. for (int i = 0; i < pt_size; ++i) {
  88. pcl::PointXYZ pt = { pPoint[i*step], pPoint[i*step + 1] , pPoint[i*step + 2] };
  89. m_raw_cloud->push_back(pt);
  90. }
  91. rst = m_raw_cloud->width * m_raw_cloud->height;
  92. if (m_pLogger) {
  93. stringstream buff;
  94. buff << m_pcdId <<": load raw point cloud " << rst << " data points";
  95. m_pLogger->INFO(buff.str());
  96. }
  97. }
  98. else if (fn != 0) {
  99. //read file
  100. rst = this->read_ply_file(fn);
  101. }
  102. else {//error
  103. if (m_pLogger) {
  104. m_pLogger->ERRORINFO(m_pcdId + ": no valid input");
  105. return (-1);
  106. }
  107. }
  108. if (m_ppImgSaver && *m_ppImgSaver) {
  109. (*m_ppImgSaver)->saveBinPly(m_raw_cloud, m_pcdId);
  110. }
  111. if (m_cparam.image_show) {
  112. pcl::visualization::PCLVisualizer viewer(m_pcdId + std::string(": raw point cloud"));
  113. viewer.setBackgroundColor(0.35, 0.35, 0.35);
  114. viewer.addCoordinateSystem();
  115. viewer.addPointCloud<pcl::PointXYZ>(m_raw_cloud, "raw_cloud");
  116. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "raw_cloud");
  117. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "raw_cloud");
  118. float xmin, ymin, zmin, xmax, ymax, zmax;
  119. xmin = m_cparam.rs_grab_xmin;
  120. ymin = m_cparam.rs_grab_ymin;
  121. zmin = m_cparam.rs_grab_zmin;
  122. xmax = m_cparam.rs_grab_xmax;
  123. ymax = m_cparam.rs_grab_ymax;
  124. zmax = m_cparam.rs_grab_zmax;
  125. if (m_dtype == 0) {
  126. xmin = m_cparam.sc_grab_xmin;
  127. ymin = m_cparam.sc_grab_ymin;
  128. zmin = m_cparam.sc_grab_zmin;
  129. xmax = m_cparam.sc_grab_xmax;
  130. ymax = m_cparam.sc_grab_ymax;
  131. zmax = m_cparam.sc_grab_zmax;
  132. }
  133. viewer.addCube(xmin, xmax,ymin, ymax, zmin, zmax, 0.75, 0.0, 0.0, "AABB_");
  134. viewer.setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_REPRESENTATION,
  135. pcl::visualization::PCL_VISUALIZER_REPRESENTATION_WIREFRAME, "AABB_");
  136. pcl::PointXYZ px0, px1, py1, pz1;
  137. px0.x = 0;
  138. px0.y = 0;
  139. px0.z = 0;
  140. px1.x = 10.0;
  141. px1.y = 0;
  142. px1.z = 0;
  143. py1.x = 0;
  144. py1.y = 10.0;
  145. py1.z = 0;
  146. pz1.x = 0;
  147. pz1.y = 0;
  148. pz1.z = 10.0;
  149. viewer.addLine(px0, px1, 255, 0, 0, "x");
  150. viewer.addLine(px0, py1, 0, 255, 0, "y");
  151. viewer.addLine(px0, pz1, 0, 0, 255, "z");
  152. while (!viewer.wasStopped()) {
  153. viewer.spinOnce(100);
  154. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  155. }
  156. }
  157. if (m_pStemInfos==0) {
  158. double seedling_distance = m_cparam.rs_grab_seedling_dist;
  159. int holes_number = m_cparam.rs_grab_holes_number;
  160. double x_min = m_cparam.rs_grab_xmin;
  161. double x_max = m_cparam.rs_grab_xmax;
  162. double z_min = m_cparam.rs_grab_zmin;
  163. double z_max = m_cparam.rs_grab_zmax;
  164. if (m_dtype == 0) {
  165. seedling_distance = m_cparam.sc_grab_seedling_dist;
  166. holes_number = m_cparam.sc_grab_holes_number;
  167. x_min = m_cparam.sc_grab_xmin;
  168. x_max = m_cparam.sc_grab_xmax;
  169. z_min = m_cparam.sc_grab_zmin;
  170. z_max = m_cparam.sc_grab_zmax;
  171. }
  172. m_pStemInfos = new CStemResultInfos(
  173. seedling_distance, holes_number,
  174. x_min, x_max, z_min, z_max,
  175. m_pcdId, m_pLogger);
  176. }
  177. return rst;
  178. }
  179. int CRootStockGrabPoint::grab_point_detect(
  180. PositionInfo& posinfo
  181. )
  182. {
  183. // 抓取位置识别入口函数,实现整个抓取位置识别功能,返回抓取位置信息
  184. // dtype == 0, scion
  185. // dtype != 0, rootstock
  186. // 输入,穗苗--0, 砧木--1
  187. if (m_raw_cloud->width * m_raw_cloud->height < 1) {
  188. if (m_pLogger) {
  189. stringstream buff;
  190. buff << m_pcdId << ": m_raw_cloud point cloud " << m_raw_cloud->width * m_raw_cloud->height << " data points";
  191. m_pLogger->ERRORINFO(buff.str());
  192. }
  193. return 1;
  194. }
  195. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  196. //1 crop filter,剪裁需要的部分点云
  197. if (m_pLogger) {
  198. m_pLogger->INFO(m_pcdId + ": begin crop");
  199. }
  200. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  201. pcl::CropBox<pcl::PointXYZ> box_filter;
  202. if (m_dtype != 0) {//rootstock
  203. box_filter.setMin(Eigen::Vector4f(m_cparam.rs_grab_xmin, m_cparam.rs_grab_ymin, m_cparam.rs_grab_zmin, 1));
  204. box_filter.setMax(Eigen::Vector4f(m_cparam.rs_grab_xmax, m_cparam.rs_grab_ymax, m_cparam.rs_grab_zmax, 1));
  205. }
  206. else {//scion
  207. box_filter.setMin(Eigen::Vector4f(m_cparam.sc_grab_xmin, m_cparam.sc_grab_ymin, m_cparam.sc_grab_zmin, 1));
  208. box_filter.setMax(Eigen::Vector4f(m_cparam.sc_grab_xmax, m_cparam.sc_grab_ymax, m_cparam.sc_grab_zmax, 1));
  209. }
  210. box_filter.setNegative(false);
  211. box_filter.setInputCloud(m_raw_cloud);
  212. box_filter.filter(*cloud_inbox);
  213. if (m_pLogger) {
  214. stringstream buff;
  215. buff << m_pcdId << ": CropBox croped point cloud " << cloud_inbox->width * cloud_inbox->height << " data points";
  216. m_pLogger->INFO(buff.str());
  217. }
  218. if (cloud_inbox->width * cloud_inbox->height < 1) {
  219. return 1;
  220. }
  221. if (m_cparam.image_show) {
  222. viewer_cloud(cloud_inbox, std::string("cloud_inbox"));
  223. }
  224. if (m_pLogger) {
  225. m_pLogger->INFO(m_pcdId + ": end crop");
  226. }
  227. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  228. //2 filtler with radius remove, 去除孤立点
  229. if (m_pLogger) {
  230. m_pLogger->INFO(m_pcdId + ": begin ror");
  231. }
  232. int nb_points = 20;
  233. double stem_radius = m_cparam.rs_grab_stem_diameter / 2.0;
  234. if(m_dtype == 0){
  235. stem_radius = m_cparam.sc_grab_stem_diameter / 2.0;
  236. }
  237. double remove_radius = stem_radius * 2.0;
  238. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_ror(new pcl::PointCloud<pcl::PointXYZ>);
  239. pcl::RadiusOutlierRemoval<pcl::PointXYZ> ror(false);
  240. ror.setInputCloud(cloud_inbox);
  241. ror.setRadiusSearch(remove_radius);
  242. ror.setMinNeighborsInRadius(nb_points);
  243. ror.filter(*cloud_ror);
  244. if (m_pLogger) {
  245. stringstream buff;
  246. buff << m_pcdId <<": RadiusOutlierRemoval filtered point cloud "
  247. << cloud_ror->width * cloud_ror->height
  248. << " data points. param stem_radius="<<
  249. stem_radius<<", nb_points="<< nb_points<< "(if < 50, return)";
  250. m_pLogger->INFO(buff.str());
  251. }
  252. if (cloud_ror->width * cloud_ror->height < 50) {
  253. return 1;
  254. }
  255. if (m_cparam.image_show) {
  256. viewer_cloud(cloud_ror, std::string("cloud_ror"));
  257. }
  258. if (m_pLogger) {
  259. m_pLogger->INFO(m_pcdId + ": end ror");
  260. }
  261. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  262. //?? 用截取到的点云,初步找到植株的位置中心
  263. m_root_centers.clear();
  264. m_root_center_with_seedling.clear();
  265. if (m_pStemInfos) {
  266. m_pStemInfos->get_root_centers(m_root_centers);
  267. if (m_pLogger) {
  268. stringstream buff;
  269. buff << m_pcdId << ": root positions "<< m_root_centers.size()<<" centers\n";
  270. for (auto& sr : m_root_centers) {
  271. buff << "\tx=" << sr.root_x << ", y=" << sr.root_y << ", z=" << sr.root_z
  272. << ", cnt=" << sr.root_count
  273. << ", pcd_size=" << sr.stem_size << "\n";
  274. }
  275. m_pLogger->INFO(buff.str());
  276. }
  277. }
  278. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  279. //3 原来的降采样没有意义,改成统计平均距离
  280. m_cloud_mean_dist = 0.0;
  281. cloud_mean_dist(cloud_ror, m_cloud_mean_dist);
  282. if (m_pLogger) {
  283. stringstream buff;
  284. buff << m_pcdId <<": cloud_mean_dist = " << m_cloud_mean_dist;
  285. m_pLogger->INFO(buff.str());
  286. }
  287. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  288. // 4 对截取的点云进行ror滤除大面积联通区域,剔除叶片
  289. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_dowm_sampled(new pcl::PointCloud<pcl::PointXYZ>);
  290. std::vector<int> non_leaf_indices;
  291. std::vector<int> leaf_indices;
  292. leaf_filter_ror(cloud_ror, non_leaf_indices, leaf_indices);
  293. pcl::copyPointCloud(*cloud_ror, non_leaf_indices, *cloud_dowm_sampled);
  294. if (m_pLogger) {
  295. stringstream buff;
  296. buff << m_pcdId << ": after leaf_filter dowm_sampled point cloud "
  297. << cloud_dowm_sampled->width * cloud_dowm_sampled->height << " data points (if < 50, return)";
  298. m_pLogger->INFO(buff.str());
  299. }
  300. if (cloud_dowm_sampled->width * cloud_dowm_sampled->height < 50) {
  301. return 1;
  302. }
  303. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  304. //判断m_root_centers位置上是否有叶片遮挡
  305. occluded_seedling_detect_by_leaf(cloud_ror, leaf_indices);
  306. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  307. //5 seedling position,找到第一个位置上的植株
  308. std::vector<int> first_seedling_cloud_idx; //存储找到的第一个位置上植株茎上直线点的index
  309. pcl::PointXYZ xz_center; //存储找到的第一个位置上植株根部坐标
  310. pcl::ModelCoefficients::Ptr line_model; //存储找到的第一个位置上植株茎直线模型
  311. if (m_pLogger) {
  312. m_pLogger->INFO(m_pcdId + ": begin find seedling position");
  313. }
  314. int first_row_seedling_number = 0;
  315. bool fund_seedling = find_seedling_position_key(cloud_dowm_sampled, first_seedling_cloud_idx, xz_center, line_model, first_row_seedling_number);
  316. if (!fund_seedling && first_row_seedling_number == 0) {
  317. if (m_pLogger) {
  318. stringstream buff;
  319. buff << m_pcdId << ": Not found seedlings" ;
  320. m_pLogger->INFO(buff.str());
  321. }
  322. return 1;
  323. }
  324. if (m_pLogger) {
  325. stringstream buff;
  326. buff << m_pcdId <<": after find_seedling_position(), foud first seedling seeds points size " << first_seedling_cloud_idx.size();
  327. m_pLogger->INFO(buff.str());
  328. m_pLogger->INFO(m_pcdId + ": end find seedling position");
  329. }
  330. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  331. //5.1 如果没有找到茎,在root_center 中找一个位置
  332. pcl::PointXYZ selected_pt;//抓取点坐标,根据茎节点的偏移
  333. pcl::PointXYZ selected_pt_ref; //返回茎节点,作为抓取点偏的基点
  334. if (!fund_seedling) {
  335. int selected_center_idx;
  336. no_seedling_detected_post_process(
  337. first_row_seedling_number, //input
  338. selected_center_idx, //output
  339. selected_pt, //output
  340. selected_pt_ref, //output
  341. posinfo);
  342. if (selected_center_idx < 0) {
  343. if (m_pLogger) {
  344. stringstream buff;
  345. buff << m_pcdId << ": Not found seedlings after no_seedling_detected_post_process()";
  346. m_pLogger->INFO(buff.str());
  347. }
  348. return 1;
  349. }
  350. }
  351. else {
  352. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  353. //6 nearest seedling grab point selection,找到最接近指定高度的最优抓取点坐标
  354. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_seedling_seed(new pcl::PointCloud<pcl::PointXYZ>);
  355. pcl::copyPointCloud(*cloud_dowm_sampled, first_seedling_cloud_idx, *cloud_seedling_seed);
  356. //改进思路:将茎直线上cloud_seedling_seed的点,提取周围邻域xz的宽度,在相同邻域cloud_dowm_sampled点云内提取xz的宽度
  357. //对比两个宽度,变化很大的,说明周边有异物(叶柄或叶子),不适合作为抓取位置;否则就是抓取位置
  358. //在众多抓取位置上,选择离指定高度最近的点作为抓取位置
  359. //
  360. //显示位置
  361. if (m_cparam.image_show) {
  362. std::vector<pcl::PointXYZ> tmp;
  363. tmp.push_back(xz_center);
  364. pcl::PointCloud<pcl::PointXYZRGB>::Ptr tttp(new pcl::PointCloud<pcl::PointXYZRGB>);
  365. pcl::copyPointCloud(*cloud_seedling_seed, *tttp);
  366. for (auto& pt : *tttp) {
  367. pt.r = 255;
  368. pt.g = 255;
  369. pt.b = 255;
  370. }
  371. viewer_cloud_cluster(tttp, tmp, string("first"));
  372. }
  373. float stem_width_mu;
  374. float stem_deflection;
  375. get_optimal_seed_compare(
  376. cloud_dowm_sampled, //input
  377. cloud_seedling_seed, //input
  378. line_model, //input
  379. selected_pt, //output
  380. selected_pt_ref, //output
  381. stem_width_mu, //output
  382. stem_deflection //output
  383. ); //output
  384. if (selected_pt.z < 0) {
  385. int selected_center_idx;
  386. no_seedling_detected_post_process(
  387. first_row_seedling_number, //input
  388. selected_center_idx, //output
  389. selected_pt, //output
  390. selected_pt_ref, //output
  391. posinfo);
  392. if (selected_center_idx < 0) {
  393. if (m_pLogger) {
  394. stringstream buff;
  395. buff << m_pcdId << ": Not found valid fork point after no_seedling_detected_post_process()";
  396. m_pLogger->INFO(buff.str());
  397. }
  398. return 1;
  399. }
  400. }
  401. else {
  402. int selected_center_idx;
  403. had_seedling_detected_post_process(
  404. first_row_seedling_number, //input
  405. stem_width_mu, //input
  406. stem_deflection, //input
  407. selected_center_idx, //output, 选择root_center的序号
  408. selected_pt, //input-output, 检测到的目标抓取点
  409. selected_pt_ref, //input-output, 检测到的目标抓取参考点
  410. posinfo);
  411. }
  412. if (m_pLogger) {
  413. m_pLogger->INFO(m_pcdId + ": end get_optimal_seed");
  414. }
  415. }
  416. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  417. //7 保存处理结果到图片
  418. cv::Mat rst_img = cv::Mat::zeros(cv::Size(1280,640),CV_8UC3);
  419. gen_result_img(cloud_ror, cloud_dowm_sampled, selected_pt, selected_pt_ref, posinfo,rst_img);
  420. if (m_ppImgSaver && *m_ppImgSaver) {
  421. (*m_ppImgSaver)->saveImage(rst_img, m_pcdId + "_rst_0");
  422. }
  423. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  424. //8 页面显示结果
  425. this->clear_imginfo();
  426. m_pImginfoResult = mat2imginfo(rst_img);
  427. posinfo.pp_images[0] = m_pImginfoResult;
  428. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  429. //9 debug 显示结果
  430. if (m_cparam.image_show) {
  431. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_cand_demo(new pcl::PointCloud<pcl::PointXYZRGB>);
  432. pcl::copyPointCloud(*cloud_dowm_sampled, *cloud_cand_demo);
  433. for (auto& pt : *cloud_cand_demo) {
  434. pt.r = 255;
  435. pt.g = 255;
  436. pt.b = 255;
  437. }
  438. pcl::PointXYZRGB pt_grab = {0,255,0};
  439. pt_grab.x = selected_pt.x;
  440. pt_grab.y = selected_pt.y;
  441. pt_grab.z = selected_pt.z;
  442. pcl::PointXYZRGB pt_grab_ = { 0,255,120 };
  443. pt_grab_.x = selected_pt.x;
  444. pt_grab_.y = selected_pt.y+0.2;
  445. pt_grab_.z = selected_pt.z;
  446. cloud_cand_demo->push_back(pt_grab);
  447. //viewer_cloud(cloud_cand_demo, std::string("cloud_cand_demo"));
  448. viewer_cloud_debug(cloud_cand_demo, selected_pt, selected_pt_ref, xz_center, std::string("cloud_cand_demo"));
  449. imshow_wait("rst_img", rst_img);
  450. }
  451. return 0;
  452. }
  453. //没有检测到苗的后处理
  454. void CRootStockGrabPoint::no_seedling_detected_post_process(
  455. int first_row_seedling_number, //input
  456. int& selected_idx, //output, 选择root_center的序号
  457. pcl::PointXYZ& selected_pt, //output
  458. pcl::PointXYZ& selected_pt_ref, //output
  459. PositionInfo& posinfo //output
  460. )
  461. {
  462. //m_dtype == 0 找最大x位置, 否则找最小x位置
  463. selected_idx = -1;
  464. if (m_dtype == 0) {
  465. for (int i = 0; i < m_root_center_with_seedling.size(); ++i) {
  466. if (m_root_center_with_seedling[i]) { selected_idx = i; }
  467. }
  468. }
  469. else {
  470. for (int i = 0; i < m_root_center_with_seedling.size(); ++i) {
  471. if (m_root_center_with_seedling[i]) {
  472. selected_idx = i;
  473. break;
  474. }
  475. }
  476. }
  477. if (selected_idx < 0) {
  478. return;
  479. }
  480. double grab_fork_ybt = m_cparam.rs_grab_fork_ybt;
  481. double grab_offset = m_cparam.rs_grab_offset;
  482. if (m_dtype == 0) {
  483. grab_fork_ybt = m_cparam.sc_grab_fork_ybt;
  484. grab_offset = m_cparam.sc_grab_offset;
  485. }
  486. selected_pt_ref.x = m_root_centers.at(selected_idx).root_x;
  487. selected_pt_ref.y = grab_fork_ybt;
  488. selected_pt_ref.z = m_root_centers.at(selected_idx).root_z;
  489. selected_pt = selected_pt_ref;
  490. selected_pt.y += static_cast<int>(grab_offset);
  491. double stem_width_mu = 0.0;
  492. double stem_deflection = 0.0;
  493. if (m_dtype == 0) {
  494. posinfo.sc_grab_x = selected_pt.x;
  495. posinfo.sc_grab_y = selected_pt.y;
  496. posinfo.sc_grab_z = selected_pt.z;
  497. posinfo.sc_count = (double)first_row_seedling_number;
  498. posinfo.sc_width = (double)stem_width_mu;
  499. posinfo.sc_tortuosity = (double)stem_deflection;
  500. }
  501. else {
  502. posinfo.rs_grab_x = selected_pt.x;
  503. posinfo.rs_grab_y = selected_pt.y;
  504. posinfo.rs_grab_z = selected_pt.z;
  505. posinfo.rs_count = (double)first_row_seedling_number;
  506. posinfo.rs_width = (double)stem_width_mu;
  507. posinfo.rs_tortuosity = (double)stem_deflection;
  508. }
  509. }
  510. //检测到苗的情况,后处理
  511. void CRootStockGrabPoint::had_seedling_detected_post_process(
  512. int first_row_seedling_number, //input
  513. float stem_width_mu, //input
  514. float stem_deflection, //input
  515. int& selected_idx, //output, 选择root_center的序号
  516. pcl::PointXYZ& selected_pt, //input-output, 检测到的目标抓取点
  517. pcl::PointXYZ& selected_pt_ref, //input-output, 检测到的目标抓取参考点
  518. PositionInfo& posinfo //output
  519. )
  520. {
  521. //通过selected_pt找到对应的gen中心位置
  522. //0 首选识别到的位置
  523. if (m_dtype == 0) {
  524. posinfo.sc_grab_x = selected_pt.x;
  525. posinfo.sc_grab_y = selected_pt.y;
  526. posinfo.sc_grab_z = selected_pt.z;
  527. posinfo.sc_count = (double)first_row_seedling_number;
  528. posinfo.sc_width = (double)stem_width_mu;
  529. posinfo.sc_tortuosity = (double)stem_deflection;
  530. }
  531. else {
  532. posinfo.rs_grab_x = selected_pt.x;
  533. posinfo.rs_grab_y = selected_pt.y;
  534. posinfo.rs_grab_z = selected_pt.z;
  535. posinfo.rs_count = (double)first_row_seedling_number;
  536. posinfo.rs_width = (double)stem_width_mu;
  537. posinfo.rs_tortuosity = (double)stem_deflection;
  538. }
  539. // 1 找到基于rootcenter应该找的中心位置
  540. //m_dtype == 0 找最大x位置, 否则找最小x位置
  541. selected_idx = -1;
  542. if (m_dtype == 0) {
  543. for (int i = 0; i < m_root_center_with_seedling.size(); ++i) {
  544. if (m_root_center_with_seedling[i]) { selected_idx = i; }
  545. }
  546. }
  547. else {
  548. for (int i = 0; i < m_root_center_with_seedling.size(); ++i) {
  549. if (m_root_center_with_seedling[i]) {
  550. selected_idx = i;
  551. break;
  552. }
  553. }
  554. }
  555. if (selected_idx < 0) {
  556. //可能没有中心(刚开始工作)
  557. //按识别的结果
  558. return;
  559. }
  560. //2 如果识别的位置和rootcenter位置接近,就按识别位置
  561. double seedling_distance = m_cparam.rs_grab_seedling_dist;
  562. if (m_dtype == 0) {
  563. seedling_distance = m_cparam.sc_grab_seedling_dist;
  564. }
  565. double root_x = m_root_centers.at(selected_idx).root_x;
  566. if (m_dtype == 0) {
  567. //穗苗找x最大的位置
  568. if (selected_pt.x > (root_x - 0.5 * seedling_distance)) {
  569. //找到x最大的位置上的苗,识别的苗位置大于有苗的位置,就认为识别结果更可信
  570. return;
  571. }
  572. else {
  573. //否则用指定根中心的位置
  574. goto obstructed;
  575. }
  576. }
  577. else {
  578. //砧木,找x最小的位置
  579. if (selected_pt.x < (root_x + 0.5 * seedling_distance)) {
  580. //找到x最大的位置上的苗,识别的苗位置大于有苗的位置,就认为识别结果更可信
  581. return;
  582. }
  583. else {
  584. //否则用指定根中心的位置
  585. goto obstructed;
  586. }
  587. }
  588. obstructed:
  589. double grab_fork_ybt = m_cparam.rs_grab_fork_ybt;
  590. double grab_offset = m_cparam.rs_grab_offset;
  591. if (m_dtype == 0) {
  592. grab_fork_ybt = m_cparam.sc_grab_fork_ybt;
  593. grab_offset = m_cparam.sc_grab_offset;
  594. }
  595. selected_pt_ref.x = m_root_centers.at(selected_idx).root_x;
  596. selected_pt_ref.y = grab_fork_ybt;
  597. selected_pt_ref.z = m_root_centers.at(selected_idx).root_z;
  598. selected_pt = selected_pt_ref;
  599. selected_pt.y += static_cast<int>(grab_offset);
  600. double stem_width_mu_obstructed = 0.0;
  601. double stem_deflection_obstructed = 0.0;
  602. if (m_dtype == 0) {
  603. posinfo.sc_grab_x = selected_pt.x;
  604. posinfo.sc_grab_y = selected_pt.y;
  605. posinfo.sc_grab_z = selected_pt.z;
  606. posinfo.sc_count = (double)first_row_seedling_number;
  607. posinfo.sc_width = stem_width_mu_obstructed;
  608. posinfo.sc_tortuosity = stem_deflection_obstructed;
  609. }
  610. else {
  611. posinfo.rs_grab_x = selected_pt.x;
  612. posinfo.rs_grab_y = selected_pt.y;
  613. posinfo.rs_grab_z = selected_pt.z;
  614. posinfo.rs_count = (double)first_row_seedling_number;
  615. posinfo.rs_width = stem_width_mu_obstructed;
  616. posinfo.rs_tortuosity = stem_deflection_obstructed;
  617. }
  618. }
  619. //根据历史根的位置,计算对应位置点云数量,进而判断此位置是否有苗
  620. void CRootStockGrabPoint::occluded_seedling_detect_by_leaf(
  621. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //input 输入点云数据
  622. std::vector<int>& leaf_idx
  623. )
  624. {
  625. //仅用叶子点云数量,有时会出现因叶子误识别(和ror_ratio相关),造成叶子点云数量小,造成误判
  626. //所以改成整个空间点云数量做判别--2024-3-2
  627. m_root_center_with_seedling.clear();
  628. m_root_center_with_seedling.assign(m_root_centers.size(), false);
  629. std::vector<int> pcd_cnt;
  630. //int th_pcd_size = m_cparam.rs_grab_seedling_min_pts;
  631. double stem_diameter = m_cparam.rs_grab_stem_diameter;
  632. double y_min = m_cparam.rs_grab_ymin;
  633. double y_max = m_cparam.rs_grab_ymax;
  634. if (m_dtype == 0) {
  635. //th_pcd_size = m_cparam.sc_grab_seedling_min_pts;
  636. stem_diameter = m_cparam.sc_grab_stem_diameter;
  637. y_min = m_cparam.sc_grab_ymin;
  638. y_max = m_cparam.sc_grab_ymax;
  639. }
  640. int th_pcd_size = 0.333 * stem_diameter * (y_max - y_min) / m_cloud_mean_dist / m_cloud_mean_dist;
  641. std::vector<pcl::PointXYZ> aabb_mins_maxs;
  642. for (int i = 0; i < m_root_centers.size(); ++i) {
  643. CStemResult& rc = m_root_centers.at(i);
  644. pcl::PointXYZ aabb_min;
  645. pcl::PointXYZ aabb_max;
  646. int cnt = get_point_count_inbox(rc, aabb_min, aabb_max, in_cloud);
  647. pcd_cnt.push_back(cnt);
  648. bool has_seedling = cnt > th_pcd_size;
  649. m_root_center_with_seedling.at(i) = has_seedling;
  650. aabb_mins_maxs.push_back(aabb_min);
  651. aabb_mins_maxs.push_back(aabb_max);
  652. }
  653. if (m_pLogger) {
  654. stringstream buff;
  655. buff << m_pcdId << ": root positions points size " << m_root_centers.size() << " centers\n";
  656. for (int i = 0; i < m_root_centers.size();++i) {
  657. CStemResult& sr = m_root_centers.at(i);
  658. buff << "\tx=" << sr.root_x << ", z=" << sr.root_z
  659. << ", stem_pcd_size=" << sr.stem_size
  660. << ", raw_pcd_size=" << pcd_cnt.at(i) << "\n";
  661. }
  662. m_pLogger->INFO(buff.str());
  663. }
  664. // 显示
  665. if (m_cparam.image_show) {
  666. std::vector<std::vector<int> > clt_line_idx;
  667. std::vector<pcl::PointXYZ>target_root;
  668. double ymin = m_cparam.rs_grab_ymin;
  669. if (m_dtype == 0) { ymin = m_cparam.sc_grab_ymin; }
  670. for (auto&rc : m_root_centers) {
  671. pcl::PointXYZ p(rc.root_x, ymin, rc.root_z);
  672. target_root.push_back(p);
  673. }
  674. viewer_cloud_cluster_box(m_raw_cloud, target_root, aabb_mins_maxs, clt_line_idx, std::string("seedling root centers"));
  675. }
  676. }
  677. //int CRootStockGrabPoint::get_point_count_inbox(const CStemResult& sr,
  678. // pcl::PointXYZ& aabb_min,
  679. // pcl::PointXYZ& aabb_max)
  680. //{
  681. // double seedling_distance = m_cparam.rs_grab_seedling_dist;
  682. // double min_y = m_cparam.rs_grab_ymin;
  683. // if (m_dtype == 0) {
  684. // seedling_distance = m_cparam.sc_grab_seedling_dist;
  685. // min_y = m_cparam.sc_grab_ymin;
  686. // }
  687. // double min_x = sr.root_x - 0.5 * seedling_distance;
  688. // double max_x = sr.root_x + 0.5 * seedling_distance;
  689. // double min_z = sr.root_z - 0.75 * seedling_distance;
  690. // double max_z = sr.root_z + 0.25 * seedling_distance;
  691. // double max_y = min_y + 150.0; //假定苗高150mm
  692. // int count = 0;
  693. // for (auto&pt : m_raw_cloud->points) {
  694. // if(pt.y >= min_y && pt.y <= max_y &&
  695. // pt.x >=min_x && pt.x<=max_x &&
  696. // pt.z >= min_z && pt.z <= max_z)
  697. // {
  698. // count++;
  699. // }
  700. // }
  701. // aabb_min.x = min_x;
  702. // aabb_min.y = min_y;
  703. // aabb_min.z = min_z;
  704. // aabb_max.x = max_x;
  705. // aabb_max.y = max_y;
  706. // aabb_max.z = max_z;
  707. // return count;
  708. //}
  709. int CRootStockGrabPoint::get_point_count_inbox(const CStemResult& sr,
  710. pcl::PointXYZ& aabb_min,
  711. pcl::PointXYZ& aabb_max,
  712. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud //input 输入点云数据
  713. )
  714. {
  715. double seedling_distance = m_cparam.rs_grab_seedling_dist;
  716. double min_y = m_cparam.rs_grab_ymin;
  717. double max_y = m_cparam.rs_grab_ymax;
  718. if (m_dtype == 0) {
  719. seedling_distance = m_cparam.sc_grab_seedling_dist;
  720. min_y = m_cparam.sc_grab_ymin;
  721. max_y = m_cparam.sc_grab_ymax;
  722. }
  723. double min_x = sr.root_x - 0.5 * seedling_distance;
  724. double max_x = sr.root_x + 0.5 * seedling_distance;
  725. double min_z = sr.root_z - 0.75 * seedling_distance;
  726. double max_z = sr.root_z + 0.25 * seedling_distance;
  727. int count = 0;
  728. for (auto&pt : in_cloud->points) {
  729. if (pt.y >= min_y && pt.y <= max_y &&
  730. pt.x >= min_x && pt.x <= max_x &&
  731. pt.z >= min_z && pt.z <= max_z)
  732. {
  733. count++;
  734. }
  735. }
  736. aabb_min.x = min_x;
  737. aabb_min.y = min_y;
  738. aabb_min.z = min_z;
  739. aabb_max.x = max_x;
  740. aabb_max.y = max_y;
  741. aabb_max.z = max_z;
  742. return count;
  743. }
  744. //生成结果图片
  745. void CRootStockGrabPoint::gen_result_img(
  746. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud_raw,//输入,未经过滤的整体点云in_cloud_raw,
  747. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //输入,整体点云cloud_dowm_sampled,
  748. pcl::PointXYZ& selected_pt, //输入,selected_pt,
  749. pcl::PointXYZ& selected_pt_ref, //输入,selected_pt_ref,
  750. const PositionInfo& posinfo, //输入,
  751. cv::Mat& rst_img //输出,图片, 640*1280
  752. )
  753. {
  754. if (rst_img.empty()) { return; }
  755. if (rst_img.rows != 640 && rst_img.cols != 1280) { return; }
  756. int cx = 640; //图像中心点0
  757. int cy = 320; //图像中心点0
  758. float res = 0.33333; //分辨率,1个像素0.333mm
  759. //绘制坐标轴
  760. int arrow_len = 20;
  761. cv::line(rst_img, cv::Point(0, cy), cv::Point(cx, cy), cv::Scalar(128,128,128));
  762. cv::line(rst_img, cv::Point(0, cy), cv::Point(arrow_len, cy - int(arrow_len/2)), cv::Scalar(128, 128, 128));
  763. cv::line(rst_img, cv::Point(0, cy), cv::Point(arrow_len, cy + int(arrow_len/2)), cv::Scalar(128, 128, 128));
  764. cv::line(rst_img, cv::Point(cx, 0), cv::Point(cx, cy), cv::Scalar(128, 128, 128));
  765. cv::line(rst_img, cv::Point(cx, 0), cv::Point(cx - int(arrow_len/2), arrow_len), cv::Scalar(128, 128, 128));
  766. cv::line(rst_img, cv::Point(cx, 0), cv::Point(cx + int(arrow_len/2), arrow_len), cv::Scalar(128, 128, 128));
  767. cv::putText(rst_img, std::string("x"), cv::Point(20, cy-10), cv::FONT_HERSHEY_SIMPLEX, 1, cv::Scalar(128, 128, 128));
  768. cv::putText(rst_img, std::string("y"), cv::Point(cx+10, 20), cv::FONT_HERSHEY_SIMPLEX, 1, cv::Scalar(128, 128, 128));
  769. //绘制所有点
  770. int x, y;
  771. for (auto&pt : in_cloud_raw->points) {
  772. x = cx - int(pt.x / res + 0.5);
  773. y = cy - int(pt.y / res + 0.5);
  774. if (x < 0 || x >= rst_img.cols) { continue; }
  775. if (y < 0 || y >= rst_img.rows) { continue; }
  776. rst_img.at<cv::Vec3b>(y, x)[0] = 60;
  777. rst_img.at<cv::Vec3b>(y, x)[1] = 20;
  778. rst_img.at<cv::Vec3b>(y, x)[2] = 220;
  779. }
  780. for (auto&pt : in_cloud->points) {
  781. x = cx - int(pt.x / res + 0.5);
  782. y = cy - int(pt.y / res + 0.5);
  783. if (x < 0 || x >= rst_img.cols) { continue; }
  784. if (y < 0 || y >= rst_img.rows) { continue; }
  785. rst_img.at<cv::Vec3b>(y, x)[0] = 0;
  786. rst_img.at<cv::Vec3b>(y, x)[1] = 225;
  787. rst_img.at<cv::Vec3b>(y, x)[2] = 0;
  788. }
  789. //绘制抓取点坐标
  790. int grab_x = cx - int(selected_pt.x / res + 0.5);
  791. int grab_y = cy - int(selected_pt.y / res + 0.5);
  792. int radius = 30;
  793. cv::line(rst_img, cv::Point(grab_x - radius, grab_y - radius), cv::Point(grab_x + radius, grab_y + radius), cv::Scalar(0, 215, 255));
  794. cv::line(rst_img, cv::Point(grab_x - radius, grab_y + radius), cv::Point(grab_x + radius, grab_y - radius), cv::Scalar(0, 215, 255));
  795. //如果是遮挡,画一个圆
  796. if ((m_dtype == 0 && posinfo.sc_width == 0) || (m_dtype == 1 && posinfo.rs_width == 0)) {
  797. cv::circle(rst_img, cv::Point(grab_x, grab_y), radius, cv::Scalar(0, 215, 255));
  798. }
  799. //绘制茎节点坐标
  800. int fork_x = cx - int(selected_pt_ref.x / res + 0.5);
  801. int fork_y = cy - int(selected_pt_ref.y / res + 0.5);
  802. radius = 15;
  803. cv::line(rst_img, cv::Point(fork_x - radius, fork_y), cv::Point(fork_x + radius, fork_y), cv::Scalar(153, 51, 255));
  804. //在图片中写入文件名字
  805. cv::putText(rst_img, m_pcdId, cv::Point(20, 20), cv::FONT_HERSHEY_SIMPLEX, 0.5, cv::Scalar(128, 128, 128));
  806. }
  807. int CRootStockGrabPoint::read_ply_file(const char* fn)
  808. {
  809. m_raw_cloud.reset( new pcl::PointCloud<pcl::PointXYZ>);
  810. if (pcl::io::loadPLYFile<pcl::PointXYZ>(fn, *m_raw_cloud) == -1) {
  811. if (m_pLogger) {
  812. m_pLogger->ERRORINFO(m_pcdId + ": could not load file: "+std::string(fn));
  813. }
  814. return (-1);
  815. }
  816. if (m_pLogger) {
  817. stringstream buff;
  818. buff << m_pcdId <<": load raw point cloud " << m_raw_cloud->width * m_raw_cloud->height << " data points";
  819. m_pLogger->INFO(buff.str());
  820. }
  821. return m_raw_cloud->width * m_raw_cloud->height;
  822. }
  823. double CRootStockGrabPoint::compute_nearest_neighbor_distance(pcl::PointCloud<pcl::PointXYZ>::Ptr pcd)
  824. {
  825. pcl::KdTreeFLANN<pcl::PointXYZ> tree;
  826. tree.setInputCloud(pcd);
  827. int k = 2;
  828. double res = 0.0;
  829. int n_points = 0;
  830. for (size_t i = 0; i < pcd->size(); ++i) {
  831. std::vector<int> idx(k);
  832. std::vector<float> sqr_distances(k);
  833. if (tree.nearestKSearch(i, k, idx, sqr_distances) == k) {
  834. for (int id = 1; id < k; ++id) {
  835. res += sqrt(sqr_distances.at(id));
  836. ++n_points;
  837. }
  838. }
  839. }
  840. if (n_points > 0) {
  841. res /= (double)n_points;
  842. }
  843. return res;
  844. }
  845. ////////////////////////////////////////////////////////////
  846. //叶子剔除, 假设叶子和茎是分离的,用欧式聚类分割
  847. void CRootStockGrabPoint::leaf_filter(
  848. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //输入,降采样的点云
  849. std::vector<int> &stem_cloud_idx //输出,除了叶子的点云序号
  850. )
  851. {
  852. //采用欧式距离聚类,对每一类的点云进行分析,剔除叶子
  853. stem_cloud_idx.clear();
  854. std::vector<pcl::PointIndices> cluster_indices;
  855. pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>);
  856. tree->setInputCloud(in_cloud);
  857. pcl::EuclideanClusterExtraction<pcl::PointXYZ> ec;
  858. ec.setClusterTolerance(2.5);
  859. ec.setMinClusterSize(30);
  860. ec.setMaxClusterSize(20000);
  861. ec.setSearchMethod(tree);
  862. ec.setInputCloud(in_cloud);
  863. ec.extract(cluster_indices);
  864. for (std::vector<pcl::PointIndices>::const_iterator it = cluster_indices.begin();
  865. it != cluster_indices.end(); ++it) {
  866. pcl::PointCloud<pcl::PointXYZ>::Ptr cluster_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  867. pcl::copyPointCloud(*in_cloud, *it, *cluster_cloud);
  868. //计算点云的pca
  869. Eigen::Vector4f pcaCenteroid;
  870. pcl::compute3DCentroid(*cluster_cloud, pcaCenteroid);
  871. Eigen::Matrix3f covariance;
  872. pcl::computeCovarianceMatrix(*cluster_cloud, pcaCenteroid, covariance);
  873. Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigen_solver(covariance, Eigen::ComputeEigenvectors);
  874. Eigen::Vector3f eigenValuesPCA = eigen_solver.eigenvalues();
  875. float e1, e2, e3;
  876. e1 = std::sqrt(eigenValuesPCA(2));
  877. e2 = std::sqrt(eigenValuesPCA(1));
  878. e3 = std::sqrt(eigenValuesPCA(0));
  879. float a1d = (e1 - e2) / e1;
  880. float a2d = (e2 - e3) / e1;
  881. float a3d = e3 / e1;
  882. if (m_pLogger) {
  883. stringstream buff;
  884. buff << m_pcdId << ": ads=(" << a1d << "," << a2d << "," << a3d << "), pca_eigem_values=(" << e1
  885. << "," << e2 << "," << e3 << "), cluster_size=" << it->indices.size();
  886. m_pLogger->INFO(buff.str());
  887. }
  888. if (m_cparam.image_show) {
  889. if (a1d < 0.75) {
  890. viewer_cloud(cluster_cloud, string("cluster_leaf"));
  891. }
  892. else {
  893. viewer_cloud(cluster_cloud, string("cluster_stem"));
  894. }
  895. }
  896. /*
  897. 特征值归一化处理: a1d = (e1 - e2) / e1
  898. a2d = (e2 - e3) / e1
  899. a3d = e3 / e1
  900. ei = sqrt(lamda_i)
  901. 当点云线状分布时,a1d = 1, a2d = 0, a3d = 0
  902. 当点云面状分布时,a1d = 0, a2d = 1, a3d = 0
  903. 当点云球状分布时,a1d = 0, a2d = 0, a3d = 1
  904. */
  905. if (a1d < 0.75) { continue; }
  906. for (int idx : it->indices) {
  907. stem_cloud_idx.push_back(idx);
  908. }
  909. }
  910. }
  911. // 基于孤立点滤除的方法检测叶子区域
  912. void CRootStockGrabPoint::leaf_filter_ror(
  913. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //input 输入点云数据
  914. std::vector<int> &stem_cloud_idx, //output, 过滤后的点云序号
  915. std::vector<int>& leaf_idx //output, 叶子的点云序号
  916. )
  917. {
  918. //计算点云平均间距
  919. //double mean_dist = 0.0;
  920. //cloud_mean_dist(in_cloud, mean_dist);
  921. //计算点云过滤半径和点数阈值
  922. double stem_diameter = m_cparam.rs_grab_stem_diameter;
  923. if (m_dtype == 0) { stem_diameter = m_cparam.sc_grab_stem_diameter; }
  924. double remove_radius = stem_diameter;
  925. double ror_ratio = m_cparam.rs_grab_ror_ratio;
  926. if (m_dtype == 0){ror_ratio = m_cparam.sc_grab_ror_ratio;}
  927. int nb_points = ror_ratio* stem_diameter * stem_diameter * 2.0 / m_cloud_mean_dist / m_cloud_mean_dist; // (2d*d + pi *d*d) /2
  928. if (m_pLogger) {
  929. stringstream buff;
  930. buff << m_pcdId << ": ror_ratio=" << ror_ratio << ", ror filter nb_points=" << nb_points;
  931. m_pLogger->INFO(buff.str());
  932. }
  933. //获取叶片的点云
  934. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_leaf(new pcl::PointCloud<pcl::PointXYZ>);
  935. pcl::RadiusOutlierRemoval<pcl::PointXYZ> ror(true);
  936. ror.setInputCloud(in_cloud);
  937. ror.setRadiusSearch(remove_radius);
  938. ror.setMinNeighborsInRadius(nb_points);
  939. ror.filter(*cloud_leaf);
  940. //通过叶片点云寻找附近remove_radius内的近邻点,当做叶子点云序号
  941. std::set<int> leaf_idx_set;
  942. int nres;
  943. std::vector<int> indices;
  944. std::vector<float> sqr_distances;
  945. pcl::search::KdTree<pcl::PointXYZ> tree;
  946. tree.setInputCloud(in_cloud);
  947. for (size_t i = 0; i < cloud_leaf->size(); ++i)
  948. {
  949. nres = tree.radiusSearch(cloud_leaf->points.at(i), remove_radius, indices, sqr_distances);
  950. for (int& idx : indices) {
  951. leaf_idx_set.insert(idx);
  952. }
  953. }
  954. leaf_idx.clear();
  955. leaf_idx.assign(leaf_idx_set.begin(), leaf_idx_set.end());
  956. //in_img是经过开运算的图像,其中像素位置的点云为叶子区域
  957. std::vector<int>::iterator lit;
  958. for (int i = 0; i < in_cloud->points.size(); ++i) {
  959. lit = std::find(leaf_idx.begin(), leaf_idx.end(), i);
  960. if (lit == leaf_idx.end()) {
  961. stem_cloud_idx.push_back(i);
  962. }
  963. }
  964. //显示开运算的结果
  965. if (m_cparam.image_show) {
  966. pcl::PointCloud<pcl::PointXYZ>::Ptr stem_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  967. pcl::PointCloud<pcl::PointXYZ>::Ptr leaf_cloud(new pcl::PointCloud < pcl::PointXYZ>);
  968. pcl::copyPointCloud(*in_cloud, stem_cloud_idx, *stem_cloud);
  969. pcl::copyPointCloud(*in_cloud, leaf_idx, *leaf_cloud);
  970. pcl::visualization::PCLVisualizer viewer("open cloud");
  971. viewer.addCoordinateSystem();
  972. viewer.addPointCloud<pcl::PointXYZ>(stem_cloud, "stem_cloud");//????????????????
  973. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "stem_cloud");
  974. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "stem_cloud");
  975. viewer.addPointCloud<pcl::PointXYZ>(leaf_cloud, "leaf_cloud");
  976. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 0, 0, "leaf_cloud");
  977. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "leaf_cloud");
  978. while (!viewer.wasStopped()) {
  979. viewer.spinOnce(100);
  980. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  981. }
  982. }
  983. }
  984. void CRootStockGrabPoint::cloud_mean_dist(
  985. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //input 输入点云数据
  986. double& mean_dist
  987. )
  988. {
  989. mean_dist = 0.0;
  990. int n_points = 0;
  991. int nres;
  992. std::vector<int> indices(2);
  993. std::vector<float> sqr_distances(2);
  994. pcl::search::KdTree<pcl::PointXYZ> tree;
  995. tree.setInputCloud(in_cloud);
  996. std::vector<double> distances;
  997. double dist;
  998. for (size_t i = 0; i < in_cloud->size(); ++i)
  999. {
  1000. //Considering the second neighbor since the first is the point itself.
  1001. nres = tree.nearestKSearch(i, 2, indices, sqr_distances);
  1002. if (nres == 2)
  1003. {
  1004. dist = std::sqrt(sqr_distances[1]);
  1005. mean_dist += dist;
  1006. ++n_points;
  1007. distances.push_back(dist);
  1008. }
  1009. }
  1010. std::sort(distances.begin(), distances.end());
  1011. int per10 = int(distances.size() * 0.1);
  1012. mean_dist = distances.at(per10);
  1013. }
  1014. //////////////////////////////////////////////////////////////////////////////////////////
  1015. void CRootStockGrabPoint::gen_all_seedling_positions(
  1016. pcl::PointXYZ&key_center, //输入,已知的苗的坐标
  1017. std::vector<pcl::PointXYZ>&candidate_center //输出,有倾斜苗的坐标
  1018. )
  1019. {
  1020. //生成所有的植株位置,并排序
  1021. candidate_center.clear();
  1022. //找到z最小的那一株,并找出和它一排的那些
  1023. float step_harf = m_cparam.rs_grab_seedling_dist / 2.0;
  1024. float xmin = m_cparam.rs_grab_xmin;
  1025. float xmax = m_cparam.rs_grab_xmax;
  1026. float zmin = m_cparam.rs_grab_zmin;
  1027. float zmax = m_cparam.rs_grab_zmax;
  1028. int col_from = -8;
  1029. int col_to = 8;
  1030. int col_step = 1;
  1031. if (m_dtype == 0) {
  1032. step_harf = m_cparam.sc_grab_seedling_dist / 2.0;
  1033. xmin = m_cparam.sc_grab_xmin;
  1034. xmax = m_cparam.sc_grab_xmax;
  1035. zmin = m_cparam.sc_grab_zmin;
  1036. zmax = m_cparam.sc_grab_zmax;
  1037. }
  1038. // 生成所有可能的植株位置中心点
  1039. for (int row = -8; row <= 4; row += 1) {
  1040. float cz = key_center.z + row * step_harf/2.0;
  1041. if (cz < zmin || cz > zmax) { continue; }
  1042. for (int col = 4*col_from; col <= 4*col_to; col += col_step) {
  1043. float cx = key_center.x + col * step_harf/2;
  1044. if (cx < xmin || cx > xmax) { continue; }
  1045. //保存
  1046. pcl::PointXYZ pc(cx, 0.0, cz);
  1047. candidate_center.push_back(pc);
  1048. }
  1049. }
  1050. }
  1051. void CRootStockGrabPoint::nms_box(
  1052. std::vector<pcl::PointXYZ>&clt_root_v,
  1053. std::vector<float>&valid_box_cc_dist,
  1054. float hole_step_radius,
  1055. std::vector<pcl::PointXYZ>& target_toot)
  1056. {
  1057. target_toot.clear();
  1058. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
  1059. for (auto&p : clt_root_v) { cloud->points.push_back(p); }
  1060. pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;
  1061. kdtree.setInputCloud(cloud);
  1062. std::vector<int>idx;
  1063. std::vector<float>dist_sqr;
  1064. std::vector<int> tar_idx;
  1065. for (size_t i = 0; i < cloud->points.size(); ++i) {
  1066. int k = kdtree.radiusSearch(cloud->points.at(i), hole_step_radius, idx, dist_sqr);
  1067. std::vector<float> cc_tmp;
  1068. for (auto& j : idx) { cc_tmp.push_back(valid_box_cc_dist.at(j)); }
  1069. int min_id = std::min_element(cc_tmp.begin(), cc_tmp.end()) - cc_tmp.begin();
  1070. if (min_id == 0) {
  1071. tar_idx.push_back(i);
  1072. }
  1073. }
  1074. for (auto& i : tar_idx) {
  1075. target_toot.push_back(clt_root_v.at(i));
  1076. }
  1077. }
  1078. //通过指定位置内,取部分点云分析是否存在真正的茎,真茎位置保存到target_filtered
  1079. void CRootStockGrabPoint::line_filter(
  1080. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //输入点云
  1081. float radius, //输入,取点半径
  1082. float ymin,
  1083. float ymax,
  1084. std::vector<pcl::PointXYZ>&target_root, //输入,位置
  1085. std::vector<pcl::PointXYZ>&target_filtered, //输出,位置
  1086. std::vector<pcl::PointXYZ>&target_filtered_root, //输出,茎上根部点坐标
  1087. std::vector<std::vector<int>>&target_filtered_element, //输出,茎上点index
  1088. std::vector<pcl::ModelCoefficients::Ptr>& target_filtered_models//输出,茎直线模型
  1089. )
  1090. {
  1091. target_filtered.clear();
  1092. target_filtered_element.clear();
  1093. target_filtered_root.clear();
  1094. if (target_root.size() == 0) { return; }
  1095. for (auto&p : target_root) {
  1096. // 构建box,获取植株点云
  1097. pcl::PointCloud<pcl::PointXYZ>::Ptr seedling_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  1098. pcl::CropBox<pcl::PointXYZ> box_filter;
  1099. std::vector<int>inbox_idx;
  1100. pcl::PointXYZ min_point_aabb_ex;
  1101. pcl::PointXYZ max_point_aabb_ex;
  1102. min_point_aabb_ex.x = p.x - radius;
  1103. min_point_aabb_ex.y = ymin;
  1104. min_point_aabb_ex.z = p.z - radius;
  1105. max_point_aabb_ex.x = p.x + radius;
  1106. max_point_aabb_ex.y = ymax;
  1107. max_point_aabb_ex.z = p.z + radius;
  1108. box_filter.setMin(Eigen::Vector4f(min_point_aabb_ex.x, min_point_aabb_ex.y, min_point_aabb_ex.z, 1));
  1109. box_filter.setMax(Eigen::Vector4f(max_point_aabb_ex.x, max_point_aabb_ex.y, max_point_aabb_ex.z, 1));
  1110. box_filter.setNegative(false);
  1111. box_filter.setInputCloud(in_cloud);
  1112. box_filter.filter(inbox_idx);
  1113. pcl::copyPointCloud(*in_cloud, inbox_idx, *seedling_inbox);
  1114. //植株点云直线查找
  1115. //找到inbox点云中的直线
  1116. float stem_radius = m_cparam.rs_grab_stem_diameter / 2.0;
  1117. if (m_dtype == 0) {
  1118. stem_radius = m_cparam.sc_grab_stem_diameter / 2.0;
  1119. }
  1120. pcl::ModelCoefficients::Ptr coefficients(new pcl::ModelCoefficients);
  1121. pcl::PointIndices::Ptr inliers(new pcl::PointIndices);
  1122. pcl::SACSegmentation<pcl::PointXYZ> seg;
  1123. seg.setOptimizeCoefficients(true);
  1124. seg.setModelType(pcl::SACMODEL_LINE);
  1125. seg.setDistanceThreshold(stem_radius);
  1126. seg.setInputCloud(seedling_inbox);
  1127. seg.segment(*inliers, *coefficients);
  1128. pcl::PointCloud<pcl::PointXYZ>::Ptr stem_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  1129. pcl::copyPointCloud(*seedling_inbox, *inliers, *stem_cloud);
  1130. //点数过滤
  1131. int min_stem_pts = m_cparam.rs_grab_stem_min_pts;
  1132. if (m_dtype == 0) { min_stem_pts = m_cparam.sc_grab_stem_min_pts; }
  1133. if (inliers->indices.size() < int(min_stem_pts / 2)) { continue; }
  1134. //y方向分布范围滤波
  1135. float y_length_th = 10.0;
  1136. pcl::PointXYZ min_v;
  1137. pcl::PointXYZ max_v;
  1138. pcl::getMinMax3D(*stem_cloud, min_v, max_v);
  1139. float dy = max_v.y - min_v.y;
  1140. if (dy<y_length_th && dy / (ymax - ymin) < 0.25) { continue; }
  1141. //y方向分布中心滤波
  1142. float dy_c = 0.5*(max_v.y + min_v.y);
  1143. if ((dy_c - ymin) / (ymax - ymin) > 0.75) { continue; }
  1144. //倾斜角度过滤
  1145. float xz = std::sqrt(coefficients->values.at(3) * coefficients->values.at(3) +
  1146. coefficients->values.at(5) * coefficients->values.at(5));
  1147. float a = std::atan2f(coefficients->values.at(4), xz);
  1148. a = std::fabs(a) * 180.0 / 3.14159;
  1149. if (a < 45.0) { continue; }
  1150. target_filtered.push_back(p);
  1151. //有效茎长计算
  1152. std::vector<int> stem_y_count_hist;
  1153. get_line_project_hist(stem_cloud, coefficients, stem_y_count_hist);
  1154. std::vector<int> stem_y_count_hist_valid;
  1155. for (auto&v : stem_y_count_hist) {
  1156. if (v > 2) {
  1157. stem_y_count_hist_valid.push_back(v);
  1158. }
  1159. }
  1160. double valid_mu = get_hist_mean(stem_y_count_hist_valid);
  1161. int cnt_th = static_cast<int>(valid_mu/2.0);
  1162. if (cnt_th < 3) { cnt_th = 2; }
  1163. float valid_length = 0.0;
  1164. std::vector<int> stem_y_mask;
  1165. for (auto& c : stem_y_count_hist) {
  1166. if (c > cnt_th) {
  1167. valid_length += 1.0;
  1168. stem_y_mask.push_back(1);
  1169. }
  1170. else {
  1171. stem_y_mask.push_back(0);
  1172. }
  1173. }
  1174. //填补空
  1175. for (int i = 1; i < stem_y_mask.size() - 1; ++i) {
  1176. if (stem_y_mask.at(i) == 0 && stem_y_mask.at(i - 1) == 1 && stem_y_mask.at(i + 1) == 1) {
  1177. stem_y_mask.at(i) = 1;
  1178. }
  1179. }
  1180. //找出最长的连续段的长度
  1181. int longest_segment = 0;
  1182. int start = 0;
  1183. for (int i = 0; i < stem_y_mask.size(); ++i) {
  1184. if (i == 0 ) {
  1185. if (stem_y_mask.at(i) == 1) {
  1186. start = i;
  1187. }
  1188. continue;
  1189. }
  1190. if (i == stem_y_mask.size() - 1 ) {
  1191. if (stem_y_mask.at(i) == 1) {
  1192. int length = i - start + 1;
  1193. if (length > longest_segment) { longest_segment = length; }
  1194. }
  1195. continue;
  1196. }
  1197. if (stem_y_mask.at(i - 1) == 0 && stem_y_mask.at(i) == 1) {
  1198. start = i;
  1199. continue;
  1200. }
  1201. if (stem_y_mask.at(i - 1) == 1 && stem_y_mask.at(i) == 0) {
  1202. int length = i - start;
  1203. if (length > longest_segment) { longest_segment = length; }
  1204. continue;
  1205. }
  1206. }
  1207. if (longest_segment<10.0 && valid_length / (ymax - ymin) < 0.35) { continue; }
  1208. float min_y = 10000.0;
  1209. int min_y_idx = -1;
  1210. for (size_t j = 0; j < stem_cloud->points.size();++j) {
  1211. pcl::PointXYZ &spt = stem_cloud->at(j);
  1212. if (spt.y < min_y) {
  1213. min_y = spt.y;
  1214. min_y_idx = j;
  1215. }
  1216. }
  1217. pcl::PointXYZ tr_pt(stem_cloud->points.at(min_y_idx).x, stem_cloud->points.at(min_y_idx).y, stem_cloud->points.at(min_y_idx).z);
  1218. target_filtered_root.push_back(tr_pt);
  1219. int idx_raw = 0;
  1220. std::vector<int> out_idx;
  1221. for (auto& idx : inliers->indices) {
  1222. idx_raw = inbox_idx.at(idx);
  1223. out_idx.push_back(idx_raw);
  1224. }
  1225. target_filtered_element.push_back(out_idx);
  1226. target_filtered_models.push_back(coefficients);
  1227. // debug show
  1228. /*if (m_cparam.image_show) {
  1229. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_line(new pcl::PointCloud<pcl::PointXYZ>);
  1230. pcl::copyPointCloud(*seedling_inbox, *inliers, *cloud_line);
  1231. pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("line in cluster"));
  1232. viewer->addPointCloud<pcl::PointXYZ>(seedling_inbox, "cluster");
  1233. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "cluster");
  1234. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "cluster");
  1235. viewer->addPointCloud(cloud_line, "fit_line");
  1236. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0, 1, 0, "fit_line");
  1237. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, "fit_line");
  1238. while (!viewer->wasStopped()) {
  1239. viewer->spinOnce(100);
  1240. }
  1241. }*/
  1242. }
  1243. }
  1244. bool CRootStockGrabPoint::find_seedling_position_key(
  1245. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  1246. std::vector<int> &first_seedling_cloud_idx,
  1247. pcl::PointXYZ&xz_center,
  1248. pcl::ModelCoefficients::Ptr& line_model,
  1249. int& first_row_size
  1250. )
  1251. {
  1252. first_row_size = 0;
  1253. // 确定植株inbox范围
  1254. float hole_step = m_cparam.rs_grab_seedling_dist - 5.0; //穴盘中穴间距
  1255. if (m_dtype == 0) {
  1256. hole_step = m_cparam.sc_grab_seedling_dist - 5.0;
  1257. }
  1258. float hole_step_radius = hole_step / 2.0;
  1259. // 点云降维到xz平面,y=0
  1260. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2d(new pcl::PointCloud < pcl::PointXYZ>);
  1261. pcl::copyPointCloud(*in_cloud, *cloud2d);
  1262. for (auto&pt : *cloud2d) {
  1263. pt.y = 0.0;
  1264. }
  1265. // 在xz平面内统计点的密度
  1266. double radius = m_cparam.rs_grab_stem_diameter;
  1267. if (m_dtype == 0) {
  1268. radius = m_cparam.sc_grab_stem_diameter;
  1269. }
  1270. std::vector<int> counter;
  1271. pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;
  1272. kdtree.setInputCloud(cloud2d);
  1273. std::vector<int>idx;
  1274. std::vector<float>dist_sqr;
  1275. for (size_t i = 0; i < cloud2d->points.size(); ++i) {
  1276. int k = kdtree.radiusSearch(cloud2d->points.at(i), radius, idx, dist_sqr);
  1277. counter.push_back(k);
  1278. }
  1279. int max_density_idx = std::max_element(counter.begin(), counter.end()) - counter.begin();
  1280. int max_density = counter.at(max_density_idx);
  1281. if (m_cparam.image_show) {
  1282. //显示密度最大的点的位置
  1283. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_tmp(new pcl::PointCloud<pcl::PointXYZRGB>);
  1284. pcl::copyPointCloud(*in_cloud, *cloud_tmp);
  1285. for (auto& pt : *cloud_tmp) {
  1286. pt.r = 255;
  1287. pt.g = 255;
  1288. pt.b = 255;
  1289. }
  1290. std::vector<pcl::PointXYZ>cc;
  1291. cc.push_back(in_cloud->points.at(max_density_idx));
  1292. viewer_cloud_cluster(cloud_tmp, cc, string("key"));
  1293. }
  1294. // 生成植株的中心及box
  1295. std::vector<pcl::PointXYZ>clt_root;
  1296. std::vector<pcl::PointXYZ> clt_box;
  1297. float ymin = m_cparam.rs_grab_ymin;
  1298. float ymax = m_cparam.rs_grab_ymax;
  1299. if (m_dtype == 0) {
  1300. ymin = m_cparam.sc_grab_ymin;
  1301. ymax = m_cparam.sc_grab_ymax;
  1302. }
  1303. gen_all_seedling_positions(in_cloud->points.at(max_density_idx), clt_root);
  1304. // 显示生成的每一个候选框
  1305. /*if (m_cparam.image_show) {
  1306. std::vector<std::vector<int> > clt_line_idx;
  1307. std::vector<pcl::PointXYZ> clt_root_v;
  1308. std::vector<pcl::PointXYZ> clt_box_v;
  1309. for (auto&p : clt_root) {
  1310. pcl::PointXYZ p_show(p.x, ymin, p.z);
  1311. clt_root_v.push_back(p_show);
  1312. pcl::PointXYZ min_point_aabb_ex;
  1313. pcl::PointXYZ max_point_aabb_ex;
  1314. min_point_aabb_ex.x = p.x - hole_step_radius;
  1315. min_point_aabb_ex.y = ymin;
  1316. min_point_aabb_ex.z = p.z - hole_step_radius;
  1317. max_point_aabb_ex.x = p.x + hole_step_radius;
  1318. max_point_aabb_ex.y = ymax;
  1319. max_point_aabb_ex.z = p.z + hole_step_radius;
  1320. clt_box_v.push_back(min_point_aabb_ex);
  1321. clt_box_v.push_back(max_point_aabb_ex);
  1322. }
  1323. viewer_cloud_cluster_box(in_cloud, clt_root_v, clt_box_v, clt_line_idx, std::string("candidate_box"));
  1324. }*/
  1325. // 每个位置点云情况,过滤
  1326. std::vector<int> valid_index; //初步有效的矩形index
  1327. std::vector<int> valid_box_pts; //立方体内点云数量
  1328. std::vector<float>valid_box_cc_dist;//重心和矩形中心距离
  1329. for (size_t i = 0; i < clt_root.size();++i) {
  1330. pcl::PointXYZ &p = clt_root.at(i);
  1331. pcl::PointCloud<pcl::PointXYZ>::Ptr seedling_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  1332. pcl::CropBox<pcl::PointXYZ> box_filter;
  1333. std::vector<int>inbox_idx;
  1334. pcl::PointXYZ min_point_aabb_ex;
  1335. pcl::PointXYZ max_point_aabb_ex;
  1336. min_point_aabb_ex.x = p.x - hole_step_radius;
  1337. min_point_aabb_ex.y = ymin;
  1338. min_point_aabb_ex.z = p.z - hole_step_radius;
  1339. max_point_aabb_ex.x = p.x + hole_step_radius;
  1340. max_point_aabb_ex.y = ymax;
  1341. max_point_aabb_ex.z = p.z + hole_step_radius;
  1342. box_filter.setMin(Eigen::Vector4f(min_point_aabb_ex.x, min_point_aabb_ex.y, min_point_aabb_ex.z, 1));
  1343. box_filter.setMax(Eigen::Vector4f(max_point_aabb_ex.x, max_point_aabb_ex.y, max_point_aabb_ex.z, 1));
  1344. box_filter.setNegative(false);
  1345. box_filter.setInputCloud(in_cloud);
  1346. box_filter.filter(inbox_idx);
  1347. pcl::copyPointCloud(*in_cloud, inbox_idx, *seedling_inbox);
  1348. //点数过滤
  1349. if (inbox_idx.size() < int(max_density/3)) { continue; }
  1350. //y方向分布范围滤波
  1351. pcl::PointXYZ min_v;
  1352. pcl::PointXYZ max_v;
  1353. pcl::getMinMax3D(*seedling_inbox, min_v, max_v);
  1354. float dy = max_v.y - min_v.y;
  1355. if (dy / (ymax - ymin) < 0.35) { continue; }
  1356. //y方向分布中心滤波
  1357. float dy_c = 0.5*(max_v.y + min_v.y);
  1358. if ((dy_c-ymin) / (ymax - ymin) > 0.75) { continue; }
  1359. //计算中心
  1360. Eigen::Vector4f mean;
  1361. pcl::compute3DCentroid(*seedling_inbox, mean);
  1362. double cc_dist = std::sqrt((mean[0] - p.x)*(mean[0] - p.x) + (mean[2] - p.z)*(mean[2] - p.z));
  1363. valid_index.push_back(i);
  1364. valid_box_pts.push_back(inbox_idx.size());
  1365. valid_box_cc_dist.push_back((float)cc_dist);
  1366. }
  1367. // 找到局部最大值,找出真正的位置
  1368. std::vector<pcl::PointXYZ> clt_root_v;
  1369. std::vector<pcl::PointXYZ> clt_box_v;
  1370. for (auto&i : valid_index) {
  1371. pcl::PointXYZ&p = clt_root.at(i);
  1372. pcl::PointXYZ p_show(p.x, ymin, p.z);
  1373. clt_root_v.push_back(p_show);
  1374. pcl::PointXYZ min_point_aabb_ex;
  1375. pcl::PointXYZ max_point_aabb_ex;
  1376. min_point_aabb_ex.x = p.x - hole_step_radius;
  1377. min_point_aabb_ex.y = ymin;
  1378. min_point_aabb_ex.z = p.z - hole_step_radius;
  1379. max_point_aabb_ex.x = p.x + hole_step_radius;
  1380. max_point_aabb_ex.y = ymax;
  1381. max_point_aabb_ex.z = p.z + hole_step_radius;
  1382. clt_box_v.push_back(min_point_aabb_ex);
  1383. clt_box_v.push_back(max_point_aabb_ex);
  1384. }
  1385. // 显示
  1386. if (m_cparam.image_show) {
  1387. std::vector<std::vector<int> > clt_line_idx;
  1388. viewer_cloud_cluster_box(in_cloud, clt_root_v, clt_box_v, clt_line_idx, std::string("valid_box"));
  1389. }
  1390. std::vector<pcl::PointXYZ> target_root;
  1391. nms_box(clt_root_v, valid_box_cc_dist, hole_step_radius, target_root);
  1392. // 显示
  1393. if (m_cparam.image_show) {
  1394. std::vector<std::vector<int> > clt_line_idx;
  1395. std::vector<pcl::PointXYZ>clt_box_tmp;
  1396. for (auto&p : target_root) {
  1397. pcl::PointXYZ p_show(p.x, ymin, p.z);
  1398. clt_root_v.push_back(p_show);
  1399. pcl::PointXYZ min_point_aabb_ex;
  1400. pcl::PointXYZ max_point_aabb_ex;
  1401. min_point_aabb_ex.x = p.x - hole_step_radius;
  1402. min_point_aabb_ex.y = ymin;
  1403. min_point_aabb_ex.z = p.z - hole_step_radius;
  1404. max_point_aabb_ex.x = p.x + hole_step_radius;
  1405. max_point_aabb_ex.y = ymax;
  1406. max_point_aabb_ex.z = p.z + hole_step_radius;
  1407. clt_box_tmp.push_back(min_point_aabb_ex);
  1408. clt_box_tmp.push_back(max_point_aabb_ex);
  1409. }
  1410. viewer_cloud_cluster_box(in_cloud, target_root, clt_box_tmp, clt_line_idx, std::string("nms_box"));
  1411. }
  1412. ///////////////////////////////////////////////////////////////////////////////////////////////////////////
  1413. ///////////////////////////////////////////////////////////////////////////////////////////////////////////
  1414. // 根据得到的位置,直线检测,并过滤
  1415. std::vector<pcl::PointXYZ>target_filtered;
  1416. std::vector<std::vector<int>> target_member;//输出,茎上点index
  1417. std::vector<pcl::PointXYZ>target_filtered_root;
  1418. std::vector<pcl::ModelCoefficients::Ptr> target_filtered_models; //茎直线模型
  1419. line_filter(in_cloud, hole_step_radius,ymin,ymax, target_root, target_filtered,
  1420. target_filtered_root,
  1421. target_member,
  1422. target_filtered_models);
  1423. if (target_filtered_root.size() == 0) {
  1424. if (m_pLogger) {
  1425. stringstream buff;
  1426. buff << m_pcdId << ": line_filter reuslt cluster is empty" ;
  1427. m_pLogger->INFO(buff.str());
  1428. }
  1429. //统计苗数量
  1430. first_row_size = 0;
  1431. for (auto&has_seedling : m_root_center_with_seedling) {
  1432. if (has_seedling) { first_row_size += 1; }
  1433. }
  1434. return false;
  1435. }
  1436. if (m_cparam.image_show) {
  1437. std::vector<pcl::PointXYZ> final_box;
  1438. for (auto&pt : target_filtered_root) {
  1439. //扩展矩形范围
  1440. pcl::PointXYZ min_point_aabb_ex;
  1441. pcl::PointXYZ max_point_aabb_ex;
  1442. min_point_aabb_ex.x = pt.x - hole_step_radius;
  1443. min_point_aabb_ex.y = ymin;
  1444. min_point_aabb_ex.z = pt.z - hole_step_radius;
  1445. max_point_aabb_ex.x = pt.x + hole_step_radius;
  1446. max_point_aabb_ex.y = ymax;
  1447. max_point_aabb_ex.z = pt.z + hole_step_radius;
  1448. final_box.push_back(min_point_aabb_ex);
  1449. final_box.push_back(max_point_aabb_ex);
  1450. }
  1451. viewer_cloud_cluster_box(in_cloud, target_filtered_root, final_box, target_member, std::string("line_filter"));
  1452. }
  1453. //保存茎的根部坐标和点云数量
  1454. if (m_pStemInfos) {
  1455. for (size_t k = 0; k < target_filtered_root.size(); ++k) {
  1456. CStemResult sr = CStemResult(target_filtered_root.at(k).x,
  1457. target_filtered_root.at(k).y,
  1458. target_filtered_root.at(k).z,
  1459. target_member.at(k).size(),
  1460. m_pcdId,
  1461. 1);
  1462. m_pStemInfos->append(sr);
  1463. }
  1464. }
  1465. //sort cluster center, get the frontest leftest one
  1466. //找最小z,然后计算平均z
  1467. float mean_root_z_all = 0.0;
  1468. for(auto&pt : target_filtered_root) {
  1469. mean_root_z_all += pt.z;
  1470. }
  1471. mean_root_z_all /= target_filtered_root.size();
  1472. //通过最小z,保守的找到和他一排的植株
  1473. std::vector<int> first_row_index;
  1474. float mean_z = 0.0;
  1475. for (int idx_r = 0; idx_r < target_filtered_root.size();++idx_r) {
  1476. pcl::PointXYZ&pt = target_filtered_root.at(idx_r);
  1477. if (std::fabs(pt.z - mean_root_z_all) < hole_step_radius) {
  1478. first_row_index.push_back(idx_r);
  1479. mean_z += pt.z;
  1480. }
  1481. }
  1482. if (m_pLogger) {
  1483. stringstream buff;
  1484. buff << m_pcdId << ": current seedling number:"<< first_row_index.size()
  1485. <<", mean_z:"<< mean_z;
  1486. m_pLogger->INFO(buff.str());
  1487. }
  1488. //如果第一排的植株大于3,计算平均值
  1489. float first_row_num = (float)first_row_index.size();
  1490. mean_z /= first_row_num;
  1491. if (first_row_num >= 4) {
  1492. if (m_dtype == 0) {
  1493. m_1st_row_zmean_sc = mean_z;
  1494. }
  1495. else{ m_1st_row_zmean_rs = mean_z; }
  1496. }
  1497. else {
  1498. if( m_dtype == 0) {
  1499. if (m_1st_row_zmean_sc > 0) {
  1500. mean_z = m_1st_row_zmean_sc; }
  1501. }
  1502. else {
  1503. if (m_1st_row_zmean_rs > 0) {
  1504. mean_z = m_1st_row_zmean_rs; }
  1505. }
  1506. if (m_pLogger) {
  1507. stringstream buff;
  1508. buff << m_pcdId << ": update with historic mean_z:" << mean_z;
  1509. m_pLogger->INFO(buff.str());
  1510. }
  1511. }
  1512. //通过mean_z再次寻找第一排的植株
  1513. first_row_index.clear();
  1514. for (int idx_r = 0; idx_r < target_filtered_root.size(); ++idx_r) {
  1515. pcl::PointXYZ&pt = target_filtered_root.at(idx_r);
  1516. if (std::fabs(pt.z - mean_z) < hole_step_radius) {
  1517. first_row_index.push_back(idx_r);
  1518. }
  1519. }
  1520. //找第一排的植株没有满足要求的
  1521. if (first_row_index.size()==0) {
  1522. if (m_pLogger) {
  1523. stringstream buff;
  1524. buff << m_pcdId << ": NOT found valid seedlings";
  1525. m_pLogger->INFO(buff.str());
  1526. }
  1527. //统计苗数量
  1528. first_row_size = 0;
  1529. for (auto&has_seedling : m_root_center_with_seedling) {
  1530. if (has_seedling) { first_row_size += 1; }
  1531. }
  1532. return false;
  1533. }
  1534. ///////////////////////////////////////////////////////////
  1535. std::vector<float> cluster_index;
  1536. for (int i=0; i<first_row_index.size();++i){
  1537. int raw_idx = first_row_index.at(i);
  1538. pcl::PointXYZ&pt = target_filtered_root.at(raw_idx);
  1539. cluster_index.push_back(pt.x);
  1540. //update m_root_center_with_seedling
  1541. double min_dist = hole_step;
  1542. int min_i = -1;
  1543. for (int rc_i = 0; rc_i < m_root_centers.size(); ++rc_i)
  1544. {
  1545. double dist = fabs(pt.x - m_root_centers.at(rc_i).root_x);
  1546. if (dist < min_dist) {
  1547. min_dist = dist;
  1548. min_i = rc_i;
  1549. }
  1550. }
  1551. if (min_i >= 0) {
  1552. m_root_center_with_seedling.at(min_i) = true;
  1553. }
  1554. }
  1555. //统计苗数量
  1556. first_row_size = 0;
  1557. for (auto&has_seedling : m_root_center_with_seedling) {
  1558. if (has_seedling) { first_row_size += 1; }
  1559. }
  1560. int first_idx = std::min_element(cluster_index.begin(), cluster_index.end()) - cluster_index.begin();
  1561. if (m_dtype == 0) {
  1562. first_idx = std::max_element(cluster_index.begin(), cluster_index.end()) - cluster_index.begin();
  1563. }
  1564. first_idx = first_row_index.at(first_idx);
  1565. first_seedling_cloud_idx.clear();
  1566. for (auto&v : target_member.at(first_idx)) {
  1567. first_seedling_cloud_idx.push_back(v);
  1568. }
  1569. xz_center.x = target_filtered_root.at(first_idx).x;
  1570. xz_center.y = target_filtered_root.at(first_idx).y;
  1571. xz_center.z = target_filtered_root.at(first_idx).z;
  1572. line_model = target_filtered_models.at(first_idx);
  1573. if (m_pLogger) {
  1574. stringstream buff;
  1575. buff << m_pcdId << ": euclidean_clustering_ttsas, find cluster center(" << xz_center.x
  1576. <<", "<< xz_center.y<<", "<< xz_center.z<<")";
  1577. m_pLogger->INFO(buff.str());
  1578. }
  1579. return true;
  1580. }
  1581. void CRootStockGrabPoint::tilted_seedling_discover(
  1582. std::vector<pcl::PointXYZ>&cluster_center, //输入,已知的苗的坐标
  1583. std::vector<pcl::PointXYZ>&tilted_center //输出,有倾斜苗的坐标
  1584. )
  1585. {
  1586. tilted_center.clear();
  1587. if (cluster_center.size() == 0) { return; }
  1588. //找到z最小的那一株,并找出和它一排的那些
  1589. float step_harf = m_cparam.rs_grab_seedling_dist / 2.0;
  1590. float xmin = m_cparam.rs_grab_xmin;
  1591. float xmax = m_cparam.rs_grab_xmax;
  1592. float zmin = m_cparam.rs_grab_zmin;
  1593. float zmax = m_cparam.rs_grab_zmax;
  1594. if (m_dtype == 0) {
  1595. step_harf = m_cparam.sc_grab_seedling_dist / 2.0;
  1596. xmin = m_cparam.sc_grab_xmin;
  1597. xmax = m_cparam.sc_grab_xmax;
  1598. zmin = m_cparam.sc_grab_zmin;
  1599. zmax = m_cparam.sc_grab_zmax;
  1600. }
  1601. auto minz_it = std::min_element(cluster_center.begin(), cluster_center.end(),
  1602. [](const pcl::PointXYZ& p1, const pcl::PointXYZ& p2)-> bool {return p1.z < p2.z; });
  1603. float minz = minz_it->z;
  1604. float meanz = 0.0;
  1605. float cnt = 0.0;
  1606. for (auto&p : cluster_center) {
  1607. if (fabs(p.z - minz) < step_harf) { meanz += p.z; cnt += 1.0; }
  1608. }
  1609. if(cnt>0.0){meanz /= cnt;}
  1610. else { meanz = minz; }
  1611. // 生成所有可能的植株位置中心点
  1612. std::vector<pcl::PointXYZ>candidate_centers;
  1613. for (int row = -1; row <= 1; row += 1) {
  1614. for (int col = -6; col <= 6; col += 1) {
  1615. float cx = minz_it->x + col * 2.0 * step_harf;
  1616. float cz = meanz + row * 2.0 * step_harf;
  1617. if (cx < xmin || cx > xmax) { continue; }
  1618. if (cz < zmin || cz > zmax) { continue; }
  1619. // 邻域内已经存在植株的跳过
  1620. bool with_seedling = false;
  1621. for (auto&p : cluster_center) {
  1622. float d = std::sqrt((p.x - cx) * (p.x - cx) + (p.z - cz) * (p.z - cz));
  1623. if (d < 2.0*step_harf) {
  1624. with_seedling=true;
  1625. break;
  1626. }
  1627. }
  1628. if (!with_seedling) {//保存
  1629. pcl::PointXYZ pc(cx,0.0,cz);
  1630. candidate_centers.push_back(pc);
  1631. }
  1632. }
  1633. }
  1634. //针对每一个candidat_centers的邻域进行检查,检测是否有茎目标
  1635. for (auto&pc : candidate_centers) {
  1636. }
  1637. }
  1638. void CRootStockGrabPoint::find_seedling_position(
  1639. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  1640. std::vector<int> &first_seedling_cloud_idx,
  1641. pcl::PointXYZ&xz_center
  1642. )
  1643. {
  1644. // 确定植株inbox范围
  1645. float hole_step = m_cparam.rs_grab_seedling_dist - 10.0; //穴盘中穴间距
  1646. if (m_dtype == 0) {
  1647. hole_step = m_cparam.sc_grab_seedling_dist - 10.0;
  1648. }
  1649. float hole_step_radius = hole_step / 2.0;
  1650. // 点云降维到xz平面,y=0
  1651. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2d(new pcl::PointCloud < pcl::PointXYZ>);
  1652. pcl::copyPointCloud(*in_cloud, *cloud2d);
  1653. for (auto&pt : *cloud2d) {
  1654. pt.y = 0.0;
  1655. }
  1656. /*if(m_cparam.image_show){
  1657. viewer_cloud(cloud2d, std::string("cloud2d"));
  1658. } */
  1659. // 在xz平面内统计点的密度
  1660. double radius = m_cparam.rs_grab_stem_diameter;
  1661. if (m_dtype == 0) {
  1662. radius = m_cparam.sc_grab_stem_diameter;
  1663. }
  1664. std::vector<int> counter;
  1665. pcl::KdTreeFLANN<pcl::PointXYZ> kdtree;
  1666. kdtree.setInputCloud(cloud2d);
  1667. std::vector<int>idx;
  1668. std::vector<float>dist_sqr;
  1669. for (size_t i = 0; i < cloud2d->points.size(); ++i) {
  1670. int k = kdtree.radiusSearch(cloud2d->points.at(i), radius, idx, dist_sqr);
  1671. counter.push_back(k);
  1672. }
  1673. // 计算根据密度分割的自动阈值
  1674. int th = (int)(otsu(counter));
  1675. // 得到茎目标的点云序号,和对应的点云cloud2d_pos
  1676. std::vector<int> index;
  1677. for (size_t i = 0; i < counter.size(); ++i) {
  1678. if (counter.at(i) >= th) {
  1679. index.push_back(i);
  1680. }
  1681. }
  1682. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud2d_pos(new pcl::PointCloud < pcl::PointXYZ>);
  1683. pcl::copyPointCloud(*cloud2d, index, *cloud2d_pos);
  1684. if (m_pLogger) {
  1685. stringstream buff;
  1686. buff << m_pcdId <<": get 2d seedling seed point cloud " << index.size() << " data points with thrreshold " << th;
  1687. m_pLogger->INFO(buff.str());
  1688. }
  1689. if (m_cparam.image_show) {
  1690. viewer_cloud(cloud2d_pos, std::string("cloud2d_pos"));
  1691. }
  1692. //对茎的点云进行clustering,得到不同的茎的分组
  1693. double d1 = m_cparam.rs_grab_stem_diameter;
  1694. double d2 = m_cparam.rs_grab_stem_diameter * 3.0;
  1695. if (m_dtype == 0) {
  1696. d1 = m_cparam.sc_grab_stem_diameter;
  1697. d2 = m_cparam.sc_grab_stem_diameter * 3.0;
  1698. }
  1699. std::vector<pcl::PointXYZ>cluster_center;
  1700. std::vector<std::vector<int>> cluster_member;
  1701. euclidean_clustering_ttsas(cloud2d_pos, d1, d2, cluster_center, cluster_member);
  1702. if (m_pLogger) {
  1703. stringstream buff;
  1704. buff << m_pcdId <<": euclidean_clustering_ttsas: " << cluster_center.size() << " t1=" << d1
  1705. << " t2=" << d2;
  1706. m_pLogger->INFO(buff.str());
  1707. }
  1708. if (m_cparam.image_show) {
  1709. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_cluster(new pcl::PointCloud<pcl::PointXYZRGB>);
  1710. pcl::copyPointCloud(*cloud2d_pos, *cloud_cluster);
  1711. for (auto& pt : *cloud_cluster) {
  1712. pt.r = 255;
  1713. pt.g = 255;
  1714. pt.b = 255;
  1715. }
  1716. viewer_cloud_cluster(cloud_cluster, cluster_center, std::string("cloud2d_cluster_pos"));
  1717. }
  1718. // 扩展检测第一排的可能位置,找到倾斜的苗
  1719. std::vector<pcl::PointXYZ> tilted_center;
  1720. tilted_seedling_discover(cluster_center, tilted_center);
  1721. //对每个类别进行检测,剔除不是茎的类别
  1722. std::vector<pcl::PointXYZ> clt_root;
  1723. std::vector<pcl::PointXYZ> clt_box;
  1724. std::vector<std::vector<int> > clt_line_idx;
  1725. for (size_t i = 0; i < cluster_center.size(); ++i) {
  1726. pcl::PointIndices cluster_idxs;
  1727. for (auto idx_clt : cluster_member.at(i)) {
  1728. int idx_raw = index.at(idx_clt);
  1729. cluster_idxs.indices.push_back(idx_raw);
  1730. }
  1731. pcl::PointCloud<pcl::PointXYZ>::Ptr clt_cloud(new pcl::PointCloud < pcl::PointXYZ>);
  1732. pcl::copyPointCloud(*in_cloud, cluster_idxs, *clt_cloud);
  1733. //计算每一个茎的外接矩形
  1734. pcl::PointXYZ min_point_aabb;
  1735. pcl::PointXYZ max_point_aabb;
  1736. pcl::MomentOfInertiaEstimation<pcl::PointXYZ> feature_extractor;
  1737. feature_extractor.setInputCloud(clt_cloud);
  1738. feature_extractor.compute();
  1739. feature_extractor.getAABB(min_point_aabb, max_point_aabb);
  1740. //扩展矩形范围
  1741. pcl::PointXYZ min_point_aabb_ex;
  1742. pcl::PointXYZ max_point_aabb_ex;
  1743. min_point_aabb_ex.x = 0.5*(min_point_aabb.x + max_point_aabb.x) - hole_step_radius;
  1744. min_point_aabb_ex.y = m_cparam.sc_grab_ymin;
  1745. if (m_dtype != 0) { min_point_aabb_ex.y = m_cparam.rs_grab_ymin; }
  1746. min_point_aabb_ex.z = 0.5*(min_point_aabb.z + max_point_aabb.z) - hole_step_radius;
  1747. max_point_aabb_ex.x = 0.5*(min_point_aabb.x + max_point_aabb.x) + hole_step_radius;
  1748. max_point_aabb_ex.y = m_cparam.sc_grab_ymax;
  1749. if (m_dtype != 0) { max_point_aabb_ex.y = m_cparam.rs_grab_ymax; }
  1750. max_point_aabb_ex.z = 0.5*(min_point_aabb.z + max_point_aabb.z) + hole_step_radius;
  1751. /////////////////////////////////
  1752. //扩展矩形内直线检测
  1753. std::vector<int>line_idx;
  1754. find_line_in_cube(in_cloud, min_point_aabb_ex, max_point_aabb_ex, line_idx);
  1755. clt_line_idx.push_back(line_idx);
  1756. //找到line_idx中y最小的那个点的idx
  1757. int line_root_idx = -1;
  1758. float line_root_y_min = 10000.0;
  1759. for (size_t lidx = 0; lidx < line_idx.size(); ++lidx) {
  1760. int raw_idx = line_idx.at(lidx);
  1761. if (in_cloud->at(raw_idx).y < line_root_y_min) {
  1762. line_root_y_min = in_cloud->at(raw_idx).y;
  1763. line_root_idx = raw_idx;
  1764. }
  1765. }
  1766. //找到底部中心点,然后找到根部坐标
  1767. pcl::PointXYZ pt_root;
  1768. pt_root = in_cloud->at(line_root_idx);
  1769. clt_root.push_back(pt_root);
  1770. clt_box.push_back(min_point_aabb_ex);
  1771. clt_box.push_back(max_point_aabb_ex);
  1772. //viewer_cloud(clt_cloud, std::string("cluster"));
  1773. }
  1774. if (m_cparam.image_show) {
  1775. viewer_cloud_cluster_box(in_cloud, clt_root, clt_box, clt_line_idx, std::string("cloud2d_cluster_box"));
  1776. }
  1777. //sort cluster center, get the frontest leftest one
  1778. std::vector<float> cluster_index;
  1779. for (auto&pt : clt_root) {
  1780. float idx = pt.x - pt.z;
  1781. cluster_index.push_back(idx);
  1782. }
  1783. int first_idx = std::min_element(cluster_index.begin(), cluster_index.end()) - cluster_index.begin();
  1784. if (m_dtype == 0) {
  1785. first_idx = std::max_element(cluster_index.begin(), cluster_index.end()) - cluster_index.begin();
  1786. }
  1787. first_seedling_cloud_idx.clear();
  1788. for (auto&v : clt_line_idx.at(first_idx)) {
  1789. first_seedling_cloud_idx.push_back(v);
  1790. }
  1791. xz_center.x = clt_root.at(first_idx).x;
  1792. xz_center.y = clt_root.at(first_idx).y;
  1793. xz_center.z = clt_root.at(first_idx).z;
  1794. if (m_pLogger) {
  1795. stringstream buff;
  1796. buff << m_pcdId <<": euclidean_clustering_ttsas, find cluster center(" << xz_center.x
  1797. << ", " << xz_center.y << ", " << xz_center.z << ")";
  1798. m_pLogger->INFO(buff.str());
  1799. }
  1800. }
  1801. void CRootStockGrabPoint::find_line_in_cube(
  1802. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //输入,整体点云
  1803. pcl::PointXYZ& min_pt_ex, //输入,
  1804. pcl::PointXYZ& max_pt_ex, //输入,
  1805. std::vector<int>& out_idx //输出,直线上点的index, 基于输入整体点云
  1806. )
  1807. {
  1808. out_idx.clear();
  1809. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  1810. pcl::CropBox<pcl::PointXYZ> box_filter;
  1811. std::vector<int>inbox_idx;
  1812. box_filter.setMin(Eigen::Vector4f(min_pt_ex.x, min_pt_ex.y, min_pt_ex.z, 1));
  1813. box_filter.setMax(Eigen::Vector4f(max_pt_ex.x, max_pt_ex.y, max_pt_ex.z, 1));
  1814. box_filter.setNegative(false);
  1815. box_filter.setInputCloud(in_cloud);
  1816. box_filter.filter(inbox_idx);
  1817. pcl::copyPointCloud(*in_cloud, inbox_idx, *cloud_inbox);
  1818. //找到inbox点云中的直线
  1819. float stem_radius = m_cparam.rs_grab_stem_diameter / 2.0;
  1820. if (m_dtype == 0) {
  1821. stem_radius = m_cparam.sc_grab_stem_diameter / 2.0;
  1822. }
  1823. pcl::ModelCoefficients::Ptr coefficients(new pcl::ModelCoefficients);
  1824. pcl::PointIndices::Ptr inliers(new pcl::PointIndices);
  1825. pcl::SACSegmentation<pcl::PointXYZ> seg;
  1826. seg.setOptimizeCoefficients(true);
  1827. seg.setModelType(pcl::SACMODEL_LINE);
  1828. seg.setDistanceThreshold(stem_radius);
  1829. seg.setInputCloud(cloud_inbox);
  1830. seg.segment(*inliers, *coefficients);
  1831. int idx_raw = 0;
  1832. for (auto& idx : inliers->indices) {
  1833. idx_raw = inbox_idx.at(idx);
  1834. out_idx.push_back(idx_raw);
  1835. }
  1836. if (m_cparam.image_show) {
  1837. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_line(new pcl::PointCloud<pcl::PointXYZ>);
  1838. pcl::copyPointCloud(*cloud_inbox, *inliers, *cloud_line);
  1839. pcl::visualization::PCLVisualizer::Ptr viewer(new pcl::visualization::PCLVisualizer("line in cluster"));
  1840. viewer->addPointCloud<pcl::PointXYZ>(cloud_inbox, "cluster");
  1841. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "cluster");
  1842. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "cluster");
  1843. viewer->addPointCloud(cloud_line, "fit_line");
  1844. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0, 1, 0, "fit_line");
  1845. viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, "fit_line");
  1846. while (!viewer->wasStopped()) {
  1847. viewer->spinOnce(100);
  1848. }
  1849. }
  1850. }
  1851. void CRootStockGrabPoint::crop_nn_analysis(
  1852. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  1853. pcl::PointCloud<pcl::PointXYZ>::Ptr seed_cloud,
  1854. double dist_mean,
  1855. std::vector<double>&mass_indices,
  1856. std::vector<int>& candidate_idx
  1857. )
  1858. {
  1859. candidate_idx.clear();
  1860. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  1861. pcl::CropBox<pcl::PointXYZ> box_filter;
  1862. box_filter.setNegative(false);
  1863. box_filter.setInputCloud(in_cloud);
  1864. double radius = 5;
  1865. double rx = 0.8;
  1866. double ry = 1.0;
  1867. double rz = 0.8;
  1868. double cx, cy, cz;
  1869. double dz;
  1870. for (size_t i = 0; i < seed_cloud->points.size(); ++i) {
  1871. cx = seed_cloud->points.at(i).x;
  1872. cy = seed_cloud->points.at(i).y;
  1873. cz = seed_cloud->points.at(i).z;
  1874. box_filter.setMin(Eigen::Vector4f(cx - rx*radius, cy - ry*radius, cz - rz*radius, 1));
  1875. box_filter.setMax(Eigen::Vector4f(cx + rx*radius, cy + ry*radius, cz + rz*radius, 1));
  1876. box_filter.filter(*cloud_inbox);
  1877. //dz
  1878. Eigen::Vector4f min_point;
  1879. Eigen::Vector4f max_point;
  1880. pcl::getMinMax3D(*cloud_inbox, min_point, max_point);
  1881. dz = max_point(2) - min_point(2);
  1882. //project to xy-plane
  1883. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox_xy(new pcl::PointCloud<pcl::PointXYZ>);
  1884. pcl::copyPointCloud(*cloud_inbox, *cloud_inbox_xy);
  1885. for (auto&pt : *cloud_inbox_xy) { pt.z = 0.0; }
  1886. //pca
  1887. double dx_obb;
  1888. double dy_obb;
  1889. double angle_obb;
  1890. cal_obb_2d(cloud_inbox_xy, 0, dx_obb, dy_obb, angle_obb);
  1891. try {
  1892. double dx_grid = dx_obb / dist_mean;
  1893. double dy_grid = dy_obb / dist_mean;
  1894. double dz_grid = dz / dist_mean;
  1895. double fill_ratio = cloud_inbox->points.size() / dx_grid / dy_grid / dz_grid;
  1896. double y_ratio = dy_obb / dx_obb/2;
  1897. y_ratio = pow(y_ratio, 2);
  1898. double a_ratio = cos((angle_obb - 90)*PI / 180.0);
  1899. a_ratio = pow(a_ratio, 3);
  1900. double mass_index = fill_ratio*y_ratio*a_ratio;
  1901. mass_indices.push_back(mass_index);
  1902. if (m_pLogger) {
  1903. stringstream buff;
  1904. buff << i << "/" << seed_cloud->points.size() << " dx=" << dx_obb << ", dy=" << dy_obb << ", fill_ratio=" << fill_ratio
  1905. << ", y_ratio=" << y_ratio << ", a_ratio=" << a_ratio << ", mass_index=" << mass_index;
  1906. m_pLogger->INFO(buff.str());
  1907. }
  1908. }
  1909. catch (...) {
  1910. mass_indices.push_back(0);
  1911. }
  1912. }
  1913. // sort by mass_indices
  1914. std::vector<size_t> sort_idx = sort_indexes_e(mass_indices, false);
  1915. for (auto&v : sort_idx) { candidate_idx.push_back((int)v); }
  1916. }
  1917. void CRootStockGrabPoint::euclidean_clustering_ttsas(
  1918. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  1919. double d1, double d2,
  1920. std::vector<pcl::PointXYZ>&cluster_center,
  1921. std::vector<std::vector<int>> &clustr_member
  1922. )
  1923. {
  1924. if (m_pLogger) {
  1925. stringstream buff;
  1926. buff << m_pcdId <<": euclidean_clustering_ttsas() begin...";
  1927. m_pLogger->INFO(buff.str());
  1928. }
  1929. std::vector<int> cluster_weight;
  1930. std::vector<int> data_stat;
  1931. std::vector<pcl::PointXYZ>cluster_center_raw;
  1932. std::vector<std::vector<int>> clustr_member_raw;
  1933. for (size_t i = 0; i < in_cloud->points.size(); ++i) { data_stat.push_back(0); }
  1934. size_t data_len = in_cloud->points.size();
  1935. int exists_change = 0;
  1936. int prev_change = 0;
  1937. int cur_change = 0;
  1938. int m = 0;
  1939. while (std::find(data_stat.begin(), data_stat.end(), 0) != data_stat.end()) {
  1940. bool new_while_first = true;
  1941. for (size_t i = 0; i < data_len; ++i) {
  1942. if (data_stat.at(i) == 0 && new_while_first && exists_change == 0) {
  1943. new_while_first = false;
  1944. std::vector<int> idx;
  1945. idx.push_back(i);
  1946. clustr_member_raw.push_back(idx);
  1947. pcl::PointXYZ center;
  1948. center.x = in_cloud->points.at(i).x;
  1949. center.y = in_cloud->points.at(i).y;
  1950. center.z = in_cloud->points.at(i).z;
  1951. cluster_center_raw.push_back(center);
  1952. data_stat.at(i) = 1;
  1953. cur_change += 1;
  1954. cluster_weight.push_back(1);
  1955. m += 1;
  1956. }
  1957. else if (data_stat[i] == 0) {
  1958. std::vector<float> distances;
  1959. for (size_t j = 0; j < clustr_member_raw.size(); ++j) {
  1960. std::vector<float> distances_sub;
  1961. for (size_t jj = 0; jj < clustr_member_raw.at(j).size(); ++jj) {
  1962. size_t ele_idx = clustr_member_raw.at(j).at(jj);
  1963. double d = sqrt(
  1964. (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) +
  1965. (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) +
  1966. (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));
  1967. distances_sub.push_back(d);
  1968. }
  1969. double min_dist = *std::min_element(distances_sub.begin(), distances_sub.end());
  1970. distances.push_back(min_dist);
  1971. }
  1972. int min_idx = std::min_element(distances.begin(), distances.end()) - distances.begin();
  1973. if (distances.at(min_idx) < d1) {
  1974. data_stat.at(i) = 1;
  1975. double w = cluster_weight.at(min_idx);
  1976. cluster_weight.at(min_idx) += 1;
  1977. clustr_member_raw.at(min_idx).push_back(i);
  1978. cluster_center_raw.at(min_idx).x = (cluster_center_raw.at(min_idx).x * w + in_cloud->points.at(i).x) / (w + 1);
  1979. cluster_center_raw.at(min_idx).y = (cluster_center_raw.at(min_idx).y * w + in_cloud->points.at(i).y) / (w + 1);
  1980. cluster_center_raw.at(min_idx).z = (cluster_center_raw.at(min_idx).z * w + in_cloud->points.at(i).z) / (w + 1);
  1981. cur_change += 1;
  1982. }
  1983. else if (distances.at(min_idx) > d2) {
  1984. std::vector<int> idx;
  1985. idx.push_back(i);
  1986. clustr_member_raw.push_back(idx);
  1987. pcl::PointXYZ center;
  1988. center.x = in_cloud->points.at(i).x;
  1989. center.y = in_cloud->points.at(i).y;
  1990. center.z = in_cloud->points.at(i).z;
  1991. cluster_center_raw.push_back(center);
  1992. data_stat.at(i) = 1;
  1993. cur_change += 1;
  1994. cluster_weight.push_back(1);
  1995. m += 1;
  1996. }
  1997. }
  1998. else if (data_stat.at(i)== 1) {
  1999. cur_change += 1;
  2000. }
  2001. }
  2002. exists_change = fabs(cur_change - prev_change);
  2003. prev_change = cur_change;
  2004. cur_change = 0;
  2005. }
  2006. // copy result
  2007. for (size_t i = 0; i < clustr_member_raw.size(); ++i) {
  2008. if (clustr_member_raw.at(i).size() < 20) { continue; }
  2009. cluster_center.push_back(cluster_center_raw.at(i));
  2010. clustr_member.push_back(clustr_member_raw.at(i));
  2011. }
  2012. if (m_pLogger) {
  2013. stringstream buff;
  2014. buff << m_pcdId <<": euclidean_clustering_ttsas() end";
  2015. m_pLogger->INFO(buff.str());
  2016. }
  2017. }
  2018. void CRootStockGrabPoint::cal_obb_2d(
  2019. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud,
  2020. int axis, //0--xy, 1--zy
  2021. double &dx_obb,
  2022. double &dy_obb,
  2023. double &angle_obb
  2024. )
  2025. {
  2026. // asign value
  2027. Eigen::MatrixXd pts(in_cloud->points.size(), 2);
  2028. for (size_t i = 0; i < in_cloud->points.size(); ++i) {
  2029. if (axis == 0) {
  2030. pts(i, 0) = in_cloud->points.at(i).x;
  2031. }
  2032. else {
  2033. pts(i, 0) = in_cloud->points.at(i).z;
  2034. }
  2035. pts(i, 1) = in_cloud->points.at(i).y;
  2036. }
  2037. // centerlize
  2038. Eigen::MatrixXd mu = pts.colwise().mean();
  2039. Eigen::RowVectorXd mu_row = mu;
  2040. pts.rowwise() -= mu_row;
  2041. //calculate covariance
  2042. Eigen::MatrixXd C = pts.adjoint()*pts;
  2043. C = C.array() / (pts.cols() - 1);
  2044. //compute eig
  2045. Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(C);
  2046. Eigen::MatrixXd d = eig.eigenvalues();
  2047. Eigen::MatrixXd v = eig.eigenvectors();
  2048. //projection
  2049. Eigen::MatrixXd p = pts * v;
  2050. Eigen::VectorXd minv = p.colwise().minCoeff();
  2051. Eigen::VectorXd maxv = p.colwise().maxCoeff();
  2052. Eigen::VectorXd range = maxv - minv;
  2053. dy_obb = range(1);
  2054. dx_obb = range(0);
  2055. angle_obb = std::atan2(v(1, 1), v(0, 1));
  2056. if (angle_obb < 0) { angle_obb = PI + angle_obb; }
  2057. angle_obb = angle_obb * 180 / PI;
  2058. }
  2059. void CRootStockGrabPoint::get_optimal_seed_compare(
  2060. pcl::PointCloud<pcl::PointXYZ>::Ptr in_cloud, //input 全部有效点云
  2061. pcl::PointCloud<pcl::PointXYZ>::Ptr in_line_cloud, //input 茎上直线点云
  2062. pcl::ModelCoefficients::Ptr line_model, //input
  2063. pcl::PointXYZ&pt, // 返回抓取的点坐标,基于pt_ref的偏移点
  2064. pcl::PointXYZ &pt_ref, // 返回点茎节的参考点
  2065. float& stem_width_mu, // 茎宽度,直径
  2066. float& stem_deflection //茎的最大挠度,最大弯曲处的偏离直径轴心的距离,毫米
  2067. )
  2068. {
  2069. // stem_deflection 的计算方法:计算点云到拟合轴线的平均距离
  2070. float th_ratio = 1.5; //原始点云和直线点云邻域尺寸增加原来的20%,不能当做抓取点
  2071. pt.x = -1000.0;
  2072. pt.y = -1000.0;
  2073. pt.z = -1000.0;
  2074. stem_width_mu = 0.0;
  2075. stem_deflection = 0.0;
  2076. float ymin, ymax;
  2077. ymin = m_cparam.rs_grab_ymin;
  2078. ymax = m_cparam.rs_grab_ymax;
  2079. if (m_dtype == 0) {
  2080. ymin = m_cparam.sc_grab_ymin;
  2081. ymax = m_cparam.sc_grab_ymax;
  2082. }
  2083. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  2084. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox_line(new pcl::PointCloud<pcl::PointXYZ>);
  2085. std::vector<float> stem_width;
  2086. std::vector<pcl::PointXYZ>online_points;
  2087. pcl::CropBox<pcl::PointXYZ> box_filter;
  2088. box_filter.setNegative(false);
  2089. box_filter.setInputCloud(in_cloud);
  2090. pcl::CropBox<pcl::PointXYZ> box_filter_line;
  2091. box_filter_line.setNegative(false);
  2092. box_filter_line.setInputCloud(in_line_cloud);
  2093. float radius = m_cparam.rs_grab_stem_diameter;
  2094. float cx, cy, cz, t;
  2095. float rx = 1.5;
  2096. float ry = 1.5;
  2097. float rz = 1.5;
  2098. float dz, dx;
  2099. // 如果opt_y_valid==false,就是用户没有指定抓取的y高度
  2100. float max_dist_to_boundary = 0.0;
  2101. int max_idx_to_boundary = -1;
  2102. find_fork(in_line_cloud, max_dist_to_boundary, max_idx_to_boundary);
  2103. stem_deflection = static_cast<float>(calculate_deflection(in_line_cloud, line_model));
  2104. cy = ymin;
  2105. //计算茎上的直径
  2106. while(cy<ymax){
  2107. t = (cy - line_model->values.at(1)) / line_model->values.at(4);
  2108. cx = line_model->values.at(3) * t + line_model->values.at(0);
  2109. cz = line_model->values.at(5) * t + line_model->values.at(2);
  2110. //////////////////////////////////////////////////////////////////
  2111. //原始点云,获取指定区域的dx,dz
  2112. box_filter.setMin(Eigen::Vector4f(cx - rx*radius, cy - 2, cz - rz*radius, 1));
  2113. box_filter.setMax(Eigen::Vector4f(cx + rx*radius, cy + 2, cz + rz*radius, 1));
  2114. box_filter.filter(*cloud_inbox);
  2115. if (cloud_inbox->points.size() > 10) {
  2116. Eigen::Vector4f min_point;
  2117. Eigen::Vector4f max_point;
  2118. pcl::getMinMax3D(*cloud_inbox, min_point, max_point);
  2119. dx = max_point(0) - min_point(0);
  2120. stem_width.push_back(dx);
  2121. }
  2122. else {
  2123. stem_width.push_back(0);
  2124. }
  2125. pcl::PointXYZ tmp_pt(cx,cy,cz);
  2126. online_points.push_back(tmp_pt);
  2127. cy += 1.0;
  2128. }
  2129. //得到有效直径数值,并计算均值、标准差
  2130. std::vector<float>valid_stem_width;
  2131. for (auto&w : stem_width) {
  2132. if (w > 0) {
  2133. valid_stem_width.push_back(w);
  2134. }
  2135. }
  2136. float mu = get_hist_mean(valid_stem_width);
  2137. stem_width_mu = mu;
  2138. float stdv = get_hist_std(valid_stem_width, mu);
  2139. //0)计算点云到轴线的最大距离
  2140. //1) original max position
  2141. int max_pos = std::max_element(stem_width.begin(), stem_width.end()) - stem_width.begin();
  2142. int max_pos_ref = max_pos;
  2143. //3 如果fork方法得到位置信息,进行更新,并记录历史
  2144. if (max_idx_to_boundary >= 0) {
  2145. max_pos = int(in_line_cloud->points.at(max_idx_to_boundary).y + 0.5 - ymin);
  2146. max_pos_ref = max_pos;
  2147. }
  2148. //4 用茎节上下限高度值进行约束,如果超限,用最低限高度作为茎节高度
  2149. double grab_fork_yup = m_cparam.rs_grab_fork_yup;
  2150. double grab_fork_ybt = m_cparam.rs_grab_fork_ybt;
  2151. if (m_dtype == 0) {
  2152. grab_fork_yup = m_cparam.sc_grab_fork_yup;
  2153. grab_fork_ybt = m_cparam.sc_grab_fork_ybt;
  2154. }
  2155. bool out_of_range = false;
  2156. if ((max_pos + ymin) > grab_fork_yup || (max_pos + ymin) < grab_fork_ybt) {
  2157. out_of_range = true;
  2158. int original_max_pos = max_pos;
  2159. max_pos = int(grab_fork_ybt - ymin + 0.5);
  2160. max_pos_ref = max_pos;
  2161. if (m_pLogger) {
  2162. stringstream buff;
  2163. buff << m_pcdId << ": warning,self fork postiont = " << original_max_pos <<
  2164. ", USE bottom limit fork postiont " << max_pos <<
  2165. ", valid fork postiont range:[" << int(grab_fork_ybt - ymin + 0.5) <<
  2166. ", " << int(grab_fork_yup - ymin + 0.5) << "]";
  2167. m_pLogger->INFO(buff.str());
  2168. }
  2169. }
  2170. else {
  2171. if (m_pLogger) {
  2172. stringstream buff;
  2173. buff << m_pcdId << ": self fork postiont = " << max_pos <<
  2174. ", valid fork postiont range:[" << int(grab_fork_ybt - ymin + 0.5)<<
  2175. ", "<< int(grab_fork_yup - ymin + 0.5) << "]";
  2176. m_pLogger->INFO(buff.str());
  2177. }
  2178. }
  2179. //5 按指定量偏移
  2180. if(m_dtype == 0){
  2181. max_pos_ref = max_pos;
  2182. max_pos += static_cast<int>(m_cparam.sc_grab_offset);
  2183. }
  2184. else{
  2185. max_pos_ref = max_pos;
  2186. max_pos += static_cast<int>(m_cparam.rs_grab_offset);
  2187. }
  2188. if (max_pos < 0) { max_pos = 0; }
  2189. if (max_pos >= online_points.size()) { max_pos = online_points.size() - 1; }
  2190. if (out_of_range) {
  2191. max_pos_ref = max_pos;
  2192. }
  2193. /////////////////////////////////////////////////////////////////////
  2194. //直线点云,获取指定区域的dx_line,dz_line
  2195. Eigen::Vector4f min_point;
  2196. Eigen::Vector4f max_point;
  2197. cx = online_points.at(max_pos).x;
  2198. cy = online_points.at(max_pos).y;
  2199. cz = online_points.at(max_pos).z;
  2200. box_filter_line.setMin(Eigen::Vector4f(cx - rx*radius, cy - 2, cz - rz*radius, 1));
  2201. box_filter_line.setMax(Eigen::Vector4f(cx + rx*radius, cy + 2, cz + rz*radius, 1));
  2202. box_filter_line.filter(*cloud_inbox_line);
  2203. float z_mu = cz;
  2204. float x_mu = cx;
  2205. if (cloud_inbox_line->points.size() > 5) {
  2206. pcl::getMinMax3D(*cloud_inbox_line, min_point, max_point);
  2207. z_mu = 0.5 * (max_point(2) + min_point(2));
  2208. x_mu = 0.5 * (max_point(0) + min_point(0));
  2209. }
  2210. pt.x = x_mu;
  2211. pt.y = cy;
  2212. pt.z = z_mu;
  2213. pt_ref.x = online_points.at(max_pos_ref).x;
  2214. pt_ref.y = online_points.at(max_pos_ref).y;
  2215. pt_ref.z = online_points.at(max_pos_ref).z;
  2216. //显示位置
  2217. if (m_cparam.image_show) {
  2218. pcl::visualization::PCLVisualizer viewer("grab line");
  2219. viewer.addCoordinateSystem();
  2220. viewer.addPointCloud<pcl::PointXYZ>(in_line_cloud, "stem_cloud");//????????????????
  2221. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "stem_cloud");
  2222. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "stem_cloud");
  2223. viewer.addLine(online_points.front(), online_points.back(), 1.0, 0.0, 0.0);
  2224. while (!viewer.wasStopped()) {
  2225. viewer.spinOnce(100);
  2226. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2227. }
  2228. }
  2229. }
  2230. //计算茎的弯曲度:点云到直线的距离95分位
  2231. double CRootStockGrabPoint::calculate_deflection(
  2232. pcl::PointCloud<pcl::PointXYZ>::Ptr in_line_cloud, //input
  2233. pcl::ModelCoefficients::Ptr line_model //input
  2234. )
  2235. {
  2236. double deflection = 0.0;
  2237. std::vector<double> dists;
  2238. Eigen::Vector4f line_pt = { line_model->values[0],line_model->values[1],line_model->values[2],0.0 };
  2239. Eigen::Vector4f line_dir = { line_model->values[3],line_model->values[4],line_model->values[5],0.0 };
  2240. for (auto ptc : in_line_cloud->points) {
  2241. //Get the square distance from a point to a line (represented by a point and a direction)
  2242. Eigen::Vector4f pt = { ptc.x, ptc.y, ptc.z, 0.0 };
  2243. double dist = (line_dir.cross3(line_pt - pt)).squaredNorm() / line_dir.squaredNorm();
  2244. dists.push_back(sqrt(dist));
  2245. }
  2246. //deflection = get_hist_mean(dists);
  2247. //95 percentile
  2248. if (dists.size() > 0) {
  2249. int idx_95 = int(0.95 * dists.size());
  2250. std::sort(dists.begin(), dists.end());
  2251. deflection = dists.at(idx_95);
  2252. }
  2253. return deflection;
  2254. }
  2255. void CRootStockGrabPoint::find_fork(
  2256. pcl::PointCloud<pcl::PointXYZ>::Ptr in_line_cloud,
  2257. float& max_dist,
  2258. int& max_idx
  2259. )
  2260. {
  2261. //1 project to xy plane
  2262. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_line_2d(new pcl::PointCloud<pcl::PointXYZ>);
  2263. pcl::copyPointCloud(*in_line_cloud, *cloud_line_2d);
  2264. for (pcl::PointXYZ&p : cloud_line_2d->points) { p.z = 0.0; }
  2265. //2 边界检测
  2266. pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree_ptr(new pcl::search::KdTree<pcl::PointXYZ>);
  2267. pcl::PointCloud<pcl::Normal>::Ptr normal(new pcl::PointCloud<pcl::Normal>);
  2268. pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
  2269. ne.setInputCloud(cloud_line_2d);
  2270. ne.setSearchMethod(kdtree_ptr);
  2271. ne.setRadiusSearch(3.0);
  2272. ne.compute(*normal);
  2273. pcl::PointCloud<pcl::Boundary>::Ptr boundaries(new pcl::PointCloud<pcl::Boundary>);
  2274. boundaries->resize(cloud_line_2d->size());
  2275. pcl::BoundaryEstimation < pcl::PointXYZ, pcl::Normal, pcl::Boundary>boundary_estimation;
  2276. boundary_estimation.setInputCloud(cloud_line_2d);
  2277. boundary_estimation.setInputNormals(normal);
  2278. boundary_estimation.setSearchMethod(kdtree_ptr);
  2279. boundary_estimation.setKSearch(30);
  2280. boundary_estimation.setAngleThreshold(M_PI*0.6);
  2281. boundary_estimation.compute(*boundaries);
  2282. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_lines(new pcl::PointCloud<pcl::PointXYZ>);
  2283. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_visual(new pcl::PointCloud<pcl::PointXYZRGB>);
  2284. cloud_visual->resize(cloud_line_2d->size());
  2285. for (size_t i = 0; i < cloud_line_2d->size(); ++i) {
  2286. cloud_visual->points.at(i).x = cloud_line_2d->points.at(i).x;
  2287. cloud_visual->points.at(i).y = cloud_line_2d->points.at(i).y;
  2288. cloud_visual->points.at(i).z = cloud_line_2d->points.at(i).z;
  2289. if (boundaries->points.at(i).boundary_point != 0) {
  2290. cloud_visual->points.at(i).r = 255;
  2291. cloud_visual->points.at(i).g = 0;
  2292. cloud_visual->points.at(i).b = 0;
  2293. cloud_lines->points.push_back(cloud_line_2d->points.at(i));
  2294. }
  2295. else {
  2296. cloud_visual->points.at(i).r = 255;
  2297. cloud_visual->points.at(i).g = 255;
  2298. cloud_visual->points.at(i).b = 255;
  2299. }
  2300. }
  2301. /*if (m_cparam.image_show) {
  2302. viewer_cloud(cloud_visual, string("boundary"));
  2303. viewer_cloud(cloud_lines, string("cloud_lines"));
  2304. }
  2305. */
  2306. //3 find points with max distance to boundaries point
  2307. pcl::search::KdTree<pcl::PointXYZ>::Ptr bound_kdtree_ptr(new pcl::search::KdTree<pcl::PointXYZ>);
  2308. bound_kdtree_ptr->setInputCloud(cloud_lines);
  2309. max_idx = -1;
  2310. max_dist = 0.0;
  2311. int k = 1;
  2312. for (int i = 0; i < cloud_line_2d->size(); ++i) {
  2313. if (boundaries->points.at(i).boundary_point != 0) { continue; }
  2314. std::vector<int> idx(k);
  2315. std::vector<float> sqr_distances(k);
  2316. int cnt = bound_kdtree_ptr->nearestKSearch(cloud_line_2d->points.at(i), k,idx,sqr_distances);
  2317. if (cnt > 0) {
  2318. if (sqr_distances.at(0) > max_dist) {
  2319. max_dist = sqr_distances.at(0);
  2320. max_idx = i;
  2321. }
  2322. }
  2323. }
  2324. if (m_cparam.image_show) {
  2325. pcl::visualization::PCLVisualizer viewer("boundary line");
  2326. viewer.addCoordinateSystem();
  2327. viewer.addPointCloud<pcl::PointXYZRGB>(cloud_visual, "boundary");
  2328. if (max_idx >= 0) {
  2329. pcl::PointXYZ p0,p1;
  2330. p0.x = cloud_line_2d->points.at(max_idx).x-5;
  2331. p0.y = cloud_line_2d->points.at(max_idx).y;
  2332. p0.z = cloud_line_2d->points.at(max_idx).z;
  2333. p1.x = cloud_line_2d->points.at(max_idx).x + 5;
  2334. p1.y = cloud_line_2d->points.at(max_idx).y;
  2335. p1.z = cloud_line_2d->points.at(max_idx).z;
  2336. viewer.addLine(p0, p1, 1.0, 0.0, 0.0);
  2337. }
  2338. while (!viewer.wasStopped()) {
  2339. viewer.spinOnce(100);
  2340. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2341. }
  2342. }
  2343. }
  2344. void CRootStockGrabPoint::get_line_project_hist(
  2345. pcl::PointCloud<pcl::PointXYZ>::Ptr in_line_cloud, //input 茎上直线点云
  2346. pcl::ModelCoefficients::Ptr line_model, //input
  2347. std::vector<int>& count_hist // 返回有效直线1mm内点云数量
  2348. )
  2349. {
  2350. count_hist.clear();
  2351. float ymin, ymax;
  2352. ymin = m_cparam.rs_grab_ymin;
  2353. ymax = m_cparam.rs_grab_ymax;
  2354. if (m_dtype == 0) {
  2355. ymin = m_cparam.sc_grab_ymin;
  2356. ymax = m_cparam.sc_grab_ymax;
  2357. }
  2358. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_inbox(new pcl::PointCloud<pcl::PointXYZ>);
  2359. pcl::CropBox<pcl::PointXYZ> box_filter_line;
  2360. box_filter_line.setNegative(false);
  2361. box_filter_line.setInputCloud(in_line_cloud);
  2362. float radius = m_cparam.rs_grab_stem_diameter;
  2363. if (m_dtype == 0) {
  2364. radius = m_cparam.sc_grab_stem_diameter;
  2365. }
  2366. float rx = 1.5;
  2367. float rz = 1.5;
  2368. float cx, cy, cz, t;
  2369. cy = ymin;
  2370. while (cy<ymax) {
  2371. t = (cy - line_model->values.at(1)) / line_model->values.at(4);
  2372. cx = line_model->values.at(3) * t + line_model->values.at(0);
  2373. cz = line_model->values.at(5) * t + line_model->values.at(2);
  2374. box_filter_line.setMin(Eigen::Vector4f(cx - rx*radius, cy - 0.5, cz - rz*radius, 1));
  2375. box_filter_line.setMax(Eigen::Vector4f(cx + rx*radius, cy + 0.5, cz + rz*radius, 1));
  2376. box_filter_line.filter(*cloud_inbox);
  2377. count_hist.push_back(cloud_inbox->points.size());
  2378. cy += 1.0;
  2379. }
  2380. }
  2381. void CRootStockGrabPoint::viewer_cloud(pcl::PointCloud<pcl::PointXYZ>::Ptr cloud, std::string&winname)
  2382. {
  2383. pcl::visualization::CloudViewer viewer(winname);
  2384. //viewer.runOnVisualizationThreadOnce(viewerOneOff);
  2385. viewer.showCloud(cloud);
  2386. while (!viewer.wasStopped()) {
  2387. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2388. }
  2389. }
  2390. void CRootStockGrabPoint::viewer_cloud(pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud, std::string&winname)
  2391. {
  2392. pcl::visualization::CloudViewer viewer(winname);
  2393. //viewer.runOnVisualizationThreadOnce(viewerOneOff);
  2394. viewer.showCloud(cloud);
  2395. while (!viewer.wasStopped()) {
  2396. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2397. }
  2398. }
  2399. void CRootStockGrabPoint::viewer_cloud_debug(
  2400. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud,
  2401. pcl::PointXYZ &p, //抓取点
  2402. pcl::PointXYZ &p_ref,//分叉点
  2403. pcl::PointXYZ &root_pt,
  2404. std::string&winname)
  2405. {
  2406. pcl::visualization::PCLVisualizer viewer(winname);
  2407. //viewer.runOnVisualizationThreadOnce(viewerOneOff);
  2408. viewer.addPointCloud(cloud);
  2409. viewer.addCoordinateSystem();
  2410. pcl::PointXYZ p0, x1, y1,p00,x0,y0, root_px0, root_px1, root_py0,root_py1;
  2411. p0.x = p.x;
  2412. p0.y = p.y;
  2413. p0.z = p.z;
  2414. x1.x = p.x + 100.0;
  2415. x1.y = p.y;
  2416. x1.z = p.z;
  2417. y1.x = p.x;
  2418. y1.y = p.y + 20.0;
  2419. y1.z = p.z;
  2420. p00.x = 0.0;
  2421. p00.y = 0.0;
  2422. p00.z = p.z;
  2423. x0.x = 600.0;
  2424. x0.y = 0;
  2425. x0.z = p.z;
  2426. y0.x = 0.0;
  2427. y0.y = 300.0;
  2428. y0.z = p.z;
  2429. root_px0.x = root_pt.x - 5.0;
  2430. root_px0.y = root_pt.y;
  2431. root_px0.z = root_pt.z;
  2432. root_px1.x = root_pt.x + 5.0;
  2433. root_px1.y = root_pt.y;
  2434. root_px1.z = root_pt.z;
  2435. root_py0.x = root_pt.x;
  2436. root_py0.y = root_pt.y;
  2437. root_py0.z = root_pt.z - 5.0;
  2438. root_py1.x = root_pt.x;
  2439. root_py1.y = root_pt.y;
  2440. root_py1.z = root_pt.z + 5.0;
  2441. //茎节点
  2442. pcl::PointXYZ fork_px0, fork_px1, fork_py0, fork_py1;
  2443. fork_px0.x = p_ref.x - 5.0;
  2444. fork_px0.y = p_ref.y;
  2445. fork_px0.z = p_ref.z;
  2446. fork_px1.x = p_ref.x + 5.0;
  2447. fork_px1.y = p_ref.y;
  2448. fork_px1.z = p_ref.z;
  2449. fork_py0.x = p_ref.x;
  2450. fork_py0.y = p_ref.y;
  2451. fork_py0.z = p_ref.z - 5.0;
  2452. fork_py1.x = p_ref.x;
  2453. fork_py1.y = p_ref.y;
  2454. fork_py1.z = p_ref.z + 5.0;
  2455. viewer.addLine(p0, x1, 255, 0, 0, "x");
  2456. viewer.addLine(p0, y1, 0, 255, 0, "y");
  2457. viewer.addLine(p00, x0, 255, 0, 0, "x0");
  2458. viewer.addLine(p00, y0, 0, 255, 0, "y0");
  2459. viewer.addLine(root_px0, root_px1, 255, 0, 0, "rootx");
  2460. viewer.addLine(root_py0, root_py1, 0, 255, 0, "rooty");
  2461. viewer.addLine(fork_px0, fork_px1, 255, 0, 0, "forkx");
  2462. viewer.addLine(fork_py0, fork_py1, 0, 255, 0, "forky");
  2463. while (!viewer.wasStopped()) {
  2464. viewer.spinOnce(100);
  2465. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2466. }
  2467. }
  2468. void CRootStockGrabPoint::viewer_cloud_cluster(
  2469. pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud,
  2470. std::vector<pcl::PointXYZ>cluster_center,
  2471. std::string&winname)
  2472. {
  2473. pcl::visualization::PCLVisualizer viewer(winname);
  2474. viewer.addPointCloud(cloud);
  2475. viewer.addCoordinateSystem();
  2476. int cnt = 0;
  2477. char buf[8];
  2478. for (auto& p : cluster_center) {
  2479. pcl::PointXYZ px0, px1, py1, py0;
  2480. px0.x = p.x;
  2481. px0.y = p.y;
  2482. px0.z = p.z;
  2483. px1.x = p.x + 10.0;
  2484. px1.y = p.y;
  2485. px1.z = p.z;
  2486. py0.x = p.x;
  2487. py0.y = p.y;
  2488. py0.z = p.z;
  2489. py1.x = p.x;
  2490. py1.y = p.y + 10.0;
  2491. py1.z = p.z;
  2492. viewer.addLine(px0, px1, 255, 0, 0, "x" + string(_itoa(cnt, buf,10)));
  2493. viewer.addLine(py0, py1, 0, 255, 0, "y" + string(_itoa(cnt, buf, 10)));
  2494. cnt += 1;
  2495. }
  2496. while (!viewer.wasStopped()) {
  2497. viewer.spinOnce(100);
  2498. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2499. }
  2500. }
  2501. void CRootStockGrabPoint::viewer_cloud_cluster_box(
  2502. pcl::PointCloud<pcl::PointXYZ>::Ptr cloud,
  2503. std::vector<pcl::PointXYZ>&cluster_center,
  2504. std::vector<pcl::PointXYZ>&cluster_box,
  2505. std::vector<std::vector<int> >& clt_line_idx,
  2506. std::string&winname)
  2507. {
  2508. pcl::visualization::PCLVisualizer viewer(winname);
  2509. viewer.addCoordinateSystem();
  2510. viewer.addPointCloud<pcl::PointXYZ>(cloud, "all_cloud");
  2511. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 1, 1, "all_cloud");
  2512. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "all_cloud");
  2513. int cnt = 0;
  2514. char buf[8];
  2515. for (size_t i = 0; i < cluster_center.size();++i) {
  2516. pcl::PointXYZ &p = cluster_center.at(i);
  2517. pcl::PointXYZ px0, px1, py1, py0;
  2518. px0.x = p.x;
  2519. px0.y = p.y;
  2520. px0.z = p.z;
  2521. px1.x = p.x + 10.0;
  2522. px1.y = p.y;
  2523. px1.z = p.z;
  2524. py0.x = p.x;
  2525. py0.y = p.y;
  2526. py0.z = p.z;
  2527. py1.x = p.x;
  2528. py1.y = p.y + 10.0;
  2529. py1.z = p.z;
  2530. viewer.addLine(px0, px1, 255, 0, 0, "x" + string(_itoa(cnt, buf, 10)));
  2531. viewer.addLine(py0, py1, 0, 255, 0, "y" + string(_itoa(cnt, buf, 10)));
  2532. //box
  2533. pcl::PointXYZ & min_pt = cluster_box.at(2 * i);
  2534. pcl::PointXYZ & max_pt = cluster_box.at(2 * i + 1);
  2535. 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)));
  2536. viewer.setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_REPRESENTATION,
  2537. pcl::visualization::PCL_VISUALIZER_REPRESENTATION_WIREFRAME, "AABB_" + string(_itoa(cnt, buf, 10)));
  2538. if (clt_line_idx.size()>i && clt_line_idx.at(i).size() > 0) {
  2539. pcl::PointCloud<pcl::PointXYZ>::Ptr line_cloud(new pcl::PointCloud<pcl::PointXYZ>);
  2540. pcl::copyPointCloud(*cloud, clt_line_idx.at(i), *line_cloud);
  2541. viewer.addPointCloud(line_cloud, "fit_line" + string(_itoa(cnt, buf, 10)));
  2542. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0, 1, 0, "fit_line" + string(_itoa(cnt, buf, 10)));
  2543. viewer.setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 3, "fit_line" + string(_itoa(cnt, buf, 10)));
  2544. }
  2545. cnt += 1;
  2546. }
  2547. while (!viewer.wasStopped()) {
  2548. viewer.spinOnce(100);
  2549. boost::this_thread::sleep(boost::posix_time::microseconds(100000));
  2550. }
  2551. }
  2552. };