diff --git a/CMakeLists.txt b/CMakeLists.txt index 3a57f04..246253d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,6 +55,7 @@ geoflow_create_plugin( src/Raster.cpp src/heightfield_nodes.cpp src/polygon_triangulate.cpp + src/line_regulariser.cpp ) target_include_directories(gfp_buildingreconstruction PRIVATE diff --git a/src/line_regulariser.cpp b/src/line_regulariser.cpp new file mode 100644 index 0000000..10e264b --- /dev/null +++ b/src/line_regulariser.cpp @@ -0,0 +1,223 @@ +#include "line_regulariser.hpp" +#include + +namespace linereg { + + double calc_mean_angle(const std::vector& lines) { + // length-weighted mean angle + double angle_sum=0, sqlenth_sum=0; + for (auto& line : lines) { + sqlenth_sum += line->sqlength; + } + for (auto& line : lines) { + angle_sum += line->angle * line->sqlength; + } + return angle_sum/sqlenth_sum; + } + + Point_2 calc_centroid(const std::vector& lines) { + double cx=0, cy=0; + for (auto& line : lines) { + auto p = line->segment.source(); + cx += CGAL::to_double(p.x()); + cy += CGAL::to_double(p.y()); + p = line->segment.target(); + cx += CGAL::to_double(p.x()); + cy += CGAL::to_double(p.y()); + } + size_t np = 2*lines.size(); + return Point_2(cx/np, cy/np); + } + + // construct a line L with centroid and mean_angle, then project all the segments from lines on L to bound it to a segment + Segment_2 calc_segment(Point_2 centroid, double mean_angle, const std::vector& lines, double extension) { + auto lv = Vector_2(std::tan(mean_angle), 1.0); + lv = lv / CGAL::sqrt(CGAL::to_double(lv.squared_length())); + + bool setminmax=false; + Point_2 pmin, pmax; + double dmin, dmax; + for (auto& line : lines) { + auto p = line->segment.source(); + auto d = CGAL::to_double(Vector_2(p.x(),p.y())*lv); + if (!setminmax) { + setminmax=true; + dmin=dmax=d; + pmin=pmax=p; + } + if (d < dmin){ + dmin = d; + pmin = p; + } + if (d > dmax) { + dmax = d; + pmax = p; + } + + p = line->segment.target(); + d = CGAL::to_double(Vector_2(p.x(),p.y())*lv); + if (d < dmin){ + dmin = d; + pmin = p; + } + if (d > dmax) { + dmax = d; + pmax = p; + } + } + Line_2 l(centroid, lv); + pmin = l.projection(pmin) - lv*extension; + pmax = l.projection(pmax) + lv*extension; + return Segment_2(pmin, pmax); + } + + double AngleCluster::distance(Cluster* other_cluster) { + return std::fabs(value - other_cluster->value); + } + void AngleCluster::calc_mean_value() { + value = calc_mean_angle(lines); + } + + double DistCluster::distance(Cluster* other_cluster) { + return CGAL::to_double(CGAL::squared_distance(value, other_cluster->value)); + } + void DistCluster::calc_mean_value() { + // a segment through the length-weighted mean centroid and elongated to 'cover' all the segments + double mean_angle = calc_mean_angle(lines); + Point_2 centroid = calc_centroid(lines); + value = calc_segment(centroid, mean_angle, lines); + } + + template DistanceTable::DistanceTable(std::set& clusters) : clusters(clusters) { + // compute only half of the distance table, since the other half will be exactly the same + for(auto cit_a = clusters.begin(); cit_a != clusters.end(); ++cit_a) { + for(auto cit_b = std::next(cit_a); cit_b != clusters.end(); ++cit_b) { + auto cluster_a = *cit_a; + auto cluster_b = *cit_b; + // do not create entries for cluster pairs that both have an intersection line, these are not to be merged + // if (cluster_a->has_intersection_line && cluster_b->has_intersection_line) continue; + // push distance to heap + auto hhandle = distances.push(ClusterPairDist( + std::make_pair(cluster_a, cluster_b), + cluster_a->distance(cluster_b.get()) + )); + // store handle for both clusters + auto hh = std::make_shared(hhandle); + cluster_to_dist_pairs[cluster_a].push_back(hh); + cluster_to_dist_pairs[cluster_b].push_back(hh); + } + } + } + template void DistanceTable::merge(ClusterH lhs, ClusterH rhs) { + // merges two clusters, then removes one from the distances map and update the affected distances + // merge + lhs->lines.insert(lhs->lines.end(), rhs->lines.begin(), rhs->lines.end()); + // lhs->has_intersection_line = lhs->has_intersection_line || rhs->has_intersection_line; + lhs->calc_mean_value(); + + // iterate distancetable + // if rhs in pair: remove pair + // if lhs has intersection line and lhs in a pair: remove pair if the other cluster also has an intersection line. Since clusters that both have intersection line are not te be merged. + auto& rhs_hhandles = cluster_to_dist_pairs[rhs]; + for (auto& hhandle : cluster_to_dist_pairs[rhs]) { + if (hhandle.use_count()==2) { + distances.erase(*hhandle); + } + } + clusters.erase(rhs); + cluster_to_dist_pairs.erase(rhs); + + auto& lhs_hhandles = cluster_to_dist_pairs[lhs]; + auto i = lhs_hhandles.begin(); + while (i != lhs_hhandles.end()) { + // we check the use count of the shared ptr to our heap handle. There should be exactly 2 for the case where both clusters still exist, otherwise one has been merged before and the heap handle was also erased before + if (i->use_count()!=2) + { + lhs_hhandles.erase(i++); + } else { + // dereference 3 times! iterator->shared_ptr->heap_handle->heap_value Notice -> is not implemented on the heap_handle, so we must use * for the last dereference here + (***i).dist = (***i).clusters.first->distance((***i).clusters.second.get()); + distances.update(*(*i)); + ++i; + } + } + } + template typename DistanceTable::ClusterPairDist DistanceTable::get_closest_pair() { + ClusterPairDist min_pair = distances.top(); //distances.pop(); <- do not pop here; the node will be removed later in the merge function! + return min_pair; + } + + template class DistanceTable; + template class DistanceTable; + + void LineRegulariser::perform_angle_clustering() { + //make clusters + angle_clusters.clear(); + for(auto& line : lines) { + auto aclusterh = std::make_shared(); + aclusterh->value = line.angle; + // aclusterh->has_intersection_line = line.priority==2; + aclusterh->lines.push_back(&line); + angle_clusters.insert(aclusterh); + } + + if (angle_clusters.size()>1) { + // make distance table + DistanceTable adt(angle_clusters); + + auto apair = adt.get_closest_pair(); + while (apair.dist < angle_threshold) { + adt.merge(apair.clusters.first, apair.clusters.second); + if (adt.distances.size()==0) break; + apair = adt.get_closest_pair(); + } + } + + size_t id_cntr=0; + for(auto& aclusterh : angle_clusters) { + for(auto line : aclusterh->lines){ + line->angle_cluster = aclusterh; + line->angle_cluster_id = id_cntr; + } + ++id_cntr; + } + } + void LineRegulariser::perform_distance_clustering() { + dist_clusters.clear(); + + // perform distance clustering for each angle cluster + for(auto& aclusterh : angle_clusters) { + std::set dclusters; + for(auto& line : aclusterh->lines) { + auto dclusterh = std::make_shared(); + dclusterh->value = line->segment; + // dclusterh->has_intersection_line = line->priority==2; + dclusterh->lines.push_back(line); + dclusters.insert(dclusterh); + } + + if (dclusters.size()>1) { + // make distance table + DistanceTable ddt(dclusters); + + // do clustering + auto dpair = ddt.get_closest_pair(); + while (dpair.dist < dist_threshold) { + ddt.merge(dpair.clusters.first, dpair.clusters.second); + if (ddt.distances.size()==0) break; + dpair = ddt.get_closest_pair(); + } + } + dist_clusters.insert(dclusters.begin(), dclusters.end()); + } + + size_t id_cntr=0; + for(auto& dclusterh : dist_clusters) { + for(auto line : dclusterh->lines) { + line->dist_cluster = dclusterh; + line->dist_cluster_id = id_cntr; + } + ++id_cntr; + } + } +} \ No newline at end of file diff --git a/src/line_regulariser.hpp b/src/line_regulariser.hpp index 434b1bf..d81c700 100644 --- a/src/line_regulariser.hpp +++ b/src/line_regulariser.hpp @@ -1,44 +1,148 @@ +#pragma once + #include #include #include #include #include -namespace linereg { - typedef CGAL::Exact_predicates_inexact_constructions_kernel::Point_2 Point_2; - typedef CGAL::Exact_predicates_inexact_constructions_kernel::Point_3 Point_3; - typedef CGAL::Exact_predicates_inexact_constructions_kernel::Vector_2 Vector_2; +#include - typedef CGAL::Exact_predicates_exact_constructions_kernel EK; - typedef CGAL::Polygon_2 Polygon_2; - typedef CGAL::Polygon_with_holes_2 Polygon_with_holes_2; - // typedef EK::Polygon_2 Polygon_2; +namespace linereg { + typedef CGAL::Exact_predicates_exact_constructions_kernel K; + // typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + + typedef K::Point_2 Point_2; + typedef K::Point_3 Point_3; + typedef K::Vector_2 Vector_2; + typedef K::Line_2 Line_2; + typedef K::Segment_2 Segment_2; + typedef CGAL::Polygon_2 Polygon_2; + typedef CGAL::Polygon_with_holes_2 Polygon_with_holes_2; + + struct linetype; + + template struct Cluster { + T value; + bool has_intersection_line; + std::vector lines; + virtual double distance(Cluster* other_cluster)=0; + virtual void calc_mean_value()=0; + }; + struct AngleCluster : public Cluster{ + double distance(Cluster* other_cluster); + void calc_mean_value(); + }; + struct DistCluster : public Cluster{ + double distance(Cluster* other_cluster); + void calc_mean_value(); + }; - class LineRegulariser { - static constexpr double pi = 3.14159265358979323846; + template struct DistanceTable { + typedef std::pair ClusterPair; - struct ValueCluster { - Vector_2 ref_vec; - Vector_2 ref_point; - std::vector idx; + // fibonacci heap from boost + struct ClusterPairDist { + ClusterPairDist(ClusterPair p, double d) : clusters(p), dist(d){} + ClusterPair clusters; + double dist; + + bool operator<(ClusterPairDist const & rhs) const + { + return dist > rhs.dist; + } }; + typedef boost::heap::fibonacci_heap DistanceHeap; + typedef typename DistanceHeap::handle_type heap_handle; + + // define hash function such that the same hash results regardless of the order of cluster handles in the pair + // struct KeyHash { + // size_t operator()(const ClusterPair& key) const { + // if (key.first.get() < key.second.get()) + // return std::hash()(key.first) ^ + // (std::hash()(key.second) << 1); + // else + // return std::hash()(key.second) ^ + // (std::hash()(key.first) << 1); + + // } + // }; + // // True equality function is needed to deal with collisions + // struct KeyEqual { + // bool operator()(const ClusterPair& lhs, const ClusterPair& rhs) const { + // return + // ((lhs.first==rhs.first) && (lhs.second==rhs.second)) + // || + // ((lhs.first==rhs.second) && (lhs.second==rhs.first)); + // } + // }; + typedef std::unordered_map>> Cluster2DistPairMap; - typedef std::tuple linetype; - // new angle, midpoint, distance in angle cluster, priority, segment_id, sqlength - typedef std::vector vec_ek_seg; + Cluster2DistPairMap cluster_to_dist_pairs; + DistanceHeap distances; + std::set& clusters; + + DistanceTable(std::set& clusters); //computes initial distances + void merge(ClusterH lhs, ClusterH rhs); // merges two clusters, then removes one from the distances map and update the affected distances + ClusterPairDist get_closest_pair(); //returns the cluster pair with the smallest distance + }; + + // extern template class Cluster; + // extern template class Cluster; + + typedef std::shared_ptr AngleClusterH; + typedef std::shared_ptr DistClusterH; + + extern template class DistanceTable; + extern template class DistanceTable; + + double calc_mean_angle(const std::vector& lines); + Point_2 calc_centroid(const std::vector& lines); + Segment_2 calc_segment(Point_2 centroid, double mean_angle, const std::vector& lines, double extension=0); + + struct linetype { + linetype(Segment_2 segment_, double angle_, Point_2 midpoint_, double dist_in_ang_cluster_, size_t priority_, size_t segment_id_, double sqlength_) : + segment(segment_), angle(angle_), midpoint(midpoint_), priority(priority_), segment_id(segment_id_), sqlength(sqlength_) {}; + + Segment_2 segment; + double angle; + Point_2 midpoint; + Vector_2 direction; + double sqlength; + // double dist_in_ang_cluster; + size_t priority; + size_t segment_id; + + Segment_2 reg_segment; + AngleClusterH angle_cluster; + DistClusterH dist_cluster; + size_t angle_cluster_id, dist_cluster_id; + }; + + static constexpr double pi = 3.14159265358979323846; + class LineRegulariser { + + typedef std::vector SegmentVec; // geoflow::SegmentCollection& input_segments; public: std::vector lines; - // vec_ek_seg input_reg_exact; + // SegmentVec input_reg_exact; double angle_threshold, dist_threshold; - std::unordered_map segments; + std::unordered_map segments; + std::set angle_clusters; + std::set dist_clusters; LineRegulariser() {}; void add_segments(size_t priority, const Polygon_2& polygon, double offset) { - size_t i=0; + size_t i; + if(segments.find(priority) == segments.end()) { + i=0; + } else { + i=segments.size(); + } auto orientation = polygon.orientation(); for(auto edge = polygon.edges_begin(); edge != polygon.edges_end(); ++edge) { auto source = edge->source(); @@ -55,306 +159,113 @@ namespace linereg { auto l = CGAL::to_double(v.squared_length()); auto angle = std::atan2(CGAL::to_double(v.x()), CGAL::to_double(v.y())); if (angle < 0) angle += pi; - lines.push_back(std::make_tuple(angle,p,0,priority,i++,l)); - segments[priority].push_back(EK::Segment_2(target,source)); + lines.push_back(linetype(*edge, angle,p,0,priority,i++,l)); + segments[priority].push_back(Segment_2(target,source)); } } void add_segments(size_t priority, geoflow::SegmentCollection& segs) { if (segs.size()==0) return; - size_t i=0; + size_t i; + if(segments.find(priority) == segments.end()) { + i=0; + } else { + i=segments.size(); + } + for(auto& edge : segs) { - auto source = EK::Point_2(edge[0][0], edge[0][1]); - auto target = EK::Point_2(edge[1][0], edge[1][1]); + auto source = Point_2(edge[0][0], edge[0][1]); + auto target = Point_2(edge[1][0], edge[1][1]); auto v = target-source; auto p_ = source + v/2; auto p = Point_2(CGAL::to_double(p_.x()),CGAL::to_double(p_.y())); auto l = CGAL::to_double(v.squared_length()); auto angle = std::atan2(CGAL::to_double(v.x()), CGAL::to_double(v.y())); if (angle < 0) angle += pi; - lines.push_back(std::make_tuple(angle,p,0,priority,i++,l)); - segments[priority].push_back(EK::Segment_2(source, target)); + lines.push_back(linetype(Segment_2(source, target), angle,p,0,priority,i++,l)); + segments[priority].push_back(Segment_2(source, target)); } } - vec_ek_seg& get_segments(const size_t& priority) { + SegmentVec& get_segments(const size_t& priority) { return segments[priority]; } - // template void get_max(const std::vector& idx, const size_t element) { - - // for(auto& i : idx) { - // std::get(lines[i]); - - // } - // } - - void cluster(bool weight_by_len, bool angle_per_distcluster) { - // cluster by angle - std::vector edge_idx(lines.size()); - for (size_t i=0; i(lines[a]) < std::get<0>(lines[b]); - }); - - std::vector angle_clusters(1); - auto angle_sum = std::get<0>(lines[edge_idx[0]]); - angle_clusters.back().idx.push_back(edge_idx[0]); - for(size_t i=1; i < edge_idx.size(); ++i) { - auto& edge_id = edge_idx[i]; - auto& line = lines[edge_id]; - auto& angle = std::get<0>(line); - // mean angle - // auto repr_angle = angle_sum/angle_clusters.back().idx.size(); - // median angle - auto mid = angle_clusters.back().idx.size()/2; - auto repr_angle = std::get<0>(lines[angle_clusters.back().idx[mid]]); - if(std::fmod((angle - repr_angle), pi) < angle_threshold) { - angle_clusters.back().idx.push_back(edge_id); - angle_sum += angle; - } else { - angle_clusters.resize(angle_clusters.size()+1); - angle_clusters.back().idx.push_back(edge_id); - angle_sum = angle; - } - } - - // determine angle for each cluster - for(auto& cluster : angle_clusters) { - // find maximum priority in this cluster - size_t max_pr=0; - for(auto& i : cluster.idx) - max_pr = std::max(max_pr, std::get<3>(lines[i])); - - // collect all ids that have the max priority - std::vector max_idx; - for(auto& i : cluster.idx) - if (std::get<3>(lines[i]) == max_pr) - max_idx.push_back(i); - - // computed average angle weighted by segment length - double sum_len=0; - double sum_all=0; - double max_len=-1; - size_t max_len_id; - for(auto& i : max_idx) { - auto& len = std::get<5>(lines[i]); - auto& angle = std::get<0>(lines[i]); - sum_all += len * angle; - sum_len += len; - if (len > max_len) { - max_len = len; - max_len_id = i; - } - } - double angle; - if (weight_by_len) - angle = sum_all/sum_len; - else - angle = std::get<0>(lines[max_len_id]); - - Vector_2 n(-1.0, std::tan(angle)); - cluster.ref_vec = n/std::sqrt(n.squared_length()); // normalize - - // or median angle: - // size_t median_id = cluster.idx[cluster.idx.size()/2]; - // cluster.value = std::get<0>(lines[median_id]); - } - - // cluster parallel lines by distance - std::vector dist_clusters; - for(auto& cluster : angle_clusters) { - auto n = cluster.ref_vec; - // compute distances along n wrt to first line in cluster - auto p = std::get<1>(lines[cluster.idx[0]]); - for(auto& i : cluster.idx) { - auto q = std::get<1>(lines[i]); - auto v = p-q; - std::get<2>(lines[i]) = v*n; - // distances.push_back(v*n); - } - // sort by distance, ascending - auto dist_idx = cluster.idx; - std::sort(dist_idx.begin(), dist_idx.end(), [&lines=lines](size_t a, size_t b){ - return std::get<2>(lines[a]) < std::get<2>(lines[b]); - }); - // cluster nearby lines using separation threshold - double dist_sum = std::get<2>(lines[dist_idx[0]]); - dist_clusters.resize(dist_clusters.size()+1); - dist_clusters.back().ref_vec = n; - dist_clusters.back().idx.push_back(dist_idx[0]); - for(size_t i=1; i < dist_idx.size(); ++i) { - auto& edge_id = dist_idx[i]; - auto& line = lines[edge_id]; - // double repr_dist = dist_sum / dist_clusters.back().idx.size(); - auto mid = dist_clusters.back().idx.size()/2; - double repr_dist = std::get<2>(lines[dist_clusters.back().idx[mid]]); - double& dist = std::get<2>(line); - if (std::abs(repr_dist-dist) < dist_threshold) { - dist_clusters.back().idx.push_back(edge_id); - dist_sum += dist; - } else { - dist_clusters.resize(dist_clusters.size()+1); - dist_clusters.back().ref_vec = n; - dist_clusters.back().idx.push_back(edge_id); - dist_sum = dist; - } - } - } - - // find direction and center point for each cluster - for(auto& cluster : dist_clusters) { - // find max priority - size_t max_pr=0; - for(auto& i : cluster.idx) - max_pr = std::max(max_pr, std::get<3>(lines[i])); - - // collect all ids that have the max priority - std::vector max_idx; - for(auto& i : cluster.idx) - if (std::get<3>(lines[i]) == max_pr) - max_idx.push_back(i); - - // computed cluster ref point - { - double sum_len=0; - double max_len=-1; - size_t max_len_id; - Vector_2 sum_all(0,0); - for(auto& i : max_idx) { - auto& len = std::get<5>(lines[i]); - auto& q = std::get<1>(lines[i]); - sum_all += len*Vector_2(q.x(), q.y()); - sum_len += len; - if (len > max_len) { - max_len = len; - max_len_id = i; - } - } - // cluster.distance = sum/cluster.idx.size(); - if (weight_by_len) - cluster.ref_point = sum_all/sum_len; - else { - auto& p = std::get<1>(lines[max_len_id]); - cluster.ref_point = Vector_2(p.x(), p.y()); - } - } - // computed cluster angle - if (angle_per_distcluster) { - double sum_len=0; - double sum_all=0; - double max_len=-1; - size_t max_len_id; - for(auto& i : max_idx) { - auto& len = std::get<5>(lines[i]); - auto& angle = std::get<0>(lines[i]); - sum_all += len * angle; - sum_len += len; - if (len > max_len) { - max_len = len; - max_len_id = i; - } - } - double angle; - if (weight_by_len) - angle = sum_all/sum_len; - else - angle = std::get<0>(lines[max_len_id]); - - Vector_2 n(-1.0, std::tan(angle)); - cluster.ref_vec = n/std::sqrt(n.squared_length()); - } - cluster.ref_vec = Vector_2(cluster.ref_vec.y(), -cluster.ref_vec.x()); - } - - // project input segments on the cluster lines - for(auto& cluster : dist_clusters) { - auto ref_v = EK::Vector_2(cluster.ref_vec.x(), cluster.ref_vec.y()); - auto ref_p = EK::Point_2(cluster.ref_point.x(), cluster.ref_point.y()); - - EK::Line_2 ref_line(ref_p, ref_v); - - for(auto& i : cluster.idx) { - auto pri = std::get<3>(lines[i]); - if (pri==2) continue; //HACK! - auto j = std::get<4>(lines[i]); - auto& edge = segments[pri][j]; - auto s_new = ref_line.projection(edge.source()); - auto t_new = ref_line.projection(edge.target()); - segments[pri][j] = EK::Segment_2(s_new, t_new); - // input_segments[i][0] = - // {float(CGAL::to_double(s_new.x())), float(CGAL::to_double(s_new.y())), 0}; - // input_segments[i][1] = - // {float(CGAL::to_double(t_new.x())), float(CGAL::to_double(t_new.y())), 0}; - } - } - }; - + void perform_angle_clustering(); + void perform_distance_clustering(); }; - template void chain( - const typename Kernel::Segment_2& a, - const typename Kernel::Segment_2& b, - typename CGAL::Polygon_2& ring, + const typename Kernel::Plane_3& plane, + const typename Kernel::Segment_3& a, + const typename Kernel::Segment_3& b, + std::vector& ring_pts, const float& snap_threshold) { - auto l_a = a.supporting_line(); - auto l_b = b.supporting_line(); - typename Kernel::Segment_2 s(a.target(), b.source()); + typedef typename Kernel::Segment_2 Segment_2; + typedef typename Kernel::Point_2 Point_2; + typedef typename Kernel::Point_3 Point_3; + typedef typename Kernel::Plane_3 Plane_3; + + auto a_2d = Segment_2(Point_2(a.source().x(), a.source().y()), Point_2(a.target().x(), a.target().y())); + auto b_2d = Segment_2(Point_2(b.source().x(), b.source().y()), Point_2(b.target().x(), b.target().y())); + + auto l_a = a_2d.supporting_line(); + auto l_b = b_2d.supporting_line(); + Segment_2 s(a_2d.target(), b_2d.source()); auto result = CGAL::intersection(l_a, l_b); if (result) { - if (auto p = boost::get(&*result)) { + if (auto p = boost::get(&*result)) { if (CGAL::squared_distance(*p, s) < snap_threshold) { - ring.push_back(*p); + double z = -plane.a()/plane.c() * p->x() - plane.b()/plane.c()*p->y() - plane.d()/plane.c(); + ring_pts.push_back( Point_3(p->x(), p->y(), z) ); + // ring_pts.push_back( plane.to_3d(*p) ); } else { - ring.push_back(a.target()); - ring.push_back(b.source()); + ring_pts.push_back(a.target()); + ring_pts.push_back(b.source()); } } // } else if (auto l = boost::get(&*result)) { - - // } } else { // there is no intersection - ring.push_back(a.target()); - ring.push_back(b.source()); + ring_pts.push_back(a.target()); + ring_pts.push_back(b.source()); } } // void chain(Segment& a, Segment& b, LinearRing& ring, const float& snap_threshold) { - template inline void check_dist(const CGAL::Polygon_2& pos, CGAL::Polygon_2& pot, const size_t a, const size_t b) { - auto d = CGAL::squared_distance(pos.vertex(a), pos.vertex(b)); - if (d > 1E-6) pot.push_back(pos.vertex(a)); + template inline void check_dist(std::vector& iring, std::vector& aring, const size_t a, const size_t b) { + auto d = CGAL::squared_distance(iring[a], iring[b]); + if (d > 1E-6) aring.push_back(iring[a]); } - template CGAL::Polygon_2 + // template CGAL::Polygon_2 + template std::vector chain_ring( - const std::vector& idx, - const std::vector& segments, + const std::vector& idx, + const typename Kernel::Plane_3& plane, + const std::vector& segments, const float& snap_threshold) { - typename CGAL::Polygon_2 ring, fixed_ring; + std::vector ring_points, new_ring_points; if (idx.size()>1) { // we need at least 2 segments for (size_t i=1; i(segments[idx[i-1]], segments[idx[i]], ring, snap_threshold); + chain(plane, segments[idx[i-1]], segments[idx[i]], ring_points, snap_threshold); } - chain(segments[idx[idx.size()-1]], segments[idx[0]], ring, snap_threshold); + chain(plane, segments[idx[idx.size()-1]], segments[idx[0]], ring_points, snap_threshold); // get rid of segments with zero length (ie duplicate points) // check again the size, to ignore degenerate case of input ring that consists of 3 co-linear segments (would get chained to eg 0 vertices) - if (ring.size()>2) { - for (size_t i=1; i(ring, fixed_ring, i-1, i); + if (ring_points.size()>2) { + for (size_t i=1; i(ring_points, new_ring_points, i-1, i); } - check_dist(ring, fixed_ring, ring.size()-1, 0); + check_dist(ring_points, new_ring_points, ring_points.size()-1, 0); } // NOTE: at this point there can still be vertices between co-linear segments (ie 3 consecutive points on the same line) } - return fixed_ring; + return new_ring_points; } } \ No newline at end of file diff --git a/src/region_growing.cpp b/src/region_growing.cpp index 186654e..83fa9c7 100644 --- a/src/region_growing.cpp +++ b/src/region_growing.cpp @@ -57,7 +57,7 @@ geoflow::Segment LineDetector::project(const size_t i1, const size_t i2) { geoflow::arr3f{float(p2n.x()), float(p2n.y()), float(p2n.z())} ); } -SCK::Segment_2 LineDetector::project_cgal(const size_t i1, const size_t i2, float extension) { +SCK::Segment_3 LineDetector::project_cgal(const size_t i1, const size_t i2, float extension) { const auto& l = segment_shapes[point_segment_idx[i1]]; const auto& p1 = indexed_points[i1].first; const auto& p2 = indexed_points[i2].first; @@ -68,9 +68,9 @@ SCK::Segment_2 LineDetector::project_cgal(const size_t i1, const size_t i2, floa v = v/CGAL::sqrt(v.squared_length()); p1n = p1n - v*extension; p2n = p2n + v*extension; - return SCK::Segment_2( - SCK::Point_2(CGAL::to_double(p1n.x()), CGAL::to_double(p1n.y())), - SCK::Point_2(CGAL::to_double(p2n.x()), CGAL::to_double(p2n.y())) + return SCK::Segment_3( + SCK::Point_3(CGAL::to_double(p1n.x()), CGAL::to_double(p1n.y()), CGAL::to_double(p1n.z())), + SCK::Point_3(CGAL::to_double(p2n.x()), CGAL::to_double(p2n.y()), CGAL::to_double(p2n.z())) ); } diff --git a/src/region_growing.h b/src/region_growing.h index 1013748..1ee0201 100644 --- a/src/region_growing.h +++ b/src/region_growing.h @@ -54,7 +54,7 @@ namespace linedect { LineDetector(vector &points, vector> neighbours); vector get_point_indices(size_t shape_id); geoflow::Segment project(const size_t i1, const size_t i2); - SCK::Segment_2 project_cgal(const size_t i1, const size_t i2, float extension); + SCK::Segment_3 project_cgal(const size_t i1, const size_t i2, float extension); size_t get_bounded_edges(geoflow::SegmentCollection& edges); std::vector detect(); diff --git a/src/stepedge_nodes.cpp b/src/stepedge_nodes.cpp index 3644af3..b9e7256 100644 --- a/src/stepedge_nodes.cpp +++ b/src/stepedge_nodes.cpp @@ -428,6 +428,10 @@ void Arr2LinearRingsNode::process(){ // input_attr_map["building_id"] = &attr_objectid_term; int j=0; + auto& plane_a = vector_output("plane_a"); + auto& plane_b = vector_output("plane_b"); + auto& plane_c = vector_output("plane_c"); + auto& plane_d = vector_output("plane_d"); for (auto face: arr.face_handles()){ if( !(only_in_footprint && !face->data().in_footprint) @@ -437,6 +441,12 @@ void Arr2LinearRingsNode::process(){ LinearRing polygon; arrangementface_to_polygon(face, polygon); linear_rings.push_back(polygon); + + // plane paramters + plane_a.push_back(float(CGAL::to_double(face->data().plane.a()))); + plane_b.push_back(float(CGAL::to_double(face->data().plane.b()))); + plane_c.push_back(float(CGAL::to_double(face->data().plane.c()))); + plane_d.push_back(float(CGAL::to_double(face->data().plane.d()))); // attributes specific to each roofpart input_attr_map["hoogte_abs"]->push_back((float)face->data().elevation_avg); @@ -1366,7 +1376,7 @@ void BuildArrFromRingsExactNode::arr_snapclean_from_fp(Arrangement_2& arr) { output("snap_fp_to_v").set(snap_to_v); output("snap_fp_to_e").set(snap_to_e); } -void BuildArrFromRingsExactNode::arr_snapclean(Arrangement_2& arr) { +void arr_snapclean(Arrangement_2& arr, float& snap_dist, bool& snap_detect_only) { typedef Arrangement_2::Traits_2 AT; typedef std::variant Candidate; typedef std::unordered_map> CandidateMap; @@ -1487,9 +1497,9 @@ void BuildArrFromRingsExactNode::arr_snapclean(Arrangement_2& arr) { } } - output("snap_v").set(snap_v); - output("snap_to_v").set(snap_to_v); - output("snap_to_e").set(snap_to_e); + // output("snap_v").set(snap_v); + // output("snap_to_v").set(snap_to_v); + // output("snap_to_e").set(snap_to_e); } @@ -1668,7 +1678,7 @@ void BuildArrFromRingsExactNode::process() { else { auto& lr = fp_in.get(); for (auto& p : lr) { - footprint.push_back(linereg::EK::Point_2(p[0], p[1])); + footprint.push_back(linereg::Point_2(p[0], p[1])); } } @@ -1718,7 +1728,7 @@ void BuildArrFromRingsExactNode::process() { } } } - if(snap_clean) arr_snapclean(arr_base); + if(snap_clean) arr_snapclean(arr_base, snap_dist, snap_detect_only); if(snap_clean_fp) arr_snapclean_from_fp(arr_base); if(extrude_unsegmented && points_per_plane.count(-1)) { @@ -1778,7 +1788,7 @@ void BuildArrFromLinesNode::process() { auto& lr = fp_term.get(); linereg::Polygon_2 poly2; for (auto& p : lr) { - poly2.push_back(linereg::EK::Point_2(p[0], p[1])); + poly2.push_back(linereg::Point_2(p[0], p[1])); } footprint = linereg::Polygon_with_holes_2(poly2); } @@ -1804,34 +1814,58 @@ void BuildArrFromLinesNode::process() { typedef std::pair PointPair; - std::vector segments; - for(size_t i=0; i(i); - Point_2 a(s[0][0],s[0][1]); - Point_2 b(s[1][0],s[1][1]); - segments.push_back(std::make_pair(a,b)); - } int arr_complexity = lines_term.size(); output("arr_complexity").set(arr_complexity); - if (arr_complexity > max_arr_complexity) { - std::sort(segments.begin(), segments.end(), [](PointPair& a, PointPair& b) { - return CGAL::squared_distance(a.first,a.second) > CGAL::squared_distance(b.first,b.second); - }); - } + if (lines_term.is_connected_type(typeid(linereg::Segment_2))) { + for(size_t i=0; i(i); + insert(arr_base, s); + } + + } else { + + std::vector segments; + for(size_t i=0; i(i); + Point_2 a(s[0][0],s[0][1]); + Point_2 b(s[1][0],s[1][1]); + segments.push_back(std::make_pair(a,b)); + } + + if (arr_complexity > max_arr_complexity) { + std::sort(segments.begin(), segments.end(), [](PointPair& a, PointPair& b) { + return CGAL::squared_distance(a.first,a.second) > CGAL::squared_distance(b.first,b.second); + }); + } + { + std::vector lines; + size_t i=0; + for(auto& s : segments) { + // Check for 0 segment length (quick fix for linux crash) + if (i++ == max_arr_complexity) break; + if (CGAL::squared_distance(s.first,s.second) > 0.001) lines.push_back(Segment_2(s.first,s.second)); + } + insert(arr_base, lines.begin(), lines.end()); + } + } + + //remove dangling edges { - std::vector lines; - size_t i=0; - for(auto& s : segments) { - // Check for 0 segment length (quick fix for linux crash) - if (i++ == max_arr_complexity) break; - if (CGAL::squared_distance(s.first,s.second) > 0.001) - lines.push_back(Line_2(s.first,s.second)); + std::vector to_remove; + for (auto he : arr_base.edge_handles()) { + if (he->face()==he->twin()->face()) + to_remove.push_back(he); + } + for (auto he : to_remove) { + arr_base.remove_edge(he); } - insert(arr_base, lines.begin(), lines.end()); } + + // if (snap_clean) arr_snapclean(arr_base, snap_dist, snap_detect_only); + output("arrangement").set(arr_base); } @@ -1899,7 +1933,7 @@ inline void DetectLinesNode::detect_lines(linedect::LineDetector& LD) { LD.detect(); } } -inline size_t DetectLinesNode::detect_lines_ring_m2(linedect::LineDetector& LD, SegmentCollection& segments_out) { +inline size_t DetectLinesNode::detect_lines_ring_m2(linedect::LineDetector& LD, Plane& plane, SegmentCollection& segments_out) { LD.dist_thres = dist_thres * dist_thres; LD.N = k; auto& c_upper = min_cnt_range.second; @@ -1980,7 +2014,7 @@ inline size_t DetectLinesNode::detect_lines_ring_m2(linedect::LineDetector& LD, } } if (perform_chaining) { - std::vector prechain_segments; + std::vector prechain_segments; std::vector idx; size_t idcnt=0; for (auto& [i0,i1] : sorted_segments) { // segments_out.push_back( LD.project(i0, i1) ); @@ -1988,22 +2022,40 @@ inline size_t DetectLinesNode::detect_lines_ring_m2(linedect::LineDetector& LD, idx.push_back(idcnt++); } // TODO: chain the ring? for better regularisation results - auto chained_segments = linereg::chain_ring(idx, prechain_segments, snap_threshold); - - // for (auto e=prechain_segments.begin(); e!=prechain_segments.end(); ++e) { - for (auto e=chained_segments.edges_begin(); e!=chained_segments.edges_end(); ++e) { + SegmentCollection new_ring; + auto chained_ring_pts = linereg::chain_ring(idx, SCK::Plane_3(plane.a(), plane.b(), plane.c(), plane.d()), prechain_segments, snap_threshold); + + if (chained_ring_pts.size() > 2) { + auto first = chained_ring_pts.begin(); + for (auto pit=std::next(first); pit!=chained_ring_pts.end(); ++pit) { + auto p2 = *pit; + auto p1 = *std::prev(pit); + segments_out.push_back({ + arr3f{ + float(CGAL::to_double(p1.x())), + float(CGAL::to_double(p1.y())), + float(CGAL::to_double(p1.z()))}, + arr3f{ + float(CGAL::to_double(p2.x())), + float(CGAL::to_double(p2.y())), + float(CGAL::to_double(p2.z()))}, + }); + } + auto p1 = *chained_ring_pts.rbegin(); + auto p2 = *first; segments_out.push_back({ arr3f{ - float(CGAL::to_double(e->source().x())), - float(CGAL::to_double(e->source().y())), - 0}, + float(CGAL::to_double(p1.x())), + float(CGAL::to_double(p1.y())), + float(CGAL::to_double(p1.z()))}, arr3f{ - float(CGAL::to_double(e->target().x())), - float(CGAL::to_double(e->target().y())), - 0}, + float(CGAL::to_double(p2.x())), + float(CGAL::to_double(p2.y())), + float(CGAL::to_double(p2.z()))}, }); } - return chained_segments.size(); + + return segments_out.size(); } else { for (const auto& e : sorted_segments) { segments_out.push_back( @@ -2018,6 +2070,7 @@ inline size_t DetectLinesNode::detect_lines_ring_m2(linedect::LineDetector& LD, void DetectLinesNode::process(){ auto input_geom = input("edge_points"); + auto planes = input("pts_per_roofplane").get(); SegmentCollection edge_segments, lines3d; @@ -2071,7 +2124,7 @@ void DetectLinesNode::process(){ linedect::LineDetector LD(cgal_pts); // SegmentCollection ring_edges; - auto n_detected = detect_lines_ring_m2(LD, edge_segments); + auto n_detected = detect_lines_ring_m2(LD, planes[plane_id].first, edge_segments); LD.get_bounded_edges(lines3d); for (size_t j=0; j(); - auto footprint = input("footprint").get(); - - typedef CGAL::Exact_predicates_inexact_constructions_kernel::Point_2 Point_2; - typedef CGAL::Exact_predicates_inexact_constructions_kernel::Point_3 Point_3; - typedef std::array arr3f; - std::vector> input_edges; - - // build vector of all input edges - for(auto edge : edges) { - input_edges.push_back(std::make_pair(Segment(edge[0], edge[1]), false)); - } - for(size_t i=0; i linetype; - // new angle, midpoint, distance in angle cluster, elevation, is_footprint, initial angle, squared distance from midpoint to an end point, id_cntr - std::vector lines; - // add non-footprint lines - size_t id_cntr = 0; - for(auto& ie : input_edges) { - auto& is_footprint = ie.second; - auto& edge = ie.first; - auto source = Point_3(edge[0][0], edge[0][1], edge[0][2]); - auto target = Point_3(edge[1][0], edge[1][1], edge[1][2]); - auto v = target-source; - auto p_ = source + v/2; - auto p = Point_2(p_.x(),p_.y()); - auto l = std::sqrt(v.squared_length()/2); - auto angle = std::atan2(v.x(),v.y()); - if (angle < 0) angle += pi; - lines.push_back(std::make_tuple(angle,p,0,p_.z(), is_footprint, angle, l, id_cntr++)); - } - - //sort by angle, smallest on top - std::vector edge_idx(input_edges.size()); - for (size_t i=0; i(lines[a]) < std::get<0>(lines[b]); - }); - //cluster by angle difference - std::vector angle_clusters(1); - auto last_angle = std::get<0>(lines[edge_idx[0]]); - for(auto edge_id : edge_idx ) { - auto& line = lines[edge_id]; - if((std::get<0>(line) - last_angle) < angle_threshold) - angle_clusters.back().idx.push_back(edge_id); - else { - angle_clusters.resize(angle_clusters.size()+1); - angle_clusters.back().idx.push_back(edge_id); - } - last_angle=std::get<0>(line); - } + auto ints_segments = input("ints_segments").get(); + // auto footprint = input("footprint").get(); + // auto ring_id = input("ring_id").get(); + // auto ring_order = input("ring_order").get(); - // get average angle for each cluster - // vec3f directions_before, directions_after; - // vec1i angles; - for(auto& cluster : angle_clusters) { - // average angle: - double sum=0; - for(auto& i : cluster.idx) { - sum+=std::get<0>(lines[i]); - } - double angle = sum/cluster.idx.size(); - Vector_2 n(-1.0, std::tan(angle)); - cluster.ref_vec = n/std::sqrt(n.squared_length()); // normalize - - // or median angle: - // size_t median_id = cluster.idx[cluster.idx.size()/2]; - // cluster.value = std::get<0>(lines[median_id]); - } - - // vec1f distances; - // snap nearby lines that are close - std::vector dist_clusters; - for(auto& cluster : angle_clusters) { - auto n = cluster.ref_vec; - // compute distances along n wrt to first line in cluster - auto p = std::get<1>(lines[cluster.idx[0]]); - for(auto& i : cluster.idx) { - auto q = std::get<1>(lines[i]); - auto v = p-q; - std::get<2>(lines[i]) = v*n; - // distances.push_back(v*n); - } - // sort by distance, ascending - auto sorted_by_dist = cluster.idx; - std::sort(sorted_by_dist.begin(), sorted_by_dist.end(), [&lines=lines](size_t a, size_t b){ - return std::get<2>(lines[a]) < std::get<2>(lines[b]); - }); - // cluster nearby lines using separation threshold - double last_dist = std::get<2>(lines[sorted_by_dist[0]]); - dist_clusters.resize(dist_clusters.size()+1); - dist_clusters.back().ref_vec = n; - for(auto& i : sorted_by_dist) { - auto& line = lines[i]; - double dist_diff = std::get<2>(line) - last_dist; - if (dist_diff < dist_threshold) { - dist_clusters.back().idx.push_back(i); - } else { - dist_clusters.resize(dist_clusters.size()+1); - dist_clusters.back().ref_vec = n; - dist_clusters.back().idx.push_back(i); - } - last_dist = std::get<2>(line); - } - } + // linereg::Polygon_2 cgal_footprint; + // std::vector ek_holes; + // for (auto& p : footprint) { + // cgal_footprint.push_back(linereg::Point_2(p[0], p[1])); + // } + // for (auto& ring : footprint.interior_rings()) { + // linereg::Polygon_2 ek_hole; + // for (auto& p : ring) { + // ek_hole.push_back(linereg::Point_2(p[0], p[1])); + // } + // ek_holes.push_back(ek_hole); + // } - // compute one line per dist cluster => the one with the highest elevation - LineStringCollection edges_out; - // std::vector> edges_out; - - for(auto& cluster : dist_clusters) { - // find average midpoint - // double sum=0; - Vector_2 sum_p(0,0); - for(auto& i : cluster.idx) { - // sum+=std::get<2>(lines[i]); - auto& q = std::get<1>(lines[i]); - sum_p += Vector_2(q.x(), q.y()); - } - // cluster.distance = sum/cluster.idx.size(); - cluster.ref_point = sum_p/cluster.idx.size(); - cluster.ref_vec = Vector_2(cluster.ref_vec.y(), -cluster.ref_vec.x()); - - - //try to find a footprint line - linetype best_line; - double best_angle; - bool found_fp=false, found_non_fp=false; - for(auto& i : cluster.idx) { - auto& line = lines[i]; - if(std::get<4>(line)){ - found_fp=true; - best_line = line; - best_angle = std::get<5>(line); // initial angle of this fp line - } else { - found_non_fp = true; + // get clusters from line regularisation + auto LR = linereg::LineRegulariser(); + LR.add_segments(0,edges); + // LR.add_segments(1,cgal_footprint, (double) 0); + // for (auto& hole : ek_holes) { + // LR.add_segments(1,hole, (double) 0); + // } + LR.add_segments(2,ints_segments); + LR.dist_threshold = dist_threshold*dist_threshold; + LR.angle_threshold = angle_threshold; + + std::cout << "angle clustering...\n"; + LR.perform_angle_clustering(); + std::cout << "distance clustering...\n"; + LR.perform_distance_clustering(); + std::cout << "...clustering complete\n"; + + // output("exact_footprint_out").set(linereg::Polygon_with_holes_2(cgal_footprint, ek_holes.begin(), ek_holes.end())); + + auto& regularised = vector_output("regularised_edges"); + auto& regularised_exact = vector_output("exact_regularised_edges"); + SegmentCollection edges_out_; + vec1i priorities, angle_cluster_ids, dist_cluster_ids; + // we should iterate of the distance clusters and output one segment per cluster + for(auto& line : LR.lines) { + linereg::Segment_2 segment; + segment = line.segment; + auto new_seg = Segment(); + new_seg[0] = {float(CGAL::to_double(segment.source().x())), float(CGAL::to_double(segment.source().y())), 0}; + new_seg[1] = {float(CGAL::to_double(segment.target().x())), float(CGAL::to_double(segment.target().y())), 0}; + edges_out_.push_back(new_seg); + priorities.push_back(line.priority); + priorities.push_back(line.priority); + angle_cluster_ids.push_back(line.angle_cluster_id); + angle_cluster_ids.push_back(line.angle_cluster_id); + dist_cluster_ids.push_back(line.dist_cluster_id); + dist_cluster_ids.push_back(line.dist_cluster_id); + } + for(auto& dclust : LR.dist_clusters) { + size_t dclust_size = dclust->lines.size(); + Segment_2 segment; + //get lines with highest priority + size_t max_priority=0; + for (auto line : dclust->lines) { + if(line->priority > max_priority) max_priority = line->priority; + } + std::vector prio_lines; + for (auto line : dclust->lines) { + if(line->priority == max_priority) { + prio_lines.push_back(line); } } - // if there are no non-footprint lines, skip this cluster - if(!found_non_fp) continue; - // if we didn't find any footprint lines, pick the line with the highest elevation - if(!found_fp) { - double max_z=0; - linetype high_line; - for(auto& i : cluster.idx) { - auto& line = lines[i]; - auto z = std::get<3>(line); - if(z > max_z) { - max_z = z; - best_line = line; - best_angle = std::get<0>(best_line); + // TODO: split clusters using common planes if intersection segment(s) are present? + // TODO: compute distance clusters with 3D dists? + // TODO: output exact segments? quick solve of spikes? + // TODO: performance optimise clustering algo + // //skip if cluster only contains fp segments + // if(max_priority==1 && (dclust_size == prio_lines.size())) continue; + if(max_priority==2 && (prio_lines.size()>1)) { + std::vector other_lines; + for (auto line : dclust->lines) { + if(line->priority != 2) { + other_lines.push_back(line); } } - } - - // compute vec orthogonal to lines in this cluster - auto p0 = std::get<1>(best_line); - auto halfdist = std::get<6>(best_line); - // Vector_2 n(-1.0, std::tan(angle)); - // n = n/std::sqrt(n.squared_length();()); // normalize - Vector_2 l(std::tan(best_angle),1.0); - l = l/std::sqrt(l.squared_length()); // normalize - auto p_center = p0;// + average_dist*n; - auto p_begin = p_center + halfdist*l; - auto p_end = p_center - halfdist*l; - edges_out.push_back({ - {float(p_begin.x()), float(p_begin.y()), 0}, - {float(p_end.x()), float(p_end.y()), 0} - }); - } - - // std::vector line_clusters; - LineStringCollection merged_edges_out; - vec1i cluster_labels; - int i=0; - for(auto& cluster : dist_clusters) { - LineCluster final_cluster; - auto ref_v = cluster.ref_vec; - auto ref_p = Point_2(cluster.ref_point.x(), cluster.ref_point.y()); - - // IntervalList interval_list; - for(auto& i : cluster.idx) { - bool& is_footprint = std::get<4>(lines[i]); - if (!is_footprint) { - auto& edge = input_edges[i].first; - auto s = Point_2(edge[0][0], edge[0][1]); - auto t = Point_2(edge[1][0], edge[1][1]); - auto d1 = (s-ref_p)*ref_v; - auto d2 = (t-ref_p)*ref_v; - // interval_list.insert({d1,d2}); - auto source = ref_p + d1*ref_v; - auto target = ref_p + d2*ref_v; - merged_edges_out.push_back({ - {float(source.x()), float(source.y()), 0}, - {float(target.x()), float(target.y()), 0} - }); - cluster_labels.push_back(i); + for (auto line : prio_lines) { + double mean_angle = line->angle; + auto centroid = line->midpoint; + segment = calc_segment(centroid, mean_angle, other_lines, extension); + auto new_seg = Segment(); + new_seg[0] = {float(CGAL::to_double(segment.source().x())), float(CGAL::to_double(segment.source().y())), 0}; + new_seg[1] = {float(CGAL::to_double(segment.target().x())), float(CGAL::to_double(segment.target().y())), 0}; + regularised.push_back(new_seg); + regularised_exact.push_back(segment); } + } else { + //compute mean line with small extensions on both ends + double mean_angle = calc_mean_angle(prio_lines); + auto centroid = calc_centroid(prio_lines); + segment = calc_segment(centroid, mean_angle, dclust->lines, extension); + auto new_seg = Segment(); + new_seg[0] = {float(CGAL::to_double(segment.source().x())), float(CGAL::to_double(segment.source().y())), 0}; + new_seg[1] = {float(CGAL::to_double(segment.target().x())), float(CGAL::to_double(segment.target().y())), 0}; + regularised.push_back(new_seg); + regularised_exact.push_back(segment); } - // merge non-footprint segments that have an overlap - // std::cout << "size of intervallist: " << interval_list.size() << "\n"; - // auto merged_intervals = interval_list.merge_overlap(); - - // for (auto& s : segments_in_cluster) { - - // } - ++i; - // clip non footprint segments on overlapping footprint segments? - // output all resulting segments as a LineCluster? } - - output("merged_edges_out").set(merged_edges_out); - output("edges_out").set(edges_out); + output("angle_cluster_id").set(angle_cluster_ids); + output("dist_cluster_id").set(dist_cluster_ids); + output("priorities").set(priorities); + output("edges_out_").set(edges_out_); + // output("edges_out").set(new_segments); + // output("rings_out").set(new_rings); + // output("footprint_out").set(new_fp); } void chain(Segment& a, Segment& b, LinearRing& ring, const float& snap_threshold) { @@ -2532,112 +2478,119 @@ void chain(Segment& a, Segment& b, LinearRing& ring, const float& snap_threshold } } void RegulariseRingsNode::process(){ - // Set up vertex data (and buffer(s)) and attribute pointers - auto edges = input("edge_segments").get(); - auto ints_segments = input("ints_segments").get(); - auto ring_idx = input("ring_idx").get>>(); - auto footprint = input("footprint").get(); - // auto ring_id = input("ring_id").get(); - // auto ring_order = input("ring_order").get(); - - // SegmentCollection all_edges; - - // build vector of all input edges - // for(auto edge : edges) { - // all_edges.push_back(edge); + // // Set up vertex data (and buffer(s)) and attribute pointers + // auto edges = input("edge_segments").get(); + // auto ints_segments = input("ints_segments").get(); + // auto ring_idx = input("ring_idx").get>>(); + // auto footprint = input("footprint").get(); + // // auto ring_id = input("ring_id").get(); + // // auto ring_order = input("ring_order").get(); + + // // SegmentCollection all_edges; + + // // build vector of all input edges + // // for(auto edge : edges) { + // // all_edges.push_back(edge); + // // } + // // size_t fpi_begin = all_edges.size(); + // // SegmentCollection fp_edges; + // // for(size_t i=0; i ek_holes; + // for (auto& p : footprint) { + // ek_footprint.push_back(linereg::Point_2(p[0], p[1])); // } - // size_t fpi_begin = all_edges.size(); - // SegmentCollection fp_edges; - // for(size_t i=0; i fp_idx; + // for (size_t i=0; i < LR.get_segments(1).size(); ++i) { + // fp_idx.push_back(i); + // } + // auto ek_reg_fp = linereg::chain_ring(fp_idx, LR.get_segments(1), snap_threshold); + // output("exact_footprint_out").set(linereg::Polygon_with_holes_2(ek_reg_fp)); + // } else { + // output("exact_footprint_out").set(linereg::Polygon_with_holes_2(ek_footprint, ek_holes.begin(), ek_holes.end())); // } - // fp_edges.push_back( - // Segment({footprint[footprint.size()-1], footprint[0]}) - // ); - linereg::Polygon_2 ek_footprint; - std::vector ek_holes; - for (auto& p : footprint) { - ek_footprint.push_back(linereg::EK::Point_2(p[0], p[1])); - } - for (auto& ring : footprint.interior_rings()) { - linereg::Polygon_2 ek_hole; - for (auto& p : ring) { - ek_hole.push_back(linereg::EK::Point_2(p[0], p[1])); - } - ek_holes.push_back(ek_hole); - } - // size_t fpi_end = all_edges.size()-1; - - // get clusters from line regularisation - auto LR = linereg::LineRegulariser(); - LR.add_segments(0,edges); - LR.add_segments(1,ek_footprint, (double) fp_offset); - for (auto& hole : ek_holes) { - LR.add_segments(1,hole, (double) fp_offset); - } - LR.add_segments(2,ints_segments); - LR.dist_threshold = dist_threshold; - LR.angle_threshold = angle_threshold; - LR.cluster(weighted_avg, angle_per_distcluster); - - if (regularise_fp) { - std::vector fp_idx; - for (size_t i=0; i < LR.get_segments(1).size(); ++i) { - fp_idx.push_back(i); - } - auto ek_reg_fp = linereg::chain_ring(fp_idx, LR.get_segments(1), snap_threshold); - output("exact_footprint_out").set(linereg::Polygon_with_holes_2(ek_reg_fp)); - } else { - output("exact_footprint_out").set(linereg::Polygon_with_holes_2(ek_footprint, ek_holes.begin(), ek_holes.end())); - } - - if (recompute_rings) { - std::unordered_map exact_polygons; - for (auto& idx : ring_idx) { - exact_polygons[idx.first] = - linereg::chain_ring(idx.second, LR.get_segments(0), snap_threshold); - // std::cout << "ch ring size : "<< exact_polygons.back().size() << ", " << exact_polygons.back().is_simple() << "\n"; - } - // std::cout << "ch fp size : "<< exact_fp.size() << ", " << exact_fp.is_simple() << "\n"; - output("exact_rings_out").set(exact_polygons); - - - LinearRingCollection lrc; - vec1i plane_ids; - for (auto& poly : exact_polygons) { - LinearRing lr; - for (auto p=poly.second.vertices_begin(); p!=poly.second.vertices_end(); ++p) { - lr.push_back({ - float(CGAL::to_double(p->x())), - float(CGAL::to_double(p->y())), - 0 - }); - plane_ids.push_back(poly.first); - } - lrc.push_back(lr); - } - output("rings_out").set(lrc); - output("plane_id").set(plane_ids); - } - auto& new_segments = vector_output("edges_out"); - vec1i priorities; - for(auto& kv : LR.segments) { - for(auto& ekseg : kv.second) { - auto new_seg = Segment(); - new_seg[0] = {float(CGAL::to_double(ekseg.source().x())), float(CGAL::to_double(ekseg.source().y())), 0}; - new_seg[1] = {float(CGAL::to_double(ekseg.target().x())), float(CGAL::to_double(ekseg.target().y())), 0}; - new_segments.push_back(new_seg); - priorities.push_back(kv.first); - priorities.push_back(kv.first); - } - } + // if (recompute_rings) { + // std::unordered_map exact_polygons; + // for (auto& idx : ring_idx) { + // exact_polygons[idx.first] = + // linereg::chain_ring(idx.second, LR.get_segments(0), snap_threshold); + // // std::cout << "ch ring size : "<< exact_polygons.back().size() << ", " << exact_polygons.back().is_simple() << "\n"; + // } + // // std::cout << "ch fp size : "<< exact_fp.size() << ", " << exact_fp.is_simple() << "\n"; + // output("exact_rings_out").set(exact_polygons); + + + // LinearRingCollection lrc; + // vec1i plane_ids; + // for (auto& poly : exact_polygons) { + // LinearRing lr; + // for (auto p=poly.second.vertices_begin(); p!=poly.second.vertices_end(); ++p) { + // lr.push_back({ + // float(CGAL::to_double(p->x())), + // float(CGAL::to_double(p->y())), + // 0 + // }); + // plane_ids.push_back(poly.first); + // } + // lrc.push_back(lr); + // } + // output("rings_out").set(lrc); + // output("plane_id").set(plane_ids); + // } + // auto& new_segments = vector_output("edges_out"); + // SegmentCollection edges_out_; + // vec1i priorities, angle_cluster_ids, dist_cluster_ids; + // for(auto& line : LR.lines) { + // auto& ekseg = line.reg_segment; + // auto new_seg = Segment(); + // new_seg[0] = {float(CGAL::to_double(ekseg.source().x())), float(CGAL::to_double(ekseg.source().y())), 0}; + // new_seg[1] = {float(CGAL::to_double(ekseg.target().x())), float(CGAL::to_double(ekseg.target().y())), 0}; + // new_segments.push_back(new_seg); + // edges_out_.push_back(new_seg); + // priorities.push_back(line.priority); + // priorities.push_back(line.priority); + // angle_cluster_ids.push_back(line.angle_cluster_id); + // angle_cluster_ids.push_back(line.angle_cluster_id); + // dist_cluster_ids.push_back(line.dist_cluster_id); + // dist_cluster_ids.push_back(line.dist_cluster_id); + // } - output("priorities").set(priorities); - // output("edges_out").set(new_segments); - // output("rings_out").set(new_rings); - // output("footprint_out").set(new_fp); + // output("angle_cluster_id").set(angle_cluster_ids); + // output("dist_cluster_id").set(dist_cluster_ids); + // output("priorities").set(priorities); + // output("edges_out_").set(edges_out_); + // // output("edges_out").set(new_segments); + // // output("rings_out").set(new_rings); + // // output("footprint_out").set(new_fp); } LinearRing simplify_footprint(const LinearRing& polygon, float& threshold_stop_cost) { @@ -2702,6 +2655,48 @@ LinearRing simplify_footprint(const LinearRing& polygon, float& threshold_stop_c return polygon; } +bool get_line_extend( + Kernel::Line_3* l, + std::vector& points, + double& dmin, + double& dmax, + Point& pmin, + Point& pmax, + float& min_dist_to_line_sq +) { + size_t cnt = 0; + bool setminmax=false; + double sqd_min = 1 + min_dist_to_line_sq; + + auto lp = l->point(); + auto lv = l->to_vector(); + lv = lv / CGAL::sqrt(lv.squared_length()); + + for (auto p : points) { + auto sqd = CGAL::squared_distance(*l, p); + if (sqd > min_dist_to_line_sq) + continue; + auto av = p-lp; + auto d = av*lv; + if (!setminmax) { + setminmax=true; + dmin=dmax=d; + pmin=pmax=p; + } + if (d < dmin){ + dmin = d; + pmin = p; + } + if (d > dmax) { + dmax = d; + pmax = p; + } + sqd_min = std::min(sqd_min, sqd); + ++cnt; + } + return cnt > 1; +} + void PlaneIntersectNode::process() { auto pts_per_roofplane = input("pts_per_roofplane").get(); auto plane_adj = input("plane_adj").get>>(); @@ -2713,47 +2708,38 @@ void PlaneIntersectNode::process() { size_t ring_cntr=0; for(auto& [id_hi, ids_lo] : plane_adj) { auto& plane_hi = pts_per_roofplane[id_hi].first; - auto& plane_pts = pts_per_roofplane[id_hi].second; + auto& plane_pts_hi = pts_per_roofplane[id_hi].second; // auto& alpha_ring = alpha_rings[ring_cntr++]; for (auto& [id_lo, cnt] : ids_lo) { // skip plain pairs with low number of neighbouring points if (cnt < min_neighb_pts) continue; auto& plane_lo = pts_per_roofplane[id_lo].first; + auto& plane_pts_lo = pts_per_roofplane[id_lo].second; auto result = CGAL::intersection(plane_hi, plane_lo); if (result) { // bound the line to extend of one plane's inliers if (auto l = boost::get(&*result)) { - auto lp = l->point(); - auto lv = l->to_vector(); - lv = lv / CGAL::sqrt(lv.squared_length()); - double dmin, dmax; - Point pmin, pmax; - bool setminmax=false; - double sqd_min = 1 + min_dist_to_line_sq; - for (auto p : plane_pts) { - auto av = p-lp; - auto d = av*lv; - if (!setminmax) { - setminmax=true; - dmin=dmax=d; - pmin=pmax=p; - } - if (d < dmin){ - dmin = d; - pmin = p; - } - if (d > dmax) { - dmax = d; - pmax = p; - } - sqd_min = std::min(sqd_min, CGAL::squared_distance(*l, p)); - } + Point pmin_lo, pmax_lo; + Point pmin_hi, pmax_hi; + double dmin_lo, dmax_lo; + double dmin_hi, dmax_hi; + // skip this line if it is too far away any of the plane_pts - if(sqd_min > min_dist_to_line_sq) + if( !get_line_extend(l, plane_pts_hi, dmin_hi, dmax_hi, pmin_hi, pmax_hi, min_dist_to_line_sq) || + !get_line_extend(l, plane_pts_lo, dmin_lo, dmax_lo, pmin_lo, pmax_lo, min_dist_to_line_sq) ) continue; - auto ppmin = l->projection(pmin); - auto ppmax = l->projection(pmax); + // take the overlap between the two extends + Point ppmin, ppmax; + if (dmin_lo > dmin_hi) + ppmin = l->projection(pmin_lo); + else + ppmin = l->projection(pmin_hi); + + if (dmax_lo < dmax_hi) + ppmax = l->projection(pmax_lo); + else + ppmax = l->projection(pmax_hi); // Check for infinity (quick fix for linux crash) auto sx = float(CGAL::to_double(ppmin.x())); diff --git a/src/stepedge_nodes.hpp b/src/stepedge_nodes.hpp index b401b46..6f30103 100644 --- a/src/stepedge_nodes.hpp +++ b/src/stepedge_nodes.hpp @@ -141,6 +141,10 @@ namespace geoflow::nodes::stepedge { add_poly_input("attributes", {typeid(bool), typeid(int), typeid(float), typeid(std::string)}); add_poly_output("attributes", {typeid(bool), typeid(int), typeid(float), typeid(std::string)}); add_vector_output("linear_rings", typeid(LinearRing)); + add_vector_output("plane_a", typeid(float)); + add_vector_output("plane_b", typeid(float)); + add_vector_output("plane_c", typeid(float)); + add_vector_output("plane_d", typeid(float)); } void process(); }; @@ -285,7 +289,7 @@ namespace geoflow::nodes::stepedge { // ImGui::Text("Arrangement valid? %s", arr_is_valid? "yes" : "no"); // ImGui::Text("vcount: %d, ecount: %d", vcount, ecount); // } - void arr_snapclean(Arrangement_2& arr); + // void arr_snapclean(Arrangement_2& arr); void arr_snapclean_from_fp(Arrangement_2& arr); void arr_process(Arrangement_2& arr); void arr_assign_pts_to_unsegmented(Arrangement_2& arr, std::vector& points); @@ -295,17 +299,24 @@ namespace geoflow::nodes::stepedge { class BuildArrFromLinesNode:public Node { float rel_area_thres = 0.1; int max_arr_complexity = 400; + // bool snap_clean = true; + // bool snap_detect_only = false; + // float snap_dist = 1.0; public: using Node::Node; void init() { - add_vector_input("lines", typeid(Segment)); + add_vector_input("lines", {typeid(Segment), typeid(linereg::Segment_2)}); add_input("footprint", {typeid(linereg::Polygon_with_holes_2), typeid(LinearRing)}); add_output("arrangement", typeid(Arrangement_2)); add_output("arr_complexity", typeid(int)); add_param(ParamBoundedFloat(rel_area_thres, 0.01, 1, "rel_area_thres", "Preserve split ring area")); add_param(ParamInt(max_arr_complexity, "max_arr_complexity", "Maximum nr of lines")); + + // add_param(ParamBool(snap_clean, "snap_clean", "Snap")); + // add_param(ParamBool(snap_detect_only, "snap_detect_only", "snap_detect_only")); + // add_param(ParamBoundedFloat(snap_dist, 0.01, 5, "snap_dist", "Snap distance")); } void process(); }; @@ -417,6 +428,7 @@ namespace geoflow::nodes::stepedge { void init() { add_input("edge_points", {typeid(PointCollection), typeid(LinearRingCollection)}); add_input("roofplane_ids", typeid(vec1i)); + add_input("pts_per_roofplane", typeid(IndexedPlanesWithPoints )); add_output("edge_segments", typeid(SegmentCollection)); add_output("lines3d", typeid(SegmentCollection)); add_output("ring_edges", typeid(SegmentCollection)); @@ -436,7 +448,7 @@ namespace geoflow::nodes::stepedge { add_param(ParamBool(remove_overlap, "remove_overlap", "Remove overlap")); } inline void detect_lines_ring_m1(linedect::LineDetector& LD, SegmentCollection& segments_out); - inline size_t detect_lines_ring_m2(linedect::LineDetector& LD, SegmentCollection& segments_out); + inline size_t detect_lines_ring_m2(linedect::LineDetector& LD, Plane&, SegmentCollection& segments_out); inline void detect_lines(linedect::LineDetector& LD); void process(); }; @@ -581,39 +593,29 @@ namespace geoflow::nodes::stepedge { void process(); }; - struct LineCluster { - // source_id, target_id, is_footprint - typedef std::tuple SegmentTuple; - Vector_2 ref_vec; // direction of reference line for this cluster (exact_exact arrangent traits) - Point_2 ref_point; // reference point on reference line - vec1f vertices; // stored as distances from ref_point in direction of ref_vec - std::vector segments; - }; - - struct ValueCluster { - Vector_2 ref_vec; - Vector_2 ref_point; - std::vector idx; - }; - class RegulariseLinesNode:public Node { - static constexpr double pi = 3.14159265358979323846; float dist_threshold = 0.5; float angle_threshold = 0.15; + float extension = 1.0; public: using Node::Node; void init() { add_input("edge_segments", typeid(SegmentCollection)); - add_input("footprint", typeid(LinearRing)); - add_output("edges_out", typeid(LineStringCollection)); - add_output("merged_edges_out", typeid(LineStringCollection)); - add_output("cluster_labels", typeid(vec1i)); - // add_output("footprint_labels", typeid(vec1i)); - // add_output("line_clusters", TT_any); // ie a LineCluster - // add_output("tmp_vec3f", typeid(vec3f)); + add_input("ints_segments", typeid(SegmentCollection)); + // add_input("footprint", typeid(LinearRing)); + + add_vector_output("regularised_edges", typeid(Segment)); + add_vector_output("exact_regularised_edges", typeid(linereg::Segment_2)); + add_output("edges_out_", typeid(SegmentCollection)); + add_output("priorities", typeid(vec1i)); + add_output("angle_cluster_id", typeid(vec1i)); + add_output("dist_cluster_id", typeid(vec1i)); + add_output("exact_footprint_out", typeid(linereg::Polygon_with_holes_2)); + add_param(ParamFloat(dist_threshold, "dist_threshold", "Distance threshold")); - add_param(ParamBoundedFloat(angle_threshold, 0.01, pi, "angle_threshold", "Angle threshold")); + add_param(ParamFloat(extension, "extension", "Line extension after regularisation")); + add_param(ParamBoundedFloat(angle_threshold, 0.01, 3.1415, "angle_threshold", "Angle threshold")); } void process(); }; @@ -639,7 +641,10 @@ namespace geoflow::nodes::stepedge { // add_input("edge_segments", typeid(SegmentCollection)); add_input("footprint", typeid(LinearRing)); add_vector_output("edges_out", typeid(Segment)); + add_output("edges_out_", typeid(SegmentCollection)); add_output("priorities", typeid(vec1i)); + add_output("angle_cluster_id", typeid(vec1i)); + add_output("dist_cluster_id", typeid(vec1i)); // add_output("rings_out", typeid(LinearRingCollection)); // add_output("footprint_out", typeid(LinearRing)); add_output("rings_out", typeid(LinearRingCollection));