### Scipy最小平方:將正方形網格擬合到2D的實驗點上。

#### [英]Scipy leastsq: fitting a square grid to experimental points in 2D

I'm trying to use Scipy `leastsq` to find the best fit of a "square" grid for a set of measured points coordinates in 2-D (the experimental points are approximately on a square grid).

The parameters of the grid are pitch (equal for x and y), the center position (`center_x` and `center_y`) and `rotation` (in degree).

I defined an error function calculating the euclidean distance for each pairs of points (experimental vs ideal grid) and taking the mean. I want to minimize this function thorugh `leastsq` but I get an error.

Here are the function definitions:

``````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) - (shape-1)/2.)*pitch,
(np.arange(shape) - (shape-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)
``````

This are the experimental coordinates:

``````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)
``````

I try to use `leastsq` in this way:

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

but I get the following error:

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

I can't figure out what's the problem here :(

## 2 个解决方案

### #1

1

Since the leastsq function assumes that the err_function return an array of residuals docs and it is a little difficult to write the err_function in this manner why not use another scipy's function - minimize. Then you add your metric - the error function you already have and it works. However, I think there is one more typo in get_spot_grid function (y_spots vs y_spads). The complete code:

``````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) - (shape-1)/2.)*pitch,
(np.arange(shape) - (shape-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))
``````

### #2

0

The function passed to leastsq (e.g. err_func) should return an array of values of the same shape as `xe`and `ye` -- that is, one residual for each value of `xe` and `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)
``````

The call to `mean()` in `get_mean_distance` is reducing the return value to a single scalar. Keep in mind that the `xe` and `ye` passed to `err_func` are arrays not scalars.

get_mean_distance中的mean()調用將返回值降低到單個標量。請記住，xe和ye傳遞到err_func是數組而不是標量。

The error message

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

is saying the number of the parameters, 4, should not exceed the number of residuals returned by `err_func`, 1.

The program can be made runnable by changing the call to `mean()` to `mean(axis=0)` (i.e. take the mean of each column) or `mean(axis=1)` (i.e. take the mean of each row):

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

I don't really understand your code well enough to know which it should be. But the idea is that there should be one value for each "point" in `xe` and `ye`.