3

我正在尝试使用 Scipyleastsq为二维中的一组测量点坐标找到“方形”网格的最佳拟合（实验点大约在方形网格上）。

网格的参数是间距（等于 x 和 y）、中心位置 (center_x和center_y) 和rotation(度数)。

我定义了一个误差函数，计算每对点的欧几里得距离（实验网格与理想网格）并取平均值。我想彻底减少这个功能leastsq，但我得到一个错误。

以下是函数定义：

import numpy as np
from scipy.optimize import leastsq

def get_spot_grid(shape, pitch, center_x, center_y, rotation=0):
    x_spots, y_spots = np.meshgrid(
             (np.arange(shape[1]) - (shape[1]-1)/2.)*pitch, 
             (np.arange(shape[0]) - (shape[0]-1)/2.)*pitch)
    theta = rotation/180.*np.pi
    x_spots = x_spots*np.cos(theta) - y_spots*np.sin(theta) + center_x
    y_spads = x_spots*np.sin(theta) + y_spots*np.cos(theta) + center_y
    return x_spots, y_spots

def get_mean_distance(x1, y1, x2, y2):
    return np.sqrt((x1 - x2)**2 + (y1 - y2)**2).mean()

def err_func(params, xe, ye):
    pitch, center_x, center_y, rotation = params
    x_grid, y_grid = get_spot_grid(xe.shape, pitch, center_x, center_y, rotation)
    return get_mean_distance(x_grid, y_grid, xe, ye)

这是实验坐标：

xe = np.array([ -23.31,  -4.01,  15.44,  34.71, -23.39,  -4.10,  15.28,  34.60, -23.75,  -4.38,  15.07,  34.34, -23.91,  -4.53,  14.82,  34.15]).reshape(4, 4)
ye = np.array([-16.00, -15.81, -15.72, -15.49,   3.29,   3.51,   3.90,   4.02,  22.75,  22.93,  23.18,  23.43,  42.19,  42.35,  42.69,  42.87]).reshape(4, 4)

我尝试以leastsq这种方式使用：

leastsq(err_func, x0=(19, 12, 5, 0), args=(xe, ye))

但我收到以下错误：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-19-ee91cf6ce7d6> in <module>()
----> 1 leastsq(err_func, x0=(19, 12, 5, 0), args=(xe, ye))

C:\Anaconda\lib\site-packages\scipy\optimize\minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    369     m = shape[0]
    370     if n > m:
--> 371         raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
    372     if epsfcn is None:
    373         epsfcn = finfo(dtype).eps

TypeError: Improper input: N=4 must not exceed M=1

我不知道这里有什么问题:(

4

3 回答 3

1

由于 minimumsq 函数假定 err_function 返回一个残差数组，并且以这种方式编写 err_function 有点困难，为什么不使用另一个 scipy 的函数 - minimize。然后你添加你的指标 - 你已经拥有的错误函数并且它有效。但是，我认为 get_spot_grid 函数中还有一个错字（y_spots vs y_spads）。完整代码：

import numpy as np
from scipy.optimize import leastsq, minimize

def get_spot_grid(shape, pitch, center_x, center_y, rotation=0):
    x_spots, y_spots = np.meshgrid(
             (np.arange(shape[1]) - (shape[1]-1)/2.)*pitch, 
             (np.arange(shape[0]) - (shape[0]-1)/2.)*pitch)
    theta = rotation/180.*np.pi
    x_spots = x_spots*np.cos(theta) - y_spots*np.sin(theta) + center_x
    y_spots = x_spots*np.sin(theta) + y_spots*np.cos(theta) + center_y
    return x_spots, y_spots


def get_mean_distance(x1, y1, x2, y2):
    return np.sqrt((x1 - x2)**2 + (y1 - y2)**2).mean()


def err_func(params, xe, ye):
    pitch, center_x, center_y, rotation = params
    x_grid, y_grid = get_spot_grid(xe.shape, pitch, center_x, center_y, rotation)
    return get_mean_distance(x_grid, y_grid, xe, ye)

xe = np.array([-23.31,  -4.01,  15.44,  34.71, -23.39,  -4.10,  15.28,  34.60, -23.75,  -4.38,  15.07,  34.34, -23.91,  -4.53,  14.82,  34.15]).reshape(4, 4)
ye = np.array([-16.00, -15.81, -15.72, -15.49,   3.29,   3.51,   3.90,   4.02,  22.75,  22.93,  23.18,  23.43,  42.19,  42.35,  42.69,  42.87]).reshape(4, 4)

# leastsq(err_func, x0=(19, 12, 5, 0), args=(xe, ye))
minimize(err_func, x0=(19, 12, 5, 0), args=(xe, ye))
于 2014-02-06T08:24:39.057 回答
0

xe传递给 leastsq 的函数（例如 err_func）应该返回一个与and具有相同形状的值的数组ye——也就是说，对于xeand的每个值都有一个残差ye。

def err_func(params, xe, ye):
    pitch, center_x, center_y, rotation = params
    x_grid, y_grid = get_spot_grid(xe.shape, pitch, center_x, center_y, rotation)
    return get_mean_distance(x_grid, y_grid, xe, ye)

调用mean()inget_mean_distance将返回值减少为单个标量。请记住，传递给的xe和ye是err_func数组而不是标量。

错误信息

TypeError: Improper input: N=4 must not exceed M=1

是说参数的数量 4 不应超过err_func1 返回的残差数量。

mean()通过将调用更改为mean(axis=0)（即取每列的平均值）或mean(axis=1)（即取每一行的平均值） ，可以使程序可运行：

def get_mean_distance(x1, y1, x2, y2):
    return np.sqrt((x1 - x2)**2 + (y1 - y2)**2).mean(axis=1)

我不太了解您的代码，无法知道它应该是什么。xe但想法是和中的每个“点”都应该有一个值ye。

于 2014-02-06T01:45:38.943 回答
0

您的问题来自这样一个事实 hat leastsq 将“列”作为误差函数而不是二维矩阵数组。您可以使用 np.ravel() 将您想要的任何“图像”转换为平面一维数组，然后可以轻松安装。

例如拟合二维高斯：

#define gaussian function, p is parameters [wx,wy,x,y,I,offset]
def specGaussian2D(Xv,Yv,width_x, width_y, CenterX, CenterY, height=1.0, offset=0.0):
    X= (Xv - CenterX)/width_x
    Y= (Yv - CenterY)/width_y
    eX= np.exp(-0.5*(X**2))
    eY= np.exp(-0.5*(Y**2))
    eY=eY.reshape(len(eY),1)
    return offset + height*eY*eX

#define gaussian fit, use gaussian function specGaussian2D, p0 is initial parameters [wx,wy,x,y,I,offset]
def Gaussfit2D(Image,p0):
    sh=Image.shape
    Xv=np.arange(0,sh[1])
    Yv=np.arange(0,sh[0])
    errorfunction = lambda p: np.ravel(specGaussian2D(Xv,Yv,*p) -Image)
    p = optimize.leastsq(errorfunction, p0)
    return p

因此，在您的情况下，我将删除均值（否则您将拟合数据的“直方图”）并在您的 err_func 中使用 np.ravel。

于 2018-07-11T15:39:50.903 回答