Why can't I make 2D gaussian converge?

Question

I want to fit multiple spots in images to find their true intensity in a crowded area where local background correction fails.

My approach is just to pick out a local area of the image and fit that. The problem is, that the fit doesn't yield any useful results, and just defaults to the initial parameters. Adding bounds to help it along makes the fit not converge at all.

What am I doing wrong?

The code:

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

Turns out my biggest mistake was forgetting that coordinates to be fitted in a subset of the image are of course relative. In fact, just using the center is fine. I ended up using the following code, without bounds at all, as I found that fitting with initials only turned out to be a bit faster overall.

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()

This results in the following amplitudes + background (estimates directly from zoom to make sure fit is not nonsense):