Polygonal_surface_reconstruction produces open mesh with some inverted faces

[pentacular]

Issue Details

Polygonal_surface_reconstruction is producing a mesh with some face inversion, leading to a non-closed mesh.

I believe this is unexpected, since it should produce a "compact and watertight surface mesh".

I have provided an example below, which should reconstruct from the planes of a triangulated cube.

Please let me know if you see any error or misunderstanding or cannot reproduce the issue.

Source Code

https://gist.github.com/pentacular/c02a8dea2a74006716a41c618e700d4c

define CGAL_USE_GLPK

include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

include <CGAL/GLPK_mixed_integer_program_traits.h>

include <CGAL/Polygonal_surface_reconstruction.h>

include <CGAL/Surface_mesh.h>

include <CGAL/property_map.h>

int main() { typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick_kernel; typedef Epick_kernel::Point_3 Epick_point; typedef Epick_kernel::Vector_3 Epick_vector; typedef CGAL::Surface_mesh Epick_surface_mesh;

typedef boost::tuple<Epick_point, Epick_vector, int> PNI;
typedef CGAL::Nth_of_tuple_property_map<0, PNI> Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI> Normal_map; 
typedef CGAL::Nth_of_tuple_property_map<2, PNI> Plane_id_map;
typedef std::vector<PNI> Point_vector;

Point_vector points{
  PNI(Epick_point(0.5, 0.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(1.5, -1.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(-1.5, 1.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(1.5, 1.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(-0.5, -0.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(-1.5, -1.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(-1.5, 1.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(1.5, -1.5, -1.5), Epick_vector(0, 0, -1), 0),
  PNI(Epick_point(0.5, 0.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(1.5, 1.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(-1.5, 1.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(1.5, -1.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(-0.5, -0.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(1.5, -1.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(-1.5, 1.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(-1.5, -1.5, 1.5), Epick_vector(0, 0, 1), 1),
  PNI(Epick_point(1.5, 0.5, 0.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, 1.5, -1.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, 1.5, 1.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, -1.5, 1.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, -0.5, -0.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, -1.5, 1.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, -1.5, -1.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(1.5, 1.5, -1.5), Epick_vector(1, 0, 0), 2),
  PNI(Epick_point(-0.5, 1.5, 0.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(-1.5, 1.5, -1.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(-1.5, 1.5, 1.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(1.5, 1.5, 1.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(0.5, 1.5, -0.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(1.5, 1.5, 1.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(1.5, 1.5, -1.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(-1.5, 1.5, -1.5), Epick_vector(0, 1, 0), 3),
  PNI(Epick_point(-1.5, -0.5, 0.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, -1.5, -1.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, -1.5, 1.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, 1.5, 1.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, 0.5, -0.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, 1.5, 1.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, 1.5, -1.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(-1.5, -1.5, -1.5), Epick_vector(-1, 0, 0), 4),
  PNI(Epick_point(0.5, -1.5, 0.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(1.5, -1.5, -1.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(1.5, -1.5, 1.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(-1.5, -1.5, 1.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(-0.5, -1.5, -0.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(-1.5, -1.5, 1.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(-1.5, -1.5, -1.5), Epick_vector(0, -1, 0), 5),
  PNI(Epick_point(1.5, -1.5, -1.5), Epick_vector(0, -1, 0), 5)
};

typedef CGAL::Polygonal_surface_reconstruction<Epick_kernel>
    Polygonal_surface_reconstruction;
typedef CGAL::GLPK_mixed_integer_program_traits<double> MIP_Solver;

Polygonal_surface_reconstruction algo(points, Point_map(), Normal_map(),
                                      Plane_id_map());

Epick_surface_mesh mesh;
if (!algo.reconstruct<MIP_Solver>(mesh)) {
  std::cerr << " Failed: " << algo.error_message() << std::endl;
  return 1;
}
std::cout << "Reconstructed mesh is closed: " << CGAL::is_closed(mesh) << std::endl;
std::cout << "Reconstructed mesh: " << mesh << std::endl;
return 0;

}

Which produces

Reconstructed mesh is closed: 0 Reconstructed mesh: OFF 24 6 0

1.5 1.5 -1.5 -1.5 1.5 -1.5 -1.5 -1.5 -1.5 1.5 -1.5 -1.5 -1.5 -1.5 1.5 1.5 1.5 1.5 -1.5 1.5 1.5 1.5 -1.5 1.5 -1.5 -1.5 1.5 1.5 1.5 -1.5 1.5 1.5 1.5 1.5 -1.5 1.5 1.5 -1.5 -1.5 -1.5 1.5 -1.5 -1.5 -1.5 -1.5 1.5 1.5 1.5 1.5 1.5 -1.5 1.5 -1.5 1.5 1.5 -1.5 -1.5 -1.5 1.5 -1.5 -1.5 1.5 1.5 -1.5 1.5 1.5 -1.5 -1.5 1.5 -1.5 -1.5 -1.5 4 2 3 0 1 4 4 7 5 6 4 11 12 9 10 4 13 21 15 16 4 14 22 17 18 4 8 23 19 20

Visualized with meshlab, we see inverted faces

Environment

• Operating system (Linux, 64 bits):
• Compiler: g++
• Release or debug mode:
• Specific flags used (if any): -lglpk -lgmp
• CGAL version: 5.5-I-900
• Boost version:
• Other libraries versions if used (Eigen, TBB, etc.): glpk-5, gmp-6.2.1

[sloriot]

@LiangliangNan would you have to look at this case or shall @soesau have a look?

[LiangliangNan]

I will take care of it in the next few days) (now got stuck in the airport)

From: Sebastien Loriot @.> Reply-To: CGAL/cgal @.> Date: Monday, August 14, 2023 at 14:35 To: CGAL/cgal @.> Cc: Liangliang Nan @.>, Mention @.***> Subject: Re: [CGAL/cgal] Polygonal_surface_reconstruction produces open mesh with some inverted faces (Issue #7647)

@LiangliangNanhttps://github.com/LiangliangNan would you have to look at this case or shall @soesauhttps://github.com/soesau have a look?

— Reply to this email directly, view it on GitHubhttps://github.com/CGAL/cgal/issues/7647#issuecomment-1676757230, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADWOVCG26V3YTVCEFLY5XPLXVHBKLANCNFSM6AAAAAA3M6Z4TA. You are receiving this because you were mentioned.Message ID: @.***>

[soesau]

@LiangliangNan Thank you, but my pull request should already take care of this issue.

[LiangliangNan]

Oh. Thanks. Please proceed.

I think in the implementation, the algorithm ignores the normal information provided in the input, but it internally fits planes using least squares fitting, and thus the orientation of the face normals are not consistent with the input.

The easiest way to fix it is to compare the normal of a fitted plane with the provided point normals (using average), and invert it if the two are pointing to opposite directions. This should work well if the input normals are reoriented.

Another way to fix it is to apply the reorient-and-stitch function provided by CGAL. This approach should work even if the input normals are not reoriented.