X,Y coordinate flip in receptive_field_mapping.fit_2d_gaussian

[michaheilbron]

Bug description Hi! In allensdk.brain_observatory.ecephys.stimulus_analysis.receptive_field_mapping, the function fit_2d_gaussian has a bug resulting in the (x,y) coordinates of the returned parameters being flipped. Specifically, the function returns 4 x/y parameters (center_y, center_x, width_y, width_x), and in each case, x,y are the wrong way around.

Bug demonstration The bug is very apparent when we fit a gaussian to a fake highly non-uniform, non-centered RF. MWE:

import numpy as np
import matplotlib.pyplot as plt
from allensdk.brain_observatory.ecephys.stimulus_analysis.receptive_field_mapping import (
    fit_2d_gaussian)

# Simulate fake highly non-uniform RF 
# (in non-uniform dimensions, to ensure not to confuse x,y)
rf_raw=np.zeros((9,11))
rf_raw[2:8,7]=10
rf_raw[2:8,6]=3
rf_raw[2:8,8]=3

# fit, and apply fitted function 
(peak_height, cy,cx,wy,wx),success=fit_2d_gaussian(rf_raw)
y,x=np.indices((rf_raw.shape)) # rows, cols --> y,x 
rf_fitted=_gaussian_function_2d(peak_height, cy,cx,wy,wx)(x,y)

# compare 
fig,axes=plt.subplots(1,2,figsize=(4,2),sharey=True,sharex=True)
axes[0].imshow(rf_raw)
axes[0].set_title('original RF')
axes[1].imshow(rf_fitted)
axes[1].set_title('fitted RF \n (current function)')

As you can see the resulting fit is all wrong:

Bug origin and how to fix The origin of the bug is in line 421 (link). This line currently reads:

errorfunction = lambda p: np.ravel(_gaussian_function_2d(*p)(*np.indices(matrix.shape)) - matrix)

... but it should be:

errorfunction = lambda p: np.ravel(_gaussian_function_2d(*p)(*np.flip(np.indices(matrix.shape),0)) - matrix)

... because np.indices(matrix.shape) returns rows,columns i.e. y,x. However, _gaussian_function_2d(*p) expects x,y. Hence, we need to add the np.flip(...,0) to make it correct.

Demonstration of the full fixed version:

# import what we need 
from scipy.optimize import curve_fit, leastsq
from allensdk.brain_observatory.ecephys.stimulus_analysis.receptive_field_mapping import (
    fit_2d_gaussian,gaussian_moments_2d,_gaussian_function_2d)

# define fixed function 
def fit_2d_gaussian_fixed(matrix):
    """
    bug-fixed version of fit_2d_gaussian from allen-institute  

    Fits a receptive field with a 2-dimensional Gaussian distribution
    Parameters
    ----------
    matrix : numpy.ndarray
        2D matrix of spike counts
    Returns
    -------
    parameters - tuple
        peak_height : peak of distribution
        center_y : y-coordinate of distribution center
        center_x : x-coordinate of distribution center
        width_y : width of distribution along x-axis
        width_x : width of distribution along y-axis
    success - bool
        True if a fit was found, False otherwise
    """

    params = gaussian_moments_2d(matrix)
    if params is None:
        return (np.nan, np.nan, np.nan, np.nan, np.nan), False

    errorfunction = lambda p: np.ravel(_gaussian_function_2d(*p)(*np.flip(np.indices(matrix.shape),0)) - matrix)
    fit_params, ier = leastsq(errorfunction, params)
    success = True if ier < 5 else False

    return fit_params, success

# Simulate fake highly non-uniform RF 
# (in non-uniform dimensions, to ensure not to confuse x,y)
rf_raw=np.zeros((9,11))
rf_raw[2:8,7]=10
rf_raw[2:8,6]=3
rf_raw[2:8,8]=3
y,x=np.indices((rf_raw.shape)) # rows, cols --> y,x 

# apply fixed function 
(peak_height, cy,cx,wy,wx),success= fit_2d_gaussian_fixed(rf_raw)
rf_fitted_fix=_gaussian_function_2d(peak_height, cy_fix,cx_fix,wy_fix, wx_fix)(x,y)

# compare results 
fig,axes=plt.subplots(1,3,figsize=(6,2),sharey=True,sharex=True)
axes[0].imshow(rf_raw)
axes[0].set_title('original RF')

axes[1].imshow(rf_fitted)
axes[1].set_title('fitted RF \n (current function)')

axes[2]
plt.imshow(rf_fitted_fix)
axes[2].set_title('fitted RF \n (fixed function)');

... which results in:

Environment

• AllenSDK version 2.13.4

Do you want to work on this issue? I gave the instructions to fix above.