Skip to content

Commit 549ef01

Browse files
steinnymirsteinnymir
committed
gaussian fit
1 parent 5337cda commit 549ef01

File tree

4 files changed

+357
-0
lines changed

4 files changed

+357
-0
lines changed

‎Math/function.py

+5
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,11 @@ def gaussian(x, mu, sig):
1919
"""
2020
return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
2121

22+
def lorentz(p,x):
23+
return p[2] / (1.0 + (x / p[0] - 1.0)**4 * p[1]**2)
24+
25+
def errorfunc(p, x, z):
26+
return lorentz(p, x) - z
2227

2328
if __name__ == '__main__':
2429
main()

‎Processing/gauss01.png

4.25 KB
Loading

‎Processing/gaussfitter.py

+178
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,178 @@
1+
import numpy as np
2+
from scipy import optimize, io as spio, stats, misc
3+
import matplotlib.pyplot as plt
4+
from datetime import datetime
5+
6+
def main():
7+
filepath = 'E:/Software/StrahlprofiL2.1/TestProfil'
8+
data_dict = spio.loadmat(filepath)
9+
data = data_dict['Profil']
10+
# gaussFig = misc.imread('gauss01.png')
11+
# print(gaussFig.shape)
12+
# data = gaussFig.sum(axis=2)
13+
t0 = datetime.now()
14+
15+
init_p = (height, amplitude, x, y, width_x, width_y, rota) = (1,1,200,350,20,20,0)
16+
params = gaussfit(data,return_all=1)
17+
fit = twodgaussian(params[0], 0,1,1)
18+
plt.imshow(data, cmap=plt.cm.gist_earth_r)
19+
plt.contour(fit(*np.indices(data.shape)), cmap=plt.cm.copper)
20+
21+
plt.show()
22+
23+
print(params)
24+
print('it took {} seconds'.format(datetime.now()-t0))
25+
26+
def moments_(data,circle,rotate,vheight):
27+
"""Returns (height, x, y, width_x, width_y)
28+
the gaussian parameters of a 2D distribution by calculating its
29+
moments """
30+
total = data.sum()
31+
X, Y = np.indices(data.shape)
32+
x = (X * data).sum() / total
33+
y = (Y * data).sum() / total
34+
col = data[:, int(y)]
35+
width_x = np.sqrt(np.abs((np.arange(col.size) - y) ** 2 * col).sum() / col.sum())
36+
row = data[int(x), :]
37+
width_y = np.sqrt(np.abs((np.arange(row.size) - x) ** 2 * row).sum() / row.sum())
38+
height = data.max()
39+
return (height, x, y, width_x, width_y)
40+
41+
42+
def moments(data,circle,rotate,vheight):
43+
"""Returns (height, amplitude, x, y, width_x, width_y, rotation angle)
44+
the gaussian parameters of a 2D distribution by calculating its
45+
moments. Depending on the input parameters, will only output
46+
a subset of the above"""
47+
total = data.sum()
48+
X, Y = np.indices(data.shape)
49+
x = (X*data).sum()/total
50+
y = (Y*data).sum()/total
51+
col = data[:, int(y)]
52+
width_x = np.sqrt(abs((np.arange(col.size)-y)**2*col).sum()/col.sum())
53+
row = data[int(x), :]
54+
width_y = np.sqrt(abs((np.arange(row.size)-x)**2*row).sum()/row.sum())
55+
width = ( width_x + width_y ) / 2.
56+
# height = stats.mode(data.ravel())[0][0]
57+
# height = data.max()
58+
59+
amplitude = data.max()#-height
60+
mylist = [amplitude,x,y]
61+
if vheight==1:
62+
mylist = [height] + mylist
63+
if circle==0:
64+
mylist = mylist + [width_x,width_y]
65+
else:
66+
mylist = mylist + [width]
67+
if rotate==1:
68+
mylist = mylist + [0.] #rotation "moment" is just zero...
69+
return tuple(mylist)
70+
71+
def twodgaussian(inpars, circle, rotate, vheight):
72+
"""Returns a 2d gaussian function of the form:
73+
x' = cos(rota) * x - sin(rota) * y
74+
y' = sin(rota) * x + cos(rota) * y
75+
(rota should be in degrees)
76+
g = b + a exp ( - ( ((x-center_x)/width_x)**2 +
77+
((y-center_y)/width_y)**2 ) / 2 )
78+
79+
However, the above values are passed by list. The list should be:
80+
inpars = (height,amplitude,center_x,center_y,width_x,width_y,rota)
81+
82+
You can choose to ignore / neglect some of the above input parameters using the following options:
83+
circle=0 - default is an elliptical gaussian (different x, y widths), but can reduce
84+
the input by one parameter if it's a circular gaussian
85+
rotate=1 - default allows rotation of the gaussian ellipse. Can remove last parameter
86+
by setting rotate=0
87+
vheight=1 - default allows a variable height-above-zero, i.e. an additive constant
88+
for the Gaussian function. Can remove first parameter by setting this to 0
89+
"""
90+
inpars_old = inpars
91+
inpars = list(inpars)
92+
if vheight == 1:
93+
height = inpars.pop(0)
94+
height = float(height)
95+
else:
96+
height = float(0)
97+
amplitude, center_x, center_y = inpars.pop(0),inpars.pop(0),inpars.pop(0)
98+
amplitude = float(amplitude)
99+
center_x = float(center_x)
100+
center_y = float(center_y)
101+
if circle == 1:
102+
width = inpars.pop(0)
103+
width_x = float(width)
104+
width_y = float(width)
105+
else:
106+
width_x, width_y = inpars.pop(0),inpars.pop(0)
107+
width_x = float(width_x)
108+
width_y = float(width_y)
109+
if rotate == 1:
110+
rota = inpars.pop(0)
111+
rota = np.pi/180. * float(rota)
112+
rcen_x = center_x * np.cos(rota) - center_y * np.sin(rota)
113+
rcen_y = center_x * np.sin(rota) + center_y * np.cos(rota)
114+
else:
115+
rcen_x = center_x
116+
rcen_y = center_y
117+
if len(inpars) > 0:
118+
raise ValueError("There are still input parameters:" + str(inpars) + \
119+
" and you've input: " + str(inpars_old) + " circle=%d, rotate=%d, vheight=%d" % (circle,rotate,vheight) )
120+
121+
def rotgauss(x,y):
122+
if rotate==1:
123+
xp = x * np.cos(rota) - y * np.sin(rota)
124+
yp = x * np.sin(rota) + y * np.cos(rota)
125+
else:
126+
xp = x
127+
yp = y
128+
g = height+amplitude*np.exp(
129+
-(((rcen_x-xp)/width_x)**2+
130+
((rcen_y-yp)/width_y)**2)/2.)
131+
return g
132+
return rotgauss
133+
134+
def gaussfit(data,err=None,params=[],autoderiv=1,return_all=0,circle=0,rotate=1,vheight=0):
135+
"""
136+
Gaussian fitter with the ability to fit a variety of different forms of 2-dimensional gaussian.
137+
138+
Input Parameters:
139+
data - 2-dimensional data array
140+
err=None - error array with same size as data array
141+
params=[] - initial input parameters for Gaussian function.
142+
(height, amplitude, x, y, width_x, width_y, rota)
143+
autoderiv=1 - use the autoderiv provided in the lmder.f function (the alternative
144+
is to us an analytic derivative with lmdif.f: this method is less robust)
145+
return_all=0 - Default is to return only the Gaussian parameters. See below for
146+
detail on output
147+
circle=0 - default is an elliptical gaussian (different x, y widths), but can reduce
148+
the input by one parameter if it's a circular gaussian
149+
rotate=1 - default allows rotation of the gaussian ellipse. Can remove last parameter
150+
by setting rotate=0
151+
vheight=1 - default allows a variable height-above-zero, i.e. an additive constant
152+
for the Gaussian function. Can remove first parameter by setting this to 0
153+
154+
Output:
155+
Default output is a set of Gaussian parameters with the same shape as the input parameters
156+
Can also output the covariance matrix, 'infodict' that contains a lot more detail about
157+
the fit (see scipy.optimize.leastsq), and a message from leastsq telling what the exit
158+
status of the fitting routine was
159+
"""
160+
if params == []:
161+
params = (moments(data,circle,rotate,vheight))
162+
if err == None:
163+
errorfunction = lambda p: np.ravel((twodgaussian(p,circle,rotate,vheight)(*np.indices(data.shape)) - data))
164+
else:
165+
errorfunction = lambda p: np.ravel((twodgaussian(p,circle,rotate,vheight)(*np.indices(data.shape)) - data)/err)
166+
if autoderiv == 0:
167+
raise ValueError("I'm sorry, I haven't implemented this feature yet.")
168+
else:
169+
p, cov, infodict, errmsg, success = optimize.leastsq(errorfunction, params, full_output=1)
170+
if return_all == 0:
171+
return p
172+
elif return_all == 1:
173+
return p,cov,infodict,errmsg
174+
175+
176+
if __name__=="__main__":
177+
178+
main()

‎Processing/gaussfitter2.py

+174
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
# gaussfitter.py
2+
# created by Adam Ginsburg ([email protected] or [email protected]) 3/17/08)
3+
import numpy as np
4+
from scipy import optimize, io as spio
5+
from scipy import stats
6+
import matplotlib.pyplot as plt
7+
8+
from datetime import datetime
9+
def main():
10+
filepath = 'E:/Software/StrahlprofiL2.1/TestProfil'
11+
data_dict = spio.loadmat(filepath)
12+
data = data_dict['Profil']
13+
# gaussFig = misc.imread('gauss01.png')
14+
# print(gaussFig.shape)
15+
# data = gaussFig.sum(axis=2)
16+
t0 = datetime.now()
17+
18+
init_p = (height, amplitude, x, y, width_x, width_y, rota) = (1,1,200,350,20,20,0)
19+
params = gaussfit(data,return_all=1)
20+
fit = twodgaussian(params[0], 0,1,1)
21+
plt.imshow(data, cmap=plt.cm.gist_earth_r)
22+
plt.contour(fit(*np.indices(data.shape)), cmap=plt.cm.copper)
23+
24+
plt.show()
25+
26+
print(params)
27+
print('it took {} seconds'.format(datetime.now()-t0))
28+
29+
30+
def moments(data,circle,rotate,vheight):
31+
"""Returns (height, amplitude, x, y, width_x, width_y, rotation angle)
32+
the gaussian parameters of a 2D distribution by calculating its
33+
moments. Depending on the input parameters, will only output
34+
a subset of the above"""
35+
total = data.sum()
36+
X, Y = np.indices(data.shape)
37+
x = (X*data).sum()/total
38+
y = (Y*data).sum()/total
39+
col = data[:, int(y)]
40+
width_x = np.sqrt(abs((np.arange(col.size)-y)**2*col).sum()/col.sum())
41+
row = data[int(x), :]
42+
width_y = np.sqrt(abs((np.arange(row.size)-x)**2*row).sum()/row.sum())
43+
width = ( width_x + width_y ) / 2.
44+
height = stats.mode(data.ravel())[0][0]
45+
amplitude = data.max()-height
46+
mylist = [amplitude,x,y]
47+
if vheight==1:
48+
mylist = [height] + mylist
49+
if circle==0:
50+
mylist = mylist + [width_x,width_y]
51+
else:
52+
mylist = mylist + [width]
53+
if rotate==1:
54+
mylist = mylist + [0.] #rotation "moment" is just zero...
55+
return tuple(mylist)
56+
57+
def twodgaussian(inpars, circle, rotate, vheight):
58+
"""Returns a 2d gaussian function of the form:
59+
x' = cos(rota) * x - sin(rota) * y
60+
y' = sin(rota) * x + cos(rota) * y
61+
(rota should be in degrees)
62+
g = b + a exp ( - ( ((x-center_x)/width_x)**2 +
63+
((y-center_y)/width_y)**2 ) / 2 )
64+
65+
where x and y are the input parameters of the returned function,
66+
and all other parameters are specified by this function
67+
68+
However, the above values are passed by list. The list should be:
69+
inpars = (height,amplitude,center_x,center_y,width_x,width_y,rota)
70+
71+
You can choose to ignore / neglect some of the above input parameters using the following options:
72+
circle=0 - default is an elliptical gaussian (different x, y widths), but can reduce
73+
the input by one parameter if it's a circular gaussian
74+
rotate=1 - default allows rotation of the gaussian ellipse. Can remove last parameter
75+
by setting rotate=0
76+
vheight=1 - default allows a variable height-above-zero, i.e. an additive constant
77+
for the Gaussian function. Can remove first parameter by setting this to 0
78+
"""
79+
inpars_old = inpars
80+
inpars = list(inpars)
81+
if vheight == 1:
82+
height = inpars.pop(0)
83+
height = float(height)
84+
else:
85+
height = float(0)
86+
amplitude, center_x, center_y = inpars.pop(0),inpars.pop(0),inpars.pop(0)
87+
amplitude = float(amplitude)
88+
center_x = float(center_x)
89+
center_y = float(center_y)
90+
if circle == 1:
91+
width = inpars.pop(0)
92+
width_x = float(width)
93+
width_y = float(width)
94+
else:
95+
width_x, width_y = inpars.pop(0),inpars.pop(0)
96+
width_x = float(width_x)
97+
width_y = float(width_y)
98+
if rotate == 1:
99+
rota = inpars.pop(0)
100+
rota = np.pi/180. * float(rota)
101+
rcen_x = center_x * np.cos(rota) - center_y * np.sin(rota)
102+
rcen_y = center_x * np.sin(rota) + center_y * np.cos(rota)
103+
else:
104+
rcen_x = center_x
105+
rcen_y = center_y
106+
if len(inpars) > 0:
107+
raise ValueError("There are still input parameters:" + str(inpars) + \
108+
" and you've input: " + str(inpars_old) + " circle=%d, rotate=%d, vheight=%d" % (circle,rotate,vheight) )
109+
110+
def rotgauss(x,y):
111+
if rotate==1:
112+
xp = x * np.cos(rota) - y * np.sin(rota)
113+
yp = x * np.sin(rota) + y *np. cos(rota)
114+
else:
115+
xp = x
116+
yp = y
117+
g = height+amplitude*np.exp(
118+
-(((rcen_x-xp)/width_x)**2+
119+
((rcen_y-yp)/width_y)**2)/2.)
120+
return g
121+
return rotgauss
122+
123+
def gaussfit(data,err=None,params=[],autoderiv=1,return_all=0,circle=0,rotate=1,vheight=1):
124+
"""
125+
Gaussian fitter with the ability to fit a variety of different forms of 2-dimensional gaussian.
126+
127+
Input Parameters:
128+
data - 2-dimensional data array
129+
err=None - error array with same size as data array
130+
params=[] - initial input parameters for Gaussian function.
131+
(height, amplitude, x, y, width_x, width_y, rota)
132+
if not input, these will be determined from the moments of the system,
133+
assuming no rotation
134+
autoderiv=1 - use the autoderiv provided in the lmder.f function (the alternative
135+
is to us an analytic derivative with lmdif.f: this method is less robust)
136+
return_all=0 - Default is to return only the Gaussian parameters. See below for
137+
detail on output
138+
circle=0 - default is an elliptical gaussian (different x, y widths), but can reduce
139+
the input by one parameter if it's a circular gaussian
140+
rotate=1 - default allows rotation of the gaussian ellipse. Can remove last parameter
141+
by setting rotate=0
142+
vheight=1 - default allows a variable height-above-zero, i.e. an additive constant
143+
for the Gaussian function. Can remove first parameter by setting this to 0
144+
145+
Output:
146+
Default output is a set of Gaussian parameters with the same shape as the input parameters
147+
Can also output the covariance matrix, 'infodict' that contains a lot more detail about
148+
the fit (see scipy.optimize.leastsq), and a message from leastsq telling what the exit
149+
status of the fitting routine was
150+
151+
Warning: Does NOT necessarily output a rotation angle between 0 and 360 degrees.
152+
"""
153+
if params == []:
154+
params = (moments(data,circle,rotate,vheight))
155+
if err == None:
156+
errorfunction = lambda p: np.ravel((twodgaussian(p,circle,rotate,vheight)(*np.indices(data.shape)) - data))
157+
else:
158+
errorfunction = lambda p: np.ravel((twodgaussian(p,circle,rotate,vheight)(*np.indices(data.shape)) - data)/err)
159+
if autoderiv == 0:
160+
# the analytic derivative, while not terribly difficult, is less efficient and useful. I only bothered
161+
# putting it here because I was instructed to do so for a class project - please ask if you would like
162+
# this feature implemented
163+
raise ValueError("I'm sorry, I haven't implemented this feature yet.")
164+
else:
165+
p, cov, infodict, errmsg, success = optimize.leastsq(errorfunction, params, full_output=1)
166+
if return_all == 0:
167+
return p
168+
elif return_all == 1:
169+
return p,cov,infodict,errmsg
170+
171+
172+
if __name__=="__main__":
173+
174+
main()

0 commit comments

Comments
 (0)