Issues with fitting 2DGaussian

[caradoct]

Hi, I am using lmfit to fit a beam on a synchrotron beamline. My issue is that the beam is wider than it is tall and when I use lmfit the resulting fit is not a good match to the beam.

My lmfit code is (x2n and X1n are the pixels in the camera converted to 1D arrays and data is the flattened 2D greyscale image data):

model = lmfit.models.Gaussian2dModel(independent_vars=['x', 'y'])
params = model.guess(data, x2n, x1n)
params['sigmax'].set(value=35, min=0)
params['sigmay'].set(value=10, min=0)

result = model.fit(data, x=x2n, y=x1n, params=params, weights=1/error)
print(result.fit_report())

The report is:

[[Model]]
Model(gaussian2d)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 145
# data points = 179776
# variables = 5
chi-square = 161168.816
reduced chi-square = 0.89652289
Akaike info crit = -19632.1942
Bayesian info crit = -19581.6969
R-squared = 0.91153405
[[Variables]]
amplitude: 707119.250 +/- 810.808671 (0.11%) (init = 949317.8)
centerx: 191.936035 +/- 0.04079530 (0.02%) (init = 188)
centery: 217.960223 +/- 0.02992817 (0.01%) (init = 221)
sigmax: 34.6185785 +/- 0.03103234 (0.09%) (init = 35)
sigmay: 25.7437069 +/- 0.02762418 (0.11%) (init = 10)
fwhmx: 81.5205210 +/- 0.07307557 (0.09%) == '2.3548200sigmax'
fwhmy: 60.6217959 +/- 0.06504996 (0.11%) == '2.3548200sigmay'
height: 126.279357 +/- 0.22882345 (0.18%) == '0.1591549*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))'

From the plots you can see the centers in X and Y seem reasonable but the sigmaY is overestimated (and FWHM in Y).

Woudl anyone be able to help me to understand what I am doing wrong? I thought stating x and y are independent vars would work or are my starting parameters too far out?

Cheers,

Tom

[newville]

@caradoct Great to see a fellow synchrotron person using this library.

We always ask to see a complete-ish example. I would suspect that the challenge with Gaussian2D is figuring out what the shapes of the arrays should be. An example usage would be:

import numpy as np

from lmfit.models import Gaussian2dModel
from lmfit.lineshapes import gaussian2d

import wxmplot.interactive as wi

x = np.linspace(0, 15, 151)
y = np.linspace(0, 20, 201)

# make the (nx, ny) shaped arrays of x, y values, and compute z for each of these nx*ny points
xmesh, ymesh = np.meshgrid(x, y)

zmesh = gaussian2d(xmesh, ymesh, amplitude=580,
                   centerx=6.6, centery=11.5,
                   sigmax=2.2, sigmay=3.1)

noise = np.random.normal(scale=1, size=zmesh.size)
noise.shape = zmesh.shape

data = zmesh + noise

# Now use the flattened (1D) arrays from the mesh data for the fit
model = Gaussian2dModel()
params = model.guess(data.flatten(), xmesh.flatten(), ymesh.flatten())
result = model.fit(data.flatten(), x=xmesh.flatten(), y=ymesh.flatten(), params=params)

print(result.fit_report())

# reshape the flattened best fit to match the data
bfit = result.best_fit.reshape(zmesh.shape)

# I'm sure this could be done with matplotlib,
# but I didn't bother to take the time to do that
wi.imshow(data, x=x, y=y, show_axis=True, title='Data', colormap='coolwarm')
wi.imshow(bfit, x=x, y=y, show_axis=True, title='Fit', colormap='coolwarm', win=2)
wi.imshow(bfit-data, x=x, y=y, show_axis=True, title='residual', colormap='coolwarm', win=3)

If you're having trouble, I might suggest starting with this example.

[caradoct]

Hi Matt, Thanks for the help! I get the same issue when using your approach as with my original one. My beam is a reasonable approximation to a gaussian in the horizontal but in the vertical has a very different profile. It is almost a delta function and I think gaussian2d is having issues with it. The complete-ish example is: 1. I take an image from our sample camera. This is 8-bit and I attached one frame converted to a numpy array. Processing code: im = get_image_from_md3_camera(np.uint16) halfheight = im.max()/2 imagemax = im.max() X, Y = np.meshgrid(np.arange(0, im.shape[1], 1), np.arange(0, im.shape[0], 1)) Xas1D = X.flatten() Yas1D = Y.flatten() data=im.flatten() error=np.sqrt(data+1) Z = data model = lmfit.models.Gaussian2dModel(independent_vars=['x', 'y']) #params = model.guess(data, x2n, x1n) params = model.guess(data, Xas1D, Yas1D) #params = model.make_params(centerx=512,centery=612,sigmax=33,sigmay=2,amplitude=960000) params['sigmax'].set(value=35, min=0) params['sigmay'].set(value=1e-4, min=0) result = model.fit(data, x=Xas1D, y=Yas1D, params=params, weights=1/error) print(result.fit_report()) I get a result from lmfit like this: [[Model]] Model(gaussian2d) [[Fit Statistics]] # fitting method = leastsq # function evals = 139 # data points = 1253376 # variables = 5 chi-square = 149886.630 reduced chi-square = 0.11958680 Akaike info crit = -2661805.64 Bayesian info crit = -2661745.44 R-squared = 0.94099492 [[Variables]] amplitude: 962701.867 +/- 345.958157 (0.04%) (init = 5803848) centerx: 608.558070 +/- 0.01421396 (0.00%) (init = 606) centery: 543.503355 +/- 0.01221779 (0.00%) (init = 545) sigmax: 38.4560917 +/- 0.01100216 (0.03%) (init = 35) sigmay: 33.4608534 +/- 0.01100894 (0.03%) (init = 0.0001) fwhmx: 90.5571740 +/- 0.02590810 (0.03%) == '2.3548200*sigmax' fwhmy: 78.7942868 +/- 0.02592406 (0.03%) == '2.3548200*sigmay' height: 119.072006 +/- 0.06711452 (0.06%) == '0.1591549*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))' I can plot it using: from scipy.interpolate import griddata X, Y = np.meshgrid(np.arange(0, im.shape[1], 1), np.arange(0, im.shape[0], 1)) Z = griddata((x2n, x1n), data, (X, Y), method='linear', fill_value=0) fig, axs = plt.subplots(2, 2, figsize=(10, 10)) vmax = np.nanpercentile(Z, 99.9) ax = axs[0, 0] art = ax.pcolor(X, Y, Z, vmin=0, vmax=vmax, shading='auto') plt.colorbar(art, ax=ax, label='z') ax.set_title('Data') ax = axs[0, 1] fit = model.func(X, Y, **result.best_values) art = ax.pcolor(X, Y, fit, vmin=0, vmax=vmax, shading='auto') plt.colorbar(art, ax=ax, label='z') ax.set_title('Fit') ax = axs[1, 0] fit = model.func(X, Y, **result.best_values) art = ax.pcolor(X, Y, Z-fit, vmin=0, vmax=10, shading='auto') plt.colorbar(art, ax=ax, label='z') ax.set_title('Data - Fit') for ax in axs.ravel(): ax.set_xlabel('x') ax.set_ylabel('y') axs[1, 1].remove() plt.show() Looks like this: ***@***.*** What is the best way to process this type of image? Cheers, Tom OFFICIAL From: Matt Newville ***@***.***> Sent: Thursday, 3 July 2025 1:11 PM To: lmfit/lmfit-py ***@***.***> Cc: CARADOC DAVIES, Tom ***@***.***>; Mention ***@***.***> Subject: Re: [lmfit/lmfit-py] Issues with fitting 2DGaussian (Discussion #1007) You don't often get email from ***@***.******@***.***>. Learn why this is important<https://aka.ms/LearnAboutSenderIdentification> CAUTION, EXTERNAL EMAIL: This message has come from outside of ANSTO. Do not take action, click links or open attachments unless you trust the source of this message and know the content is safe. Report spam and phishing using the Report Message button or if unsure, forward this message to ***@***.******@***.***> as an attachment. @caradoct<https://github.com/caradoct> Great to see a fellow synchrotron person using this library. We always ask to see a complete-ish example. I would suspect that the challenge with Gaussian2D is figuring out what the shapes of the arrays should be. An example usage would be: import numpy as np from lmfit.models import Gaussian2dModel from lmfit.lineshapes import gaussian2d import wxmplot.interactive as wi x = np.linspace(0, 15, 151) y = np.linspace(0, 20, 201) # make the (nx, ny) shaped arrays of x, y values, and compute z for each of these nx*ny points xmesh, ymesh = np.meshgrid(x, y) zmesh = gaussian2d(xmesh, ymesh, amplitude=580, centerx=6.6, centery=11.5, sigmax=2.2, sigmay=3.1) noise = np.random.normal(scale=1, size=zmesh.size) noise.shape = zmesh.shape data = zmesh + noise # Now use the flattened (1D) arrays from the mesh data for the fit model = Gaussian2dModel() params = model.guess(data.flatten(), xmesh.flatten(), ymesh.flatten()) result = model.fit(data.flatten(), x=xmesh.flatten(), y=ymesh.flatten(), params=params) print(result.fit_report()) # reshape the flattened best fit to match the data bfit = result.best_fit.reshape(zmesh.shape) # I'm sure this could be done with matplotlib, # but I didn't bother to take the time to do that wi.imshow(data, x=x, y=y, show_axis=True, title='Data', colormap='coolwarm') wi.imshow(bfit, x=x, y=y, show_axis=True, title='Fit', colormap='coolwarm', win=2) wi.imshow(bfit-data, x=x, y=y, show_axis=True, title='residual', colormap='coolwarm', win=3) If you're having trouble, I might suggest starting with this example. - Reply to this email directly, view it on GitHub<#1007 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AM3DJX7MOKWBN37VOW4AVMD3GSNL5AVCNFSM6AAAAACAVJWE2OVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTGNRUG4YTONA>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

[newville]

@caradoct that seems like the fit worked, giving non-trivial values for center and sigma in both directions. Could it just be a display problem? Like, what if you do:

plt.imshow(result.best_fit.reshape(im.shape))

I would expect that to show a 2D Gaussian with center at X,Y = 608, 543, and sigmax, sigmay = 38,33.

[caradoct]

Hi Matt, The beam at sample is around 9x3 microns (FWHM) so I expect the vertical to be ~1/3rd the horizontal but lmfit gives close to 1:1. I have been using tungsten blades to do knife-edge scans at the sample position with a downstream photodiode to measure the beam size and it matches what we see on the scintillator. The vertical sigma seems too high but I don't know what the issue is. Cheers, Tom OFFICIAL From: Matt Newville ***@***.***> Sent: Tuesday, 22 July 2025 10:05 PM To: lmfit/lmfit-py ***@***.***> Cc: CARADOC DAVIES, Tom ***@***.***>; Mention ***@***.***> Subject: Re: [lmfit/lmfit-py] Issues with fitting 2DGaussian (Discussion #1007) You don't often get email from ***@***.******@***.***>. Learn why this is important<https://aka.ms/LearnAboutSenderIdentification> CAUTION, EXTERNAL EMAIL: This message has come from outside of ANSTO. Do not take action, click links or open attachments unless you trust the source of this message and know the content is safe. Report spam and phishing using the Report Message button or if unsure, forward this message to ***@***.******@***.***> as an attachment. @caradoct<https://github.com/caradoct> that seems like the fit worked, giving non-trivial values for center and sigma in both directions. Could it just be a display problem? Like, what if you do: plt.imshow(result.best_fit.reshape(im.shape)) I would expect that to show a 2D Gaussian with center at X,Y = 608, 543, and sigmax, sigmay = 38,33. - Reply to this email directly, view it on GitHub<#1007 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AM3DJX7TICTFI2MO64WVBST3JYSFRAVCNFSM6AAAAACAVJWE2OVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTGOBUG4ZTINY>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

[caradoct]

Hi Matt, I attached pics using: plt.imshow(result.best_fit.reshape(im2.shape)) and plt.imshow(im2) As you can see the input data is very different in Y compared to the lmfit result. Cheers, Tom OFFICIAL From: Matt Newville ***@***.***> Sent: Tuesday, 22 July 2025 10:05 PM To: lmfit/lmfit-py ***@***.***> Cc: CARADOC DAVIES, Tom ***@***.***>; Mention ***@***.***> Subject: Re: [lmfit/lmfit-py] Issues with fitting 2DGaussian (Discussion #1007) You don't often get email from ***@***.******@***.***>. Learn why this is important<https://aka.ms/LearnAboutSenderIdentification> CAUTION, EXTERNAL EMAIL: This message has come from outside of ANSTO. Do not take action, click links or open attachments unless you trust the source of this message and know the content is safe. Report spam and phishing using the Report Message button or if unsure, forward this message to ***@***.******@***.***> as an attachment. @caradoct<https://github.com/caradoct> that seems like the fit worked, giving non-trivial values for center and sigma in both directions. Could it just be a display problem? Like, what if you do: plt.imshow(result.best_fit.reshape(im.shape)) I would expect that to show a 2D Gaussian with center at X,Y = 608, 543, and sigmax, sigmay = 38,33. - Reply to this email directly, view it on GitHub<#1007 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AM3DJX7TICTFI2MO64WVBST3JYSFRAVCNFSM6AAAAACAVJWE2OVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTGOBUG4ZTINY>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

[newville]

@caradoct I do not see any images.

The fit report you showed earlier definitely worked and gave reasonable values, with uncertainties.
The code (complete and runnable) I attached definitely runs and gives reasonable results.

If you want further help, you must show what you are doing. Post actual code that runs.

[caradoct]

Hi Matt, I attached the images here.
input data:

lmfit result:

The results are not reasonable values. The vertical size (sigmaY) is too large.

[caradoct]

Apologies, I was replying in email and the code and images I attached did not come through correctly to github. I am adding them now.

[caradoct]

The working image (as a numpy array and github won't let me upload that directly):
https://drive.google.com/file/d/1r9wKWq6Kd_RSorkJM2IQV3a1lvQl6MJq/view?usp=sharing

[caradoct]

Working code (please edit the path to the numpy array on your local system):

#im = get_image_from_md3_camera(np.uint16)
#edit the below path to your local copy of the np array
im = np.load("/mnt/disk/commissioning/notebooks/scandata/md3_image_for_lmfit.npy")

halfheight = im.max()/2
imagemax = im.max()

X, Y = np.meshgrid(np.arange(0, im.shape[1], 1), np.arange(0, im.shape[0], 1))
Xas1D = X.flatten()
Yas1D = Y.flatten()

data=im.flatten()
error=np.sqrt(data+1)

Z = data

model = lmfit.models.Gaussian2dModel(independent_vars=['x', 'y'])
#params = model.guess(data, x2n, x1n)
params = model.guess(data, Xas1D, Yas1D)
#params = model.make_params(centerx=512,centery=612,sigmax=33,sigmay=2,amplitude=960000)

params['sigmax'].set(value=35, min=0)
params['sigmay'].set(value=1e-4, min=0)

result = model.fit(data, x=Xas1D, y=Yas1D, params=params, weights=1/error)
print(result.fit_report())

from scipy.interpolate import griddata

X, Y = np.meshgrid(np.arange(0, im.shape[1], 1), np.arange(0, im.shape[0], 1))
Z = griddata((x2n, x1n), data, (X, Y), method='linear', fill_value=0)

fig, axs = plt.subplots(2, 2, figsize=(10, 10))

vmax = np.nanpercentile(Z, 99.9)

ax = axs[0, 0]
art = ax.pcolor(X, Y, Z, vmin=0, vmax=vmax, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Data')

ax = axs[0, 1]
fit = model.func(X, Y, **result.best_values)
art = ax.pcolor(X, Y, fit, vmin=0, vmax=vmax, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Fit')

ax = axs[1, 0]
fit = model.func(X, Y, **result.best_values)
art = ax.pcolor(X, Y, Z-fit, vmin=0, vmax=10, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Data - Fit')

for ax in axs.ravel():
ax.set_xlabel('x')
ax.set_ylabel('y')
axs[1, 1].remove()
plt.show()

[caradoct]

As I said above lmfit runs and gives a result, the issue is that while the result is reasonable in X it is out by a factor of 3 in Y.

[newville]

@caradoct

For me, doing

import numpy as np
from lmfit.models import Gaussian2dModel
import matplotlib.pyplot as plt

img = np.load("md3_image_for_lmfit.npy")

ny, nx = img.shape
x, y = np.arange(nx), np.arange(ny)
xmesh, ymesh = np.meshgrid(x, y)

model = Gaussian2dModel()
params = model.guess(img.flatten(), xmesh.flatten(), ymesh.flatten())
result = model.fit(img.flatten(), params, x=xmesh.flatten(), y=ymesh.flatten())

imgfit = result.best_fit.reshape((ny, nx))

print(result.fit_report())

fig, axes = plt.subplots(2, 1, figsize=(5, 8))
axes[0].pcolor(x, y, img)
axes[0].set_title('data')

axes[1].pcolor(x, y, imgfit)
axes[1].set_title('fit')
plt.show()

gives

[[Model]]
    Model(gaussian2d)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 69
    # data points      = 1253376
    # variables        = 5
    chi-square         = 14241415.5
    reduced chi-square = 11.3624900
    Akaike info crit   = 3046106.73
    Bayesian info crit = 3046166.93
    R-squared          = 0.95449183
[[Variables]]
    amplitude:  1966184.54 +/- 538.129872 (0.03%) (init = 8862164)
    centerx:    612.590150 +/- 0.01079401 (0.00%) (init = 594)
    centery:    510.640449 +/- 0.00703724 (0.00%) (init = 493)
    sigmax:     39.4384984 +/- 0.01079412 (0.03%) (init = 203.8333)
    sigmay:     25.7123721 +/- 0.00703731 (0.03%) (init = 170.5)
    fwhmx:      92.8705648 +/- 0.02541820 (0.03%) == '2.3548200*sigmax'
    fwhmy:      60.5480081 +/- 0.01657160 (0.03%) == '2.3548200*sigmay'
    height:     308.589948 +/- 0.08445896 (0.03%) == '0.1591549*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))'
[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, sigmay) = +0.5000
    C(amplitude, sigmax) = +0.5000

and

that seems like it is working to me. sigmax is 1.5 times sigmay.

I'm not sure if what you are seeing is a difference in image display, or maybe from the weighting you are using -- I would tend to think that the uncertainty in a 16-bit image would not be "1+sqrt(im"). And the data does have stronger tails in the vertical than in the horizontal, so if you weight by 1/sqrt(im), you might be prone to getting a "rounder" (more round? English!) result.