Remove small clusters from a CGAL::Point_set_3

small_cluster_removal.cpp

template < typename Iterable_t >
static std::map<int, int>
compute_cluster_size(const Iterable_t& cluster_map) {
    std::map<int, int> cluster_size;

    for ( const auto& idx : cluster_map ) {
        if ( cluster_size.find(idx) == cluster_size.end() ) {
            // New idx.
            cluster_size[idx] = 1;
        } else {
            // Existing cluster.
            cluster_size[idx]++;
        }
    }

    return cluster_size;
}

// sr::PointSet_t is CGAL::Point_set_3.
void sr::compute_and_remove_small_clusters( 
    sr::PointSet_t& cloud, 
    double avg_spacing, 
    int size_limit, 
    std::ostream& os ) {
    sr::PointSet_t::Property_map<int> cluster_map;

    // Reset or create a new property map.
    if ( cloud.has_property_map<int>("cluster") ) {
        boost::tie( cluster_map, boost::tuples::ignore ) = 
            cloud.property_map<int>("cluster");

        cloud.remove_property_map(cluster_map);
    }

    bool success = false;
    boost::tie( cluster_map, success ) = 
        cloud.add_property_map<int>("cluster", -1);

    if ( !success ) {
        os << "Cannot add cluster property map. \n";
        return;
    }

    // Clustering.
    CGAL::Real_timer rt; rt.start();
    std::size_t n_clusters = 
        CGAL::cluster_point_set(
            cloud,
            cluster_map,
            cloud.parameters().neighbor_radius(avg_spacing)
        );
    rt.stop();
    os << "Find " << n_clusters << " clusters in " << rt.time() << " seoncds. \n";

    // Compute the cluster size.
    std::map<int, int> cluster_size = compute_cluster_size( cluster_map );

    // Remove small clusters.
    for ( auto it = cloud.begin(); it != cloud.end(); ++it ) {
        // Get the cluster index.
        const auto c_idx = cluster_map[*it];

        // Get the cluster size.
        const auto c_size = cluster_size[c_idx];

        if ( c_size <= size_limit ) cloud.remove(it);
    }

    cloud.collect_garbage();
}

Thank you for posting this. It is exactly what I was looking for.
I think there is a small bug at line 66. When the point is remove, the loop counter "it" needs to be decremented by one. Otherwise points are skipped. Adding this one line of code made this work very well.

Hi jaffarrell,

Thank you for pointing that out.

I was under the impression that remove() only marks an entity as to be removed and before I call collect_garbage() the index/iterator of a point won't change. But after giving the documentation a second look, it seems anytime the user calls remove() the end() iterator also changes its value if we dereference it. This supports what you are talking about the bug on line 66.

It has been too long for me to reproduce the code, if possible, can you help me to double-check this is the case where after a call to remove(), subsequent iterator dereference will be pointing to different entities?

Thank you!

I misread the documentation as well, so I thought the code would work as written when downloaded. It did not. Once I added the decrement to the "it" loop counter it worked great. I embedded your call into my program, so I cannot just push a commit. I have pasted the corrected code below. It shows my counters to check that everything added up correctly while debugging, but really only one line of your code (i.e., it--;) changed.

num_pnts = 0;
size_t num_rmvd = 0,num_kept = 0;
for ( auto it = wall_points.begin(); it != wall_points.end(); ++it ) 
{
    num_pnts++;
    if ( rmv_point.at(*it) ) 
    {
        wall_points.remove(it);
        it--;
        num_rmvd ++;
    }
    else
       num_kept ++;
}    
std::cout << num_pnts << " were processed with "<<num_rmvd<<" removed and "<<num_kept<<" kept\n\n"<<std::endl;
wall_points.collect_garbage();