Is there an easy way to calculate the gaussian fit for data that has been cut off by pixel saturation?

Question

I have image data for a flash of light where a few of the pixels have reached 255 at the centre of the flash. What would be the best way to go about fitting a 2D gaussian to the data that isn't at 255 in order to model what the hypothetical pixel value of the pixels would be.

I have been playing around with astropy.convolution, and the interpolate_replace_nans function, but this always returns values below 255 rather than the expected value of above 255.

The surface plot of the data

The surface plot of the data clearly follows a single gaussian shape, with the top cut off; is there a simple way to work out the parameters of the gaussian to reconstruct what the values should be?

Answer 1

If you mask out the pixels that are saturated, you can then fit a 2D Gaussian model to the data.

The masking can be as simple as not passing the bad data to the fitter, since the fitter expects to receive x,y coordinates and a corresponding z value.

(the convolution approach is not correct here; it can never return a value larger than was in the original data)

Answer 2

You can model the saturated gaussian, fit the data with the model, and recover the original gaussian with the fitting parameter.

Here is an example:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from numba import njit

@njit
def get_model(mu_x, mu_y, sig_x, sig_y, amp, cutoff, image_size):
    """
    simulate a 2D gaussian were the value higher than cutoff
        were saturated.
    """
    dom_x = -1 / 2 / sig_x ** 2
    dom_y = -1 / 2 / sig_y ** 2
    result = np.empty(image_size)
    for x in range(image_size[0]):
        for y in range(image_size[1]):
            gx = (x - mu_x)**2 * dom_x
            gy = (y - mu_y)**2 * dom_y
            g = amp * np.exp(gx + gy)
            g = min(g, cutoff)
            result[x, y] = g
    return result

def cost(args, data):
    model = get_model(*args, data.shape)
    return np.power(model - data, 2).sum()

img = get_model(25, 25, 10, 5, 300, 255, (70, 70))

result = minimize(
    cost, x0=(25, 25, 1, 1, 255, 255),
    args=(img)
)

fit = get_model(*result.x, img.shape)
gauss = get_model(
    *result.x[:-1],
    cutoff=np.inf,  # set cutoff to inf == no saturation
    image_size=img.shape
)


plt.figure(figsize=(10, 4))
plt.subplot(131).imshow(img, vmin=0, vmax=300)
plt.title("data")
plt.subplot(132).imshow(fit, vmin=0, vmax=300)
plt.title("fitting result")
plt.subplot(133).imshow(gauss, vmin=0, vmax=300)
plt.title("gaussian from fit")
plt.tight_layout()
plt.show()