Scipy leastsq: fitting a square grid to experimental points in 2D

Question

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

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

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

Answer 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[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))

Answer 2

The function passed to leastsq (e.g. err_func) should return an array of values of the same shape as xeand 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.

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.

Answer 3

Your problem comes from the fact hat leastsq take a 'column' as errorfunction and not a 2D matrix array. You can transform any 'image' you want to a flat 1D array with np.ravel(), that can then easily be fitted.

for example for fitting a 2D gaussian:

#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

So in your case I would remove mean (or else you will fit a 'histogram' of your data) and use np.ravel in your err_func.