/usr/lib/python3/dist-packages/photutils/isophote/fitter.py is in python3-photutils 0.4-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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[0] or
sample.geometry.y0 < 1. or
sample.geometry.y0 > sample.image.shape[1]):
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)
|