/
fitter.py
389 lines (312 loc) · 15.4 KB
/
fitter.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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from astropy import log
import math
import numpy as np
from .harmonics import (fit_first_and_second_harmonics,
first_and_second_harmonic_function)
from .isophote import Isophote, CentralPixel
from .sample import EllipseSample
__all__ = ['EllipseFitter']
__doctest_skip__ = ['EllipseFitter.fit']
PI2 = np.pi / 2
MAX_EPS = 0.95
MIN_EPS = 0.05
DEFAULT_CONVERGENCE = 0.05
DEFAULT_MINIT = 10
DEFAULT_MAXIT = 50
DEFAULT_FFLAG = 0.7
DEFAULT_MAXGERR = 0.5
class EllipseFitter(object):
"""
Class to fit ellipses.
Parameters
----------
sample : `~photutils.isophote.EllipseSample` instance
The sample data to be fitted.
"""
def __init__(self, sample):
self._sample = sample
def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
maxit=DEFAULT_MAXIT, fflag=DEFAULT_FFLAG, maxgerr=DEFAULT_MAXGERR,
going_inwards=False):
"""
Fit an elliptical isophote.
Parameters
----------
conver : float, optional
The main convergence criterion. Iterations stop when the
largest harmonic amplitude becomes smaller (in absolute
value) than ``conver`` times the harmonic fit rms. The
default is 0.05.
minit : int, optional
The minimum number of iterations to perform. A minimum of 10
(the default) iterations guarantees that, on average, 2
iterations will be available for fitting each independent
parameter (the four harmonic amplitudes and the intensity
level). For the first isophote, the minimum number of
iterations is 2 * ``minit`` to ensure that, even departing
from not-so-good initial values, the algorithm has a better
chance to converge to a sensible solution.
maxit : int, optional
The maximum number of iterations to perform. The default is
50.
fflag : float, optional
The acceptable fraction of flagged data points in the
sample. If the actual fraction of valid data points is
smaller than this, the iterations will stop and the current
`~photutils.isophote.Isophote` will be returned. Flagged
data points are points that either lie outside the image
frame, are masked, or were rejected by sigma-clipping. The
default is 0.7.
maxgerr : float, optional
The maximum acceptable relative error in the local radial
intensity gradient. This is the main control for preventing
ellipses to grow to regions of too low signal-to-noise
ratio. It specifies the maximum acceptable relative error
in the local radial intensity gradient. `Busko (1996; ASPC
101, 139)
<http://adsabs.harvard.edu/abs/1996ASPC..101..139B>`_ showed
that the fitting precision relates to that relative error.
The usual behavior of the gradient relative error is to
increase with semimajor axis, being larger in outer, fainter
regions of a galaxy image. In the current implementation,
the ``maxgerr`` criterion is triggered only when two
consecutive isophotes exceed the value specified by the
parameter. This prevents premature stopping caused by
contamination such as stars and HII regions.
A number of actions may happen when the gradient error
exceeds ``maxgerr`` (or becomes non-significant and is set
to `None`). If the maximum semimajor axis specified by
``maxsma`` is set to `None`, semimajor axis growth is
stopped and the algorithm proceeds inwards to the galaxy
center. If ``maxsma`` is set to some finite value, and this
value is larger than the current semimajor axis length, the
algorithm enters non-iterative mode and proceeds outwards
until reaching ``maxsma``. The default is 0.5.
going_inwards : bool, optional
Parameter to define the sense of SMA growth. When fitting
just one isophote, this parameter is used only by the code
that defines the details of how elliptical arc segments
("sectors") are extracted from the image, when using area
extraction modes (see the ``integrmode`` parameter in the
`~photutils.isophote.EllipseSample` class). The default is
`False`.
Returns
-------
result : `~photutils.isophote.Isophote` instance
The fitted isophote, which also contains fit status
information.
Examples
--------
>>> from photutils.isophote import EllipseSample, EllipseFitter
>>> sample = EllipseSample(data, sma=10.)
>>> fitter = EllipseFitter(sample)
>>> isophote = fitter.fit()
"""
sample = self._sample
# this flag signals that limiting gradient error (`maxgerr`)
# wasn't exceeded yet.
lexceed = False
# here we keep track of the sample that caused the minimum harmonic
# amplitude(in absolute value). This will eventually be used to
# build the resulting Isophote in cases where iterations run to
# the maximum allowed (maxit), or the maximum number of flagged
# data points (fflag) is reached.
minimum_amplitude_value = np.Inf
minimum_amplitude_sample = None
for iter in range(maxit):
# Force the sample to compute its gradient and associated values.
sample.update()
# The extract() method returns sampled values as a 2-d numpy array
# with the following structure:
# values[0] = 1-d array with angles
# values[1] = 1-d array with radii
# values[2] = 1-d array with intensity
values = sample.extract()
# Fit harmonic coefficients. Failure in fitting is
# a fatal error; terminate immediately with sample
# marked as invalid.
try:
coeffs = fit_first_and_second_harmonics(values[0], values[2])
except Exception as e:
log.info(e)
return Isophote(sample, iter+1, False, 3)
coeffs = coeffs[0]
# largest harmonic in absolute value drives the correction.
largest_harmonic_index = np.argmax(np.abs(coeffs[1:]))
largest_harmonic = coeffs[1:][largest_harmonic_index]
# see if the amplitude decreased; if yes, keep the
# corresponding sample for eventual later use.
if abs(largest_harmonic) < minimum_amplitude_value:
minimum_amplitude_value = abs(largest_harmonic)
minimum_amplitude_sample = sample
# check if converged
model = first_and_second_harmonic_function(values[0], coeffs)
residual = values[2] - model
if ((conver * sample.sector_area * np.std(residual))
> np.abs(largest_harmonic)):
# Got a valid solution. But before returning, ensure
# that a minimum of iterations has run.
if iter >= minit-1:
sample.update()
return Isophote(sample, iter+1, True, 0)
# it may not have converged yet, but the sample contains too
# many invalid data points: return.
if sample.actual_points < (sample.total_points * fflag):
# when too many data points were flagged, return the
# best fit sample instead of the current one.
minimum_amplitude_sample.update()
return Isophote(minimum_amplitude_sample, iter+1, True, 1)
# pick appropriate corrector code.
corrector = _correctors[largest_harmonic_index]
# generate *NEW* EllipseSample instance with corrected
# parameter. Note that this instance is still devoid of other
# information besides its geometry. It needs to be explicitly
# updated for computations to proceed. We have to build a new
# EllipseSample instance every time because of the lazy
# extraction process used by EllipseSample code. To minimize
# the number of calls to the area integrators, we pay a
# (hopefully smaller) price here, by having multiple calls to
# the EllipseSample constructor.
sample = corrector.correct(sample, largest_harmonic)
sample.update()
# see if any abnormal (or unusual) conditions warrant
# the change to non-iterative mode, or go-inwards mode.
proceed, lexceed = self._check_conditions(
sample, maxgerr, going_inwards, lexceed)
if not proceed:
sample.update()
return Isophote(sample, iter+1, True, -1)
# Got to the maximum number of iterations. Return with
# code 2, and handle it as a valid isophote. Use the
# best fit sample instead of the current one.
minimum_amplitude_sample.update()
return Isophote(minimum_amplitude_sample, maxit, True, 2)
def _check_conditions(self, sample, maxgerr, going_inwards, lexceed):
proceed = True
# If center wandered more than allowed, put it back
# in place and signal the end of iterative mode.
# if wander:
# if abs(dx) > WANDER(al)) or abs(dy) > WANDER(al):
# sample.geometry.x0 -= dx
# sample.geometry.y0 -= dy
# STOP(al) = ST_NONITERATE
# proceed = False
# check if an acceptable gradient value could be computed.
if sample.gradient_error:
if (not going_inwards and
(sample.gradient_relative_error > maxgerr or
sample.gradient >= 0.)):
if lexceed:
proceed = False
else:
lexceed = True
else:
proceed = False
# check if ellipse geometry diverged.
if abs(sample.geometry.eps > MAX_EPS):
proceed = False
if (sample.geometry.x0 < 1. or
sample.geometry.x0 > sample.image.shape[1] or
sample.geometry.y0 < 1. or
sample.geometry.y0 > sample.image.shape[0]):
proceed = False
# See if eps == 0 (round isophote) was crossed.
# If so, fix it but still proceed
if sample.geometry.eps < 0.:
sample.geometry.eps = min(-sample.geometry.eps, MAX_EPS)
if sample.geometry.pa < PI2:
sample.geometry.pa += PI2
else:
sample.geometry.pa -= PI2
# If ellipse is an exact circle, computations will diverge.
# Make it slightly flat, but still proceed
if sample.geometry.eps == 0.0:
sample.geometry.eps = MIN_EPS
return proceed, lexceed
class _ParameterCorrector(object):
def correct(self, sample, harmonic):
raise NotImplementedError
class _PositionCorrector(_ParameterCorrector):
def finalize_correction(self, dx, dy, sample):
new_x0 = sample.geometry.x0 + dx
new_y0 = sample.geometry.y0 + dy
return EllipseSample(sample.image, sample.geometry.sma, x0=new_x0,
y0=new_y0, astep=sample.geometry.astep,
sclip=sample.sclip, nclip=sample.nclip,
eps=sample.geometry.eps,
position_angle=sample.geometry.pa,
linear_growth=sample.geometry.linear_growth,
integrmode=sample.integrmode)
class _PositionCorrector_0(_PositionCorrector):
def correct(self, sample, harmonic):
aux = -harmonic * (1. - sample.geometry.eps) / sample.gradient
dx = -aux * math.sin(sample.geometry.pa)
dy = aux * math.cos(sample.geometry.pa)
return self.finalize_correction(dx, dy, sample)
class _PositionCorrector_1(_PositionCorrector):
def correct(self, sample, harmonic):
aux = -harmonic / sample.gradient
dx = aux * math.cos(sample.geometry.pa)
dy = aux * math.sin(sample.geometry.pa)
return self.finalize_correction(dx, dy, sample)
class _AngleCorrector(_ParameterCorrector):
def correct(self, sample, harmonic):
eps = sample.geometry.eps
sma = sample.geometry.sma
gradient = sample.gradient
correction = (harmonic * 2. * (1. - eps) / sma / gradient /
((1. - eps)**2 - 1.))
# '% np.pi' to make angle lie between 0 and np.pi radians
new_pa = (sample.geometry.pa + correction) % np.pi
return EllipseSample(sample.image, sample.geometry.sma,
x0=sample.geometry.x0, y0=sample.geometry.y0,
astep=sample.geometry.astep, sclip=sample.sclip,
nclip=sample.nclip, eps=sample.geometry.eps,
position_angle=new_pa,
linear_growth=sample.geometry.linear_growth,
integrmode=sample.integrmode)
class _EllipticityCorrector(_ParameterCorrector):
def correct(self, sample, harmonic):
eps = sample.geometry.eps
sma = sample.geometry.sma
gradient = sample.gradient
correction = harmonic * 2. * (1. - eps) / sma / gradient
new_eps = min((sample.geometry.eps - correction), MAX_EPS)
return EllipseSample(sample.image, sample.geometry.sma,
x0=sample.geometry.x0, y0=sample.geometry.y0,
astep=sample.geometry.astep, sclip=sample.sclip,
nclip=sample.nclip, eps=new_eps,
position_angle=sample.geometry.pa,
linear_growth=sample.geometry.linear_growth,
integrmode=sample.integrmode)
# instances of corrector code live here:
_correctors = [_PositionCorrector_0(), _PositionCorrector_1(),
_AngleCorrector(), _EllipticityCorrector()]
class CentralEllipseFitter(EllipseFitter):
"""
A special Fitter class to handle the case of the central pixel in
the galaxy image.
"""
def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
maxit=DEFAULT_MAXIT, fflag=DEFAULT_FFLAG, maxgerr=DEFAULT_MAXGERR,
going_inwards=False):
"""
Perform just a simple 1-pixel extraction at the current (x0, y0)
position using bilinear interpolation.
The input parameters are ignored, but included simple to match
the calling signature of the parent class.
Returns
-------
result : `~photutils.isophote.CentralEllipsePixel` instance
The central pixel value. For convenience, the
`~photutils.isophote.CentralEllipsePixel` class inherits
from the `~photutils.isophote.Isophote` class, although it's
not really a true isophote but just a single intensity value
at the central position. Thus, most of its attributes are
hardcoded to `None` or other default value when appropriate.
"""
self._sample.update()
return CentralPixel(self._sample)