forked from bblonder/flir_thermal_control
-
Notifications
You must be signed in to change notification settings - Fork 0
/
RadiometricData.py
executable file
·215 lines (165 loc) · 5.07 KB
/
RadiometricData.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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
import math
import numpy
ASY_SAFEGUARD = 1.0002
EXP_SAFEGUARD = 709.78
H2O_K1 = +1.5587e+0
H2O_K2 = +6.9390e-2
H2O_K3 = -2.7816e-4
H2O_K4 = +6.8455e-7
TAO_TATM_MIN = -30.0
TAO_TATM_MAX = 90.0
TAO_SQRTH2OMAX = 6.2365
TAO_COMP_MIN = 0.400
TAO_COMP_MAX = 1.000
class RadiometricData:
"""TODO: What does this do?"""
lPixval = None # the radiometric data
Tkelvin = None # the output temperature data
m_RelHum = 0.5
m_AtmTemp = 293.15
m_ObjectDistance = 1
m_X = 1.9
m_alpha1 = 0.006569
m_beta1 = -0.002276
m_beta2 = -0.00667
m_alpha2 = 0.01262
m_AtmTao = None
m_Emissivity = 0.95
m_ExtOptTransm = 1
m_K1 = None
m_K2 = None
m_ExtOptTemp = 293.15
m_AmbTemp = 293.15
m_J0 = 4214
m_J1 = 69.6245
m_R = 16671.9043
m_F = 1
m_B = 1430.1
def __init__(self):
pass
def doCalcAtmTao(self):
#double tao, dtao
#double H, T, sqrtD, X, a1, b1, a2, b2
#double sqrtH2O
#double TT
#double a1b1sqH2O, a2b2sqH2O, exp1, exp2
H = self.m_RelHum
C = self.m_AtmTemp
T = (C - 273.15)
#print("t atm = %.2f" % T)
#T = C.Value() # We need Celsius to use constants defined above
sqrtD = math.sqrt(self.m_ObjectDistance)
X = self.m_X
a1 = self.m_alpha1
b1 = self.m_beta1
a2 = self.m_alpha2
b2 = self.m_beta2
if (T < TAO_TATM_MIN):
T = TAO_TATM_MIN
elif (T > TAO_TATM_MAX):
T = TAO_TATM_MAX
TT = T*T
sqrtH2O = math.sqrt(H*math.exp(H2O_K1 + H2O_K2*T + H2O_K3*TT + H2O_K4*TT*T))
if (sqrtH2O > TAO_SQRTH2OMAX):
sqrtH2O = TAO_SQRTH2OMAX
a1b1sqH2O = (a1+b1*sqrtH2O)
a2b2sqH2O = (a2+b2*sqrtH2O)
exp1 = math.exp(-sqrtD*a1b1sqH2O)
exp2 = math.exp(-sqrtD*a2b2sqH2O)
tao = X*exp1 + (1-X)*exp2
dtao = -(a1b1sqH2O*X*exp1+a2b2sqH2O*(1-X)*exp2)
# The real D-derivative is also divided by 2 and sqrtD.
# Here we only want the sign of the slope!
if (tao < TAO_COMP_MIN):
tao = TAO_COMP_MIN # below min value, clip
elif (tao > TAO_COMP_MAX):
# check tao at 1_000_000 m dist
tao = X*math.exp(-(1.0E3)*a1b1sqH2O)+(1.0-X)*math.exp(-(1.0E3)*a2b2sqH2O)
# above max, staying up, assume \/-shape
if (tao > 1.0):
tao = TAO_COMP_MIN
else:
tao = TAO_COMP_MAX # above max, going down, assume /\-shape
elif ( dtao > 0.0 and self.m_ObjectDistance > 0.0):
tao = TAO_COMP_MIN # beween max & min, going up, assume \/
# else between max & min, going down => OK as it is, -)
return( tao)
def doCalcK1(self):
#dblVal = 1.0
dblVal = self.m_AtmTao * self.m_Emissivity * self.m_ExtOptTransm
if (dblVal > 0.0):
dblVal = 1/dblVal
return (dblVal)
def doCalcK2(self, dAmbObjSig, dAtmObjSig, dExtOptTempObjSig):
#double emi
temp1 = 0.0
temp2 = 0.0
temp3 = 0.0
emi = self.m_Emissivity
if (emi > 0.0):
temp1 = (1.0 - emi)/emi * dAmbObjSig
if (self.m_AtmTao > 0.0):
temp2 = (1.0 - self.m_AtmTao)/(emi*self.m_AtmTao)* dAtmObjSig
if (self.m_ExtOptTransm > 0.0 and self.m_ExtOptTransm < 1.0):
temp3 = (1.0 - self.m_ExtOptTransm) / (emi*self.m_AtmTao*self.m_ExtOptTransm)* dExtOptTempObjSig
return (temp1 + temp2 + temp3)
def tempToObjSig(self, dblKelvin):
objSign = 0.0
#objSign = R / (exp(B/T) - F)
dbl_reg = dblKelvin
if (dbl_reg > 0.0):
dbl_reg = self.m_B / dbl_reg
if (dbl_reg < EXP_SAFEGUARD):
dbl_reg = math.exp(dbl_reg)
if (self.m_F <= 1.0):
if ( dbl_reg < ASY_SAFEGUARD ):
dbl_reg = ASY_SAFEGUARD # Don't get above a R/(1-F) (horizontal) asymptote
else:
# F > 1.0
if ( dbl_reg < self.m_F*ASY_SAFEGUARD ):
dbl_reg = self.m_F*ASY_SAFEGUARD
# Don't get too close to a B/ln(F) (vertical) asymptote
objSign = self.m_R/(dbl_reg - self.m_F)
return(objSign)
def doUpdateCalcConst(self):
self.m_AtmTao = self.doCalcAtmTao()
self.m_K1 = self.doCalcK1()
self.m_K2 = self.doCalcK2(self.tempToObjSig(self.m_AmbTemp),self.tempToObjSig(self.m_AtmTemp),self.tempToObjSig(self.m_ExtOptTemp))
def imgToPow(self, lPixval):
pow = (lPixval - self.m_J0) / self.m_J1
return (pow)
def powToObjSig(self, dPow):
return (self.m_K1 * dPow - self.m_K2)
def objSigToTemp(self, dObjSig):
Tkelvin = 0.0
# Tkelvin = B /log(R / objSign + F)
#print(dObjSig)
dbl_reg = self.m_R / dObjSig + self.m_F
if (self.m_F <= 1.0):
if (dbl_reg < ASY_SAFEGUARD):
dbl_reg = ASY_SAFEGUARD # Don't get above a R/(1-F) (horizontal) asymptote
else:
tmp = self.m_F * ASY_SAFEGUARD
if (dbl_reg < tmp):
dbl_reg = tmp
# Don't get too close to a B/ln(F) (vertical) asymptote
Tkelvin = self.m_B / math.log(dbl_reg) # changed from quicklog to log
return (Tkelvin)
def getTempFast(self):
self.doUpdateCalcConst()
dPow = (self.lPixval - self.m_J0) / self.m_J1
dSig = self.m_K1 * dPow - self.m_K2
dbl_reg = self.m_R / dSig + self.m_F
if (self.m_F <= 1.0):
dbl_reg[dbl_reg < ASY_SAFEGUARD] = ASY_SAFEGUARD
else:
tmp = m_F * ASY_SAFEGUARD
dbl_reg[dbl_reg < tmp] = tmp
self.Tkelvin = self.m_B / numpy.log(dbl_reg)
return(self.Tkelvin)
d = RadiometricData()
d.lPixval = numpy.ones((640, 480)) * 14000
#print(d.getTemp())
print(d.getTempFast())
d.lPixval = numpy.ones((640, 480)) * 13000
print(d.getTempFast())