如何用Python通过高斯函数从x、y坐标估算z值（两点位于峰值）

核心思路 #

二维高斯函数的基本形式为：

( G(x,y) = A \cdot e^{{-\frac{(x-x_0)}2 + (y-y_0)^2}{2\sigma2}} )
其中( A )是峰值高度，( (x_0,y_0) )是峰值中心，( \sigma )是标准差（控制高斯分布的扩散程度）。

要让两个指定点都处于峰值，我们可以叠加两个独立的高斯函数，将这两个点分别作为两个高斯的中心。这样两个目标点的z值就是两个高斯峰值的和（若设定每个高斯的峰值为( A )，则这两个点的z值为( 2A )，其他点的z值为两个高斯在该点的取值之和）。

具体实现步骤 #

1. 定义双峰值高斯函数 #

import numpy as np

def dual_gaussian(x, y, peak1, peak2, A=1.0, sigma=1.0):
    # peak1: 第一个峰值点的(x,y)坐标
    # peak2: 第二个峰值点的(x,y)坐标
    x0_1, y0_1 = peak1
    x0_2, y0_2 = peak2
    
    # 计算两个独立高斯函数的取值
    gauss1 = A * np.exp(-((x - x0_1)**2 + (y - y0_1)**2) / (2 * sigma**2))
    gauss2 = A * np.exp(-((x - x0_2)**2 + (y - y0_2)**2) / (2 * sigma**2))
    
    # 返回叠加后的结果
    return gauss1 + gauss2

2. 输入坐标数组并计算z值 #

将你的6个点坐标替换到示例数组中，指定两个峰值点后即可计算所有点的高度：

# 替换成你实际的6个(x,y)坐标
points = np.array([
    [1, 2],   # 第一个峰值点
    [3, 4],   # 第二个峰值点
    [1.5, 2.5],
    [2.5, 3.5],
    [0, 0],
    [4, 4]
])

# 从坐标数组中指定两个峰值点（示例选索引0和1的点）
peak1 = points[0]
peak2 = points[1]

# 批量计算每个点的z值
z_values = [dual_gaussian(x, y, peak1, peak2, A=1.0, sigma=0.8) for x, y in points]

# 输出结果
for (x,y), z in zip(points, z_values):
    print(f"点({x},{y})的高度z值: {z:.4f}")

3. 可选：可视化结果 #

如果需要直观查看高斯分布和点的位置，可使用matplotlib绘制3D图：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 生成网格用于绘制完整高斯曲面
x_grid = np.linspace(min(points[:,0])-1, max(points[:,0])+1, 50)
y_grid = np.linspace(min(points[:,1])-1, max(points[:,1])+1, 50)
X, Y = np.meshgrid(x_grid, y_grid)
Z = dual_gaussian(X, Y, peak1, peak2)

# 绘制3D可视化图
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, alpha=0.7, cmap='viridis')
# 标记输入的6个点
ax.scatter(points[:,0], points[:,1], z_values, color='red', s=50, label='输入点')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z (高度)')
ax.legend()
plt.show()

参数调整说明 #

• A：控制单个高斯的峰值高度，两个峰值点的z值为2*A（若两个峰值点距离较远，另一个高斯在该点的取值趋近于0，峰值点z值近似等于各自的A值）。
• sigma：值越小，高斯分布越集中，离峰值点越近的点z值下降越快；值越大，分布越平缓。
• 若需要两个峰值高度不同，可给两个高斯设置不同的A值（如A1=1.5、A2=1.0），此时两个峰值点的z值会分别接近A1和A2。

内容的提问来源于stack exchange，提问作者nooshin boroumand