From df58e7c4f893b0794e5a2a9602f9acc944aa4813 Mon Sep 17 00:00:00 2001 From: Laurent Saboret Date: Tue, 28 Oct 2008 08:41:41 +0000 Subject: [PATCH] New normals orientation method: radial_normals_orientation_3() does a radial orientation of the normals of a point set. Normals are oriented towards exterior of the point set. This very simple and very fast method is intended to convex objects. --- .../poisson/Poisson.rc | 9 +- .../poisson/PoissonDoc.cpp | 131 +++++++++++++++--- .../poisson/PoissonDoc.h | 8 ++ .../poisson/Poisson_vc80.vcproj | 4 + .../poisson/resource.h | 6 +- .../APSS_reconstruction.cpp | 4 +- .../poisson_reconstruction.cpp | 4 +- .../Implicit_fct_delaunay_triangulation_3.h | 12 +- .../include/CGAL/Poisson_implicit_function.h | 84 +++-------- .../orient_normals_minimum_spanning_tree_3.h | 2 +- .../CGAL/radial_normals_orientation_3.h | 112 +++++++++++++++ .../APSS_reconstruction_test.cpp | 4 +- .../poisson_reconstruction_test.cpp | 4 +- 13 files changed, 279 insertions(+), 105 deletions(-) create mode 100644 Surface_reconstruction_3/include/CGAL/radial_normals_orientation_3.h diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson.rc b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson.rc index cc283c5119f..87cc0e90384 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson.rc +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson.rc @@ -140,16 +140,19 @@ BEGIN MENUITEM SEPARATOR MENUITEM "Orient Normals with Minimum Spanning Tree", ID_ALGORITHMS_ORIENTNORMALSMST MENUITEM "Orient Normals w.r.t. Cameras (Gyroviz)", ID_ALGORITHMS_ORIENTNORMALSCAMERAS + MENUITEM "Radial Normals Orientation", ID_RADIAL_NORMALS_ORIENTATION + MENUITEM "Flip Normals", ID_FLIP_NORMALS END POPUP "&Reconstruction" BEGIN MENUITEM "1-step Poisson reconstruction\tP", ID_ONE_STEP_POISSON_RECONSTRUCTION - MENUITEM "1-step Poisson reconstruction with Normalized Divergence", 32838 + MENUITEM "1-step Poisson reconstruction with Normalized Divergence", ID_RECONSTRUCTION_ONE_STEP_NORMALIZED POPUP "Step-by-step Poisson reconstruction" BEGIN MENUITEM "Create Poisson Delaunay Triangulation\t1", ID_CREATE_POISSON_TRIANGULATION MENUITEM SEPARATOR MENUITEM "Calculate average spacing", ID_STEP_CALCULATEAVERAGESPACING + MENUITEM SEPARATOR MENUITEM "Delaunay refinement\t2", ID_RECONSTRUCTION_DELAUNAYREFINEMENT MENUITEM "Delaunay refinement in shell", ID_ALGORITHMS_REFINEINSHELL MENUITEM SEPARATOR @@ -157,9 +160,9 @@ BEGIN MENUITEM "Extrapolate Normals using gaussian kernels", ID_STEP_EXTRAPOLATENORMALSUSINGGAUSSIANKERNELS MENUITEM SEPARATOR MENUITEM "Solve Poisson Equation\t4", ID_RECONSTRUCTION_POISSON - MENUITEM "Solve Poisson Equation With normalized Divergence", ID_RECONSTRUCTION_POISSON_NORMALIZED - MENUITEM SEPARATOR + MENUITEM "Solve Poisson Equation With Normalized Divergence", ID_RECONSTRUCTION_POISSON_NORMALIZED MENUITEM "Poisson Statistics", ID_ALGORITHMS_POISSONSTATISTICS + MENUITEM SEPARATOR MENUITEM "Surface Meshing\t5", ID_RECONSTRUCTION_POISSON_SURFACE_MESHING MENUITEM "Marching Tet Contouring\t6", ID_ALGORITHMS_MARCHINGTETCONTOURING END diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp index 98fb67d5dda..f15adef3a8d 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.cpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -77,6 +78,7 @@ BEGIN_MESSAGE_MAP(CPoissonDoc, CDocument) ON_COMMAND(ID_RECONSTRUCTION_DELAUNAYREFINEMENT, OnReconstructionDelaunayRefinement) ON_COMMAND(ID_RECONSTRUCTION_POISSON, OnReconstructionPoisson) ON_COMMAND(ID_RECONSTRUCTION_POISSON_NORMALIZED, OnReconstructionPoissonNormalized) + ON_UPDATE_COMMAND_UI(ID_RECONSTRUCTION_POISSON_NORMALIZED, OnUpdateReconstructionPoissonNormalized) ON_COMMAND(ID_ALGORITHMS_REFINEINSHELL, OnAlgorithmsRefineInShell) ON_COMMAND(ID_RECONSTRUCTION_POISSON_SURFACE_MESHING, OnReconstructionPoissonSurfaceMeshing) ON_UPDATE_COMMAND_UI(ID_RECONSTRUCTION_POISSON, OnUpdateReconstructionPoisson) @@ -110,21 +112,28 @@ BEGIN_MESSAGE_MAP(CPoissonDoc, CDocument) ON_COMMAND(ID_ALGORITHMS_OUTLIERS_REMOVAL_WRT_AVG_KNN_SQ_DIST, OnOutliersRemovalWrtAvgKnnSqDist) ON_UPDATE_COMMAND_UI(ID_ALGORITHMS_OUTLIERS_REMOVAL_WRT_AVG_KNN_SQ_DIST, OnUpdateOutliersRemovalWrtAvgKnnSqDist) ON_COMMAND(ID_ANALYSIS_AVERAGE_SPACING, OnAnalysisAverageSpacing) - ON_COMMAND(ID_ONE_STEP_POISSON_RECONSTRUCTION, OnOneStepPoissonReconstruction) - ON_COMMAND(ID_RECONSTRUCTION_ONE_STEP_NORMALIZED, OnOneStepPoissonReconstructionWithNormalizedDivergence) - ON_UPDATE_COMMAND_UI(ID_ONE_STEP_POISSON_RECONSTRUCTION, OnUpdateOneStepPoissonReconstruction) ON_UPDATE_COMMAND_UI(ID_ANALYSIS_AVERAGE_SPACING, OnUpdateAnalysisAverageSpacing) + ON_COMMAND(ID_ONE_STEP_POISSON_RECONSTRUCTION, OnOneStepPoissonReconstruction) + ON_UPDATE_COMMAND_UI(ID_ONE_STEP_POISSON_RECONSTRUCTION, OnUpdateOneStepPoissonReconstruction) + ON_COMMAND(ID_RECONSTRUCTION_ONE_STEP_NORMALIZED, OnOneStepPoissonReconstructionWithNormalizedDivergence) ON_COMMAND(ID_RECONSTRUCTION_APSS_RECONSTRUCTION, OnReconstructionApssReconstruction) ON_UPDATE_COMMAND_UI(ID_RECONSTRUCTION_APSS_RECONSTRUCTION, OnUpdateReconstructionApssReconstruction) ON_COMMAND(ID_MODE_APSS, OnModeAPSS) ON_UPDATE_COMMAND_UI(ID_MODE_APSS, OnUpdateModeAPSS) - ON_COMMAND(ID_RECONSTRUCTION_SAVEAS, &CPoissonDoc::OnReconstructionSaveas) + ON_COMMAND(ID_RECONSTRUCTION_SAVEAS, OnReconstructionSaveas) ON_COMMAND(ID_STEP_CALCULATEAVERAGESPACING, OnCalculateAverageSpacing) ON_COMMAND(ID_STEP_EXTRAPOLATENORMALSUSINGGAUSSIANKERNELS,OnExtrapolateNormalsUsingGaussianKernel) - ON_COMMAND(ID_POINT_CLOUD_SIMPLIFICATION_BY_CLUSTERING, &CPoissonDoc::OnPointCloudSimplificationByClustering) - ON_UPDATE_COMMAND_UI(ID_POINT_CLOUD_SIMPLIFICATION_BY_CLUSTERING, &CPoissonDoc::OnUpdatePointCloudSimplificationByClustering) - ON_COMMAND(ID_POINT_CLOUD_SIMPLIFICATION_RANDOM, &CPoissonDoc::OnPointCloudSimplificationRandom) - ON_UPDATE_COMMAND_UI(ID_POINT_CLOUD_SIMPLIFICATION_RANDOM, &CPoissonDoc::OnUpdatePointCloudSimplificationRandom) + ON_UPDATE_COMMAND_UI(ID_STEP_EXTRAPOLATENORMALSUSINGGAUSSIANKERNELS, OnUpdateStepExtrapolatenormalsusinggaussiankernels) + ON_COMMAND(ID_POINT_CLOUD_SIMPLIFICATION_BY_CLUSTERING, OnPointCloudSimplificationByClustering) + ON_UPDATE_COMMAND_UI(ID_POINT_CLOUD_SIMPLIFICATION_BY_CLUSTERING, OnUpdatePointCloudSimplificationByClustering) + ON_COMMAND(ID_POINT_CLOUD_SIMPLIFICATION_RANDOM, OnPointCloudSimplificationRandom) + ON_UPDATE_COMMAND_UI(ID_POINT_CLOUD_SIMPLIFICATION_RANDOM, OnUpdatePointCloudSimplificationRandom) + ON_COMMAND(ID_RADIAL_NORMALS_ORIENTATION, OnRadialNormalsOrientation) + ON_UPDATE_COMMAND_UI(ID_RADIAL_NORMALS_ORIENTATION, OnUpdateRadialNormalsOrientation) + ON_COMMAND(ID_FLIP_NORMALS, OnFlipNormals) + ON_UPDATE_COMMAND_UI(ID_FLIP_NORMALS, OnUpdateFlipNormals) + ON_UPDATE_COMMAND_UI(ID_RECONSTRUCTION_ONE_STEP_NORMALIZED, &CPoissonDoc::OnUpdateReconstructionOneStepNormalized) + ON_UPDATE_COMMAND_UI(ID_STEP_CALCULATEAVERAGESPACING, &CPoissonDoc::OnUpdateStepCalculateaveragespacing) END_MESSAGE_MAP() @@ -374,6 +383,10 @@ void CPoissonDoc::OnFileSaveAs() path = path.Left(path.ReverseFind('\\')); SetCurrentDirectory(path); + // Save normals? + assert(m_points.begin() != m_points.end()); + bool points_have_normals = (m_points.begin()->normal() != CGAL::NULL_VECTOR); + // if .pwn extension if(extension.CompareNoCase(".pwn") == 0) { @@ -388,7 +401,8 @@ void CPoissonDoc::OnFileSaveAs() else if(extension.CompareNoCase(".xyz") == 0) { if( ! CGAL::surface_reconstruction_write_xyz(dlgExport.m_ofn.lpstrFile, - m_points.begin(), m_points.end()) ) + m_points.begin(), m_points.end(), + points_have_normals) ) { prompt_message("Unable to save file"); return; @@ -867,6 +881,12 @@ void CPoissonDoc::OnReconstructionPoisson() EndWaitCursor(); } +// Enable "Solve Poisson equation" if Delaunay refinement is applied +void CPoissonDoc::OnUpdateReconstructionPoisson(CCmdUI *pCmdUI) +{ + pCmdUI->Enable(m_edit_mode == POISSON && m_triangulation_refined); +} + void CPoissonDoc::OnReconstructionPoissonNormalized() { assert(m_poisson_dt != NULL); @@ -897,8 +917,7 @@ void CPoissonDoc::OnReconstructionPoissonNormalized() EndWaitCursor(); } -// Enable "Solve Poisson equation" if Delaunay refinement is applied -void CPoissonDoc::OnUpdateReconstructionPoisson(CCmdUI *pCmdUI) +void CPoissonDoc::OnUpdateReconstructionPoissonNormalized(CCmdUI *pCmdUI) { pCmdUI->Enable(m_edit_mode == POISSON && m_triangulation_refined); } @@ -954,10 +973,10 @@ void CPoissonDoc::OnReconstructionPoissonSurfaceMeshing() << " angle="< triangles; @@ -1148,8 +1167,10 @@ void CPoissonDoc::OnOneStepPoissonReconstructionWithNormalizedDivergence() //OnCalculateAverageSpacing(); OnReconstructionDelaunayRefinement(); //OnExtrapolateNormalsUsingGaussianKernel(); // TO TEST - OnReconstructionPoissonNormalized(); - OnReconstructionPoissonSurfaceMeshing(); + if (m_triangulation_refined) + OnReconstructionPoissonNormalized(); + if (m_poisson_solved) + OnReconstructionPoissonSurfaceMeshing(); status_message("1-step Poisson reconstruction with normalized divergence...done (%.2lf s)", task_timer.time()); update_status(); @@ -1157,6 +1178,11 @@ void CPoissonDoc::OnOneStepPoissonReconstructionWithNormalizedDivergence() EndWaitCursor(); } +void CPoissonDoc::OnUpdateReconstructionOneStepNormalized(CCmdUI *pCmdUI) +{ + OnUpdateCreatePoissonTriangulation(pCmdUI); +} + // Remove vertices / cameras cone's angle is low void CPoissonDoc::OnAlgorithmsOutliersRemovalWrtCamerasConeAngle() { @@ -1165,6 +1191,7 @@ void CPoissonDoc::OnAlgorithmsOutliersRemovalWrtCamerasConeAngle() CGAL::Timer task_timer; task_timer.start(); // Select points to delete + m_points.select(m_points.begin(), m_points.end(), false); Point_set::iterator first_iterator_to_remove = remove_outliers_wrt_camera_cone_angle_3( m_points.begin(), m_points.end(), @@ -1194,6 +1221,7 @@ void CPoissonDoc::OnOutliersRemovalWrtAvgKnnSqDist() CGAL::Timer task_timer; task_timer.start(); // Select points to delete + m_points.select(m_points.begin(), m_points.end(), false); Point_set::iterator first_iterator_to_remove = CGAL::remove_outliers_wrt_avg_knn_sq_distance_3( m_points.begin(), m_points.end(), @@ -1284,10 +1312,10 @@ void CPoissonDoc::OnReconstructionApssReconstruction() << " angle="< triangles; @@ -1342,12 +1370,18 @@ void CPoissonDoc::OnCalculateAverageSpacing() Point_set output; m_poisson_function->average_spacing_avg_knn_sq_distance_3(); + status_message("Average spacing calculated...took %f seconds", task_timer.time()); update_status(); UpdateAllViews(NULL); EndWaitCursor(); } +void CPoissonDoc::OnUpdateStepCalculateaveragespacing(CCmdUI *pCmdUI) +{ + pCmdUI->Enable(m_edit_mode == POISSON); +} + void CPoissonDoc::OnExtrapolateNormalsUsingGaussianKernel() { BeginWaitCursor(); @@ -1355,13 +1389,18 @@ void CPoissonDoc::OnExtrapolateNormalsUsingGaussianKernel() CGAL::Timer task_timer; task_timer.start(); int a = m_poisson_function->extrapolate_normals_using_gaussian_kernel(); - status_message("Extrapolation of normals done... %d normals inserted",a); + status_message("Extrapolation of normals done... %d normals inserted",a); update_status(); UpdateAllViews(NULL); EndWaitCursor(); } +void CPoissonDoc::OnUpdateStepExtrapolatenormalsusinggaussiankernels(CCmdUI *pCmdUI) +{ + pCmdUI->Enable(m_edit_mode == POISSON && m_triangulation_refined); +} + void CPoissonDoc::OnEditDelete() { BeginWaitCursor(); @@ -1407,6 +1446,7 @@ void CPoissonDoc::OnPointCloudSimplificationByClustering() FT size = sqrt(bounding_sphere.squared_radius()); // Select points to delete + m_points.select(m_points.begin(), m_points.end(), false); Point_set::iterator first_iterator_to_remove = CGAL::merge_epsilon_nearest_points_3( m_points.begin(), m_points.end(), @@ -1431,6 +1471,7 @@ void CPoissonDoc::OnPointCloudSimplificationRandom() CGAL::Timer task_timer; task_timer.start(); // Select points to delete + m_points.select(m_points.begin(), m_points.end(), false); Point_set::iterator first_iterator_to_remove = CGAL::random_simplification_points_3( m_points.begin(), m_points.end(), @@ -1447,3 +1488,57 @@ void CPoissonDoc::OnUpdatePointCloudSimplificationRandom(CCmdUI *pCmdUI) { pCmdUI->Enable(m_edit_mode == POINT_SET); } + +void CPoissonDoc::OnRadialNormalsOrientation() +{ + BeginWaitCursor(); + status_message("Radial Normals Orientation..."); + CGAL::Timer task_timer; task_timer.start(); + + CGAL::radial_normals_orientation_3(m_points.begin(), m_points.end(), + get(boost::vertex_index, m_points), + get(CGAL::vertex_point, m_points), + get(boost::vertex_normal, m_points)); + + // Select non-oriented normals + m_points.select(m_points.begin(), m_points.end(), false); + for (int i=0; inormal() != CGAL::NULL_VECTOR); + pCmdUI->Enable(m_edit_mode == POINT_SET && points_have_normals); +} + +void CPoissonDoc::OnFlipNormals() +{ + BeginWaitCursor(); + status_message("Flip Normals..."); + CGAL::Timer task_timer; task_timer.start(); + + // Flip normals + for (int i=0; inormal() != CGAL::NULL_VECTOR); + pCmdUI->Enable(m_edit_mode == POINT_SET && points_have_normals); +} diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.h b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.h index fba2eb7602c..ec860d8619f 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.h +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/PoissonDoc.h @@ -247,6 +247,14 @@ public: afx_msg void OnPointCloudSimplificationRandom(); afx_msg void OnUpdatePointCloudSimplificationByClustering(CCmdUI *pCmdUI); afx_msg void OnUpdatePointCloudSimplificationRandom(CCmdUI *pCmdUI); + afx_msg void OnRadialNormalsOrientation(); + afx_msg void OnUpdateRadialNormalsOrientation(CCmdUI *pCmdUI); + afx_msg void OnFlipNormals(); + afx_msg void OnUpdateFlipNormals(CCmdUI *pCmdUI); + afx_msg void OnUpdateStepExtrapolatenormalsusinggaussiankernels(CCmdUI *pCmdUI); + afx_msg void OnUpdateReconstructionPoissonNormalized(CCmdUI *pCmdUI); + afx_msg void OnUpdateReconstructionOneStepNormalized(CCmdUI *pCmdUI); + afx_msg void OnUpdateStepCalculateaveragespacing(CCmdUI *pCmdUI); }; diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj index 10ff437f6a1..814377e2c94 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/Poisson_vc80.vcproj @@ -444,6 +444,10 @@ RelativePath="..\..\..\include\CGAL\poisson_refinement_3.h" > + + diff --git a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/resource.h b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/resource.h index 9ffab47a8bd..d839861596b 100644 --- a/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/resource.h +++ b/Surface_reconstruction_3/demo/Surface_reconstruction_3/poisson/resource.h @@ -73,13 +73,17 @@ #define ID_POINT_CLOUD_SIMPLIFICATION_RANDOM 32845 #define ID_EDIT_DELETE 32867 #define ID_EDIT_RESET_SELECTION 32869 +#define ID_PROCESSING_RADIALNORMALSORIENTATION 32870 +#define ID_RADIAL_NORMALS_ORIENTATION 32871 +#define ID_PROCESSING_FLIPNORMALS 32872 +#define ID_FLIP_NORMALS 32873 // Next default values for new objects // #ifdef APSTUDIO_INVOKED #ifndef APSTUDIO_READONLY_SYMBOLS #define _APS_NEXT_RESOURCE_VALUE 134 -#define _APS_NEXT_COMMAND_VALUE 32870 +#define _APS_NEXT_COMMAND_VALUE 32874 #define _APS_NEXT_CONTROL_VALUE 1015 #define _APS_NEXT_SYMED_VALUE 101 #endif diff --git a/Surface_reconstruction_3/examples/Surface_reconstruction_3/APSS_reconstruction.cpp b/Surface_reconstruction_3/examples/Surface_reconstruction_3/APSS_reconstruction.cpp index 606929342da..dfa49cd68b7 100644 --- a/Surface_reconstruction_3/examples/Surface_reconstruction_3/APSS_reconstruction.cpp +++ b/Surface_reconstruction_3/examples/Surface_reconstruction_3/APSS_reconstruction.cpp @@ -272,10 +272,10 @@ int main(int argc, char * argv[]) << " angle="< vertices; // use to walk in 3D Delaunay - // vertices.push(v); std::vector v_neighbors; m_dt.incident_vertices(v,std::back_inserter(v_neighbors)); typename std::vector::iterator it; - for(it = v_neighbors.begin(); - it != v_neighbors.end(); - it++) + for(it = v_neighbors.begin(); it != v_neighbors.end(); it++) { sq_distance = sq_distance + distance(*it,v)*distance(*it,v); counter++; } - - /*while(!vertices.empty() && counter < 0.001 * number_vertices ) - { - Vertex_handle v_cur = vertices.top(); - vertices.pop(); - // get incident_vertices - std::vector v_neighbors; - sq_distance = sq_distance + distance(v_cur,v)* distance(v_cur,v); - counter++; - m_dt.incident_vertices(v_cur,std::back_inserter(v_neighbors)); - std::vector::iterator it; - for(it = v_neighbors.begin(); - it != v_neighbors.end(); - it++) - { - Vertex_handle nv = *it; - int tag = nv->tag(); - int index = v_cur->index(); - if (tag != index) - { - vertices.push(nv); - nv->tag() = index; - } - } - */ - - v->average_spacing() = std::sqrt(sq_distance/counter); - - //} - /*for(v = m_dt.finite_vertices_begin(); - v != m_dt.finite_vertices_end(); - v++) - v->tag() = -1;*/ + + v->average_spacing() = std::sqrt(sq_distance/counter); } } @@ -527,9 +484,7 @@ public: // Compute extrapolated normals and store them in extrapolated_normals[] std::map extrapolated_normals; // vector + orientation Finite_vertices_iterator v; - for(v = m_dt.finite_vertices_begin(); - v != m_dt.finite_vertices_end(); - v++) + for(v = m_dt.finite_vertices_begin(); v != m_dt.finite_vertices_end(); v++) { if(v->normal() != CGAL::NULL_VECTOR) continue; @@ -554,9 +509,7 @@ public: } // set normals - for(v = m_dt.finite_vertices_begin(); - v != m_dt.finite_vertices_end(); - v++) + for(v = m_dt.finite_vertices_begin(); v != m_dt.finite_vertices_end(); v++) { if(v->normal() != CGAL::NULL_VECTOR) continue; @@ -567,21 +520,19 @@ public: } } - FT gaussian_function( FT sigma , FT distance) - { - FT answer = (1 / std::sqrt(2 * 3.14)) * std::exp(-1 * distance * distance /(2 * sigma * sigma)); - return answer; - } + FT gaussian_function( FT sigma , FT distance) + { + FT answer = (1 / std::sqrt(2 * 3.14)) * std::exp(-1 * distance * distance /(2 * sigma * sigma)); + return answer; + } - - /// Extrapolate the normals field + /// Extrapolate the normals field. + /// Return the number of normals computed. int extrapolate_normals_using_gaussian_kernel() { int counter = 0; Finite_vertices_iterator v; - for(v = m_dt.finite_vertices_begin(); - v != m_dt.finite_vertices_end(); - v++) + for(v = m_dt.finite_vertices_begin(); v != m_dt.finite_vertices_end(); v++) { if(v->type() == Triangulation::INPUT) { @@ -618,12 +569,9 @@ public: } } } - } - for(v = m_dt.finite_vertices_begin(); - v != m_dt.finite_vertices_end(); - v++) + for(v = m_dt.finite_vertices_begin(); v != m_dt.finite_vertices_end(); v++) { if(v->type() != Triangulation::INPUT ) { @@ -633,9 +581,9 @@ public: v->normal() = v->normal() / sq_norm; counter++; } - //v->type() = Triangulation::INPUT; } } + return counter; } diff --git a/Surface_reconstruction_3/include/CGAL/orient_normals_minimum_spanning_tree_3.h b/Surface_reconstruction_3/include/CGAL/orient_normals_minimum_spanning_tree_3.h index e7bbb47a965..0885948311b 100644 --- a/Surface_reconstruction_3/include/CGAL/orient_normals_minimum_spanning_tree_3.h +++ b/Surface_reconstruction_3/include/CGAL/orient_normals_minimum_spanning_tree_3.h @@ -187,7 +187,7 @@ struct propagate_normal // Is orientation robust? //normal2 = Normal(vec2, true /* oriented */); bool oriented = normal1.is_oriented() && - (std::abs(dot) > 1.0/std::sqrt(2.0)); // oriented iff angle < 45� + (std::abs(dot) > 1.0/std::sqrt(2.0)); // oriented iff angle < 45 degrees normal2 = Normal(vec2, oriented); } }; diff --git a/Surface_reconstruction_3/include/CGAL/radial_normals_orientation_3.h b/Surface_reconstruction_3/include/CGAL/radial_normals_orientation_3.h new file mode 100644 index 00000000000..329d77141b0 --- /dev/null +++ b/Surface_reconstruction_3/include/CGAL/radial_normals_orientation_3.h @@ -0,0 +1,112 @@ +// Copyright (c) 2007-08 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// 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) : Laurent Saboret + +#ifndef CGAL_RADIAL_NORMALS_ORIENTATION_3_H +#define CGAL_RADIAL_NORMALS_ORIENTATION_3_H + +#include +#include +#include + +#include +#ifndef M_PI + #define M_PI 3.14159265358979323846 +#endif + +CGAL_BEGIN_NAMESPACE + + +/// Radial orientation of the normals of a point set. +/// Normals are oriented towards exterior of the point set. +/// This very fast method is intended to convex objects. +/// +/// Preconditions: +/// - VertexIterator is a model of ForwardIterator. +/// - VertexIndexMap is a model of boost::readable_property_map. +/// - VertexPointMap is a model of boost::readable_property_map. +/// - VertexNormalMap is a model of boost::lvalue_property_map. +/// - Normals must be unit vectors. + +template +void +radial_normals_orientation_3(VertexIterator first, ///< first input vertex + VertexIterator beyond, ///< past-the-end input vertex + VertexIndexMap vertex_index_map, ///< property map VertexIterator -> index + VertexPointMap vertex_point_map, ///< property map VertexIterator -> Point_3 + VertexNormalMap vertex_normal_map) ///< property map VertexIterator -> Normal (in and out) +{ + CGAL_TRACE("Call radial_normals_orientation_3()\n"); + + // Input mesh's types + typedef typename boost::property_traits::value_type Point; + typedef typename boost::property_traits::value_type Normal; + typedef typename Normal::Vector Vector; + typedef typename Normal::Geom_traits Geom_traits; // Kernel + typedef typename Geom_traits::FT FT; + + // Precondition: at least one element in the container. + CGAL_surface_reconstruction_precondition(first != beyond); + + // Find points barycenter. + // Note: We should use CGAL::centroid() from PCA component. + // Unfortunately, it is not compatible with property maps. + Vector sum = CGAL::NULL_VECTOR; + int nb_points = 0; + for (VertexIterator it = first; it != beyond; it++) + { + Point point = get(vertex_point_map, it); + sum = sum + (point - CGAL::ORIGIN); + nb_points++; + } + Point barycenter = CGAL::ORIGIN + sum / (FT)nb_points; + + // Radial orientation of the normals of a point set. + // Normals are oriented towards exterior of the point set. + for (VertexIterator it = first; it != beyond; it++) + { + Point point = get(vertex_point_map, it); + + // Radial vector towards exterior of the point set + Vector vec1 = point - barycenter; + + // Point's normal + Normal& normal2 = vertex_normal_map[it]; + Vector vec2 = normal2; + + double dot = vec1 * vec2; + + // -> -> + // Orient vec2 parallel to vec1 + if (dot < 0) + vec2 = -vec2; + + // Is orientation robust? + //normal2 = Normal(vec2, true /* oriented */); + bool oriented = (std::abs(dot) > std::cos(80.*M_PI/180.)); // oriented iff angle < 80 degrees + normal2 = Normal(vec2, oriented); + } + + long memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); + CGAL_TRACE("End of radial_normals_orientation_3()\n"); +} + + +CGAL_END_NAMESPACE + +#endif // CGAL_RADIAL_NORMALS_ORIENTATION_3_H + diff --git a/Surface_reconstruction_3/test/Surface_reconstruction_3/APSS_reconstruction_test.cpp b/Surface_reconstruction_3/test/Surface_reconstruction_3/APSS_reconstruction_test.cpp index 047dfe96860..59ae1b134f7 100644 --- a/Surface_reconstruction_3/test/Surface_reconstruction_3/APSS_reconstruction_test.cpp +++ b/Surface_reconstruction_3/test/Surface_reconstruction_3/APSS_reconstruction_test.cpp @@ -264,10 +264,10 @@ int main(int argc, char * argv[]) << " angle="<