Пересечение сплайна с плоскостью в python scipy

Question

У меня есть несколько точек трехмерных данных, и я подгоняю поверхность через них, используя шлицы тонких пластин следующим образом:

import numpy as np
import scipy as sp
import scipy.interpolate

# x, y, z are the 3D point coordinates

spline = sp.interpolate.Rbf(x, y, z, function='thin_plate', smooth=5, episilon=5)
x_grid = np.linspace(0, 512, 1024)
y_grid = np.linspace(0, 512, 1024)
B1, B2 = np.meshgrid(x_grid, y_grid, indexing='xy')
Z = spline(B1, B2)

Это соответствует поверхности по желанию, как показано на прикрепленном изображении.

введите описание изображения здесь

Теперь я хочу узнать, где этот сплайн пересекает заданную плоскость.

Итак, учитывая эту подобранную поверхность, как я могу узнать, что (x, y) указывает, что эта поверхность разрезает плоскость (z = 25), например.

Итак, приведенный выше код подходит:

z = f(x, y)

и теперь, когда f установлен, мне интересно, подходит ли он можно сделать обратный поиск, т.е. я хочу сделать f^{-1}(z)

Answer 1

Трехмерный контурный график прекрасно интерполирует контур на желаемой высоте:

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import scipy as sp
import scipy.interpolate

N = 10
x = np.random.uniform(100, 400, N)
y = np.random.uniform(100, 400, N)
z = np.random.uniform(0, 100, N)
# x, y, z are the 3D point coordinates
spline = sp.interpolate.Rbf(x, y, z, function='thin_plate', smooth=5, episilon=5)
x_grid = np.linspace(0, 512, 1024)
y_grid = np.linspace(0, 512, 1024)
B1, B2 = np.meshgrid(x_grid, y_grid, indexing='xy')
Z = spline(B1, B2)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.contour(B1, B2, Z, levels=[25], offset=25, colors=['red'])
ax.plot_surface(B1, B2, Z, cmap='autumn_r', lw=1, rstride=10, cstride=10, alpha=0.5)

plt.show()

контур

PS: Если вам нужны координаты xy кривой (ей), они сохраняются внутри контура в виде списка списков 2d координат

contour = ax.contour(B1, B2, Z, levels=[25], offset=25, colors=['red'])
for segments in contour.allsegs:
    for segment in segments:
        print("X:", segment[:,0])
        print("Y:", segment[:,1])

Answer 2

Не уверен, что этого достаточно для вашей конечной цели, но можно использовать функцию numpy .isclose :

import numpy as np

z_target = 25
msk = np.isclose(Z, z_target)

x_target = B1[msk]
y_target = B2[msk]

Обратите внимание, что вы можете настроить допуск выровняйте по своему усмотрению в np.isclose.

Тогда вы можете ожидать, что Z_target = spline(x_target, y_target) - это допуск, отличный от z_target.