I have been trying to fit signals to a 2D Gaussian function, and while I have bene able to use sciKit-image's
curve_fit function to find the covariance matrix for the parameters, I am at a loss as to how to estimate the error of the measured Gaussian, as well as the error in the total integrated intensity.
The current parameters I have for my 2D Gaussian are:
amplitude- The scaling factor,
mu_y, the centres of the Gaussian along the x- and y-axes,
sig_v, the width of the Gaussian as defined along and perpendicular to its diagonal,
theta, the degree of rotation, and
offset, the offset from 0.
Based on my reading, I have tried a couple of methods to quantify the error in the Gaussian (fitted over a 64 x 64 pixel grid, for reference), and would like advice as to which is more appopriate (or if none of them are):
1. Squared Difference Between Raw and Fitted Signals
To start with, I have taken the variance of the Gaussian to be the squared difference between the raw signal and the fitted signal for each pixel. The equation used is as follows:
# Variance for each pixel var_g = (signal - model)**2 # returns a 64 x 64 array with variances at each point
Then, in order to find the total integrated intensity, I sum the model array:
# Calculating the total integrated intensity int_sum = model.sum() # This is summed over 64 x 64 = 4096 pixels int_sum = g + g + g + ... + g
From my understanding, the formula for the combination of errors would then be:
var_int = diff(I, g)**2 * var_g + diff(I, g)**2 * var_g + ...
diff(I, g) means differentiating the intensity function with respect to each individual point along the Gaussian function. Since it's a linear function, my understanding is that
diff(I, g) is just one, and
var_int is therefore:
var_int = var_g.sum()
2. Combination of Errors
The other method I tried involved taking the parameters from the 2D Gaussian and using them to calculate the variance the Gaussian for each pixel, again using the law of combination of errors:
# Variance at each pixel var_g = diff(g, amp)**2 * var_amp + diff(g, mu_x)**2 * mu_x + ... + diff(g, offset)**2 * offset
Again, because the total integrated signal is a linear equation in terms of the Gaussian function, the variance of the total integrated signal is then just:
var_int = var_g.sum()
Which method would be more appropriate in this case, and if not, what suggestions/changes would you recommend I make?
UPDATE: The code I used to generate the 2D Gaussian is as follows.
import numpy as np def Gaussian2D_v5(coords=None, # x and y coordinates for each image. amplitude=1, # Highest intensity in image. xo=0, # x-coordinate of peak centre. yo=0, # y-coordinate of peak centre. sig_u=1, # Standard deviation in u. sig_v=1, # Standard deviation in v. theta=0, # Orientation of the diagonal matrix, in degrees offset=0): # Offset from zero (background radiation). # Set x and y coordinates x, y = coords # Convert numbers to floats xo = float(xo) yo = float(yo) # Convert theta to radians theta = np.deg2rad(theta) # Define diagonal matrix mat_diag = np.asarray([[sig_u**2, 0], [0, sig_v**2]]) # mat_diag = np.asarray(mat_diag) # Define rotation matrix mat_rot = np.asarray([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) # mat_rot = np.asarray(mat_rot) # Define covariance matrix mat_cov = np.matmul(np.matmul(mat_rot, mat_diag), np.linalg.inv(mat_rot)) # Set the coordinates and stack them along the last axis mat_coords = np.stack((x - xo, y - yo), axis=-1) mat_coords_row = mat_coords[:, :, np.newaxis, :] mat_coords_col = mat_coords[:, :, :, np.newaxis] # Perform calculation across the meshgrid G = amplitude * np.exp(-0.5 * np.matmul(np.matmul(mat_coords_row, np.linalg.inv(mat_cov)), mat_coords_col)) + offset # Reduce the dimensions of the answer return G.squeeze().ravel() # Generating the noisy Gaussian coords = np.meshgrid(np.arange(0, 64), np.arange(0, 64)) model = Gaussian2D_v5(coords, amplitude=20, xo=32, yo=32, sig_u=6, sig_v=3, theta=20, offset=20).reshape(64, 64) noise = np.random.normal(0, # mean 10, # standard deviation (64, 64)) # array shape signal = model + noise # Curve fitting # Set initial guess p0 = [20, # amplitude 32, # xo 32, # yo 2, # sig_u 2, # sig_v 0, # theta 0] # offset # Set upper and lower bounds bound_lo = [0, 0, 0, 0, 0, 0, 0] bound_hi = [100, 64, 64, np.inf, np.inf, 90, np.inf] bounds = (bound_lo, bound_hi) # Fit to 2D Gaussian equation popt, pcov = curve_fit(Gaussian2D_v5, coords, signal.ravel(), p0=p0, bounds=bounds) # Generated fitted curve model = Gaussian2D_v5(coords, *popt).reshape(64, 64)