Почему я не могу заставить 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()

Это приводит к следующим амплитудам + фон (оценивается напрямуюот увеличения, чтобы убедиться, что подгонка не чепуха):