为什么我不能使 2D 高斯收敛？

Question

我想在图像中拟合多个点，以在局部背景校正失败的拥挤区域中找到它们的真实强度。

我的方法只是选择图像的一个局部区域并进行拟合。问题是，拟合不会产生任何有用的结果，只是默认为初始参数。添加边界来帮助它使拟合根本不收敛。

我做错了什么？

代码：

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


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
    a = (np.cos(rotation) ** 2) / (2 * sigma_x ** 2) + (
        np.sin(rotation) ** 2
    ) / (2 * sigma_y ** 2)
    b = -(np.sin(2 * rotation)) / (4 * sigma_x ** 2) + (
        np.sin(2 * rotation)
    ) / (4 * sigma_y ** 2)
    c = (np.sin(rotation) ** 2) / (2 * sigma_x ** 2) + (
        np.cos(rotation) ** 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

    initial_guess = Params(
        amp=1, x=x_guess, y=y_guess, sigma_x=1, sigma_y=1, rotation=0, offset=0
    )

    # Bounds makes it even harder to converge
    min_bounds = Params(
        amp=eps,
        x=x_guess * 0.5,
        y=y_guess * 0.5,
        sigma_x=eps,
        sigma_y=eps,
        rotation=-np.inf,
        offset=eps,
    )

    max_bounds = Params(
        amp=np.max(array),
        x=x_guess * 1.5,
        y=y_guess * 1.5,
        sigma_x=np.inf,
        sigma_y=np.inf,
        rotation=2 * np.pi,
        offset=np.max(array),
    )

    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),
        )
        popt = Params(*np.round(popt))

    except ValueError:
        print("fit didn't converge!")
        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):
    """
    Evaluates gaussian in coordinate positions.
    NOTE: this is not necessary for extracting intensity,
    as the pure signal is fitted as the amplitude.
    """
    z = gaussian_2d((x, y), *popt)
    return z


if __name__ == "__main__":
    # Create x and y indices
    np.random.seed(4)
    h, w = 200, 200
    x, y = create_grid(h=h, w=w)

    # create data
    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,
        )
        # d = d + np.random.normal(0, 5, d.shape) # add noise
        # d += 100  # add offset
        img.append(d)
    img = np.sum(img, axis=0)
    img = img.reshape(h, w)
    print("max intensity: {:.2f}".format(img.max()))

    # Detect soem possible spots first
    spots = skimage.feature.peak_local_max(img, num_peaks=20, min_distance=10)
    fig, ax = plt.subplots(ncols=2)

    h, w = img.shape
    local_area = 20
    fit = []

    # skimage returns rows, columns (y,x) while matplotlib operates in (x,y)
    for idx, (pre_y, pre_x) in enumerate(spots):
        # Fit gaussian in local area
        popt, pcov = fit_gaussian_spots(
            x_guess=pre_x,
            y_guess=pre_y,
            # Avoid falling off the edge of the image
            array=img[
                max(pre_y - local_area, 0) : pre_y + local_area,
                max(pre_x - local_area, 0) : pre_x + local_area,
            ],
        )
        if popt is None:
            continue
        print(popt)

        ax[0].add_patch(
            plt.Circle(
                (pre_x, pre_y), 5, linewidth=0.5, fill=False, color="red"
            )
        )
        ax[1].add_patch(
            plt.Rectangle(
                (pre_x - local_area, pre_y - local_area),
                width=local_area * 2,
                height=local_area * 2,
                fill=False,
                color="yellow",
            )
        )

        fit.append(evaluate_gaussian(x, y, popt))
    fit = np.sum(fit, axis=0)

    ax[0].set_title("true")
    ax[0].imshow(
        img, origin="bottom", extent=(x.min(), x.max(), y.min(), y.max())
    )
    ax[1].set_title("predicted")
    ax[1].imshow(
        fit.reshape(img.shape),
        origin="bottom",
        extent=(x.min(), x.max(), y.min(), y.max()),
    )

    plt.show()

Answer 1

原来我最大的错误是忘记了要安装在图像子集中的坐标当然是相对的。事实上，只使用中心就可以了。我最终使用了以下代码，完全没有限制，因为我发现使用首字母进行拟合总体上会更快一些。

import scipy.optimize as opt
import numpy as np
import matplotlib.pyplot as plt
import skimage.feature
from collections import namedtuple
import skimage.io
import matplotlib.patches
import skimage.filters
import warnings
from scipy.optimize import OptimizeWarning


def zoom_array(array, xy, square_radius):
    """
    Return a zoomed array at location
    """
    x, y = xy
    minix = int(max(x - square_radius, 0))
    miniy = int(max(y - square_radius, 0))
    maxix = int(x + square_radius)
    maxiy = int(y + square_radius)
    return array[miniy:maxiy, minix:maxix]


def gaussian_2d(
    xy_array, amplitude, pos_x, pos_y, sigma_x, sigma_y, angle, offset
):
    """
    Expression for a 2D gaussian function with variance in both x and y
    """
    x, y = xy_array

    a = (np.cos(angle) ** 2) / (2 * sigma_x ** 2) + (np.sin(angle) ** 2) / (
        2 * sigma_y ** 2
    )
    b = -(np.sin(2 * angle)) / (4 * sigma_x ** 2) + (np.sin(2 * angle)) / (
        4 * sigma_y ** 2
    )
    c = (np.sin(angle) ** 2) / (2 * sigma_x ** 2) + (np.cos(angle) ** 2) / (
        2 * sigma_y ** 2
    )

    g = offset + amplitude * np.exp(
        -(
            a * ((x - pos_x) ** 2)
            + 2 * b * (x - pos_x) * (y - pos_y)
            + c * ((y - pos_y) ** 2)
        )
    )
    return g.ravel()


def fit_gaussian_spots(x_guess, y_guess, array):
    Params = namedtuple(
        "Parameters", "amp, x, y, sigma_x, sigma_y, angle, offset"
    )

    initial_guess = Params(
        amp=np.max(array),
        x=x_guess,
        y=y_guess,
        sigma_x=1,
        sigma_y=1,
        angle=0,
        offset=np.abs(np.min(array)),
    )

    with warnings.catch_warnings():
        warnings.simplefilter("error", OptimizeWarning)
        try:
            X, Y = create_grid(*array.shape)
            popt, pcov = opt.curve_fit(
                f=gaussian_2d,
                xdata=(X, Y),
                ydata=array.ravel(),
                p0=initial_guess,
                maxfev=200,
                method="lm"
                # constraints make it slower. Better to time out bad fits
                # bounds=(min_bounds, max_bounds),
            )
            popt = Params(*np.round(popt))
        except (OptimizeWarning, ValueError, RuntimeError):
            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):
    """
    Evaluates gaussian in coordinate positions.
    NOTE: this is not necessary for extracting intensity,
    as the pure signal is fitted as the amplitude.
    """
    z = gaussian_2d((x, y), *popt)
    return z


def _simulate_data():
    """Create data"""
    img = []
    for _ in range(N_SPOTS):
        POSX = np.random.randint(0, W)
        POSY = np.random.randint(0, H)
        AMP = 100
        g = gaussian_2d(
            xy_array=(Xi, Yi),
            amplitude=AMP,
            pos_x=POSX,
            pos_y=POSY,
            sigma_x=2,
            sigma_y=2,
            angle=0,
            offset=0,
        )
        img.append(g)
    img = np.sum(img, axis=0)
    img = img.reshape(H, W)
    img = img + np.random.normal(5, 5, len(img.ravel())).reshape(img.shape)
    return img


if __name__ == "__main__":
    # Create x and y indices
    np.random.seed(4)
    PLOT_ROWS = 3
    PLOT_COLS = 3

    H, W = 200, 400
    N_SPOTS = 50
    Xi, Yi = create_grid(h=H, w=W)
    img = _simulate_data()

    # Detect a generous amount of spots to be safe
    spots = skimage.feature.peak_local_max(img, num_peaks=300, min_distance=5)

    figimg, aximg = plt.subplots()
    aximg.imshow(
        img, origin="bottom", extent=(Xi.min(), Xi.max(), Yi.min(), Yi.max())
    )

    figzoom, axzoom = plt.subplots(nrows=PLOT_ROWS, ncols=PLOT_COLS)
    axzoom = axzoom.ravel()

    zoom_ctr = 6
    # skimage returns rows, columns (y,x) while matplotlib operates in (x,y)
    idx = 0
    for guessy, guessx in spots:
        # Plot on the full image
        # Initial
        aximg.add_patch(
            plt.Circle(
                (guessx, guessy),
                3,
                linewidth=0.5,
                fill=False,
                alpha=0.5,
                color="yellow",
            )
        )

        # Fit
        local_arr = zoom_array(img, (guessx, guessy), square_radius=zoom_ctr)
        popt, pcov = fit_gaussian_spots(
            x_guess=zoom_ctr, y_guess=zoom_ctr, array=local_arr
        )

        # Throw away bad fits
        if popt is None or popt.sigma_x < 1 or popt.sigma_y < 1:
            continue

        predx = guessx + popt.x - zoom_ctr
        predy = guessy + popt.y - zoom_ctr

        # Plot on each of zooms
        # Predicted
        try:
            axzoom[idx].imshow(local_arr, origin="bottom")
            axzoom[idx].add_patch(
                matplotlib.patches.Ellipse(
                    (popt.x, popt.y),
                    width=popt.sigma_x * 3,
                    height=popt.sigma_y * 3,
                    angle=popt.angle,
                    color="red",
                    fill=False,
                )
            )
            axzoom[idx].set_title(
                "fit: {:.1f} + {:.1f}\n"
                "est: {:.1f} + {:.1f}".format(
                    popt.amp, popt.offset, np.max(local_arr), np.min(local_arr)
                )
            )
        except IndexError:
            pass

        # Predicted
        aximg.add_patch(
            plt.Circle((predx, predy), 4, linewidth=2, fill=False, color="red")
        )

        idx += 1

    plt.tight_layout()
    plt.show()

这会产生以下幅度 + 背景（直接从缩放中估计以确保拟合不是废话）：