Scipy наименьший квадрат: подгонка квадратной сетки к экспериментальным точкам в 2D

Поскольку функция наименьшего квадрата предполагает, что функция err_function возвращает массив остатков docs, и немного сложно написать функцию err_function таким образом, почему бы не использовать другую функцию scipy - свернуть. Затем вы добавляете свою метрику — функция ошибок у вас уже есть, и она работает. Однако я думаю, что в функции 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))

Функция, переданная наименьшему квадрату (например, err_func), должна возвращать массив значений той же формы, что и xeи ye, то есть один остаток для каждого значения xe и 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() в get_mean_distance уменьшает возвращаемое значение до одного скаляра. Имейте в виду, что xe и ye, переданные в err_func, являются массивами, а не скалярами.

Сообщение об ошибке

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

говорит, что количество параметров, 4, не должно превышать количество остатков, возвращаемых err_func, 1.

Программу можно сделать работоспособной, изменив вызов с 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 должно быть одно значение.

Ваша проблема связана с тем, что наименьший квадрат принимает «столбец» в качестве функции ошибки, а не массив двумерной матрицы. Вы можете преобразовать любое «изображение», которое хотите, в плоский одномерный массив с помощью 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

Так что в вашем случае я бы удалил среднее значение (иначе вы будете соответствовать «гистограмме» ваших данных) и использовал np.ravel в своем err_func.