Boolean diff produces incorrect mesh

[GlennWSo]

Hi

Im been testing boolean operations, to see if i can use this lib for my next project.

I set up my test to avoid coplanar faces between the base mesh and the tool mesh, but i still ran into issues :-(

Depending the resolution i choose when generating RandomHills, The boolean opeations will sometimes create mesh with some anomalies. When run my test case: boolean op starts fail on randomhills with 30x30 resoultion

what im trying to diff

The anomalies

At first glance the result seems fine but the mesh is no longer a closed envelope.
And with when looking carefully at the mesh, there is a long edge that should be divided more.

code:

import madcad as cad
from madcad import dvec3, Axis

import pyvista as pv
import numpy as np


def poly2mad(poly: pv.PolyData) -> cad.Mesh:
    """
    Helper to convert pyvista.PolyData to madcad
    """
    assert poly.n_faces > 0
    tri = poly.triangulate()  # avoid mutation and make all faces tri
    faces = tri.faces.reshape((-1, 4))[:, 1:].tolist()
    points = tri.points.tolist()

    mesh = cad.Mesh(points, faces)
    mesh.check()
    assert mesh.isvalid()
    assert mesh.issurface()
    return mesh


def mad2poly(mesh: cad.Mesh) -> pv.PolyData:
    """
    helper to convert madcad.Mesh to pyvista.PolyData
    """
    face_arr = np.array([tuple(v) for v in mesh.faces])
    faces = np.pad(
        face_arr,
        pad_width=((0, 0), (1, 0)),
        constant_values=3,
    ).ravel()

    points = np.array([tuple(v) for v in mesh.points])

    poly = pv.PolyData(points, faces)
    return poly


def randsurf(res=30, seed=1, **kwargs) -> cad.Mesh:
    """
    Creates a surface with random hills
    """
    poly = pv.ParametricRandomHills(
        u_res=res, v_res=res, w_res=res, randomseed=seed, **kwargs
    )
    return poly2mad(poly)


def inspect(mesh: cad.Mesh, **kwargs):
    poly = mad2poly(mesh)
    edges = poly.extract_feature_edges(
        feature_edges=False, non_manifold_edges=False, manifold_edges=False
    )
    pmesh = pv.PolyData(edges.points)
    p = pv.Plotter()
    p.add_mesh(poly)
    p.add_mesh(edges, color="red", line_width=2)
    p.add_mesh(pmesh, color="blue", point_size=9.0)
    p.show(**kwargs)
    print(poly)
    print("poly n open edges", poly.n_open_edges)
    print("poly is manifold", poly.is_manifold)
    print("diff1 is envelope?", mesh.isenvelope())


def test_boolean(res=30, dbg=False):
    """
    Test boolean operations robustness with a operations on a random generated Mesh
    """
    surf1 = randsurf(res=res, seed=1)
    axis1 = Axis(dvec3(0, 10, -1), dvec3(0, 0, 1))
    plane1 = cad.square(axis1, 20)
    ex1 = cad.extrusion(dvec3(0, 0, -10), surf1)

    assert ex1.isvalid()
    assert ex1.issurface()
    assert ex1.isenvelope()
    if dbg:
        cad.show([ex1, plane1])

    diff1 = cad.difference(ex1, plane1)
    if dbg:
        cad.show([diff1])
        inspect(diff1, title=f"Boolean diff with RandomHills resolution {res} * {res}")

    assert diff1.isvalid()
    assert diff1.issurface()
    assert diff1.isenvelope()


if __name__ == "__main__":
    # booleans work on low res
    for n in range(5, 55, 5):

        try:
            test_boolean(res=n)
            print(f"resolution: {n} worked")
        except AssertionError:
            test_boolean(res=n, dbg=True)  # set dbg to inspect what is going on
            print(f"resolution: {n} failed")

[GlennWSo]

update: boolean works fine as long as the base in madcad.difference(base, tool) is not tangent with the global y axis.

testing res: 10, ymin: -1.0, ok! :)
testing res: 10, ymin: -0.5, ok! :)
testing res: 10, ymin: 0.0, ok! :)
testing res: 10, ymin: 0.5, ok! :)
testing res: 10, ymin: 1.0, ok! :)
testing res: 20, ymin: -1.0, ok! :)
testing res: 20, ymin: -0.5, ok! :)
testing res: 20, ymin: 0.0, ok! :)
testing res: 20, ymin: 0.5, ok! :)
testing res: 20, ymin: 1.0, ok! :)
testing res: 30, ymin: -1.0, ok! :)
testing res: 30, ymin: -0.5, ok! :)
testing res: 30, ymin: 0.0, ok! :)
testing res: 30, ymin: 0.5, ok! :)
testing res: 30, ymin: 1.0, ok! :)
testing res: 40, ymin: -1.0, ok! :)
testing res: 40, ymin: -0.5, ok! :)
testing res: 40, ymin: 0.0, ok! :)
testing res: 40, ymin: 0.5, ok! :)
testing res: 40, ymin: 1.0, ok! :)
testing res: 50, ymin: -1.0, ok! :)
testing res: 50, ymin: -0.5, ok! :)
testing res: 50, ymin: 0.0, resolution: 50 failed
base_min: dvec3( -10, -4.48437e-07, -9.85842 )
base_max: dvec3( 10, 20, 7.47731 )

[GlennWSo]

with a small translation :testing res: 10, ymin: -1.0, ok! :)
testing res: 20, ymin: -1.0, ok! :)
testing res: 30, ymin: -1.0, ok! :)
testing res: 40, ymin: -1.0, ok! :)
testing res: 50, ymin: -1.0, ok! :)
testing res: 60, ymin: -1.0, ok! :)
testing res: 70, ymin: -1.0, ok! :)
testing res: 80, ymin: -1.0, ok! :)
testing res: 90, ymin: -1.0, ok! :)
testing res: 100, ymin: -1.0, ok! :)
base_min: dvec3( -10, -1, -9.85842 )
base_max: dvec3( 10, 19, 7.48309 )
pv.poly n open edges 0
pv.poly is manifold True
cad.mesh is envelope? True

[GlennWSo]

code:

import madcad as cad
from madcad import typedlist, dvec3, Axis

import pyvista as pv
import numpy as np

from typing import Tuple


def poly2mad(poly: pv.PolyData) -> cad.Mesh:
    """
    Helper to convert pyvista.PolyData to madcad
    """
    assert poly.n_faces > 0
    tri = poly.triangulate()  # avoid mutation and make all faces tri
    faces = tri.faces.reshape((-1, 4))[:, 1:].tolist()
    points = tri.points.tolist()

    mesh = cad.Mesh(points, faces)
    mesh.check()
    assert mesh.isvalid()
    assert mesh.issurface()
    return mesh


def mad2poly(mesh: cad.Mesh) -> pv.PolyData:
    """
    helper to convert madcad.Mesh to pyvista.PolyData
    """
    face_arr = np.array([tuple(v) for v in mesh.faces])
    faces = np.pad(
        face_arr,
        pad_width=((0, 0), (1, 0)),
        constant_values=3,
    ).ravel()

    points = np.array([tuple(v) for v in mesh.points])

    poly = pv.PolyData(points, faces)
    return poly


def randsurf(res=30, seed=1, ydist=0, **kwargs) -> cad.Mesh:
    """
    Creates a surface with random hills
    """
    poly = pv.ParametricRandomHills(
        u_res=res, v_res=res, w_res=res, randomseed=seed, **kwargs
    )
    poly.translate((0, ydist, 0), inplace=True)
    return poly2mad(poly)


def plot_normals(mesh: cad.Mesh):
    mad2poly(mesh).plot_normals()


def inspect_open_edges(mesh: cad.Mesh, **kwargs):
    poly = mad2poly(mesh)

    edges = poly.extract_feature_edges(
        feature_edges=False, non_manifold_edges=False, manifold_edges=False
    )
    pmesh = pv.PolyData(edges.points)
    print("pv.poly n open edges", poly.n_open_edges)
    print("pv.poly is manifold", poly.is_manifold)
    print("cad.mesh is envelope?", mesh.isenvelope())
    if edges.n_points > 0:
        p = pv.Plotter()
        p.add_mesh(poly)
        p.add_mesh(edges, color="red", line_width=2)
        p.add_mesh(pmesh, color="blue", point_size=9.0)
        p.show(**kwargs)
        print(poly)


def random_block(res=20, seed=1, ydist=0) -> Tuple[cad.Mesh, ...]:
    """A block shape that has been cut with randomhills"""
    surf1 = randsurf(res=res, seed=1, ydist=ydist)
    center = surf1.box().center

    base = cad.extrusion(dvec3(0, 0, -10), surf1)
    axis1 = Axis(dvec3(center[0], center[1], -2), dvec3(0, 0, 1))
    tool = cad.square(axis1, 20)
    # tool = randsurf(res=res, seed=1)
    # for i, p in enumerate(tool.points):
    #     tool.points[i] = dvec3(p[0] * 1.2, p[1] * 1.2, -2.0)

    res = cad.boolean.boolean(base, tool, (False, True))
    return base, tool, res


def random_block1(res=20, seed=1) -> cad.Mesh:

    top = randsurf(res=res, seed=1)

    skirt = cad.extrusion(dvec3(0, 0, -3), top.outlines())

    n = len(skirt.points) // 2
    skirt.points[n:] = typedlist(
        [[p[0], p[1], -2.0] for p in skirt.points[n:]],
        dtype=dvec3,
    )

    bot_boundary = skirt.outlines().islands()[1]
    cad.show([bot_boundary])
    # bot = cad.flatsurface(bot_boundary)

    # cad.show([top, skirt, bot])


def debug_boolean(base, tool, res):
    print("base_min:", base.box().min)
    print("base_max:", base.box().max)
    cad.show([base, tool])
    plot_normals(base)
    plot_normals(tool)

    cad.show([res])
    plot_normals(res)
    inspect_open_edges(res, title=f"Boolean op  {res} * {res}")


def test_boolean(base: cad.Mesh, tool: cad.Mesh, res: cad.Mesh, dbg=False):
    assert base.isvalid()
    assert base.issurface()
    assert base.isenvelope()

    assert tool.isvalid()
    assert tool.issurface()

    assert res.isvalid()
    assert res.issurface()
    assert res.isenvelope()


if __name__ == "__main__":
    # booleans work on low res
    # test_boolean(res=100, dbg=True)  # set dbg to inspect what is going on
    for n in range(10, 110, 10):
        for y in np.linspace(-1, 1, 1):
            parts = random_block(res=n, seed=1, ydist=y)

            try:
                print(f"testing res: {n}, ymin: {y}", end=", ")
                test_boolean(*parts)

                print("ok! :)")
            except AssertionError as e:
                print(f"resolution: {n} failed")
                debug_boolean(*parts)
                raise e

    debug_boolean(*parts)

[jimy-byerley]

I didn't tested your example yet, I have to
I can see on the wireframe that the bottom flat face has none of the subdivisions of the side extruded surface, which means the boolean operation noticed the points were aligned and so could be simplified. What is odd is that the side extruded surface didn't simplified those points every time.
It might be a floating-point precision issue ... does pyvista uses 32 bit floats ?

[GlennWSo]

It might be a floating-point precision issue ... does pyvista uses 32 bit floats ?
Yea i agree, precision.

Pyvista uses 32 bit floats by default but it possible to use others. But when i construct madcad. Mesh i use python std floats

edit:
just to be certain i updated my madcad converter to:
points = typedlist(tri.points.tolist(), dtype=dvec3)

to issue persists.

Ill try to conjure up a simplier senario with same issue

[jimy-byerley]

This is a floating point issue
I can avoid it by adding these lines after your randsurf()

surf1 = randsurf(res=res, seed=1)

for i, p in enumerate(surf1.points):
	if abs(p[1]) < 1e-5:
		p[1] = 0
		surf1.points[i] = p

I then have as a result

resolution: 5 worked
resolution: 10 worked
resolution: 15 worked
resolution: 20 worked
resolution: 25 worked
resolution: 30 worked
resolution: 35 worked
resolution: 40 worked
resolution: 45 worked
resolution: 50 worked

This issue is because the coordinates returned by pyvista appear flat for the bottom surface, and not flat for the side surface. That both do not have the same flatness analysis of the same edges is odd and is indeed a bug. I must check the implementation in madcad.mesh.line_simplification()

[jimy-byerley]

Ah, I finally found the cause of this oversimplification of the bottom face: it comes from the ngon retriangulation after intersection
triangulation_outline() is processing triplets of points in the outline by priority, and in some case where the line is almost flat (but not flat according to the precision at 1e-11 expected from float64) like this one, the low priority brings the function to make very thin triangles along this almost flat edge.
And next steps in the boolean operations are removing flat triangles, maybe a bit too agressively.

This problem could be avoided by chaning the triangulation algorithm ... or reducing the flat triangles filtering
I have to figure out.