arrays - ValueError：形状（2,210）和（2,210）未对齐：210（dim 1）！= 2（dim 0）

Question

我想使用 2D 图像数据拟合强度分布函数，scipy.optimize.curve_fit但无法在我的代码中找到错误：

# Define doughnut beam intensity distribution function
def doughnut(x, y, x0, y0, A, FWHM):
    '''2D intensity distribution function of doughnut beams (DOI: 10.1126/science.aak9913,
    https://science.sciencemag.org/content/sci/suppl/2016/12/21/science.aak9913.DC1/Balzarotti_SM.pdf).

    Parameters
    ----------
    x, y : float
        X and Y coordinates, orthogonal to beam axis
    x0 : float
        X offset
    y0 : float
        Y offset
    A : float
        Peak intensity
    FWHM : float
        Full width at half maximum
    '''
    return A*np.exp(1)*4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2*np.exp(-4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2)

# Read image file names
pathname = '/home/user/doughnut_beam/'
filenameList = [filename for filename in os.listdir(pathname)
                if filename.endswith('.tif')]

# Open image files, fit doughnut beam intensity distribution function
for filename in filenameList:
    img = Image.open(pathname + filename)

    X, Y = img.size
    xRange = np.arange(1, X+1)
    yRange = np.arange(1, Y+1)
    xGrid, yGrid = np.meshgrid(xRange, yRange)
    xyGrid = np.vstack((xGrid.ravel(), yGrid.ravel()))    # scipy.optimize.curve_fit requires 2xN-array 

    imgArray = np.array(img)
    imgArrayFlat = imgArray.ravel()    # Flatten 2D pixel data into 1D array for scipy.optimize.curve_fit
    
    params_opt, params_cov = curve_fit(doughnut, xyGrid, imgArrayFlat)

这是 Jupyter Notebook 的输出：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-44-eaa3ebdb6469> in <module>()
     17     imgArrayFlat = imgArray.ravel()    # Flatten 2D pixel data into 1D array for scipy.optimize.curve_fit
     18 
---> 19     params_opt, params_cov = curve_fit(doughnut, xyGrid, imgArrayFlat)

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    749         # Remove full_output from kwargs, otherwise we're passing it in twice.
    750         return_full = kwargs.pop('full_output', False)
--> 751         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    752         popt, pcov, infodict, errmsg, ier = res
    753         cost = np.sum(infodict['fvec'] ** 2)

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    381     if not isinstance(args, tuple):
    382         args = (args,)
--> 383     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    384     m = shape[0]
    385     if n > m:

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     25 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     26                 output_shape=None):
---> 27     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     28     if (output_shape is not None) and (shape(res) != output_shape):
     29         if (output_shape[0] != 1):

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in func_wrapped(params)
    461     if transform is None:
    462         def func_wrapped(params):
--> 463             return func(xdata, *params) - ydata
    464     elif transform.ndim == 1:
    465         def func_wrapped(params):

<ipython-input-43-3e0adae6fbe0> in doughnut(x, y, x0, y0, A, FWHM)
     17         Full width at half maximum
     18     '''
---> 19     return A*np.exp(1)*4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2*np.exp(-4*np.log(2)*(np.dot(x+x0,x+x0) + np.dot(y+y0,y+y0))/FWHM**2)

ValueError: shapes (2,210) and (2,210) not aligned: 210 (dim 1) != 2 (dim 0)

更新：由于某种原因，使用numpy.dot平方（偏移）变量x+x0和y+y0函数定义不起作用。只需更改为**运算符即可得到正确的绘图：

# UPDATED: Define doughnut beam intensity distribution function
def doughnut(x, y, x0, y0, A, FWHM):
    '''2D intensity distribution function of doughnut beams (DOI: 10.1126/science.aak9913,
    https://science.sciencemag.org/content/sci/suppl/2016/12/21/science.aak9913.DC1/Balzarotti_SM.pdf).

    Parameters
    ----------
    x, y : float
        X and Y coordinates, orthogonal to beam axis
    x0 : float
        X offset
    y0 : float
        Y offset
    A : float
        Peak intensity
    FWHM : float
        Full width at half maximum
    '''
    return A*np.exp(1)*4*np.log(2)*((x+x0)**2 + (y+y0)**2)/FWHM**2*np.exp(-4*np.log(2)*((x+x0)**2 + (y+y0)**2)/FWHM**2)


fig = plt.figure()
ax = fig.gca(projection='3d')

# Make data
X = np.arange(-10, 10, 0.25)
Y = np.arange(-10, 10, 0.25)
X, Y = np.meshgrid(X, Y)
Z = doughnut(X, Y, x0=0, y0=0, A=1.5, FWHM=7)

# Plot the surface
surf = ax.plot_surface(X, Y, Z)
plt.show()

=>情节

但是：现在我在尝试拟合数据时遇到了一个新错误：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-61-eaa3ebdb6469> in <module>()
     17     imgArrayFlat = imgArray.ravel()    # Flatten 2D pixel data into 1D array for scipy.optimize.curve_fit
     18 
---> 19     params_opt, params_cov = curve_fit(doughnut, xyGrid, imgArrayFlat)

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    749         # Remove full_output from kwargs, otherwise we're passing it in twice.
    750         return_full = kwargs.pop('full_output', False)
--> 751         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    752         popt, pcov, infodict, errmsg, ier = res
    753         cost = np.sum(infodict['fvec'] ** 2)

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    384     m = shape[0]
    385     if n > m:
--> 386         raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
    387     if epsfcn is None:
    388         epsfcn = finfo(dtype).eps

TypeError: Improper input: N=5 must not exceed M=2

标签： arraysscipycurve-fittingshapesscipy-optimize

Answer 1

这应该可以解决问题。查看该for_fitting函数，了解如何以curve_fit可接受的方式打包所有内容。

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit


def doughnut(y, x, y0, x0, A, FWHM):
    """2D intensity distribution function of doughnut beams (DOI: 10.1126/science.aak9913,
    https://science.sciencemag.org/content/sci/suppl/2016/12/21/science.aak9913.DC1/Balzarotti_SM.pdf).

    Parameters
    ----------
    y, x : float
        X and Y coordinates, orthogonal to beam axis
    y0 : float
        Y offset
    x0 : float
        X offset
    A : float
        Peak intensity
    FWHM : float
        Full width at half maximum
    """
    return (
        A
        * np.e
        * 4
        * np.log(2)
        * ((x + x0) ** 2 + (y + y0) ** 2)
        / FWHM ** 2
        * np.exp(-4 * np.log(2) * ((x + x0) ** 2 + (y + y0) ** 2) / FWHM ** 2)
    )


fig0, (ax0, ax1, ax2) = plt.subplots(1, 3, sharex=True, sharey=True)

# Make data
X = np.arange(-10, 10, 0.25)
Y = np.arange(-10, 10, 0.25)
X, Y = np.meshgrid(X, Y)

true_params = (0, 0, 100, 7)

Z = doughnut(Y, X, *true_params)

# Plot the surface
ax0.matshow(Z, extent=(-10, 10, 10, -10))
ax0.set_title("Ground Truth")


def for_fitting(xdata, y0, x0, A, FWHM):
    yy, xx = xdata
    return doughnut(yy, xx, y0, x0, A, FWHM).ravel()


noisy_data = np.random.poisson(Z) + np.random.randn(*Z.shape)

ax1.matshow(noisy_data, extent=(-10, 10, 10, -10))
ax1.set_title("Noisy Data")

opt_params, cov = curve_fit(for_fitting, (Y, X), noisy_data.ravel(), p0=(0, 0, 10, 1))

print(opt_params)

fit_Z = doughnut(Y, X, *opt_params)

ax2.matshow(fit_Z, extent=(-10, 10, 10, -10))
ax2.set_title("Fit")

fig0.tight_layout()

plt.show()