二维高斯拟合不收敛排查：密集区域光斑强度拟合求助

我瞅了你的代码和问题，发现几个关键的坑导致拟合失败，咱们一个个捋：

1. 局部区域拟合时的坐标系统错误 #

这是最致命的问题！你传入fit_gaussian_spots的x_guess和y_guess是原图的全局坐标，但你拟合的是裁剪后的局部子图。比如你从原图中裁剪了pre_y - local_area到pre_y + local_area的区域，那子图里的x/y坐标范围是0到2*local_area，光斑在子图里的位置应该是相对于子图左上角的偏移，而不是原图的全局坐标。用全局坐标当初始猜测，相当于让算法在子图里找一个不存在的位置，肯定收敛不了！

修正方法：计算局部区域内的相对坐标：

# 在主循环调用拟合函数前，转换坐标
y_start = max(pre_y - local_area, 0)
x_start = max(pre_x - local_area, 0)
local_x = pre_x - x_start
local_y = pre_y - y_start

2. 初始猜测参数太粗糙，和真实值差距过大 #

你当前的初始猜测initial_guess里，amp=1、sigma_x=1、sigma_y=1，但模拟数据里amp是100，sigma是3/9，差距非常大。curve_fit的Levenberg-Marquardt算法对初始值很敏感，初始值离最优解太远的话，很容易陷入局部最小值或者直接不收敛。

修正方法：基于局部区域的统计值生成合理初始值：

local_max = np.max(array)
local_min = np.min(array)
initial_guess = Params(
    amp=local_max - local_min,  # 信号幅度=峰值-背景
    x=x_guess, 
    y=y_guess, 
    sigma_x=array.shape[1]/4,  # 基于局部宽度的初始sigma
    sigma_y=array.shape[0]/4,
    rotation=0, 
    offset=local_min  # 背景初始值设为局部最小值
)

3. 参数边界设置不合理 #

你当前的边界有几个明显问题：

• x/y的边界用x_guess*0.5和x_guess*1.5，如果x_guess是全局坐标，完全超出局部区域范围；
• rotation设为-inf到2π，旋转角是周期的，[-π, π]足够，没必要设无穷大；
• offset上限设为全局最大值，背景应该低于局部区域的平均值，而非最大值。

修正后的边界：

min_bounds = Params(
    amp=1e-8, 
    x=0,  # x不能小于局部区域左边界
    y=0,  # y不能小于局部区域上边界
    sigma_x=1e-8, 
    sigma_y=1e-8, 
    rotation=-np.pi,
    offset=local_min * 0.5
)
max_bounds = Params(
    amp=(local_max - local_min)*1.2,
    x=array.shape[1]-1,  # x不能超过局部区域宽度-1
    y=array.shape[0]-1,  # y不能超过局部区域高度-1
    sigma_x=array.shape[1]/2,
    sigma_y=array.shape[0]/2,
    rotation=np.pi,
    offset=local_min + (local_max - local_min)*0.2
)

4. 拟合后直接round参数丢失精度 #

你在拟合后做了popt = Params(*np.round(popt))，这会把所有参数取整，比如sigma本来是3.2，直接变成3，导致后续评估的高斯和真实光斑偏差极大，看起来像拟合失败。必须去掉这个round操作，保留浮点精度。

修改后的完整代码 #

把所有修正整合后的代码如下：

import scipy.optimize as opt
import numpy as np
import matplotlib.pyplot as plt
import skimage.feature
from collections import namedtuple

def gaussian_2d(xy_array, amplitude, pos_x, pos_y, sigma_x, sigma_y, rotation, offset):
    """ Expression for a 2D gaussian function with variance in both x and y """
    x, y = xy_array
    cos_rot = np.cos(rotation)
    sin_rot = np.sin(rotation)
    a = (cos_rot ** 2) / (2 * sigma_x ** 2) + (sin_rot ** 2) / (2 * sigma_y ** 2)
    b = -(np.sin(2 * rotation)) / (4 * sigma_x ** 2) + (np.sin(2 * rotation)) / (4 * sigma_y ** 2)
    c = (sin_rot ** 2) / (2 * sigma_x ** 2) + (cos_rot ** 2) / (2 * sigma_y ** 2)
    g = amplitude * np.exp(-(a * ((x - pos_x) ** 2) + 2 * b * (x - pos_x) * (y - pos_y) + c * ((y - pos_y) ** 2)))
    g += offset
    return g.ravel()

def fit_gaussian_spots(x_guess, y_guess, array):
    Params = namedtuple("Parameters", "amp, x, y, sigma_x, sigma_y, rotation, offset")
    eps = 1e-8
    local_max = np.max(array)
    local_min = np.min(array)
    
    # 更合理的初始猜测
    initial_guess = Params(
        amp=local_max - local_min,
        x=x_guess,
        y=y_guess,
        sigma_x=array.shape[1]/4,
        sigma_y=array.shape[0]/4,
        rotation=0,
        offset=local_min
    )
    
    # 修正后的边界
    min_bounds = Params(
        amp=eps,
        x=0,
        y=0,
        sigma_x=eps,
        sigma_y=eps,
        rotation=-np.pi,
        offset=local_min * 0.5
    )
    max_bounds = Params(
        amp=(local_max - local_min) * 1.2,
        x=array.shape[1]-1,
        y=array.shape[0]-1,
        sigma_x=array.shape[1]/2,
        sigma_y=array.shape[0]/2,
        rotation=np.pi,
        offset=local_min + (local_max - local_min)*0.2
    )
    
    try:
        X, Y = create_grid(*array.shape)
        popt, pcov = opt.curve_fit(
            f=gaussian_2d,
            xdata=(X, Y),
            ydata=array.ravel(),
            p0=initial_guess,
            bounds=(min_bounds, max_bounds),
            maxfev=10000  # 增加最大迭代次数
        )
        popt = Params(*popt)
    except (ValueError, RuntimeError) as e:
        print(f"Fit failed for spot at ({x_guess}, {y_guess}): {str(e)}")
        popt, pcov = None, None
    return popt, pcov

def create_grid(h, w):
    """ Creates a grid of x and y points to fit and evaluate over """
    x = np.arange(0, w, 1)
    y = np.arange(0, h, 1)
    x, y = np.meshgrid(x, y)
    return x, y

def evaluate_gaussian(x, y, popt, local_x_start=0, local_y_start=0):
    """ 把局部坐标还原为全局坐标后评估高斯 """
    global_popt = (popt.amp, popt.x + local_x_start, popt.y + local_y_start,
                   popt.sigma_x, popt.sigma_y, popt.rotation, popt.offset)
    z = gaussian_2d((x, y), *global_popt)
    return z

if __name__ == "__main__":
    np.random.seed(4)
    h, w = 200, 200
    x, y = create_grid(h=h, w=w)
    # 创建模拟数据
    img = []
    for _ in range(10):
        randx = np.random.randint(10, w - 10)
        randy = np.random.randint(10, h - 10)
        amp = 100
        d = gaussian_2d(
            xy_array=(x, y),
            amplitude=amp,
            pos_x=randx,
            pos_y=randy,
            sigma_x=9,
            sigma_y=3,
            rotation=3,
            offset=0
        )
        img.append(d)
    img = np.sum(img, axis=0).reshape(h, w)
    print("max intensity: {:.2f}".format(img.max()))
    
    # 检测光斑
    spots = skimage.feature.peak_local_max(img, num_peaks=20, min_distance=10)
    fig, ax = plt.subplots(ncols=2, figsize=(12,6))
    local_area = 20
    fit = []
    
    for idx, (pre_y, pre_x) in enumerate(spots):
        # 计算局部区域边界
        y_start = max(pre_y - local_area, 0)
        y_end = pre_y + local_area
        x_start = max(pre_x - local_area, 0)
        x_end = pre_x + local_area
        local_img = img[y_start:y_end, x_start:x_end]
        
        # 转换为局部坐标
        local_x = pre_x - x_start
        local_y = pre_y - y_start
        
        # 拟合高斯
        popt, pcov = fit_gaussian_spots(x_guess=local_x, y_guess=local_y, array=local_img)
        if popt is None:
            continue
        
        # 打印全局坐标参数
        global_popt = (popt.amp, popt.x + x_start, popt.y + y_start,
                       popt.sigma_x, popt.sigma_y, popt.rotation, popt.offset)
        print(f"Spot {idx}: amp={global_popt[0]:.2f}, x={global_popt[1]:.2f}, y={global_popt[2]:.2f}, sigma_x={global_popt[3]:.2f}, sigma_y={global_popt[4]:.2f}")
        
        # 绘图
        ax[0].add_patch(plt.Circle((pre_x, pre_y), 5, linewidth=0.5, fill=False, color="red"))
        ax[1].add_patch(plt.Rectangle((x_start, y_start), width=x_end-x_start, height=y_end-y_start, fill=False, color="yellow"))
        
        # 评估全局高斯并保存
        fit.append(evaluate_gaussian(x, y, popt, local_x_start=x_start, local_y_start=y_start))
    
    fit = np.sum(fit, axis=0)
    ax[0].set_title("Original Image")
    ax[0].imshow(img, origin="bottom", extent=(x.min(), x.max(), y.min(), y.max()))
    ax[1].set_title("Fitted Spots")
    ax[1].imshow(fit.reshape(img.shape), origin="bottom", extent=(x.min(), x.max(), y.min(), y.max()))
    plt.tight_layout()
    plt.show()

额外建议 #

• 如果光斑重叠严重，单光斑拟合效果有限，可以考虑多高斯联合拟合（同时拟合多个光斑参数），但复杂度会高很多；
• 先对图像做简单背景扣除（比如局部均值拟合），减少offset参数的拟合压力；
• 保留maxfev参数，给算法足够的迭代次数收敛。

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