如何使用mgrid在矩形和圆形之间进行插值
How to use mgrid to interpolate between a rectangle and a circle
我正在尝试创建一个 3D 表面，其中 1/4 矩形用于外部，1/4 圆用于内部。之前我曾帮助过以椭圆作为外部创建 3D 表面，但由于某种原因我无法对矩形执行此操作。我已经手工完成了有意义的数学运算，但我的代码没有。如果能提供任何帮助，我将不胜感激。
import numpy as np
import pyvista as pv
# parameters for the waveguide
# diameter of the inner circle
waveguide_throat = 30
# axes of the outer ellipse
ellipse_x = 250
ellipse_y = 170
# shape parameters for the z profile
depth_factor = 4
angle_factor = 40
# number of grid points in radial and angular direction
array_length = 100
phase_plug = 0
phase_plug_dia = 20
plug_offset = 5
dome_dia = 28
# theta is angle where x and y intersect
theta = np.arctan(ellipse_x / ellipse_y)
# chi is for x direction and lhi is for y direction
chi = np.linspace(0, theta, 100)
lhi = np.linspace(theta, np.pi/2, 100)
# mgrid to create structured grid
r, phi = np.mgrid[0:1:array_length*1j, 0:np.pi/2:array_length*1j]
# Rectangle exterior, circle interior
x = (ellipse_y * np.tan(chi)) * r + ((waveguide_throat / 2 * (1 - r)) * np.cos(phi))
y = (ellipse_x / np.tan(lhi)) * r + ((waveguide_throat / 2 * (1 - r)) * np.sin(phi))
# compute z profile
angle_factor = angle_factor / 10000
z = (ellipse_x / 2 * r / angle_factor) ** (1 / depth_factor)
plotter = pv.Plotter()
waveguide_mesh = pv.StructuredGrid(x, y, z)
plotter.add_mesh(waveguide_mesh)
plotter.show()
您尝试使用的线性插值法是一种通用工具，应该可以使用（有一个小警告）。所以问题首先出在你的矩形边缘上。
这是绘制内部和外部线条的健全性检查：
# debugging: plot interior and exterior
exterior_points = np.array([
ellipse_y * np.tan(chi),
ellipse_x / np.tan(lhi),
np.zeros_like(chi)
]).T
phi_aux = np.linspace(0, np.pi/2, array_length)
interior_points = np.array([
waveguide_throat / 2 * np.cos(phi_aux),
waveguide_throat / 2 * np.sin(phi_aux),
np.zeros_like(phi_aux)
]).T
plotter = pv.Plotter()
plotter.add_mesh(pv.wrap(exterior_points))
plotter.add_mesh(pv.wrap(interior_points))
plotter.show()
左下角是你的内圈，看起来不错。右上角应该是一个矩形，但实际上不是。
要了解为什么您的原始曲面看起来是这样，我们还必须注意一件事（这是我提到的小警告）：您的曲线方向也是相反的。这意味着您将内部曲线的“顶部”（在屏幕截图中）点与外部曲线的“底部”点进行插值。这解释了奇怪的扇形。
所以你需要固定外曲线，并确保两条边的方向相同。请注意，您可以只为两条边创建两个一维数组，然后对它们进行插值。您不必提出插入插值步骤的符号公式。如果你有相同形状的 1d 数组 x_interior, y_interior, x_exterior, y_exterior 那么你可以做 x_exterior * r + x_interior * (1 - r) 并且对 y 做同样的事情。这意味着删除 mgrid 调用，仅使用形状为 (n, 1) 的数组 r，并利用数组广播进行插值。这意味着做 r = np.linspace(0, 1, array_length)[:, None].
所以问题是如何定义矩形。矩形曲线上的点数需要与圆上的点数相同（我强烈建议在任何地方都使用 array_length 参数来确保这一点！）。因为你想跨越整个矩形，我相信你必须选择一个数组索引（即圆弧中的某个角度），它将映射到矩形的角。然后这是一个简单的问题，只改变 y 直到该索引的点，以及 x 其余的（反之亦然）。
这就是我的意思：您知道矩形的角在您的代码中成 theta 角（尽管我认为如果我们假设 x 和 y 混淆了“x”、“y”和角度正切之间的约定关系）。由于 theta 从 0 变为 pi/2，并且您的 phi 值也从 0 变为 pi/2，您应该选择索引 (array_length * (2*theta/np.pi)).round().astype(int) - 1（或类似的东西）将映射到矩形的角。如果你有一个正方形，这会给你 theta = pi/4，因此 (array_length / 2).round().astype(int) - 1。对于 array_length = 3，这是索引 (2 - 1) == 1，它是 3 长度数组的中间索引。 （边缘上的点越多，如果您在这里犯 off-by-one 错误，问题就越小。）
唯一剩下的复杂问题是我们必须显式地将 1d z 数组广播到公共形状。我们可以使用相同的数学方法来获得角度等距的矩形边。
您的代码已根据此建议修正（请注意，我已将 1 添加到角索引，因为我将其用作 right-exclusive 范围索引）：
import numpy as np
import pyvista as pv
# parameters for the waveguide
# diameter of the inner circle
waveguide_throat = 30
# axes of the outer ellipse
ellipse_x = 250
ellipse_y = 170
# shape parameters for the z profile
depth_factor = 4
angle_factor = 40
# number of grid points in radial and angular direction
array_length = 100
# quarter circle interior line
phi = np.linspace(0, np.pi/2, array_length)
x_interior = waveguide_throat / 2 * np.cos(phi)
y_interior = waveguide_throat / 2 * np.sin(phi)
# theta is angle where x and y intersect
theta = np.arctan2(ellipse_y, ellipse_x)
# find array index which maps to the corner of the rectangle
corner_index = (array_length * (2*theta/np.pi)).round().astype(int)
# construct rectangular coordinates manually
x_exterior = np.zeros_like(x_interior)
y_exterior = x_exterior.copy()
phi_aux = np.linspace(0, theta, corner_index)
x_exterior[:corner_index] = ellipse_x
y_exterior[:corner_index] = ellipse_x * np.tan(phi_aux)
phi_aux = np.linspace(np.pi/2, theta, array_length - corner_index, endpoint=False)[::-1] # mind the reverse!
x_exterior[corner_index:] = ellipse_y / np.tan(phi_aux)
y_exterior[corner_index:] = ellipse_y
# interpolate between two curves
r = np.linspace(0, 1, array_length)[:, None] # shape (array_length, 1) for broadcasting
x = x_exterior * r + x_interior * (1 - r)
y = y_exterior * r + y_interior * (1 - r)
# debugging: plot interior and exterior
exterior_points = np.array([
x_exterior,
y_exterior,
np.zeros_like(x_exterior),
]).T
interior_points = np.array([
x_interior,
y_interior,
np.zeros_like(x_interior),
]).T
plotter = pv.Plotter()
plotter.add_mesh(pv.wrap(exterior_points))
plotter.add_mesh(pv.wrap(interior_points))
plotter.show()
# compute z profile
angle_factor = angle_factor / 10000
z = (ellipse_x / 2 * r / angle_factor) ** (1 / depth_factor)
# explicitly broadcast to the shape of x and y
z = np.broadcast_to(z, x.shape)
plotter = pv.Plotter()
waveguide_mesh = pv.StructuredGrid(x, y, z)
plotter.add_mesh(waveguide_mesh, style='wireframe')
plotter.show()
曲线看起来很合理：
与插值表面一样：
我正在尝试创建一个 3D 表面，其中 1/4 矩形用于外部，1/4 圆用于内部。之前我曾帮助过以椭圆作为外部创建 3D 表面，但由于某种原因我无法对矩形执行此操作。我已经手工完成了有意义的数学运算，但我的代码没有。如果能提供任何帮助，我将不胜感激。
import numpy as np
import pyvista as pv
# parameters for the waveguide
# diameter of the inner circle
waveguide_throat = 30
# axes of the outer ellipse
ellipse_x = 250
ellipse_y = 170
# shape parameters for the z profile
depth_factor = 4
angle_factor = 40
# number of grid points in radial and angular direction
array_length = 100
phase_plug = 0
phase_plug_dia = 20
plug_offset = 5
dome_dia = 28
# theta is angle where x and y intersect
theta = np.arctan(ellipse_x / ellipse_y)
# chi is for x direction and lhi is for y direction
chi = np.linspace(0, theta, 100)
lhi = np.linspace(theta, np.pi/2, 100)
# mgrid to create structured grid
r, phi = np.mgrid[0:1:array_length*1j, 0:np.pi/2:array_length*1j]
# Rectangle exterior, circle interior
x = (ellipse_y * np.tan(chi)) * r + ((waveguide_throat / 2 * (1 - r)) * np.cos(phi))
y = (ellipse_x / np.tan(lhi)) * r + ((waveguide_throat / 2 * (1 - r)) * np.sin(phi))
# compute z profile
angle_factor = angle_factor / 10000
z = (ellipse_x / 2 * r / angle_factor) ** (1 / depth_factor)
plotter = pv.Plotter()
waveguide_mesh = pv.StructuredGrid(x, y, z)
plotter.add_mesh(waveguide_mesh)
plotter.show()
您尝试使用的线性插值法是一种通用工具，应该可以使用（有一个小警告）。所以问题首先出在你的矩形边缘上。
这是绘制内部和外部线条的健全性检查：
# debugging: plot interior and exterior
exterior_points = np.array([
ellipse_y * np.tan(chi),
ellipse_x / np.tan(lhi),
np.zeros_like(chi)
]).T
phi_aux = np.linspace(0, np.pi/2, array_length)
interior_points = np.array([
waveguide_throat / 2 * np.cos(phi_aux),
waveguide_throat / 2 * np.sin(phi_aux),
np.zeros_like(phi_aux)
]).T
plotter = pv.Plotter()
plotter.add_mesh(pv.wrap(exterior_points))
plotter.add_mesh(pv.wrap(interior_points))
plotter.show()
左下角是你的内圈，看起来不错。右上角应该是一个矩形，但实际上不是。
要了解为什么您的原始曲面看起来是这样，我们还必须注意一件事（这是我提到的小警告）：您的曲线方向也是相反的。这意味着您将内部曲线的“顶部”（在屏幕截图中）点与外部曲线的“底部”点进行插值。这解释了奇怪的扇形。
所以你需要固定外曲线，并确保两条边的方向相同。请注意，您可以只为两条边创建两个一维数组，然后对它们进行插值。您不必提出插入插值步骤的符号公式。如果你有相同形状的 1d 数组 x_interior, y_interior, x_exterior, y_exterior 那么你可以做 x_exterior * r + x_interior * (1 - r) 并且对 y 做同样的事情。这意味着删除 mgrid 调用，仅使用形状为 (n, 1) 的数组 r，并利用数组广播进行插值。这意味着做 r = np.linspace(0, 1, array_length)[:, None].
所以问题是如何定义矩形。矩形曲线上的点数需要与圆上的点数相同（我强烈建议在任何地方都使用 array_length 参数来确保这一点！）。因为你想跨越整个矩形，我相信你必须选择一个数组索引（即圆弧中的某个角度），它将映射到矩形的角。然后这是一个简单的问题，只改变 y 直到该索引的点，以及 x 其余的（反之亦然）。
这就是我的意思：您知道矩形的角在您的代码中成 theta 角（尽管我认为如果我们假设 x 和 y 混淆了“x”、“y”和角度正切之间的约定关系）。由于 theta 从 0 变为 pi/2，并且您的 phi 值也从 0 变为 pi/2，您应该选择索引 (array_length * (2*theta/np.pi)).round().astype(int) - 1（或类似的东西）将映射到矩形的角。如果你有一个正方形，这会给你 theta = pi/4，因此 (array_length / 2).round().astype(int) - 1。对于 array_length = 3，这是索引 (2 - 1) == 1，它是 3 长度数组的中间索引。 （边缘上的点越多，如果您在这里犯 off-by-one 错误，问题就越小。）
唯一剩下的复杂问题是我们必须显式地将 1d z 数组广播到公共形状。我们可以使用相同的数学方法来获得角度等距的矩形边。
您的代码已根据此建议修正（请注意，我已将 1 添加到角索引，因为我将其用作 right-exclusive 范围索引）：
import numpy as np
import pyvista as pv
# parameters for the waveguide
# diameter of the inner circle
waveguide_throat = 30
# axes of the outer ellipse
ellipse_x = 250
ellipse_y = 170
# shape parameters for the z profile
depth_factor = 4
angle_factor = 40
# number of grid points in radial and angular direction
array_length = 100
# quarter circle interior line
phi = np.linspace(0, np.pi/2, array_length)
x_interior = waveguide_throat / 2 * np.cos(phi)
y_interior = waveguide_throat / 2 * np.sin(phi)
# theta is angle where x and y intersect
theta = np.arctan2(ellipse_y, ellipse_x)
# find array index which maps to the corner of the rectangle
corner_index = (array_length * (2*theta/np.pi)).round().astype(int)
# construct rectangular coordinates manually
x_exterior = np.zeros_like(x_interior)
y_exterior = x_exterior.copy()
phi_aux = np.linspace(0, theta, corner_index)
x_exterior[:corner_index] = ellipse_x
y_exterior[:corner_index] = ellipse_x * np.tan(phi_aux)
phi_aux = np.linspace(np.pi/2, theta, array_length - corner_index, endpoint=False)[::-1] # mind the reverse!
x_exterior[corner_index:] = ellipse_y / np.tan(phi_aux)
y_exterior[corner_index:] = ellipse_y
# interpolate between two curves
r = np.linspace(0, 1, array_length)[:, None] # shape (array_length, 1) for broadcasting
x = x_exterior * r + x_interior * (1 - r)
y = y_exterior * r + y_interior * (1 - r)
# debugging: plot interior and exterior
exterior_points = np.array([
x_exterior,
y_exterior,
np.zeros_like(x_exterior),
]).T
interior_points = np.array([
x_interior,
y_interior,
np.zeros_like(x_interior),
]).T
plotter = pv.Plotter()
plotter.add_mesh(pv.wrap(exterior_points))
plotter.add_mesh(pv.wrap(interior_points))
plotter.show()
# compute z profile
angle_factor = angle_factor / 10000
z = (ellipse_x / 2 * r / angle_factor) ** (1 / depth_factor)
# explicitly broadcast to the shape of x and y
z = np.broadcast_to(z, x.shape)
plotter = pv.Plotter()
waveguide_mesh = pv.StructuredGrid(x, y, z)
plotter.add_mesh(waveguide_mesh, style='wireframe')
plotter.show()
曲线看起来很合理：
与插值表面一样：