ValueError：形状2,210和2,210不对齐：210dim 1！= 2dim 0

问题描述

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

# 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,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,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,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,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,n)
    384     m = shape[0]
    385     if n > m:

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py in _check_func(checker,argname,thefunc,numinputs,output_shape)
     25 def _check_func(checker,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,FWHM)
     17         Full width at half maximum
     18     '''
---> 19     return A*np.exp(1)*4*np.log(2)*(np.dot(x+x0,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,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,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,Z)
plt.show()

=> Plot

但是：现在，在尝试拟合数据时出现了一个新错误：

---------------------------------------------------------------------------
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,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

解决方法

这应该可以解决问题。看一下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,0.25)
X,Y = np.meshgrid(X,Y)

true_params = (0,100,7)

Z = doughnut(Y,X,*true_params)

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


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


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

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

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

print(opt_params)

fit_Z = doughnut(Y,*opt_params)

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

fig0.tight_layout()

plt.show()