Boolean union two meshes only works for ellipsoid

Question

I want to create a pothole (ellipsoid) in a street (plane). After some testing, I realized that the boolean union is probably the one to use (as I want to combine both of the meshes however do not want the street to cover up the pothole).
When I try this using a Sphere it works perfectly fine.
However, I am using a ParameticEllipsoid. There I constantly get the error "vtkMath::Jacobi: Error extracting eigenfunctions" and the wanted effect only works with the boolean difference instead of the union. I suppose that it has to do something with the vtk/py and its geometric object creation.

While taking the boolean difference brings the wanted result it takes a lot of time (which I wouldn't mind too much) and only works rarely (and only if the ellipsoid isn't meshed with other objects).

Is there a way to avoid this? My idea was to extract the points of the meshes (NumPy array) and calculate the union with them however I wasn't able to.

def add_pothole():
    street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=100, j_resolution=100)

    points = street_mesh.points
    for i in range(4):
        pothole_mesh = pv.ParametricEllipsoid(10,5,3)
        # alternatively to get the half ellipsoid: (but boolean difference seemed to only work with the closed one
        # pothole_mesh = pv.ParametricEllipsoid(10,5,3, max_v=pi /2)
        pothole_mesh.translate(choice(points))
    
        # "add" the pothole to the street
        street_mesh = pothole_mesh.boolean_difference(street_mesh.triangulate())
    
    street_mesh.plot()

Alternatively, I thought about defining the ellipsoid-shaped depressions directly, using height as a function of in-plane position. Is it possible do this without using boolean operations?

Answer 1

To address your immediate issue, there are two issues with your current code:

1. Unlike Sphere, ParametricEllipsoid has a lot of stray points due to the way VTK constructs it, so the end result is not technically a closed surface. This is a known issue that I plan to fix in pyvista. For now you can fix this by calling clean(tolerance=1e-5) on your ellipsoid. This should stop those vtkMath errors from popping up.
2. subtracting the street from your pothole will reverse the normals of the plane after each call to boolean_difference, which will give you bumps on alternating sides of the plane, which is probably not what you want. You can fix this by calling street_mesh.flip_normals() at the end of your loop.

With these two changes your code will run, albeit very slowly (I don't really know why, but I also don't know how boolean operations work under the hood):

from random import choice
import numpy as np
import pyvista as pv

def add_pothole():
    street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=100, j_resolution=100)
    pothole_template = pv.ParametricEllipsoid(10, 5, 3).clean(tolerance=1e-5).triangulate(

    points = street_mesh.points
    for i in range(4):
        pothole_mesh = pothole_template.copy()
        pothole_mesh.translate(choice(points))

        # "add" the pothole to the street
        street_mesh = pothole_mesh.boolean_difference(street_mesh.triangulate())
        street_mesh.flip_normals()

    street_mesh.plot()

add_pothole()

---

But to address your underlying issue, you don't have to do this in the first place if your street is a z(x, y) surface to begin with. Considering Perlin noise, you probably start from a flat plane with scalars sampled from the Perlin noise, and then call warp_by_scalar() to implement the noise modulation. Well, you can just add ellipsoid-shaped bumps to the same scalars, and then do the warping in one step at the end (exactly the same way as you would superimpose multiple Perlin noise samples).

For this you have to compute the z(x, y) parametrization of an ellipsoid. This is not trivial, but also not very difficult. It's probably more well-known that the equation of a sphere with radius R around the origin is

x^2 + y^2 + z^2= R^2.

This is better (and canonical) when divided by R^2:

(x/R)^2 + (y/R)^2 + (z/R)^2 = 1.

The way you get an ellipsoid is by linearly transforming each axis, changing their respective radius (so they are no longer all R). Here's the implicit equation for an ellipsoid:

(x/a)^2 + (y/b)^2 + (z/c)^2 = 1.

If you want to have a non-trivial center for the ellipsoid, you have to translate each coordinate:

(x - x0)^2/a^2 + (y - y0)^2/b^2 + (z - z0)^2/c^2 = 1.

This is the general form of an ellipsoid in 3 dimensions, assuming that its axes are oriented with the axes of the Cartesian coordinate system.

Well, this is great, because we can rearrange it to get

(z - z0)^2/c^2 = 1 - (x - x0)^2/a^2 - (y - y0)^2/b^2
(z - z0)^2 = c^2 - (c/a)^2 (x - x0)^2 - (c/b)^2 (y - y0)^2
z = z0 +- sqrt(c^2 - (c/a)^2 (x - x0)^2 - (c/b)^2 (y - y0)^2)

The plus-minus part corresponds to the fact that an ellipsoid is not a z(x, y) function, since for each (x, y) point you have two possible values. But you can construct an ellipsoid from two functions: a top and a bottom surface.

Coming back to your problem at last, you can just choose a random (x0, y0, z0) point, and choose the bottom surface of the above ellipsoid (the one that goes z = z0 - sqrt(...)). The only thing you have to worry about is that for (x, y) points where there's no ellipsoid, you'd get an imaginary square root. So you must first filter out the points that are inside the ellipsoid. The simplest way to do this is compute the square root anyway, and discard NaNs:

import numpy as np
import pyvista as pv

# denser plane for sharper ellipsoid
street_mesh = pv.Plane(i_size=100, j_size=100, i_resolution=1000, j_resolution=1000)
street_mesh['Normals'] *= -1  # make normals point up for warping later
# add Perlin noise for example
freq = [0.03, 0.03, 0.03]
noise = pv.perlin_noise(5, freq, [0, 0, 0])
street_mesh['height'] = [noise.EvaluateFunction(point) for point in street_mesh.points]

# the relevant part: ellipsoid function
x0, y0, z0 = street_mesh.points[street_mesh.points.size//6, :]  # or whatever random point you want
x, y, _ = street_mesh.points.T  # two 1d arrays of coordinates
a, b, c = 10, 5, 3  # semi-axes
ellipsoid_fun = z0 - np.sqrt(c**2 - (c/a)**2 * (x - x0)**2 - (c/b)**2 *(y - y0)**2)  # RuntimeWarning
keep_inds = ~np.isnan(ellipsoid_fun)
street_mesh['height'][keep_inds] += ellipsoid_fun[keep_inds]
street_mesh.set_active_scalars('height')

# warp by 'height' and do a quick plot
street_mesh.warp_by_scalar().plot()

The above will emit a warning about the imaginary square roots (leading to NaNs in the data). If this bothers you you can explicitly silence it. Alternatively, we could follow LBYL and check the values ourselves before computing the square root:

arg = c**2 - (c/a)**2 * (x - x0)**2 - (c/b)**2 *(y - y0)**2
keep_inds = arg >= 0
ellipsoid = z0 - np.sqrt(arg[keep_inds])
street_mesh['height'][keep_inds] = ellipsoid

Here's the result, using the increased plane density for a sharper ellipsoid:

(I originally suggested that you could compute every pothole in one sitting using numpy broadcasting, but that's not actually easy, or maybe even possible, due to the NaN filtering which breaks this.)