diff --git a/Point_set_processing_3/examples/Point_set_processing_3/regularize_and_simplify_point_set_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/regularize_and_simplify_point_set_example.cpp index 22ba632bd19..07ae63d4782 100644 --- a/Point_set_processing_3/examples/Point_set_processing_3/regularize_and_simplify_point_set_example.cpp +++ b/Point_set_processing_3/examples/Point_set_processing_3/regularize_and_simplify_point_set_example.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -38,8 +39,9 @@ int main(void) std::vector points_sampled; points_sampled.resize(points.size() * (retain_percentage / 100.)); - int starttime, stoptime, timeused; - starttime = clock(); + CGAL::Timer task_timer; + task_timer.start(); + std::cout << "Run algorithm example: " << std::endl; // Run algorithm and copy results to sample points using kdtree //std::copy(CGAL::regularize_and_simplify_point_set( @@ -52,7 +54,7 @@ int main(void) // points.end(), // points_sampled.begin()); - // Run algorithm + // Run algorithm using balltree std::vector::const_iterator sample_points_begin = CGAL::regularize_and_simplify_point_set_using_balltree( points.begin(), @@ -61,15 +63,16 @@ int main(void) neighbor_radius, iter_number, need_compute_density); - // Copy results to sample points using balltree + // Copy results to sample points std::copy(sample_points_begin, static_cast::const_iterator>(points.end()), points_sampled.begin()); - stoptime = clock(); - timeused = stoptime - starttime; - std::cout << "##" << " time used: " << timeused / double(CLOCKS_PER_SEC) << " seconds." << std::endl; - + + long memory = CGAL::Memory_sizer().virtual_size(); + std::cout << "done: " << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated" << std::endl; + task_timer.stop(); // Saves point set. // Note: write_xyz_points_and_normals() requires an output iterator @@ -91,40 +94,3 @@ int main(void) - - - - -template -void find_original_neighbors( - typename IteratorType::iterator starta, - typename IteratorType::iterator enda, - typename IteratorType::iterator startb, - typename IteratorType::iterator endb, - const typename Kernel::FT radius) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - FT radius2 = radius*radius; - FT iradius16 = -4/radius2; - const FT PI = 3.1415926; - - for(CGrid::iterator dest = starta; dest != enda; dest++) - { - RichPoint &v = *(*dest); - - Point &p = v.pt; - for(CGrid::iterator origin = startb; origin != endb; origin++) - { - RichPoint &t = *(*origin); - Point &q = t.pt; - - double dist2 = CGAL::squared_distance(p, q); - - if(dist2 < radius2) - { - v.neighbors.push_back((*origin)->index); - } - } - } -} \ No newline at end of file diff --git a/Point_set_processing_3/include/CGAL/Rich_grid.h b/Point_set_processing_3/include/CGAL/Rich_grid.h new file mode 100644 index 00000000000..cafd74103ac --- /dev/null +++ b/Point_set_processing_3/include/CGAL/Rich_grid.h @@ -0,0 +1,516 @@ +// Copyright (c) 2013-06 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Shihao Wu, Cl¨Śment Jamin + +#ifndef CGAL_RICH_GRID_H +#define CGAL_RICH_GRID_H + +#include +#include + +#include +#include +#include +#include + + +//#include +//#include + +namespace CGAL { + +/// \cond SKIP_IN_MANUAL + +// ---------------------------------------------------------------------------- +// Rich Grid section +// ---------------------------------------------------------------------------- +//namespace rich_grid_internel{ + +namespace rich_grid_internel{ + +/// The Rich_point class represents a 3D point with inedxes of neighbor points; +/// - a position, +/// - a index. +/// - self point set neighbors. +/// - other point set neighbors. +/// +/// @heading Parameters: +/// @param Kernel Geometric traits class. +template +class Rich_point +{ +public: + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + +public: + Rich_point(){} + Rich_point(const Point& p, const int& i):pt(p),index(i){} + + Point pt; + unsigned int index; + std::vector neighbors; + std::vector original_neighbors; +}; + +template +class Rich_box +{ +public: + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + + Rich_box(){init();} + + void init() + { + FT max_double = (FT)(std::numeric_limits::max)(); + min_x = min_y = min_z = max_double; + max_x = max_y = max_z = -max_double; + } + + void add_point(const Point& p) + { + if (p.x() < min_x) min_x = p.x(); + if (p.y() < min_y) min_y = p.y(); + if (p.z() < min_z) min_z = p.z(); + if (p.x() > max_x) max_x = p.x(); + if (p.y() > max_y) max_y = p.y(); + if (p.z() > max_z) max_z = p.z(); + } + + Point get_min(){return Point(min_x, min_y, min_z);} + Point get_max(){return Point(max_x, max_y, max_z);} + +private: + FT min_x, min_y, min_z, max_x, max_y, max_z; +}; + + +template +class Rich_grid +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + typedef std::vector*> Point_Pointer_vector; + typedef typename Point_Pointer_vector::iterator iterator; + +public: + + Rich_grid() {} + + void init(std::vector > &vert, + Rich_box& box, + const FT _radius); + + // Travel for the point set itself + void travel_itself(void (*self)(iterator starta, iterator enda, FT radius), + void (*other)(iterator starta, iterator enda, + iterator startb, iterator endb, FT radius)); + + // Travel other self between two point set(original and samples) + void travel_others(Rich_grid &points, + void (*travel_others)(iterator starta, iterator enda, + iterator startb, iterator endb, FT radius)); + + + void static __cdecl find_original_neighbors(iterator starta, iterator enda, + iterator startb, iterator endb, FT radius); + void static __cdecl find_self_neighbors(iterator start, + iterator end, FT radius); + void static __cdecl find_other_neighbors(iterator starta, iterator enda, + iterator startb, iterator endb, FT radius); + +private: + + std::vector*> rich_points; + std::vector index; + int xside, yside, zside; + FT radius; + + int cell(int x, int y, int z) { return x + xside*(y + yside*z); } + bool isEmpty(int cell) { return index[cell+1] == index[cell]; } + iterator startV(int origin) { return rich_points.begin() + index[origin]; } + iterator endV(int origin) { return rich_points.begin() + index[origin+1]; } +}; + + +template +class XSort { +public: + bool operator()(const Rich_point *a, const Rich_point *b) { + return a->pt.x() < b->pt.x(); + } +}; + +template +class YSort { +public: + bool operator()(const Rich_point *a, const Rich_point *b) { + return a->pt.y() < b->pt.y(); + } +}; + +template +class ZSort { +public: + bool operator()(const Rich_point *a, const Rich_point *b) { + return a->pt.z() < b->pt.z(); + } +}; + + +// divide spoints into some grids +// and each grid has their points index in the index vector of sample. +template +void Rich_grid::init(std::vector > &vert, + Rich_box& box, const typename Kernel::FT _radius) +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + + radius = _radius; + + Point min = box.get_min(); + Point max = box.get_max(); + + rich_points.resize(vert.size()); + for(int i = 0; i < rich_points.size(); i++) + { + Point& pt = vert[i].pt; + rich_points[i] = &vert[i]; + } + + xside = (int)ceil((max.x() - min.x())/radius); + yside = (int)ceil((max.y() - min.y())/radius); + zside = (int)ceil((max.z() - min.z())/radius); + + xside = (xside > 0) ? xside : 1; + yside = (yside > 0) ? yside : 1; + zside = (zside > 0) ? zside : 1; + + index.resize(xside*yside*zside + 1, -1); + + std::sort(rich_points.begin(), rich_points.end(), ZSort()); + + unsigned int startz = 0; + for(unsigned int z = 0; z < zside; z++) + { + int endz = startz; + FT maxz = min.z() + (z+1)*radius; + while(endz < rich_points.size() && rich_points[endz]->pt.z()< maxz) + ++endz; + + sort(rich_points.begin()+startz, + rich_points.begin()+endz, YSort()); + + int starty = startz; + for(int y = 0; y < yside; y++) + { + int endy = starty; + FT maxy = min.y() + (y+1)*radius; + while(endy < endz && rich_points[endy]->pt.y() < maxy) + ++endy; + + sort(rich_points.begin()+starty, + rich_points.begin()+endy, XSort()); + + int startx = starty; + for(int x = 0; x < xside; x++) + { + int endx = startx; + index[x + xside*y + xside*yside*z] = endx; + FT maxx = min.x() + (x+1)*radius; + while(endx < endy && rich_points[endx]->pt.x() < maxx) + ++endx; + + startx = endx; + } + starty = endy; + } + startz = endz; + } + + //compute the last grid's range + index[xside*yside*zside] = startz; +} + +/// define how to travel in the same gird +template +void Rich_grid::travel_itself( + void (*self)(iterator starta, iterator enda, + const typename Kernel::FT radius), + void (*other)(iterator starta, iterator enda, + iterator startb, iterator endb, FT radius) +) +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + + static int corner[8*3] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, + 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1 }; + + static int diagonals[14*2] = { 0, 0, //remove this to avoid self intesextion + 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, + 2, 3, 1, 3, 1, 2, + 1, 4, 2, 5, 3, 6 }; + + for(int z = 0; z < zside; z++) { + for(int y = 0; y < yside; y++) { + for(int x = 0; x < xside; x++) { + int origin = cell(x, y, z); + self(startV(origin), endV(origin), radius); + // compute between other girds + for(int d = 2; d < 28; d += 2) { // skipping self + int *cs = corner + 3*diagonals[d]; + int *ce = corner + 3*diagonals[d+1]; + if((x + cs[0] < xside) && (y + cs[1] < yside) && (z + cs[2] < zside) + && (x + ce[0] < xside) && (y + ce[1] < yside) && (z + ce[2] < zside)) + { + origin = cell(x+cs[0], y+cs[1], z+cs[2]); + int dest = cell(x+ce[0], y+ce[1], z+ce[2]); + other(startV(origin), endV(origin), + startV(dest), endV(dest), radius); + } + } + } + } + } +} + +/// define how to travel in other gird +template +void Rich_grid::travel_others( + Rich_grid &points, + void (*travel_others)(iterator starta, iterator enda, + iterator startb, iterator endb, + const typename Kernel::FT radius) +) +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + + static int corner[8*3] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, + 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1 }; + + static int diagonals[14*2] = { 0, 0, //remove this to avoid self intesextion + 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, + 2, 3, 1, 3, 1, 2, + 1, 4, 2, 5, 3, 6 }; + + for(int z = 0; z < zside; z++) { + for(int y = 0; y < yside; y++) { + for(int x = 0; x < xside; x++) { + int origin = cell(x, y, z); + + + if(!isEmpty(origin) && !points.isEmpty(origin)) + travel_others(startV(origin), endV(origin), + points.startV(origin), points.endV(origin), radius); + + + for(int d = 2; d < 28; d += 2) { //skipping self + int *cs = corner + 3*diagonals[d]; + int *ce = corner + 3*diagonals[d+1]; + if((x+cs[0] < xside) && (y+cs[1] < yside) && (z+cs[2] < zside) && + (x+ce[0] < xside) && (y+ce[1] < yside) && (z+ce[2] < zside)) { + + origin = cell(x+cs[0], y+cs[1], z+cs[2]); + + int dest = cell(x+ce[0], y+ce[1], z+ce[2]); + + if(!isEmpty(origin) && !points.isEmpty(dest)) // Locally + travel_others(startV(origin), endV(origin), + points.startV(dest), points.endV(dest), radius); + + if(!isEmpty(dest) && !points.isEmpty(origin)) + travel_others(startV(dest), endV(dest), + points.startV(origin), points.endV(origin), radius); + } + } + } + } + } +} + +/// grid travel function to find the neighbors in the original point set +template +void Rich_grid::find_original_neighbors( + iterator starta, + iterator enda, + iterator startb, + iterator endb, + FT radius +) +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + FT radius2 = radius*radius; + + for(Rich_grid::iterator dest = starta; dest != enda; dest++) + { + Rich_point &v = *(*dest); + + Point &p = v.pt; + for(Rich_grid::iterator origin = startb; origin != endb; origin++) + { + Rich_point &t = *(*origin); + Point &q = t.pt; + + FT dist2 = CGAL::squared_distance(p, q); + + if(dist2 < radius2) + { + v.original_neighbors.push_back((*origin)->index); + } + } + } +} + +/// grid travel function to find the neighbors in the same point set +template +void Rich_grid::find_self_neighbors( + iterator start, iterator end, FT radius) +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + FT radius2 = radius*radius; + for(iterator dest = start; dest != end; dest++) + { + Rich_point &v = *(*dest); + Point &p = v.pt; + + for(iterator origin = dest+1; origin != end; origin++) + { + Rich_point &t = *(*origin); + Point &q = t.pt; + + FT dist2 = CGAL::squared_distance(p, q); + + if(dist2 < radius2) + { + v.neighbors.push_back((*origin)->index); + t.neighbors.push_back((*dest)->index); + } + } + } +} + +/// grid travel function to find the neighbors in the same point set +template +void Rich_grid::find_other_neighbors( + iterator starta, iterator enda, + iterator startb, iterator endb, FT radius) +{ + typedef typename Kernel::Point_3 Point; + typedef typename Kernel::FT FT; + FT radius2 = radius*radius; + for(iterator dest = starta; dest != enda; dest++) + { + Rich_point &v = *(*dest); + Point &p = v.pt; + + for(iterator origin = startb; origin != endb; origin++) + { + Rich_point &t = *(*origin); + Point &q = t.pt; + + FT dist2 = CGAL::squared_distance(p, q); + + if(dist2 < radius2) + { + v.neighbors.push_back((*origin)->index); + t.neighbors.push_back((*dest)->index); + } + } + } +} + + + + +/// Compute ball neighbors for each point in the same point set. +/// +/// \pre `radius > 0` +/// +/// @tparam Kernel Geometric traits class. +/// +/// @return +template +void compute_ball_neighbors_one_self( + std::vector >& points, + Rich_box& box, + const typename Kernel::FT radius) +{ + typedef typename Kernel::FT FT; + CGAL_point_set_processing_precondition(radius > 0); + + for (unsigned int i = 0; i < points.size(); i++) + { + points[i].neighbors.clear(); + } + + Rich_grid points_grid; + points_grid.init(points, box, radius); + points_grid.travel_itself(Rich_grid::find_self_neighbors, + Rich_grid::find_other_neighbors); +} + +/// Compute ball neighbors for each (sample)points in the other point set +/// +/// \pre `radius > 0` +/// +/// @tparam Kernel Geometric traits class. +/// +/// @return +template +void compute_ball_neighbors_one_to_another( + std::vector >& samples, ///< sample point set + std::vector >& original,///< original point set + Rich_box& box, ///< bounding box + const typename Kernel::FT radius ///< neighbor radius +) +{ + typedef typename Kernel::FT FT; + if (radius < FT(0.0)) + { + return; + } + + for (unsigned int i = 0; i < samples.size(); i++) + { + samples[i].original_neighbors.clear(); + } + + Rich_grid samples_grid; + samples_grid.init(samples, box, radius); + + // here can be optimized by initial the original grid just one time. + Rich_grid original_grid; + original_grid.init(original, box, radius); + + samples_grid.travel_others(original_grid, + Rich_grid::find_original_neighbors); +} + + +} //namespace internal + +} //namespace CGAL + +#endif // CGAL_RICH_GRID_H diff --git a/Point_set_processing_3/include/CGAL/denoise_point_set.h b/Point_set_processing_3/include/CGAL/denoise_point_set.h index 1a121358f18..bb062a45963 100644 --- a/Point_set_processing_3/include/CGAL/denoise_point_set.h +++ b/Point_set_processing_3/include/CGAL/denoise_point_set.h @@ -15,7 +15,7 @@ // $URL$ // $Id$ // -// Author(s) : Shihao Wu, Clement Jamin +// Author(s) : Shihao Wu, Cl¨Śment Jamin #ifndef CGAL_DENOSISE_POINTS_WITH_NORMALS_H #define CGAL_DENOSISE_POINTS_WITH_NORMALS_H diff --git a/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set.h b/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set.h index 686c7f91c67..66cb9bed6f6 100644 --- a/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set.h +++ b/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set.h @@ -24,6 +24,8 @@ #include #include #include +#include +#include #include #include @@ -36,35 +38,6 @@ namespace CGAL { -/// \cond SKIP_IN_MANUAL - -class Timer -{ -public: - - void start(const std::string& str) - { - std::cout << std::endl; - starttime = clock(); - mid_start = clock(); - // std::cout << "@@@@@ Time Count Strat For: " << str << std::endl; - - _str = str; - } - - void end() - { - stoptime = clock(); - timeused = stoptime - starttime; - std::cout << /*endl <<*/ "@@@@ finish " << _str << " time used: " << timeused / double(CLOCKS_PER_SEC) << " seconds." << std::endl; - std::cout << std::endl; - } - -private: - int starttime, mid_start, mid_end, stoptime, timeused; - std::string _str; -}; - // ---------------------------------------------------------------------------- // Private section @@ -571,7 +544,7 @@ regularize_and_simplify_point_set( ) { CGAL_point_set_processing_precondition(k > 1); - Timer time; + Timer task_timer; // basic geometric types typedef typename Kernel::Point_3 Point; @@ -612,7 +585,9 @@ regularize_and_simplify_point_set( sample_points[i] = get(point_pmap, it); // Initiate a KD-tree search for original points - time.start("Build Original Neighbor Tree"); + task_timer.start(); + std::cout << "Initialization / Compute Density For Original" << endl; + std::vector original_treeElements; for (it = first_original_point, i=0 ; it != beyond ; ++it, ++i) { @@ -620,10 +595,8 @@ regularize_and_simplify_point_set( original_treeElements.push_back(KdTreeElement(p0,i)); } Tree original_tree(original_treeElements.begin(), original_treeElements.end()); - time.end(); // Guess spacing - time.start("Guess Neighborhood Radiuse"); FT guess_neighbor_radius = (FT)(std::numeric_limits::max)(); // Or a better max number: (numeric_limits::max)()? for(it = first_original_point; it != beyond ; ++it) { @@ -631,11 +604,10 @@ regularize_and_simplify_point_set( guess_neighbor_radius = max_spacing < guess_neighbor_radius ? max_spacing : guess_neighbor_radius; } guess_neighbor_radius *= 0.95; - time.end(); std::cout << "Guess Neighborhood Radius:" << guess_neighbor_radius << std::endl; // Compute original density weight for original points if user needed - time.start("Compute Density For Original"); + task_timer.start("Compute Density For Original"); std::vector original_density_weight_set; if (need_compute_density) { @@ -645,10 +617,8 @@ regularize_and_simplify_point_set( original_density_weight_set.push_back(density); } } - time.end(); // Compute guess KNN set - time.start("Compute guess KNN set"); std::vector guess_KNN_set; for (i=0 ; i < sample_points.size(); i++) { @@ -656,14 +626,19 @@ regularize_and_simplify_point_set( unsigned int guess_knn = regularize_and_simplify_internal::guess_KNN_number_for_original_set(p0, original_tree, k, guess_neighbor_radius); guess_KNN_set.push_back(guess_knn); } - time.end(); + + long memory = CGAL::Memory_sizer().virtual_size(); + std::cout << "done: " << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated" << std::endl << std::endl; + task_timer.stop(); for (unsigned int iter_n = 0; iter_n < iter_number; iter_n++) { + task_timer.start(); + std::cout << "Compute average term and repulsion term " << endl; + // Initiate a KD-tree search for sample points - time.start("Build Sample Neighbor Tree"); - std::vector sample_treeElements; unsigned int k_for_sample = 30; // Or it can be conducted by the "guess_neighbor_radius" for (i=0 ; i < sample_points.size(); i++) @@ -672,11 +647,10 @@ regularize_and_simplify_point_set( sample_treeElements.push_back(KdTreeElement(p0,i)); } Tree sample_tree(sample_treeElements.begin(), sample_treeElements.end()); - time.end(); // Compute sample density weight for sample points if user needed std::vector sample_density_weight_set; - time.start("Compute Density For Sample"); + // task_timer.start("Compute Density For Sample"); if (need_compute_density) { for (i=0 ; i < sample_points.size(); i++) @@ -685,29 +659,24 @@ regularize_and_simplify_point_set( sample_density_weight_set.push_back(density); } } - time.end(); // Compute average term and repulsion term for each sample points separately, // then update each sample points std::vector average_set(nb_points_sample); std::vector repulsion_set(nb_points_sample); - time.start("Compute Average Term"); + //task_timer.start("Compute Average Term"); for (i = 0; i < sample_points.size(); i++) { Point& p = sample_points[i]; - average_set[i] = regularize_and_simplify_internal::compute_average_term(p, original_tree, k, guess_neighbor_radius, original_density_weight_set); // Before speed up - //average_set[i] = regularize_and_simplify_internal::compute_average_term(p, original_tree, guess_KNN_set[i], guess_neighbor_radius, original_density_weight_set); - + average_set[i] = regularize_and_simplify_internal::compute_average_term(p, original_tree, guess_KNN_set[i], guess_neighbor_radius, original_density_weight_set); } - time.end(); - time.start("Compute Repulsion Term"); + //task_timer.start("Compute Repulsion Term"); for (i = 0; i < sample_points.size(); i++) { Point& p = sample_points[i]; repulsion_set[i] = regularize_and_simplify_internal::compute_repulsion_term(p, sample_tree, k_for_sample, guess_neighbor_radius, sample_density_weight_set); } - time.end(); for (i = 0; i < sample_points.size(); i++) { @@ -715,7 +684,12 @@ regularize_and_simplify_point_set( p = CGAL::ORIGIN + average_set[i] + (FT)0.5 * repulsion_set[i]; } - std::cout << "iterate: " << iter_n + 1 << " "<< std::endl; + long memory = CGAL::Memory_sizer().virtual_size(); + std::cout << "done: " << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated" << std::endl; + task_timer.stop(); + + std::cout << "iterate: " << iter_n + 1 << " "<< std::endl << std::endl;; } //Copy back modified sample points to original points for output diff --git a/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set_using_balltree.h b/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set_using_balltree.h index a1ff89f3e2a..67cfb3dbed3 100644 --- a/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set_using_balltree.h +++ b/Point_set_processing_3/include/CGAL/regularize_and_simplify_point_set_using_balltree.h @@ -22,12 +22,13 @@ #include #include +#include +#include +#include #include -#include #include #include -#include //#include @@ -36,498 +37,13 @@ namespace CGAL { /// \cond SKIP_IN_MANUAL -// ---------------------------------------------------------------------------- -// Testing code section -// ---------------------------------------------------------------------------- - - - -// ---------------------------------------------------------------------------- -// Rich Grid section -// ---------------------------------------------------------------------------- -//namespace rich_grid_internel -//{ -//} // ---------------------------------------------------------------------------- // Private section // ---------------------------------------------------------------------------- namespace regularize_and_simplify_internal{ -class Timer -{ -public: - void start(const std::string& str) - { - std::cout << std::endl; - starttime = clock(); - //std::cout << "@@@@@ Time Count Strat For: " << str << std::endl; - _str = str; - } - - void end() - { - stoptime = clock(); - timeused = stoptime - starttime; - std::cout << /*endl <<*/ "@@@@ finish " << _str << " time used: " << timeused / double(CLOCKS_PER_SEC) << " seconds." << std::endl; - std::cout << std::endl; - } - -private: - int starttime, mid_start, mid_end, stoptime, timeused; - std::string _str; -}; - - -// ---------------------------------------------------------------------------- -// Ball Tree section -// ---------------------------------------------------------------------------- -template -class Rich_point -{ -public: - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - -public: - Rich_point(){} - Rich_point(const Point& p, const int& i):pt(p),index(i){} - - Point pt; - unsigned int index; - std::vector neighbors; - std::vector original_neighbors;//need more memory, but make the code a littel easier to read. -}; - -template -class Rich_box -{ -public: - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - - Rich_box(){init();} - - void init() - { - FT max_double = std::numeric_limits::max(); - min_x = min_y = min_z = max_double; - max_x = max_y = max_z = -max_double; - } - - void add_point(const Point& p) - { - if (p.x() < min_x) min_x = p.x(); - if (p.y() < min_y) min_y = p.y(); - if (p.z() < min_z) min_z = p.z(); - if (p.x() > max_x) max_x = p.x(); - if (p.y() > max_y) max_y = p.y(); - if (p.z() > max_z) max_z = p.z(); - } - - Point get_min(){return Point(min_x, min_y, min_z);} - Point get_max(){return Point(max_x, max_y, max_z);} - -private: - FT min_x, min_y, min_z, max_x, max_y, max_z; -}; - - -template -class Rich_grid -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - typedef std::vector*> Point_Pointer_vector; - typedef typename Point_Pointer_vector::iterator iterator; - -public: - - Rich_grid() {} - - void init(std::vector > &vert, Rich_box& box, const FT _radius); - - // Travel for the point set itself - void travel_itself(void (*self)(iterator starta, iterator enda, FT radius), - void (*other)(iterator starta, iterator enda, - iterator startb, iterator endb, FT radius)); - - // Travel other self between two point set(original and samples) - void travel_others(Rich_grid &points, - void (*travel_others)(iterator starta, iterator enda, - iterator startb, iterator endb, FT radius)); - - - void static __cdecl find_original_neighbors(iterator starta, iterator enda, - iterator startb, iterator endb, FT radius); - void static __cdecl find_self_neighbors(iterator start, iterator end, FT radius); - void static __cdecl find_other_neighbors(iterator starta, iterator enda, - iterator startb, iterator endb, FT radius); - -private: - - std::vector*> rich_points; - std::vector index; // the start index of each grid in the sample points which is order by Zsort - int xside, yside, zside; - FT radius; - - int cell(int x, int y, int z) { return x + xside*(y + yside*z); } - bool isEmpty(int cell) { return index[cell+1] == index[cell]; } - iterator startV(int origin) { return rich_points.begin() + index[origin]; } - iterator endV(int origin) { return rich_points.begin() + index[origin+1]; } -}; - - -template -class XSort { -public: - bool operator()(const Rich_point *a, const Rich_point *b) { - return a->pt.x() < b->pt.x(); - } -}; - -template -class YSort { -public: - bool operator()(const Rich_point *a, const Rich_point *b) { - return a->pt.y() < b->pt.y(); - } -}; - -template -class ZSort { -public: - bool operator()(const Rich_point *a, const Rich_point *b) { - return a->pt.z() < b->pt.z(); - } -}; - - -// divide spoints into some grids -// and each grid has their points index in the index vector of sample. -template -void Rich_grid::init(std::vector > &vert, Rich_box& box, const typename Kernel::FT _radius) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - - radius = _radius; - - Point min = box.get_min(); - Point max = box.get_max(); - - rich_points.resize(vert.size()); - for(int i = 0; i < rich_points.size(); i++) - { - Point& pt = vert[i].pt; - rich_points[i] = &vert[i]; - } - - xside = (int)ceil((max.x() - min.x())/radius); - yside = (int)ceil((max.y() - min.y())/radius); - zside = (int)ceil((max.z() - min.z())/radius); - - xside = (xside > 0) ? xside : 1; - yside = (yside > 0) ? yside : 1; - zside = (zside > 0) ? zside : 1; - - index.resize(xside*yside*zside + 1, -1); - - std::sort(rich_points.begin(), rich_points.end(), ZSort()); //this would be very slow - - unsigned int startz = 0; - for(unsigned int z = 0; z < zside; z++) - { - int endz = startz; - FT maxz = min.z() + (z+1)*radius; - while(endz < rich_points.size() && rich_points[endz]->pt.z()< maxz) - ++endz; - - sort(rich_points.begin()+startz, rich_points.begin()+endz, YSort()); - - int starty = startz; - for(int y = 0; y < yside; y++) - { - int endy = starty; - FT maxy = min.y() + (y+1)*radius; - while(endy < endz && rich_points[endy]->pt.y() < maxy) - ++endy; - - sort(rich_points.begin()+starty, rich_points.begin()+endy, XSort()); - - int startx = starty; - for(int x = 0; x < xside; x++) - { - int endx = startx; - index[x + xside*y + xside*yside*z] = endx; - FT maxx = min.x() + (x+1)*radius; - while(endx < endy && rich_points[endx]->pt.x() < maxx) - ++endx; - - startx = endx; - } - starty = endy; - } - startz = endz; - } - index[xside*yside*zside] = startz; // in order to compute the last grid's range -} - -/// define how to travel in the same gird -template -void Rich_grid::travel_itself( - void (*self)(iterator starta, iterator enda, const typename Kernel::FT radius), - void (*other)(iterator starta, iterator enda, - iterator startb, iterator endb, FT radius) -) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - - static int corner[8*3] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, - 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1 }; - - static int diagonals[14*2] = { 0, 0, //remove this line to avoid self intesextion - 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, - 2, 3, 1, 3, 1, 2, - 1, 4, 2, 5, 3, 6 }; - - for(int z = 0; z < zside; z++) { - for(int y = 0; y < yside; y++) { - for(int x = 0; x < xside; x++) { - int origin = cell(x, y, z); - self(startV(origin), endV(origin), radius); - // compute between other girds - for(int d = 2; d < 28; d += 2) { // skipping self - int *cs = corner + 3*diagonals[d]; - int *ce = corner + 3*diagonals[d+1]; - if((x + cs[0] < xside) && (y + cs[1] < yside) && (z + cs[2] < zside) && - (x + ce[0] < xside) && (y + ce[1] < yside) && (z + ce[2] < zside)) - { - origin = cell(x+cs[0], y+cs[1], z+cs[2]); - int dest = cell(x+ce[0], y+ce[1], z+ce[2]); - other(startV(origin), endV(origin), - startV(dest), endV(dest), radius); - } - } - } - } - } -} - -/// define how to travel in other gird -template -void Rich_grid::travel_others( - Rich_grid &points, - void (*travel_others)(iterator starta, iterator enda, - iterator startb, iterator endb, - const typename Kernel::FT radius) -) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - - static int corner[8*3] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, - 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1 }; - - static int diagonals[14*2] = { 0, 0, //remove this line to avoid self intesextion - 0, 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, 6, 0, 7, - 2, 3, 1, 3, 1, 2, - 1, 4, 2, 5, 3, 6 }; - - for(int z = 0; z < zside; z++) { - for(int y = 0; y < yside; y++) { - for(int x = 0; x < xside; x++) { - int origin = cell(x, y, z); - - - if(!isEmpty(origin) && !points.isEmpty(origin)) - travel_others(startV(origin), endV(origin), - points.startV(origin), points.endV(origin), radius); - - - for(int d = 2; d < 28; d += 2) { //skipping self - int *cs = corner + 3*diagonals[d]; - int *ce = corner + 3*diagonals[d+1]; - if((x+cs[0] < xside) && (y+cs[1] < yside) && (z+cs[2] < zside) && - (x+ce[0] < xside) && (y+ce[1] < yside) && (z+ce[2] < zside)) { - - origin = cell(x+cs[0], y+cs[1], z+cs[2]); - - int dest = cell(x+ce[0], y+ce[1], z+ce[2]); - - if(!isEmpty(origin) && !points.isEmpty(dest)) // Locally - travel_others(startV(origin), endV(origin), - points.startV(dest), points.endV(dest), radius); - - if(!isEmpty(dest) && !points.isEmpty(origin)) - travel_others(startV(dest), endV(dest), - points.startV(origin), points.endV(origin), radius); - } - } - } - } - } -} - -/// grid travel function to find the neighbors in the original point set -template -void Rich_grid::find_original_neighbors( - iterator starta, - iterator enda, - iterator startb, - iterator endb, - FT radius -) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - FT radius2 = radius*radius; - - for(Rich_grid::iterator dest = starta; dest != enda; dest++) - { - Rich_point &v = *(*dest); - - Point &p = v.pt; - for(Rich_grid::iterator origin = startb; origin != endb; origin++) - { - Rich_point &t = *(*origin); - Point &q = t.pt; - - FT dist2 = CGAL::squared_distance(p, q); - - if(dist2 < radius2) - { - v.original_neighbors.push_back((*origin)->index); - } - } - } -} - -/// grid travel function to find the neighbors in the same point set -template -void Rich_grid::find_self_neighbors( - iterator start, iterator end, FT radius) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - FT radius2 = radius*radius; - for(iterator dest = start; dest != end; dest++) - { - Rich_point &v = *(*dest); - Point &p = v.pt; - - for(iterator origin = dest+1; origin != end; origin++) - { - Rich_point &t = *(*origin); - Point &q = t.pt; - - FT dist2 = CGAL::squared_distance(p, q); - - if(dist2 < radius2) - { - v.neighbors.push_back((*origin)->index); - t.neighbors.push_back((*dest)->index); - } - } - } -} - -/// grid travel function to find the neighbors in the same point set -template -void Rich_grid::find_other_neighbors( - iterator starta, iterator enda, - iterator startb, iterator endb, FT radius) -{ - typedef typename Kernel::Point_3 Point; - typedef typename Kernel::FT FT; - FT radius2 = radius*radius; - for(iterator dest = starta; dest != enda; dest++) - { - Rich_point &v = *(*dest); - Point &p = v.pt; - - for(iterator origin = startb; origin != endb; origin++) - { - Rich_point &t = *(*origin); - Point &q = t.pt; - - FT dist2 = CGAL::squared_distance(p, q); - - if(dist2 < radius2) - { - v.neighbors.push_back((*origin)->index); - t.neighbors.push_back((*dest)->index); - } - } - } -} - - - - -/// Compute ball neighbors for each point in the same point set(sample or original points) -/// -/// \pre `radius > 0` -/// -/// @tparam Kernel Geometric traits class. -/// -/// @return -template -void compute_ball_neighbors_one_self( - std::vector >& points, - Rich_box& box, - const typename Kernel::FT radius) -{ - typedef typename Kernel::FT FT; - CGAL_point_set_processing_precondition(radius > 0); - - for (unsigned int i = 0; i < points.size(); i++) - { - points[i].neighbors.clear(); - } - - Rich_grid points_grid; - points_grid.init(points, box, radius); - points_grid.travel_itself(Rich_grid::find_self_neighbors, Rich_grid::find_other_neighbors); -} - -/// Compute ball neighbors for each point(sample points) in the other point set(original points) -/// -/// \pre `radius > 0` -/// -/// @tparam Kernel Geometric traits class. -/// -/// @return -template -void compute_ball_neighbors_one_to_another( - std::vector >& samples, ///< sample point set - std::vector >& original,///< original point set - Rich_box& box, ///< bounding box - const typename Kernel::FT radius ///< neighbor radius -) -{ - typedef typename Kernel::FT FT; - if (radius < FT(0.0)) - { - return; - } - - for (unsigned int i = 0; i < samples.size(); i++) - { - samples[i].original_neighbors.clear(); - } - - Rich_grid samples_grid; - samples_grid.init(samples, box, radius); - - // can be speed up by initial the original grid just one time(as paramter of this function) - Rich_grid original_grid; - original_grid.init(original, box, radius); - - samples_grid.travel_others(original_grid, Rich_grid::find_original_neighbors); -} // ---------------------------------------------------------------------------- // WLOP algorithm section @@ -544,9 +60,9 @@ template typename Kernel::Vector_3 compute_average_term( const typename Kernel::Point_3& query, ///< 3D point to project - const std::vector > neighbor_original_points, ///< neighbor original points + const std::vector > neighbor_original_points, const typename Kernel::FT radius, ///& density_weight_set ///& density_weight_set ///< ) { CGAL_point_set_processing_precondition(radius > 0); @@ -597,7 +113,7 @@ template typename Kernel::Vector_3 compute_repulsion_term( const typename Kernel::Point_3& query, ///< 3D point to project - const std::vector > neighbor_sample_points, ///< neighbor sample points + const std::vector > neighbor_sample_points, ///< neighbor sample points const typename Kernel::FT radius, ///& density_weight_set ///< if user need density ) @@ -763,13 +279,13 @@ regularize_and_simplify_point_set_using_balltree( ) { CGAL_point_set_processing_precondition(neighbor_radius > 0); - regularize_and_simplify_internal::Timer time; + Timer task_timer; // basic geometric types typedef typename Kernel::Point_3 Point; typedef typename Kernel::Vector_3 Vector; typedef typename Kernel::FT FT; - typedef typename regularize_and_simplify_internal::Rich_point Rich_point; + typedef typename rich_grid_internel::Rich_point Rich_point; // precondition: at least one element in the container. // to fix: should have at least three distinct points @@ -807,7 +323,7 @@ regularize_and_simplify_point_set_using_balltree( std::vector original_rich_point_set(nb_points_original); std::vector sample_rich_point_set(nb_points_sample); - regularize_and_simplify_internal::Rich_box box; + rich_grid_internel::Rich_box box; for (it = first_original_point, i = 0; it != beyond ; ++it, i++) { Point& p0 = get(point_pmap,it); @@ -820,11 +336,11 @@ regularize_and_simplify_point_set_using_balltree( std::vector original_density_weight_set; if (need_compute_density) { - time.start("Buile Ball Tree For Original"); - regularize_and_simplify_internal::compute_ball_neighbors_one_self(original_rich_point_set, box, neighbor_radius); - time.end(); - - time.start("Compute Density For Original"); + task_timer.start(); + std::cout << "Initialization / Compute Density For Original" << std::endl; + + rich_grid_internel::compute_ball_neighbors_one_self(original_rich_point_set, box, neighbor_radius); + for (it = first_original_point, i = 0; it != beyond ; ++it, i++) { std::vector original_neighbors; @@ -840,26 +356,30 @@ regularize_and_simplify_point_set_using_balltree( original_density_weight_set.push_back(density); original_rich_point_set[i].neighbors.clear(); } - time.end(); + + long memory = CGAL::Memory_sizer().virtual_size(); + std::cout << "done: " << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated" << std::endl << std::endl; + task_timer.stop(); } for (unsigned int iter_n = 0; iter_n < iter_number; iter_n++) { + task_timer.start(); + std::cout << "Compute average term and repulsion term " << std::endl; + // Build Ball Tree For Sample Neighbor - time.start("Build Ball Tree For Sample Neighbor"); for (i=0 ; i < sample_points.size(); i++) { Point& p0 = sample_points[i]; Rich_point rp(p0, i); sample_rich_point_set[i] = rp; } - regularize_and_simplify_internal::compute_ball_neighbors_one_self(sample_rich_point_set, box, neighbor_radius); - time.end(); + rich_grid_internel::compute_ball_neighbors_one_self(sample_rich_point_set, box, neighbor_radius); // Compute sample density weight for sample points if user needed std::vector sample_density_weight_set; - time.start("Compute Density For Sample"); if (need_compute_density) { for (i=0 ; i < sample_points.size(); i++) @@ -876,19 +396,17 @@ regularize_and_simplify_point_set_using_balltree( sample_density_weight_set.push_back(density); } } - time.end(); // Build Ball Tree For Sample-Original Neighbor - time.start("Build Ball Tree For Sample-Original Neighbor"); - regularize_and_simplify_internal::compute_ball_neighbors_one_to_another(sample_rich_point_set, + rich_grid_internel::compute_ball_neighbors_one_to_another(sample_rich_point_set, original_rich_point_set, box, neighbor_radius); - time.end(); // Compute average term and repulsion term for each sample points separately, // then update each sample points std::vector average_set(nb_points_sample); std::vector repulsion_set(nb_points_sample); - time.start("Compute Average Term"); + + // average term for (i = 0; i < sample_points.size(); i++) { Point& p = sample_points[i]; @@ -910,9 +428,8 @@ regularize_and_simplify_point_set_using_balltree( average_set[i] = regularize_and_simplify_internal::compute_average_term(p, rich_original_neighbors, neighbor_radius, original_density_weight_set); } - time.end(); - time.start("Compute Repulsion Term"); + //repulsion term for (i = 0; i < sample_points.size(); i++) { std::vector rich_sample_neighbors; @@ -934,15 +451,20 @@ regularize_and_simplify_point_set_using_balltree( Point& p = sample_points[i]; repulsion_set[i] = regularize_and_simplify_internal::compute_repulsion_term(p, rich_sample_neighbors, neighbor_radius, sample_density_weight_set); } - time.end(); - + + // update points positions according to average and repulsion term for (i = 0; i < sample_points.size(); i++) { Point& p = sample_points[i]; p = CGAL::ORIGIN + average_set[i] + (FT)0.5 * repulsion_set[i]; } - std::cout << "iterate: " << iter_n + 1 << " "<< std::endl; + long memory = CGAL::Memory_sizer().virtual_size(); + std::cout << "done: " << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated" << std::endl; + task_timer.stop(); + + std::cout << "iterate: " << iter_n + 1 << " "<< std::endl << std::endl;; } //Copy back modified sample points to original points for output