From 728b4a2f5f71d7b27f6ef6e5a9e1d4c5aa68e1f3 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 23 Jun 2015 14:59:50 +0200 Subject: [PATCH] Replace incremental search with a range query and a k-neighbor search --- ...cale_space_surface_reconstruction_3_impl.h | 61 ++++++++++++- .../Scale_space_surface_reconstruction_3.h | 3 +- .../rangequery-incremental.cpp | 45 ++++++++++ Spatial_searching/include/CGAL/Kd_tree_node.h | 6 +- .../test/Spatial_searching/CMakeLists.txt | 86 +++++++++++++++++++ .../Spatial_searching/rangeVSincremental.cpp | 77 +++++++++++++++++ 6 files changed, 270 insertions(+), 8 deletions(-) create mode 100644 Spatial_searching/benchmark/Spatial_searching/rangequery-incremental.cpp create mode 100644 Spatial_searching/test/Spatial_searching/CMakeLists.txt create mode 100644 Spatial_searching/test/Spatial_searching/rangeVSincremental.cpp diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h index 8da661da552..7cb7b0002a4 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h @@ -24,6 +24,10 @@ #include "tbb/parallel_for.h" #endif // CGAL_LINKED_WITH_TBB +#include +#include +#include + namespace CGAL { template < class Gt, class FS, class Sh, class wA, class Ct > @@ -71,6 +75,21 @@ public: }; // In_surface_tester +struct Inc { + unsigned int * i; + + Inc(unsigned int& i) + : i(&i) + {} + + template + void operator()(const T&) const + { + ++(*i); + } + +}; + // Compute the number of neighbors of a point that lie within a fixed radius. template < class Gt, class FS, class Sh, class wA, class Ct > class @@ -85,13 +104,19 @@ private: const Pointset& _pts; const Search_tree& _tree; - const FT& _sq_rd; + const FT _sq_rd; CountVec& _nn; public: ComputeNN(const Pointset& points, const Search_tree& tree, const FT& sq_radius, CountVec& nn) - : _pts(points), _tree(tree), _sq_rd(sq_radius), _nn(nn) {} + : _pts(points), _tree(tree), _sq_rd( +#if 1 + sqrt(sq_radius) +#else +sq_radius +#endif +), _nn(nn) {} #ifdef CGAL_LINKED_WITH_TBB void operator()( const tbb::blocked_range< std::size_t >& range ) const { @@ -100,12 +125,26 @@ public: } #endif // CGAL_LINKED_WITH_TBB void operator()( const std::size_t& i ) const { +#if 1 + Sphere sp(_pts[i], _sq_rd, 0); + + Inc inc(_nn[i]); + _tree.search(boost::make_function_output_iterator(inc),sp); +#endif +#if 0 // Iterate over the neighbors until the first one is found that is too far. Dynamic_search search( _tree, _pts[i] ); for( typename Dynamic_search::iterator nit = search.begin(); nit != search.end() && compare( _pts[i], nit->first, _sq_rd ) != LARGER; ++nit ) ++_nn[i]; + /* + if(count != _nn[i]){ + std::string s = boost::lexical_cast(count) + " != " + boost::lexical_cast(_nn[i]) + "\n"; + std::cerr << s; + } + */ +#endif } }; // class ComputeNN @@ -141,7 +180,16 @@ public: // If the neighborhood is too small, the vertex is not moved. if( _nn[i] < 4 ) return; - +#if 1 + Static_search search(_tree, _pts[i], _nn[i]); + Approximation pca( _nn[i] ); + unsigned int column = 0; + for( typename Static_search::iterator nit = search.begin(); + nit != search.end() && column < _nn[i]; + ++nit, ++column ) { + pca.set_point( column, nit->first, 1.0 / _nn[ _ind.at( nit->first ) ] ); + } +#else // Collect the vertices within the ball and their weights. Dynamic_search search( _tree, _pts[i] ); Approximation pca( _nn[i] ); @@ -152,6 +200,7 @@ public: pca.set_point( column, nit->first, 1.0 / _nn[ _ind.at( nit->first ) ] ); } CGAL_assertion( column == _nn[i] ); +#endif // Compute the weighted least-squares planar approximation of the point set. if( !pca.compute() ) @@ -409,7 +458,11 @@ increase_scale( unsigned int iterations ) { // This can be done concurrently. CountVec neighbors( _tree.size(), 0 ); try_parallel( ComputeNN( _points, _tree, _squared_radius, neighbors ), 0, _tree.size() ); - + /* + for(int i =0; i < neighbors.size(); i++){ + std::cerr << neighbors[i] << std::endl; + } + */ // Construct a mapping from each point to its index. PIMap indices; std::size_t index = 0; diff --git a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h index f0e01b8a081..35a5b8c0426 100644 --- a/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h +++ b/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -109,7 +110,7 @@ private: typedef Orthogonal_incremental_neighbor_search< Search_traits > Dynamic_search; typedef typename Dynamic_search::Tree Search_tree; - + typedef Fuzzy_sphere< Search_traits > Sphere; typedef CGAL::Random Random; // Constructing the surface. diff --git a/Spatial_searching/benchmark/Spatial_searching/rangequery-incremental.cpp b/Spatial_searching/benchmark/Spatial_searching/rangequery-incremental.cpp new file mode 100644 index 00000000000..4630e09c977 --- /dev/null +++ b/Spatial_searching/benchmark/Spatial_searching/rangequery-incremental.cpp @@ -0,0 +1,45 @@ +#include +#include +#include +#include + +const int D = 3; + +typedef CGAL::Simple_cartesian K; +typedef K::Point_3 Point_d; +typedef CGAL::Search_traits_d > Traits; +typedef CGAL::Random_points_in_cube_d Random_points_iterator; +typedef CGAL::Counting_iterator N_Random_points_iterator; +typedef CGAL::Kd_tree Tree; +typedef CGAL::Fuzzy_sphere Fuzzy_sphere; +typedef CGAL::Fuzzy_iso_box Fuzzy_iso_box; + +int main() { + + const int N = 1000; + // generator for random data points in the square ( (-1000,-1000), (1000,1000) ) + Random_points_iterator rpit(4, 1000.0); + + // Insert N points in the tree + Tree tree(N_Random_points_iterator(rpit,0), + N_Random_points_iterator(rpit,N)); + + // define range query objects + double pcoord[D] = { 300, 300, 300, 300 }; + double qcoord[D] = { 900.0, 900.0, 900.0, 900.0 }; + Point_d p(D, pcoord+0, pcoord+D); + Point_d q(D, qcoord+0, qcoord+D); + Fuzzy_sphere fs(p, 700.0, 100.0); + Fuzzy_iso_box fib(p, q, 100.0); + + std::cout << "points approximately in fuzzy spherical range query" << std::endl; + std::cout << "with center (300, 300, 300, 300)" << std::endl; + std::cout << "and fuzzy radius [600, 800] are:" << std::endl; + tree.search(std::ostream_iterator(std::cout, "\n"), fs); + + std::cout << "points approximately in fuzzy rectangular range query "; + std::cout << "[[200, 400], [800,1000]]^4 are:" << std::endl; + + tree.search(std::ostream_iterator(std::cout, "\n"), fib); + return 0; +} diff --git a/Spatial_searching/include/CGAL/Kd_tree_node.h b/Spatial_searching/include/CGAL/Kd_tree_node.h index c8554942845..f0b2f785fbe 100644 --- a/Spatial_searching/include/CGAL/Kd_tree_node.h +++ b/Spatial_searching/include/CGAL/Kd_tree_node.h @@ -120,8 +120,8 @@ namespace CGAL { else { Internal_node_const_handle node = static_cast(this); - it=node->lower()->tree_items(it); - it=node->upper()->tree_items(it); + it=node->lower()->tree_items(it); + it=node->upper()->tree_items(it); } return it; } @@ -173,7 +173,7 @@ namespace CGAL { if (node->size()>0) for (iterator i=node->begin(); i != node->end(); i++) if (q.contains(*i)) - {*it=*i; ++it;} + {*it++=*i;} } else { Internal_node_const_handle node = diff --git a/Spatial_searching/test/Spatial_searching/CMakeLists.txt b/Spatial_searching/test/Spatial_searching/CMakeLists.txt new file mode 100644 index 00000000000..f82bf6e23ca --- /dev/null +++ b/Spatial_searching/test/Spatial_searching/CMakeLists.txt @@ -0,0 +1,86 @@ +# Created by the script cgal_create_CMakeLists +# This is the CMake script for compiling a set of CGAL applications. + +project( Spatial_searching ) + + +cmake_minimum_required(VERSION 2.6.2) +if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" VERSION_GREATER 2.6) + if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}.${CMAKE_PATCH_VERSION}" VERSION_GREATER 2.8.3) + cmake_policy(VERSION 2.8.4) + else() + cmake_policy(VERSION 2.6) + endif() +endif() + +set( CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true ) + +if ( COMMAND cmake_policy ) + + cmake_policy( SET CMP0003 NEW ) + +endif() + +# CGAL and its components +find_package( CGAL QUIET COMPONENTS ) + +if ( NOT CGAL_FOUND ) + + message(STATUS "This project requires the CGAL library, and will not be compiled.") + return() + +endif() + +# include helper file +include( ${CGAL_USE_FILE} ) + + +# Boost and its components +find_package( Boost REQUIRED ) + +if ( NOT Boost_FOUND ) + + message(STATUS "This project requires the Boost library, and will not be compiled.") + + return() + +endif() + +# include for local directory + +# include for local package +include_directories( BEFORE ../../include ) + + +# Creating entries for all .cpp/.C files with "main" routine +# ########################################################## + +include( CGAL_CreateSingleSourceCGALProgram ) + +create_single_source_cgal_program( "Building_kd_tree_with_ddim_points.cpp" ) + +create_single_source_cgal_program( "Building_kd_tree_with_own_pointtype.cpp" ) + +create_single_source_cgal_program( "Circular_query.cpp" ) + +create_single_source_cgal_program( "Compare_methods.cpp" ) + +create_single_source_cgal_program( "Iso_rectangle_2_query.cpp" ) + +create_single_source_cgal_program( "iso_rectangle_2_query_2.cpp" ) + +create_single_source_cgal_program( "K_neighbor_search_manhattan_distance_isobox_point.cpp" ) + +create_single_source_cgal_program( "K_neighbor_search_with_circle.cpp" ) + +create_single_source_cgal_program( "Orthogonal_incremental_neighbor_search.cpp" ) + +create_single_source_cgal_program( "Orthogonal_k_neighbor_search.cpp" ) + +create_single_source_cgal_program( "Range_searching.cpp" ) + +create_single_source_cgal_program( "rangeVSincremental.cpp" ) + +create_single_source_cgal_program( "Splitters.cpp" ) + + diff --git a/Spatial_searching/test/Spatial_searching/rangeVSincremental.cpp b/Spatial_searching/test/Spatial_searching/rangeVSincremental.cpp new file mode 100644 index 00000000000..cd413a66df1 --- /dev/null +++ b/Spatial_searching/test/Spatial_searching/rangeVSincremental.cpp @@ -0,0 +1,77 @@ +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian K; +typedef CGAL::Search_traits_3 TreeTraits; +typedef CGAL::Orthogonal_incremental_neighbor_search NN_incremental_search; +typedef NN_incremental_search::iterator NN_iterator; +typedef NN_incremental_search::Tree Tree; +#include + + +typedef CGAL::Simple_cartesian K; +typedef K::Point_3 Point; + +struct Inc { + unsigned int * i; + + Inc(unsigned int& i) + : i(&i) + {} + + // template + void operator()(const Point& t) const + { + std::cout << t << std::endl; + ++(*i); + } + +}; + + +int main(int argc, char* argv[]) { + + std::cout.precision(17); + std::ifstream in(argv[1]); + std::vector points; + Point p; + while(in >> p){ + points.push_back(p); + } + + Tree tree; + tree.insert(points.begin(), points.end()); + + Point query(5,2,1); + + NN_incremental_search NN(tree, query); + + double sd = 0; + std::cout << "The first 20 nearest neighbours with positive x-coord are: " << std::endl; + NN_iterator it = NN.begin(); + for (int i=0; i<20; i++){ + sd = (*it).second; + std::cout << (*it).first << " at squared distance = " << (*it).second << std::endl; + ++it; + } + + CGAL::Fuzzy_sphere fs(query,sqrt(sd)); + std::vector result; + tree.search(std::back_inserter(result), fs); + for(int i=0; i < result.size(); i++){ + std::cout << result[i] << " " << squared_distance(query, result[i]) << std::endl; + } + + + unsigned int count = 0; + Inc inc(count); + tree.search(boost::make_function_output_iterator(inc),fs); + std::cout << count << std::endl; + + + return 0; +}