Get line intersection between interpolated plane ( CloughTocher2DInterpolator ) and flat plane

Question

I have a number of points from which I want to interpolate a plane which I do as followed:

rng = np.random.default_rng()
x = rng.random(10) - 0.5
y = rng.random(10) - 0.5
z = np.hypot(x, y)
X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
interp = CloughTocher2DInterpolator(list(zip(x, y)), z)

I defined x, y and z as random for this example. x, y and z can also be sparse.

Now I want to be able to get an intersection from ths interpolation with a plane I've defined as x*a + y*b + c = z (with a, b and c variables)

I couldn't find anything to do this in scipy. I've tried getting all the points, getting a triangulation and using this for intersection (see: How to visualize 3D delaunay triangulation in Python?) but the triangulation I get is really bad. It doesn't connect at all to the closest points.

I've also tried using curve_fit but didn't have any luck since the interpolator might return nan values when you request values outside of the convex hull and then the optimizer seems to get stuck:

from scipy import optimize

def func(x, y):
    z = uslope[0] * x + uslope[1] * y + uslope[2]
    dist = interp(x, y) - z
    dist[np.isnan(dist)] = np.Infinity
    return dist

xdata = np.linspace(0, 1000, 50)
ydata = np.zeros(xdata.shape)
ydata.fill(500)


popt, pcov = optimize.curve_fit(func, xdata, ydata, p0=500, check_finite=True, method="dogbox")
print(popt)
print(pcov)

If it's complicated to do with CloughTocher2DInterpolator I can also use another 2D interpolator that is continuous.

score 0 · Answer 1 · answered Apr 25 '23 at 13:08

You can use the plt.contour function to find the xy-coordinates of the polygonal line that approximates the intersection curves of the interpolation surface with the plane.

The following code starts from your example and then finds the zero-contour of the function Z_interp - Z_plane, which gives the xy-polylines of the intersection curve(s) (there can be more than one curve branch in the general case).

np.random.seed(7)  # so you can reconstruct my results :)
x = np.random.random(10) - 0.5
y = np.random.random(10) - 0.5
z = np.hypot(x, y)
X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
interp = CloughTocher2DInterpolator(list(zip(x, y)), z)
Z = interp(X, Y)

plt.pcolormesh(X, Y, Z, shading='auto')
plt.plot(x, y, "ok", label="input point")
plt.axis("equal")

(up to here just your code..)

# plane coefficients
a = 0
b = 0
c = 0.03

Z_plane = a*X + b*Y + c

cc = plt.contour(X, Y, Z-Z_pl, levels=[0]) 
polys = cc.allsegs[0]  # always one level [0] in our case

The figure below shows the result of this code:

When setting a=1; b=1; c=0.01, we get the following figure.

If you want to use the coordinate values, you can access them through the polys variable (which is just a list of lists of coordinate lists). Here is an example of extracting the x, y coordinates of each polyline:

for poly in polys:
    xs, ys = zip(*poly)  
    # do something with the coordinates...

This method is actually similar to the function you tried to implement for your curve fit, but the implementation uses the discrete array for contour extraction.

From here, it depends on what your application needs. If you need the z-coordinates, you can plug the x, y values back into the plane (or the surface) to get them. If you need a smooth curve (rather than a polygonal line), you can fit a spline (or any other smooth fitting method) to the resulting polyline.

Get line intersection between interpolated plane ( CloughTocher2DInterpolator ) and flat plane

1 Answers1