forked from ProfessorBrunner/sgl-model-detect
-
Notifications
You must be signed in to change notification settings - Fork 0
/
fittest.py
73 lines (54 loc) · 1.89 KB
/
fittest.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import scipy.optimize as opt
import numpy as np
import pylab as plt
from scipy import misc
from scipy import misc
import math
'''
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Fitting for two dimentional gaussians and one dimentional
sersics.
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
'''
def twoD_Gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
xo = float(xo)
yo = float(yo)
a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
g = offset + amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo)
+ c*((y-yo)**2)))
return g.ravel()
def sersic((x,y),kx,nx,x0,ky,ny,y0,s):
S = s * math.e ** (-1.0 * kx * ( (x - x0) ** (1.0/nx) ) ) * math.e ** (-1.0 * ky * ( (y - y0) ** (1.0/ny) ) )
return S.ravel()
z = misc.imread('arc.JPG',True)
m , n = z.shape
x, y = np.mgrid[:m, :n]
Max = 0
Max_x = 0
Max_y = 0
Val_sum = 0
for i in range (0,m-1):
for j in range (0,n-1):
Val_sum += z[i][j] #Find Sum to calculate Mean
if(z[i][j]>Max):
Max = z[i][j]
Max_x = i
Max_y = j
Mean = Val_sum / (float)( m * n )
# add some noise to the data and try to fit the data generated beforehand
#initial_guess = (Max,Max_x,Max_y,1,1,0,0)
initial_guess = (100,4,Max_x,100,4,Max_y,0)
#popt, pcov = opt.curve_fit(twoD_Gaussian, (x, y), z.ravel(), p0=initial_guess)
#data_fitted = twoD_Gaussian((x, y), *popt)
popt, pcov = opt.curve_fit(sersic, (x, y), z.ravel(), p0=initial_guess)
data_fitted = sersic((x, y), *popt)
plt.figure()
#plt.imshow(z, cmap=plt.cm.jet, extent=(x.min(), x.max(), y.min(), y.max()))
plt.imshow(z)
plt.figure()
plt.imshow(data_fitted.reshape(m, n))
plt.figure()
plt.imshow(z-data_fitted.reshape(m, n))
plt.show()