Bad results after MST normal orientation and Poisson reconstruction

DamonsJ commented 6 years ago

Please use the following template to help us solving your issue.

Issue Details

I want to estimate normals of point cloud using pca_estimate_normals, but i found it generates very bad result compare with pcl using same parameters(k neighbor )

test file is attached in the attachment:please change the file extension to ply

Source Code

CGAL::pca_estimate_normals(points.begin(), points.end(), CGAL::First_of_pair_property_map(), CGAL::Second_of_pair_property_map(), nb_neighbors_pca_normals);

Environment

Operating system (Windows/Mac/Linux, 32/64 bits):windows 10
Compiler: visual studio 2017
Release or debug mode: release
Specific flags used (if any):
CGAL version:4.11.1
Boost version:1.66
Other libraries versions if used (Eigen, TBB, etc.):eigen tbb qt

sgiraudot commented 6 years ago

Are the normals in your point set the result of the PCL algorithm?
What is the value of nb_neighbors_pca_normals?
Can you also share the source code of the PCL version so we can compare?

DamonsJ commented 6 years ago

the normals in my point set is not the result of pcl or cgal,it is the origin.
nb_neighbors_pca_normals is 18or 6or 10 ,every parameter generates bad result. when i say bad result you need to perform mst_orient_normals function after pca_estimate_normals
i can not find the source code of pcl currently,but you can try the open source softe named CloudCompare :http://www.cloudcompare.org/ there is a function named estimate normals and curvature ,it use pcl's api,you can try.

by the way ,I use cloudcompare follow this flow : 1.estimate normals and curvature

edit->normals->orient normals->with minimum spanning tree
possion surface reconstruction the result is much more good than cgal ,i can not believe,may be it's my fault.

here is the main code i used 👍

//1. k-scale calculating m_nKSize = CGAL::estimate_global_k_neighbor_scale(m_ptSet.begin(), m_ptSet.end(), CGAL::First_of_pair_property_map()); //2. Estimate scale of the point set with average spacing m_dAvgSpacing = CGAL::compute_average_spacing (m_ptSet.begin(), m_ptSet.end(), CGAL::First_of_pair_property_map(), m_nKSize); // I don't know the ratio of outliers present in the point set PointList::iterator first_to_remove = CGAL::remove_outliers(m_ptSet.begin(), m_ptSet.end(), CGAL::First_of_pair_property_map(), m_nKSize, 100., // No limit on the number of outliers to remove

- m_dAvgSpacing); // Point with distance above 2*average_spacing are considered outliers m_ptSet.erase(first_to_remove,m_ptSet.end());
// Estimates normals direction. bool success = run_pca_estimate_normals(m_ptSet, m_nKSize, original_normals); //bool success = run_jet_estimate_normals(m_ptSet, m_nKSize, original_normals); //if (!success) // return false; // set error and continue

// Orients normals. success = run_mst_orient_normals(m_ptSet, m_nKSize, original_normals,_bIsEraseUnOriented); if (!success) return false; // set error and continue

FT sm_angle = 20.0; // Min triangle angle in degrees. FT sm_radius = 30; // Max triangle size w.r.t. point set average spacing. FT sm_distance = 0.1*m_dAvgSpacing; // Surface Approximation error w.r.t. point set average spacing.

// Creates implicit function from the read points using the default solver.
// Note: this method requires an iterator over points
// + property maps to access each point's position and normal.
// The position property map can be omitted here as we use iterators over Point_3 elements.
Poisson_reconstruction_function function(m_ptSet.begin(), m_ptSet.end(),
    CGAL::First_of_pair_property_map<PointVectorPair>(), CGAL::Second_of_pair_property_map<PointVectorPair>());
// Computes the Poisson indicator function f()
// at each vertex of the triangulation.
if (!function.compute_implicit_function())
    return ;
// Computes average spacing
FT average_spacing = m_dAvgSpacing;
// Gets one point inside the implicit surface
// and computes implicit function bounding sphere radius.
Point inner_point = function.get_inner_point();
Sphere bsphere = function.bounding_sphere();
FT radius = std::sqrt(bsphere.squared_radius());
// Defines the implicit surface: requires defining a
// conservative bounding sphere centered at inner point.
FT sm_sphere_radius = 5.0 * radius;
FT sm_dichotomy_error = sm_distance * average_spacing / 1000.0; // Dichotomy error must be << sm_distance
Surface_3 surface(function,
    Sphere(inner_point, sm_sphere_radius*sm_sphere_radius),
    sm_dichotomy_error / sm_sphere_radius);
// Defines surface mesh generation criteria
CGAL::Surface_mesh_default_criteria_3<STr> criteria(sm_angle,  // Min triangle angle (degrees)
    sm_radius*average_spacing,  // Max triangle size
    sm_distance*average_spacing); // Approximation error
                                  // Generates surface mesh with manifold option
STr tr; // 3D Delaunay triangulation for surface mesh generation
C2t3 c2t3(tr); // 2D complex in 3D Delaunay triangulation
CGAL::make_surface_mesh(c2t3,                                 // reconstructed mesh
    surface,                              // implicit surface
    criteria,                             // meshing criteria
    CGAL::Manifold_with_boundary_tag() /*CGAL::Manifold_tag()*/);  // require manifold mesh
if (tr.number_of_vertices() == 0)
    return ;

CGAL::output_surface_facets_to_polyhedron(c2t3, output_mesh);

////////////////////////////////////////////////////////////////////////////// like this flow: m_cgalPro.ReadPointCloudFromFile(filename); m_cgalPro.RemoveOutlierPoints(); m_cgalPro.Simplification(); m_cgalPro.Smoothing(); m_cgalPro.NormalEstimate();

m_cgalPro.PossionSurfaceReconstruction();

sgiraudot commented 6 years ago

@DamonsJ commented on 9 avr. 2018 à 09:42 UTC+2:

the normals in my point set is not the result of pcl or cgal,it is the origin.

nb_neighbors_pca_normals is 18or 6or 10 ,every parameter generates bad result. when i say bad result you need to perform mst_orient_normals function after pca_estimate_normals

i can not find the source code of pcl currently,but you can try the open source softe named CloudCompare :http://www.cloudcompare.org/ there is a function named estimate normals and curvature ,it use pcl's api,you can try.

by the way ,I use cloudcompare follow this flow : 1.estimate normals and curvature

edit->normals->orient normals->with minimum spanning tree

possion surface reconstruction the result is much more good than cgal ,i can not believe,may be it's my fault.

Could you also share the parameters for these functions? I can't get a good reconstruction using these steps on CloudCompare, but I believe this is strongly dependent on the parameters…

here is the main code i used 👍

//1. k-scale calculating m_nKSize = CGAL::estimate_global_k_neighbor_scale(m_ptSet.begin(), m_ptSet.end(), CGAL::First_of_pair_property_map()); (...) CGAL::output_surface_facets_to_polyhedron(c2t3, output_mesh);

////////////////////////////////////////////////////////////////////////////// like this flow: m_cgalPro.ReadPointCloudFromFile(filename); m_cgalPro.RemoveOutlierPoints(); m_cgalPro.Simplification(); m_cgalPro.Smoothing(); m_cgalPro.NormalEstimate();

m_cgalPro.PossionSurfaceReconstruction();

I'm not sure I understand: do you do simplification and smoothing or not? It is not on your code, but it is listed in your "flow" at the end.

Using outlier removal + simplification + smoothing + normal estimation + Poisson, I manage to get a pretty decent mesh, see here: https://cgal.geometryfactory.com/~sgiraudot/keypoints_mesh.ply (I cropped it by hand as Poisson creates an airbag to close the shape.)

DamonsJ commented 6 years ago

sorry for replying late. here is all the step i used in cloudcompare and the result shown in the attchment. all the code i used also attached in the attachment . please read readme.txt first. very thanks for checking this problem.

step1 step2 step3

CGALPointCloudProcess.cpp.txt CGALPointCloudProcess.h.txt readme.txt testmain.cpp.txt

sgiraudot commented 6 years ago

I'm sorry but I can't seem to reproduce your results with CloudCompare. Using normal estimation + normal orientation, both with 18 neighbors using exactly the algorithms you use in your screenshots, I get quite a bad result with CloudCompare, see left image here (the result using the same parameters with PCA + MST in CGAL is given as a comparison on the right): normals

About the results that you get with CGAL, they indeed seem not good. Are you using the library Eigen? If not, that may be a reason why PCA is giving you bad results: without Eigen, CGAL uses an internal fallback mode that doesn't work as well as Eigen. If you're not using it, I would strongly advise trying again using Eigen (it is a header-only library that is quite easy to integrate, see for example: https://github.com/CGAL/cgal/blob/master/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt#L77 ).

(Note also that we usually prefer jet_estimate_normals() in CGAL which produces better results than PCA.)

DamonsJ commented 6 years ago

thanks for giving me hint,would please share your code which can generate result on the right? by the way I used Eigen in CGAL,in the attachment,I have given my code. Thanks very much

sgiraudot commented 6 years ago

@DamonsJ commented on 17 avr. 2018 à 03:42 UTC+2:

thanks for giving me hint,would please share your code which can generate result on the right?

I'm using the following piece of code:

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>

#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>

// Concurrency
#ifdef CGAL_LINKED_WITH_TBB
typedef CGAL::Parallel_tag Concurrency_tag;
#else
typedef CGAL::Sequential_tag Concurrency_tag;
#endif

typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
typedef CGAL::Point_set_3<Point> Point_set;

int main (int argc, char ** argv)
{
  std::ifstream in (argv[1], std::ios_base::binary);
  Point_set points;
  in >> points;

  CGAL::pca_estimate_normals<Concurrency_tag> (points, 18);
  CGAL::mst_orient_normals (points, 18);

  std::ofstream out ("out.ply", std::ios_base::binary);
  CGAL::set_binary_mode (out);
  out << points;

  return EXIT_SUCCESS;
}

I'm using CGAL::Point_set_3 for convenience, but it shouldn't change anything that you are using pairs of points + vectors.

I compiled using the following CMakeLists.txt:

project( normal_estimation )

cmake_minimum_required(VERSION 2.8.11)

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( ${CGAL_USE_FILE} )

find_package( Boost REQUIRED )
if ( NOT Boost_FOUND )
  message(STATUS "This project requires the Boost library, and will not be compiled.")
  return()  
endif()

find_package( TBB QUIET )

find_package(Eigen3 3.1.0)
if (EIGEN3_FOUND)
  include( ${EIGEN3_USE_FILE} )
endif()

include( CGAL_CreateSingleSourceCGALProgram )

create_single_source_cgal_program( "normal_estimation.cpp" CXX_FEATURES cxx_rvalue_references cxx_variadic_templates)    
if(TBB_FOUND)
  CGAL_target_use_TBB("normal_estimation")
endif()

by the way I used Eigen in CGAL,in the attachment,I have given my code.

Eigen does not appear in your code, which is normal as the inclusion of Eigen is done at compilation (in my case, using this CMakeLists.txt).

DamonsJ commented 6 years ago

@sgiraudot thanks very much I have used Eigen in my cgal and I get the same result as yours. Thanks very much! By the way,I use the point set after normal estimate(right result) to reconstruct surface(using possion surface reconstruct) and can not get a good result, how can I set the parameters to get a good result? The code I used is in the last attachment. Thanks here is what I generate using possion surface reconstruction:

CGAL / cgal