Created
July 31, 2013 15:00
-
-
Save andrewgiessel/6122739 to your computer and use it in GitHub Desktop.
fit 2d gaussian with numpy and scipy, including rotation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
def gaussian(self, height, center_x, center_y, width_x, width_y, rotation): | |
"""Returns a gaussian function with the given parameters""" | |
width_x = float(width_x) | |
width_y = float(width_y) | |
rotation = np.deg2rad(rotation) | |
center_x = center_x * np.cos(rotation) - center_y * np.sin(rotation) | |
center_y = center_x * np.sin(rotation) + center_y * np.cos(rotation) | |
def rotgauss(x,y): | |
xp = x * np.cos(rotation) - y * np.sin(rotation) | |
yp = x * np.sin(rotation) + y * np.cos(rotation) | |
g = height*np.exp( | |
-(((center_x-xp)/width_x)**2+ | |
((center_y-yp)/width_y)**2)/2.) | |
return g | |
return rotgauss | |
def moments(self, data): | |
"""Returns (height, x, y, width_x, width_y) | |
the gaussian parameters of a 2D distribution by calculating its | |
moments """ | |
total = data.sum() | |
X, Y = np.indices(data.shape) | |
x = (X*data).sum()/total | |
y = (Y*data).sum()/total | |
col = data[:, int(y)] | |
width_x = np.sqrt(abs((np.arange(col.size)-y)**2*col).sum()/col.sum()) | |
row = data[int(x), :] | |
width_y = np.sqrt(abs((np.arange(row.size)-x)**2*row).sum()/row.sum()) | |
height = data.max() | |
return height, x, y, width_x, width_y, 0.0 | |
def fitgaussian(self, data): | |
"""Returns (height, x, y, width_x, width_y) | |
the gaussian parameters of a 2D distribution found by a fit""" | |
params = self.moments(data) | |
errorfunction = lambda p: np.ravel(self.gaussian(*p)(*np.indices(data.shape)) - data) | |
p, success = scipy.optimize.leastsq(errorfunction, params) | |
return p | |
There's a mistake on line 8, center_y' needs to be calculated with the old value of
center_x' but as it is it uses the rotated value
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Yes, 0.0 is the rotation parameter which is just passed into the gaussian function. However, it is then adjusted when called for a fit where p returns all the params of the function - height, x, y, width_x, width_y, rotation. You would then know the best parameters to fit the function so 0 is not always the value assigned to rotation I believe...