00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include <pcl/common/common.h>
00041 #include <pcl/console/time.h>
00042 #include <pcl/common/angles.h>
00043 #include <pcl/console/print.h>
00044 #include <pcl/console/parse.h>
00045
00046 #include <pcl/io/pcd_io.h>
00047
00048 #include <pcl/features/integral_image_normal.h>
00049 #include <pcl/features/normal_3d.h>
00050
00051 #include <pcl/filters/extract_indices.h>
00052 #include <pcl/filters/project_inliers.h>
00053
00054 #include <pcl/segmentation/organized_multi_plane_segmentation.h>
00055 #include <pcl/segmentation/euclidean_cluster_comparator.h>
00056 #include <pcl/segmentation/organized_connected_component_segmentation.h>
00057 #include <pcl/segmentation/edge_aware_plane_comparator.h>
00058 #include <pcl/segmentation/sac_segmentation.h>
00059 #include <pcl/segmentation/extract_polygonal_prism_data.h>
00060 #include <pcl/segmentation/extract_clusters.h>
00061
00062 #include <pcl/surface/convex_hull.h>
00063
00064 #include <pcl/geometry/polygon_operations.h>
00065
00066 #include <pcl/visualization/point_cloud_handlers.h>
00067 #include <pcl/visualization/pcl_visualizer.h>
00068 #include <pcl/visualization/image_viewer.h>
00069
00070 using namespace pcl;
00071 using namespace pcl::console;
00072 using namespace std;
00073
00075 template <typename PointT>
00076 class ObjectSelection
00077 {
00078 public:
00079 ObjectSelection ()
00080 : plane_comparator_ (new EdgeAwarePlaneComparator<PointT, Normal>)
00081 , rgb_data_ ()
00082 {
00083
00084 plane_comparator_->setDistanceThreshold (0.01f, false);
00085 }
00086
00088 virtual ~ObjectSelection ()
00089 {
00090 if (rgb_data_)
00091 delete[] rgb_data_;
00092 }
00093
00095 void
00096 estimateNormals (const typename PointCloud<PointT>::ConstPtr &input, PointCloud<Normal> &normals)
00097 {
00098 if (input->isOrganized ())
00099 {
00100 IntegralImageNormalEstimation<PointT, Normal> ne;
00101
00102 ne.setNormalEstimationMethod (ne.COVARIANCE_MATRIX);
00103 ne.setMaxDepthChangeFactor (0.02f);
00104 ne.setNormalSmoothingSize (20.0f);
00105
00106 ne.setInputCloud (input);
00107 ne.compute (normals);
00108
00109
00110 float *map=ne.getDistanceMap ();
00111 distance_map_.assign(map, map+input->size() );
00112 plane_comparator_->setDistanceMap(distance_map_.data());
00113 }
00114 else
00115 {
00116 NormalEstimation<PointT, Normal> ne;
00117 ne.setInputCloud (input);
00118 ne.setRadiusSearch (0.02f);
00119 ne.setSearchMethod (search_);
00120 ne.compute (normals);
00121 }
00122 }
00123
00125 void
00126 keyboard_callback (const visualization::KeyboardEvent&, void*)
00127 {
00128
00129
00130
00131
00132
00133
00134
00135
00136 }
00137
00139 void
00140 mouse_callback (const visualization::MouseEvent& mouse_event, void*)
00141 {
00142 if (mouse_event.getType() == visualization::MouseEvent::MouseButtonPress && mouse_event.getButton() == visualization::MouseEvent::LeftButton)
00143 {
00144 cout << "left button pressed @ " << mouse_event.getX () << " , " << mouse_event.getY () << endl;
00145 }
00146 }
00147
00149
00157 void
00158 segmentObject (int picked_idx,
00159 const typename PointCloud<PointT>::ConstPtr &cloud,
00160 const PointIndices::Ptr &plane_indices,
00161 PointCloud<PointT> &object)
00162 {
00163 typename PointCloud<PointT>::Ptr plane_hull (new PointCloud<PointT>);
00164
00165
00166 ConvexHull<PointT> chull;
00167 chull.setDimension (2);
00168 chull.setInputCloud (cloud);
00169 chull.setIndices (plane_indices);
00170 chull.reconstruct (*plane_hull);
00171
00172
00173 typename PointCloud<PointT>::Ptr plane (new PointCloud<PointT>);
00174 ExtractIndices<PointT> extract (true);
00175 extract.setInputCloud (cloud);
00176 extract.setIndices (plane_indices);
00177 extract.setNegative (false);
00178 extract.filter (*plane);
00179 PointIndices::Ptr indices_but_the_plane (new PointIndices);
00180 extract.getRemovedIndices (*indices_but_the_plane);
00181
00182
00183 PointIndices::Ptr points_above_plane (new PointIndices);
00184 ExtractPolygonalPrismData<PointT> exppd;
00185 exppd.setInputCloud (cloud);
00186 exppd.setIndices (indices_but_the_plane);
00187 exppd.setInputPlanarHull (plane_hull);
00188 exppd.setViewPoint (cloud->points[picked_idx].x, cloud->points[picked_idx].y, cloud->points[picked_idx].z);
00189 exppd.setHeightLimits (0.001, 0.5);
00190 exppd.segment (*points_above_plane);
00191
00192 vector<PointIndices> euclidean_label_indices;
00193
00194 if (cloud_->isOrganized ())
00195 {
00196
00197 typename EuclideanClusterComparator<PointT, Normal, Label>::Ptr euclidean_cluster_comparator (new EuclideanClusterComparator<PointT, Normal, Label>);
00198 euclidean_cluster_comparator->setInputCloud (cloud);
00199 euclidean_cluster_comparator->setDistanceThreshold (0.03f, false);
00200
00201 Label l; l.label = 0;
00202 PointCloud<Label>::Ptr scene (new PointCloud<Label> (cloud->width, cloud->height, l));
00203
00204 for (int i = 0; i < static_cast<int> (points_above_plane->indices.size ()); ++i)
00205 scene->points[points_above_plane->indices[i]].label = 1;
00206 euclidean_cluster_comparator->setLabels (scene);
00207
00208 vector<bool> exclude_labels (2); exclude_labels[0] = true; exclude_labels[1] = false;
00209 euclidean_cluster_comparator->setExcludeLabels (exclude_labels);
00210
00211 OrganizedConnectedComponentSegmentation<PointT, Label> euclidean_segmentation (euclidean_cluster_comparator);
00212 euclidean_segmentation.setInputCloud (cloud);
00213
00214 PointCloud<Label> euclidean_labels;
00215 euclidean_segmentation.segment (euclidean_labels, euclidean_label_indices);
00216 }
00217 else
00218 {
00219 print_highlight (stderr, "Extracting individual clusters from the points above the reference plane ");
00220 TicToc tt; tt.tic ();
00221
00222 EuclideanClusterExtraction<PointT> ec;
00223 ec.setClusterTolerance (0.02);
00224 ec.setMinClusterSize (100);
00225 ec.setSearchMethod (search_);
00226 ec.setInputCloud (cloud);
00227 ec.setIndices (points_above_plane);
00228 ec.extract (euclidean_label_indices);
00229
00230 print_info ("[done, "); print_value ("%g", tt.toc ()); print_info (" ms : "); print_value ("%zu", euclidean_label_indices.size ()); print_info (" clusters]\n");
00231 }
00232
00233
00234 bool cluster_found = false;
00235 for (size_t i = 0; i < euclidean_label_indices.size (); i++)
00236 {
00237 if (cluster_found)
00238 break;
00239
00240 for (size_t j = 0; j < euclidean_label_indices[i].indices.size (); ++j)
00241 {
00242 if (picked_idx != euclidean_label_indices[i].indices[j])
00243 continue;
00244 copyPointCloud (*cloud, euclidean_label_indices[i].indices, object);
00245 cluster_found = true;
00246 break;
00247 }
00248 }
00249 }
00250
00252 void
00253 segment (const PointT &picked_point,
00254 int picked_idx,
00255 PlanarRegion<PointT> ®ion,
00256 typename PointCloud<PointT>::Ptr &object)
00257 {
00258 object.reset ();
00259
00260
00261 vector<ModelCoefficients> model_coefficients;
00262 vector<PointIndices> inlier_indices, boundary_indices;
00263 vector<PlanarRegion<PointT>, Eigen::aligned_allocator<PlanarRegion<PointT> > > regions;
00264
00265
00266 if (cloud_->isOrganized ())
00267 {
00268 print_highlight (stderr, "Estimating normals ");
00269 TicToc tt; tt.tic ();
00270
00271 PointCloud<Normal>::Ptr normal_cloud (new PointCloud<Normal>);
00272 estimateNormals (cloud_, *normal_cloud);
00273 print_info ("[done, "); print_value ("%g", tt.toc ()); print_info (" ms : "); print_value ("%zu", normal_cloud->size ()); print_info (" points]\n");
00274
00275 OrganizedMultiPlaneSegmentation<PointT, Normal, Label> mps;
00276 mps.setMinInliers (1000);
00277 mps.setAngularThreshold (deg2rad (3.0));
00278 mps.setDistanceThreshold (0.03);
00279 mps.setMaximumCurvature (0.001);
00280 mps.setProjectPoints (true);
00281 mps.setComparator (plane_comparator_);
00282 mps.setInputNormals (normal_cloud);
00283 mps.setInputCloud (cloud_);
00284
00285
00286 PointCloud<Label>::Ptr labels (new PointCloud<Label>);
00287 vector<PointIndices> label_indices;
00288 mps.segmentAndRefine (regions, model_coefficients, inlier_indices, labels, label_indices, boundary_indices);
00289 }
00290 else
00291 {
00292 SACSegmentation<PointT> seg;
00293 seg.setOptimizeCoefficients (true);
00294 seg.setModelType (SACMODEL_PLANE);
00295 seg.setMethodType (SAC_RANSAC);
00296 seg.setMaxIterations (10000);
00297 seg.setDistanceThreshold (0.005);
00298
00299
00300 typename PointCloud<PointT>::Ptr cloud_segmented (new PointCloud<PointT> (*cloud_));
00301 typename PointCloud<PointT>::Ptr cloud_remaining (new PointCloud<PointT>);
00302
00303 ModelCoefficients coefficients;
00304 ExtractIndices<PointT> extract;
00305 PointIndices::Ptr inliers (new PointIndices ());
00306
00307
00308 int i = 1;
00309 while (double (cloud_segmented->size ()) > 0.3 * double (cloud_->size ()))
00310 {
00311 seg.setInputCloud (cloud_segmented);
00312
00313 print_highlight (stderr, "Searching for the largest plane (%2.0d) ", i++);
00314 TicToc tt; tt.tic ();
00315 seg.segment (*inliers, coefficients);
00316 print_info ("[done, "); print_value ("%g", tt.toc ()); print_info (" ms : "); print_value ("%zu", inliers->indices.size ()); print_info (" points]\n");
00317
00318
00319 if (inliers->indices.empty ())
00320 break;
00321
00322
00323 PlanarRegion<PointT> region;
00324 region.setCoefficients (coefficients);
00325 regions.push_back (region);
00326
00327 inlier_indices.push_back (*inliers);
00328 model_coefficients.push_back (coefficients);
00329
00330
00331 extract.setInputCloud (cloud_segmented);
00332 extract.setIndices (inliers);
00333 extract.setNegative (true);
00334 extract.filter (*cloud_remaining);
00335 cloud_segmented.swap (cloud_remaining);
00336 }
00337 }
00338 print_highlight ("Number of planar regions detected: %zu for a cloud of %zu points\n", regions.size (), cloud_->size ());
00339
00340 double max_dist = numeric_limits<double>::max ();
00341
00342 int idx = -1;
00343 for (size_t i = 0; i < regions.size (); ++i)
00344 {
00345 double dist = pointToPlaneDistance (picked_point, regions[i].getCoefficients ());
00346 if (dist < max_dist)
00347 {
00348 max_dist = dist;
00349 idx = static_cast<int> (i);
00350 }
00351 }
00352
00353
00354 if (idx != -1)
00355 {
00356 plane_indices_.reset (new PointIndices (inlier_indices[idx]));
00357
00358 if (cloud_->isOrganized ())
00359 {
00360 approximatePolygon (regions[idx], region, 0.01f, false, true);
00361 print_highlight ("Planar region: %zu points initial, %zu points after refinement.\n", regions[idx].getContour ().size (), region.getContour ().size ());
00362 }
00363 else
00364 {
00365
00366 region = regions[idx];
00367
00368 print_highlight (stderr, "Obtaining the boundary points for the region ");
00369 TicToc tt; tt.tic ();
00370
00371 typename PointCloud<PointT>::Ptr cloud_projected (new PointCloud<PointT>);
00372 ModelCoefficients::Ptr coefficients (new ModelCoefficients (model_coefficients[idx]));
00373 ProjectInliers<PointT> proj;
00374 proj.setModelType (SACMODEL_PLANE);
00375 proj.setInputCloud (cloud_);
00376 proj.setIndices (plane_indices_);
00377 proj.setModelCoefficients (coefficients);
00378 proj.filter (*cloud_projected);
00379
00380
00381 ConvexHull<PointT> chull;
00382 chull.setDimension (2);
00383 chull.setInputCloud (cloud_projected);
00384 PointCloud<PointT> plane_hull;
00385 chull.reconstruct (plane_hull);
00386 region.setContour (plane_hull);
00387 print_info ("[done, "); print_value ("%g", tt.toc ()); print_info (" ms : "); print_value ("%zu", plane_hull.size ()); print_info (" points]\n");
00388 }
00389
00390 }
00391
00392
00393 if (plane_indices_ && !plane_indices_->indices.empty ())
00394 {
00395 plane_.reset (new PointCloud<PointT>);
00396 copyPointCloud (*cloud_, plane_indices_->indices, *plane_);
00397 object.reset (new PointCloud<PointT>);
00398 segmentObject (picked_idx, cloud_, plane_indices_, *object);
00399 }
00400 }
00401
00403
00408 void
00409 pp_callback (const visualization::PointPickingEvent& event, void*)
00410 {
00411
00412 int idx = event.getPointIndex ();
00413 if (idx == -1)
00414 return;
00415
00416 vector<int> indices (1);
00417 vector<float> distances (1);
00418
00419
00420 PointT picked_pt;
00421 event.getPoint (picked_pt.x, picked_pt.y, picked_pt.z);
00422
00423 print_info (stderr, "Picked point with index %d, and coordinates %f, %f, %f.\n", idx, picked_pt.x, picked_pt.y, picked_pt.z);
00424
00425
00426 stringstream ss;
00427 ss << "sphere_" << idx;
00428 cloud_viewer_->addSphere (picked_pt, 0.01, 1.0, 0.0, 0.0, ss.str ());
00429
00430
00431 search_->nearestKSearch (picked_pt, 1, indices, distances);
00432
00433
00434 if (image_viewer_)
00435 {
00436
00437 uint32_t width = search_->getInputCloud ()->width,
00438 height = search_->getInputCloud ()->height;
00439 int v = height - indices[0] / width,
00440 u = indices[0] % width;
00441
00442 image_viewer_->addCircle (u, v, 5, 1.0, 0.0, 0.0, "circles", 1.0);
00443 image_viewer_->addFilledRectangle (u-5, u+5, v-5, v+5, 0.0, 1.0, 0.0, "boxes", 0.5);
00444 image_viewer_->markPoint (u, v, visualization::red_color, visualization::blue_color, 10);
00445 }
00446
00447
00448
00449
00450 PlanarRegion<PointT> region;
00451 segment (picked_pt, indices[0], region, object_);
00452
00453
00454 if (region.getContour ().empty ())
00455 {
00456 PCL_ERROR ("No planar region detected. Please select another point or relax the thresholds and continue.\n");
00457 return;
00458 }
00459
00460 else
00461 {
00462 cloud_viewer_->addPolygon (region, 0.0, 0.0, 1.0, "region");
00463 cloud_viewer_->setShapeRenderingProperties (visualization::PCL_VISUALIZER_LINE_WIDTH, 10, "region");
00464
00465
00466 if (image_viewer_)
00467 {
00468 image_viewer_->addPlanarPolygon (search_->getInputCloud (), region, 0.0, 0.0, 1.0, "refined_region", 1.0);
00469 }
00470 }
00471
00472
00473 if (!object_)
00474 {
00475 PCL_ERROR ("No object detected. Please select another point or relax the thresholds and continue.\n");
00476 return;
00477 }
00478 else
00479 {
00480
00481 visualization::PointCloudColorHandlerCustom<PointT> red (object_, 255, 0, 0);
00482 if (!cloud_viewer_->updatePointCloud (object_, red, "object"))
00483 cloud_viewer_->addPointCloud (object_, red, "object");
00484
00485 if (image_viewer_)
00486 {
00487 image_viewer_->removeLayer ("object");
00488 image_viewer_->addMask (search_->getInputCloud (), *object_, "object");
00489 }
00490
00491
00492 if (image_viewer_)
00493 image_viewer_->addRectangle (search_->getInputCloud (), *object_);
00494 }
00495 }
00496
00498 void
00499 compute ()
00500 {
00501
00502 while (!cloud_viewer_->wasStopped ())
00503 {
00504
00505
00506
00507
00508
00509
00510
00511
00512 cloud_viewer_->spinOnce ();
00513 if (image_viewer_)
00514 {
00515 image_viewer_->spinOnce ();
00516 if (image_viewer_->wasStopped ())
00517 break;
00518 }
00519 boost::this_thread::sleep (boost::posix_time::microseconds (100));
00520 }
00521 }
00522
00524 void
00525 initGUI ()
00526 {
00527 cloud_viewer_.reset (new visualization::PCLVisualizer ("PointCloud"));
00528
00529 if (cloud_->isOrganized ())
00530 {
00531
00532 vector<pcl::PCLPointField> fields;
00533 int rgba_index = -1;
00534 rgba_index = getFieldIndex (*cloud_, "rgba", fields);
00535
00536 if (rgba_index >= 0)
00537 {
00538 image_viewer_.reset (new visualization::ImageViewer ("RGB PCLImage"));
00539
00540 image_viewer_->registerMouseCallback (&ObjectSelection::mouse_callback, *this);
00541 image_viewer_->registerKeyboardCallback(&ObjectSelection::keyboard_callback, *this);
00542 image_viewer_->setPosition (cloud_->width, 0);
00543 image_viewer_->setSize (cloud_->width, cloud_->height);
00544
00545 int poff = fields[rgba_index].offset;
00546
00547 rgb_data_ = new unsigned char [cloud_->width * cloud_->height * 3];
00548 for (uint32_t i = 0; i < cloud_->width * cloud_->height; ++i)
00549 {
00550 RGB rgb;
00551 memcpy (&rgb, reinterpret_cast<unsigned char*> (&cloud_->points[i]) + poff, sizeof (rgb));
00552
00553 rgb_data_[i * 3 + 0] = rgb.r;
00554 rgb_data_[i * 3 + 1] = rgb.g;
00555 rgb_data_[i * 3 + 2] = rgb.b;
00556 }
00557 image_viewer_->showRGBImage (rgb_data_, cloud_->width, cloud_->height);
00558 }
00559 cloud_viewer_->setSize (cloud_->width, cloud_->height);
00560 }
00561
00562 cloud_viewer_->registerMouseCallback (&ObjectSelection::mouse_callback, *this);
00563 cloud_viewer_->registerKeyboardCallback(&ObjectSelection::keyboard_callback, *this);
00564 cloud_viewer_->registerPointPickingCallback (&ObjectSelection::pp_callback, *this);
00565 cloud_viewer_->setPosition (0, 0);
00566
00567 cloud_viewer_->addPointCloud (cloud_, "scene");
00568 cloud_viewer_->resetCameraViewpoint ("scene");
00569 cloud_viewer_->addCoordinateSystem (0.1, 0, 0, 0);
00570 }
00571
00573 bool
00574 load (const std::string &file)
00575 {
00576
00577 TicToc tt; tt.tic ();
00578 print_highlight (stderr, "Loading ");
00579 print_value (stderr, "%s ", file.c_str ());
00580 cloud_.reset (new PointCloud<PointT>);
00581 if (io::loadPCDFile (file, *cloud_) < 0)
00582 {
00583 print_error (stderr, "[error]\n");
00584 return (false);
00585 }
00586 print_info ("[done, "); print_value ("%g", tt.toc ());
00587 print_info (" ms : "); print_value ("%zu", cloud_->size ()); print_info (" points]\n");
00588
00589 if (cloud_->isOrganized ())
00590 search_.reset (new search::OrganizedNeighbor<PointT>);
00591 else
00592 search_.reset (new search::KdTree<PointT>);
00593
00594 search_->setInputCloud (cloud_);
00595
00596 return (true);
00597 }
00598
00600 void
00601 save (const std::string &object_file, const std::string &plane_file)
00602 {
00603 PCDWriter w;
00604 if (object_ && !object_->empty ())
00605 {
00606 w.writeBinaryCompressed (object_file, *object_);
00607 w.writeBinaryCompressed (plane_file, *plane_);
00608 print_highlight ("Object succesfully segmented. Saving results in: %s, and %s.\n", object_file.c_str (), plane_file.c_str ());
00609 }
00610 }
00611
00612 boost::shared_ptr<visualization::PCLVisualizer> cloud_viewer_;
00613 boost::shared_ptr<visualization::ImageViewer> image_viewer_;
00614
00615 typename PointCloud<PointT>::Ptr cloud_;
00616 typename search::Search<PointT>::Ptr search_;
00617 private:
00618
00619 typename EdgeAwarePlaneComparator<PointT, Normal>::Ptr plane_comparator_;
00620 PointIndices::Ptr plane_indices_;
00621 unsigned char* rgb_data_;
00622 std::vector<float> distance_map_;
00623
00624
00625 typename PointCloud<PointT>::Ptr plane_;
00626 typename PointCloud<PointT>::Ptr object_;
00627 };
00628
00629
00630 int
00631 main (int argc, char** argv)
00632 {
00633
00634 std::vector<int> p_file_indices;
00635 p_file_indices = parse_file_extension_argument (argc, argv, ".pcd");
00636 if (p_file_indices.empty ())
00637 {
00638 print_error (" Need at least an input PCD file (e.g. scene.pcd) to continue!\n\n");
00639 print_info ("Ideally, need an input file, and three output PCD files, e.g., object.pcd, plane.pcd, rest.pcd\n");
00640 return (-1);
00641 }
00642
00643 string object_file = "object.pcd", plane_file = "plane.pcd", rest_file = "rest.pcd";
00644 if (p_file_indices.size () >= 4)
00645 rest_file = argv[p_file_indices[3]];
00646 if (p_file_indices.size () >= 3)
00647 plane_file = argv[p_file_indices[2]];
00648 if (p_file_indices.size () >= 2)
00649 object_file = argv[p_file_indices[1]];
00650
00651
00652 PCDReader reader;
00653
00654 pcl::PCLPointCloud2 dummy;
00655 reader.readHeader (argv[p_file_indices[0]], dummy);
00656 if (dummy.height != 1 && getFieldIndex (dummy, "rgba") != -1)
00657 {
00658 print_highlight ("Enabling 2D image viewer mode.\n");
00659 ObjectSelection<PointXYZRGBA> s;
00660 if (!s.load (argv[p_file_indices[0]])) return (-1);
00661 s.initGUI ();
00662 s.compute ();
00663 s.save (object_file, plane_file);
00664 }
00665 else
00666 {
00667 ObjectSelection<PointXYZ> s;
00668 if (!s.load (argv[p_file_indices[0]])) return (-1);
00669 s.initGUI ();
00670 s.compute ();
00671 s.save (object_file, plane_file);
00672 }
00673
00674 return (0);
00675 }
00676