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_x
andmu_y
, the centres of the Gaussian along the x- and y-axes,sig_u
andsig_v
, the width of the Gaussian as defined along and perpendicular to its diagonal,theta
, the degree of rotation, andoffset
, 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[0][0] + g[0][1] + g[0][2] + ... + g[63][63]
From my understanding, the formula for the combination of errors would then be:
var_int = diff(I, g)**2 * var_g[0][0] + diff(I, g)**2 * var_g[0][1] + ...
Where 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?
Thanks!
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)