MeshSurfaceIntersection: use initial geometry as much as possible.

[BrunoLevy]

There is room for optimization in RadialSort: we can probably compute all the necessary quantities directly from the input points instead of the generated points.

[BrunoLevy]

• The orientations of the triangles created by the CDT2d and the orientation of the initial triangle are the same
• The attribute original_triangle_id gives me the original triangle (in mesh_copy) associated with a new triangle
• Now I need to find when the normal is flipped. I thought that normal was flipped only for the triangles coming from duplicating the interfaces, and that I could determine that by comparing a triangle t with facets.nb()/2 but it did not work (and I don't understand why). Still investigating ...

Got it: in epilogue(), mesh_remove_bad_facets_no_check() flips some facets. Need to keep track of facet orientation...

Found out why my previous fix did not work in the first place: I was restoring initial facet orientations before identifying the sequences of identical facets in the sorted vector.

[BrunoLevy]

It will be also probably possible to also optimize ray-tracing queries for determining inclusions.

[BrunoLevy]

Let us remind how RadialSort works. It uses two predicates:

• RadialSort::h_orient(h1,h2) = orient3d(p1,p2,q1,q2) where (p1,p2) is the radial edge and q1,q2 are the two opposite vertices of the two triangles.
• RadialSort::h_refNorient(h2) = sign(dot(normal(h2),normal(h_ref)) where normal(h) denotes the normal to the triangle incident to h.

RadialSort::init(index_t h_ref) stores h_ref and computes normal(h_ref) = U x V for h_refNorient().

The degree of h_refNorient() can be considerably reduced by computing the normals to the input triangles that contain the considered triangles, that involve only input point coordinates.

Attribute<index_t> original_facet_id(
       mesh_.facets.attributes(), "original_facet_id"
);
index_t f =  mesh_.halfedges_.facet(h);
bool flipped = (f >= mesh_.facets.nb()/2);
index_t orig_f = original_facet_id[f];
vec3 p1 =   mesh_copy_.facets.point(orig_f,0);
vec3 p2 =   mesh_copy_.facets.point(orig_f,1);
vec3 p3 =   mesh_copy_.facets.point(orig_f,2);
if(flipped) {
   std::swap(p1,p3);
} 
exact::vec3 N = cross(exact::make_vec(p1,p2), exact::make_vec(p1,p3));

For h_orient() there is no trivial optimization. The way orient3d() is written could be maybe improved a bit.

[BrunoLevy]

OK, now we have an algorithm that is

• much easier to understand
• that uses initial mesh vertices directly in most cases

Two things to do:

• do we still need a predicate cache for hrefNorient() ? (to be measured with Expansions and geogram+) -> not needed anymore
• understand what's going on with example022.scad with the Expansions kernel
  • does not come from radial sort
  • underflows / overflows in expansions may come from:
    • three_planes_intersection() / line_plane_intersection()
    • CDT 2D and in_circle() predicate

[BrunoLevy]

Working on understanding what's going on with example022.scad and expansions kernel:

• some co-planar zones are not triangulated the same
• ... some vertices in these zones do not have exactly the same Z (but I did the test using the double-precision coordinates, not the exact ones). The geogram+-based CDT does not see the intersections (and it is normal, if the Z's are different, there is no intersection).
• when snapping the Z's, the exansion-based CDT manages to triangulate it