How to get correct parameters for 2D Gaussian fit to an image with noise

Question

I have multiple images of objects where the object of interest (one per image) is near the center of said image. However, it's not necessarily the brightest source in the image, because some of the images also include sources that I am not interested in, but happen to be brighter/more intense. There is also usually a considerable amount of noise.

I would like to fit Gaussians to these 2D numpy image arrays, but am unsure how to effectively do so with these bright sources I don't want. In the end, I'll be stacking all the images (median) with the source of interest centered, so these other bright sources will just disappear. Below is my attempted code so far, where data is a list of multiple 2D arrays (images).

import numpy as np
import scipy

def Gaussian2D(x, y, x_0, y_0, theta, sigma_x, sigma_y, amp):
    a = np.cos(theta)**2/(2*sigma_x**2) + np.sin(theta)**2/(2*sigma_y**2)
    b = -np.sin(2*theta)/(4*sigma_x**2) + np.sin(2*theta)/(4*sigma_y**2)
    c = np.sin(theta)**2/(2*sigma_x**2) + np.cos(theta)**2/(2*sigma_y**2)
    exp_term = a * (x-x_0)**2
    exp_term += 2*b*(x-x_0)*(y-y_0)
    exp_term += c * (y-y_0)**2
    return amp * np.exp(-exp_term)

def GaussianFit(data):
    for data_set in data:
        y_0, x_0 = (np.shape(data_set)[0]//2, np.shape(data_set)[1]//2)
        sigma_x, sigma_y = np.std(data_set, axis=1), np.std(data_set, axis=0)
        fit = scipy.optimize.curve_fit(Gaussian2D(x, y, x_0, y_0, 0, sigma_x, sigma_y, amp), data_set)
    return fit

I've never done function fitting in a code, so I feel pretty lost. My specific questions are:

1. How can I define my parameters correctly? Do I need to flatten by array to get the sigma parameters? Also, I noticed in some example code that people made x and y arrays with linspace, so I'm not sure if I need to do that, and I'm also not sure what to put for the amplitude.

2. How would I handle the fact that I have multiple bright sources per image but only want to fit for the one closest to the center? Can I somehow specify to look near the center of the image?

3. I will also need the coordinates of the center source after fitting. How can I do this will ensuring it doesn't give me coordinates of other sources instead?

Any other help or advice is also appreciated. Thank you!

Answer 1

You can do this using a Gaussian Mixture Model. I don't think there is a function in SciPy, but there is one in scikit-learn

Here is a tutorial on this.

(from my answer to this question)

Then just remove the unwanted distribution from the image and fit to it.

Or there is skimage's blob detection.

On fitting a 2d Gaussian, read here. To use this you have to flatten the array as scipy's curve_fit only takes a 1d array. But it works fine.

Another approach is described here. A fit function with already three Gaussians in it is used. This would work if you know that there are always three (or in your case two) peaks on the image.