AllenInstitute / AllenSDK

code for reading and processing Allen Institute for Brain Science data
https://allensdk.readthedocs.io/en/latest/
Other
333 stars 149 forks source link

Fix gaussian_function_2d flipped axes #2678

Closed mikejhuang closed 1 year ago

mikejhuang commented 1 year ago

Summary

https://github.com/AllenInstitute/AllenSDK/issues/2655

_gaussian_function_2d expects an input of (x, y), while the fitting function is assuming _gaussian_function_2d takes an input of (y, x).

This causes the x, y values to be flipped for the respected x, y parameters (e.g. center_y, center_x, width_y, width_x)

Micha opened a ticket to report this issue. He suggested a fix to flip the input of the estimated gaussian in the error function after it takes in the incorrect (y, x) indices.

A simpler fix would be to correct the _gaussian_function_2d to take in (y, x) input.

Test result of fix:

Before the fix:

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

# 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)')

image

After the fix and updating rf_fitted to abide by the input specification for _gaussian_function_2d

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

# 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)(y, x)

# 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)')

image