如何在Python中拟合带指定法向量的3D点集曲线？

云原生阿岚

2026-8-14

3D空间带法向量的三点曲面拟合问题

问题背景

我在3D空间中有三个点，每个点都带有对应的法向量：A、B、C为点的位置，N_A、N_B、N_C为各自的法向量：

A = np.array([  348.92065834, -1402.3305998,     32.69313966])
N_A = np.array([-0.86925426,  0.02836434, -0.49355091])

B =  np.array([282.19332067,  82.52027998,  -5.92595371])
N_B = np.array([-0.82339849,  0.43041935,  0.3698028])

C = np.array([247.37475615,  -3.70129865, -22.10494737])
N_C = np.array([-0.83989222, 0.23796899, 0.48780305])

这三个点几乎共面，但距离最近的B和C之间存在轻微方向变化。从B、C到A点，我推测X-Y坐标平面上存在曲率，而X-Z坐标平面上呈抛物线形态。

若曲率更明显，圆柱面可拟合这三个点；结合X、Z坐标的法向量来看，B和C的法向量朝向下方，A的法向量朝向上方，因此整体可能为抛物面。问题在于如何结合法向量进行拟合；若无法实现，如何拟合出带有X、Y方向曲率的曲面？

以下是绘图代码：

import numpy as np
import matplotlib.pyplot as plt
def set_axes_radius(ax, origin, radius):

    ax.set_xlim3d([origin[0] - radius, origin[0] + radius])
    ax.set_ylim3d([origin[1] - radius, origin[1] + radius])
    ax.set_zlim3d([origin[2] - radius, origin[2] + radius])


def set_axes_equal(ax, zoom=1.):
    '''
        Make axes of 3D plot have equal scale so that spheres appear as spheres,
        cubes as cubes, etc..  This is one possible solution to Matplotlib's
        ax.set_aspect("equal") and ax.axis("equal") not working for 3D.
        input:
          ax:   a matplotlib axis, e.g., as output from plt.gca().

    '''

    limits = np.array([
        ax.get_xlim3d(),
        ax.get_ylim3d(),
        ax.get_zlim3d(),
    ])

    origin = np.mean(limits, axis=1)
    radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0])) / zoom
    set_axes_radius(ax, origin, radius)
%matplotlib qt 


# positions and their respective normal vectors

A = np.array([  348.92065834, -1402.3305998,     32.69313966])
N_A = np.array([-0.86925426,  0.02836434, -0.49355091])

B =  np.array([282.19332067,  82.52027998,  -5.92595371])
N_B = np.array([-0.82339849,  0.43041935,  0.3698028])

C = np.array([247.37475615,  -3.70129865, -22.10494737])
N_C = np.array([-0.83989222, 0.23796899, 0.48780305])

# A plane is given by
# a*x + b*y + c*z + d = 0
# where (a, b, c) is the normal.
# If the point (x, y, z) lies on the plane, then solving for d yield:
# d = -(a*x + b*y + c*z)
d_A = -np.sum(N_A * A)
d_B = -np.sum(N_B * B)
d_C = -np.sum(N_C * C)

# Create a meshgrid:
delta = 200

xlim_A = A[0] - delta, A[0] + delta
ylim_A = A[1] - delta, A[1] + delta
xx_A, yy_A = np.meshgrid(np.arange(*xlim_A), np.arange(*ylim_A))

xlim_B = B[0] - delta, B[0] + delta
ylim_B = B[1] - delta, B[1] + delta
xx_B, yy_B = np.meshgrid(np.arange(*xlim_B), np.arange(*ylim_B))

xlim_C = C[0] - delta, C[0] + delta
ylim_C = C[1] - delta, C[1] + delta
xx_C, yy_C = np.meshgrid(np.arange(*xlim_C), np.arange(*ylim_C))

# Solving the equation above for z:
# z = -(a*x + b*y +d) / c
zz_A = -(N_A[0] * xx_A + N_A[1] * yy_A + d_A) / N_A[2]
zz_B = -(N_B[0] * xx_B + N_B[1] * yy_B + d_B) / N_B[2]
zz_C = -(N_C[0] * xx_C + N_C[1] * yy_C + d_C) / N_C[2]


fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(xx_A, yy_A, zz_A, alpha=0.5, color='g')
ax.plot_surface(xx_B, yy_B, zz_B, alpha=0.5, color='cyan')
ax.plot_surface(xx_C, yy_C, zz_C, alpha=0.5, color='crimson')





# Plot point.


x_A, y_A, z_A = A
x_B, y_B, z_B = B
x_C, y_C, z_C = C

Plane_A, = ax.plot(x_A, y_A, z_A, marker='o', markersize=5, color='g')
Plane_A.set_label('A position')
ax.legend()

Plane_B, = ax.plot(x_B, y_B, z_B, marker='o', markersize=5, color='cyan')
Plane_B.set_label('B position')
ax.legend()

Plane_C, = ax.plot(x_C, y_C, z_C, marker='o', markersize=5, color='crimson')
Plane_C.set_label('C position')
ax.legend()

# Plot normal.


dx_A, dy_A, dz_A = delta * N_A
ax.quiver(x_A, y_A, z_A, dx_A, dy_A, dz_A, arrow_length_ratio=0.15, linewidth=3, color='g')

dx_B, dy_B, dz_B = delta * N_B
ax.quiver(x_B, y_B, z_B, dx_B, dy_B, dz_B, arrow_length_ratio=0.15, linewidth=3, color='cyan')

dx_C, dy_C, dz_C = delta * N_C
ax.quiver(x_B, y_C, z_C, dx_C, dy_C, dz_C, arrow_length_ratio=0.15, linewidth=3, color='crimson')

# Enforce equal axis aspects so that the normal also appears to be normal.



ax.set_xlim(xmax=1500,xmin=-1500)
ax.set_ylim(ymax=400, ymin=-400)

zlim = max(A[2], B[2], C[2]) - delta, max(A[2], B[2], C[2]) + delta
ax.set_zlim(*zlim)

ax = plt.gca()

#ax.set_box_aspect([1,1,1])
set_axes_equal(ax)
ax.set_xlabel('X', fontsize=20)
ax.set_ylabel('Y', fontsize=20)
ax.set_zlabel('Z', fontsize=20)

plt.show()

解决方案

一、结合法向量拟合抛物面

针对3个点+对应法向量的场景，优先选择二次抛物面模型，贴合你推测的X-Z抛物线、X-Y曲率特征，简化后的模型形式为：
$$z = ax^2 + by^2 + dx + ey + f$$

约束条件构建

每个点需满足两个核心约束：

点在曲面上：将点坐标代入曲面方程，得到第一个等式；
法向量与曲面梯度垂直：曲面在该点的梯度为$(\frac{\partial z}{\partial x}, \frac{\partial z}{\partial y}, -1)$，需与法向量$N=(n_x,n_y,n_z)$点积为0，即：
$$n_x \cdot (2ax + d) + n_y \cdot (2by + e) - n_z = 0$$

三个点共生成6个方程，模型有5个待求参数，属于超定方程组，用最小二乘法求解即可。

代码实现示例

import numpy as np

# 点和法向量数据
points = np.array([A, B, C])
normals = np.array([N_A, N_B, N_C])

# 构建方程组矩阵M和向量b
M = []
b = []
for p, n in zip(points, normals):
    x, y, z = p
    nx, ny, nz = n
    # 曲面约束方程
    M.append([x**2, y**2, x, y, 1])
    b.append(z)
    # 法向量约束方程
    M.append([2 * nx * x, 2 * ny * y, nx, ny, 0])
    b.append(nz)

M = np.array(M)
b = np.array(b)

# 最小二乘求解参数
params, residuals, _, _ = np.linalg.lstsq(M, b, rcond=None)
a, b_coeff, d, e, f = params

# 定义拟合的曲面函数
def fit_surface(x, y):
    return a * x**2 + b_coeff * y**2 + d * x + e * y + f

二、拟合带双曲率的曲面（不考虑法向量）

如果暂时不需要结合法向量，仅基于三个点拟合带X、Y方向曲率的曲面，依然可以用上述二次曲面模型。此时三个点仅能生成3个方程，参数自由度较高，可加入正则化项避免过拟合，或选择圆柱面模型。

圆柱面拟合思路

圆柱面的一般方程为：
$$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - ( (x - x_0)u + (y - y_0)v + (z - z_0)w )^2 = r^2$$
其中$(u,v,w)$是圆柱轴线的单位方向向量，$(x_0,y_0,z_0)$是轴线上一点，$r$是半径。

可先通过PCA分析点集的主成分确定轴线方向，再代入点坐标用最小二乘法求解半径和轴线位置，降低计算复杂度。

内容的提问来源于stack exchange，提问作者The first man

如何在Python中拟合带指定法向量的3D点集曲线？

问题背景 #

解决方案 #

一、结合法向量拟合抛物面 #