Merge pull request #922 from ibusko/fixed_parameters

Adds ability to individually fix the fitting parameters in 'ellipse' (#893)

CHANGES.rst

@@ -42,6 +42,9 @@ New Features
   - Significantly improved the performance (~5 times faster) of
     ellipse fitting. [#826]
 
+  - Added the ability to individually fix the ellipse-fitting
+    parameters. [#922]
+
 - ``photutils.psf``
 
   - Added new centroiding function ``centroid_epsf``. [#816]

photutils/isophote/ellipse.py

@@ -203,7 +203,8 @@ def fit_image(self, sma0=None, minsma=0., maxsma=None, step=0.1,
                   conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
                   maxit=DEFAULT_MAXIT, fflag=DEFAULT_FFLAG,
                   maxgerr=DEFAULT_MAXGERR, sclip=3., nclip=0,
-                  integrmode=BILINEAR, linear=None, maxrit=None):
+                  integrmode=BILINEAR, linear=None, maxrit=None,
+                  fix_center=False, fix_pa=False, fix_eps=False):
         # This parameter list is quite large and should in principle be
         # simplified by re-distributing these controls to somewhere else.
         # We keep this design though because it better mimics the flat
@@ -216,6 +217,15 @@ def fit_image(self, sma0=None, minsma=0., maxsma=None, step=0.1,
         isophote at each sma.  The entire set of isophotes is returned
         in an `~photutils.isophote.IsophoteList` instance.
 
+        Note that the fix_XXX parameters act in unison. Meaning, if one
+        of them is set via this call, the others will assume their default
+        (False) values. This effectively overrides any settings that are
+        present in the internal `~photutils.isophote.EllipseGeometry`
+        instance that is carried along as a property of this class. If an
+        instance of `~photutils.isophote.EllipseGeometry` was passed to this
+        class' constructor, that instance will be effectively overriden by
+        the fix_XXX parameters in this call.
+
         Parameters
         ----------
         sma0 : float, optional
@@ -334,6 +344,13 @@ def fit_image(self, sma0=None, minsma=0., maxsma=None, step=0.1,
             whenever the ellipticity exceeds 1.0 or the ellipse center
             crosses the image boundaries.  If `None` (default), then no
             maximum value is used.
+        fix_center : bool, optional
+            Keep center of ellipse fixed during fit? The default is False.
+        fix_pa : bool, optional
+            Keep position angle of semi-major axis of ellipse fixed during fit?
+            The default is False.
+        fix_eps : bool, optional
+            Keep ellipticity of ellipse fixed during fit? The default is False.
 
         Returns
         -------
@@ -360,6 +377,13 @@ def fit_image(self, sma0=None, minsma=0., maxsma=None, step=0.1,
             self._geometry.linear_growth = linear
         else:
             linear = self._geometry.linear_growth
+        if fix_center and fix_pa and fix_eps:
+            warnings.warn(': Everything is fixed. Fit not possible.',
+                          AstropyUserWarning)
+            return IsophoteList([])
+        if fix_center or fix_pa or fix_eps:
+            # Note that this overrides the geometry instance for good.
+            self._geometry.fix = np.array([fix_center, fix_center, fix_pa, fix_eps])
 
         # first, go from initial sma outwards until
         # hitting one of several stopping criteria.
@@ -612,7 +636,7 @@ def _non_iterative(self, sma, step, linear, geometry, sclip, nclip,
         sample = EllipseSample(self.image, sma, astep=step, sclip=sclip,
                                nclip=nclip, linear_growth=linear,
                                geometry=geometry, integrmode=integrmode)
-        sample.update()
+        sample.update(geometry.fix)
 
         # build isophote without iterating with an EllipseFitter
         isophote = Isophote(sample, 0, True, stop_code=4)
@@ -632,7 +656,7 @@ def _fix_last_isophote(isophote_list, index):
             # force new extraction of raw data, since
             # geometry changed.
             isophote.sample.values = None
-            isophote.sample.update()
+            isophote.sample.update(isophote.sample.geometry.fix)
 
             # we take the opportunity to change an eventual
             # negative stop code to its' positive equivalent.

photutils/isophote/fitter.py

@@ -7,6 +7,7 @@
 
 from astropy import log
 import numpy as np
+import numpy.ma as ma
 
 from .harmonics import (first_and_second_harmonic_function,
                         fit_first_and_second_harmonics)
@@ -138,9 +139,12 @@ def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
         minimum_amplitude_value = np.Inf
         minimum_amplitude_sample = None
 
+        # these must be passed throughout the execution chain.
+        fixed_parameters = self._sample.geometry.fix
+
         for i in range(maxit):
             # Force the sample to compute its gradient and associated values.
-            sample.update()
+            sample.update(fixed_parameters)
 
             # The extract() method returns sampled values as a 2-d numpy array
             # with the following structure:
@@ -154,15 +158,18 @@ def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
             # marked as invalid.
             try:
                 coeffs = fit_first_and_second_harmonics(values[0], values[2])
+                coeffs = coeffs[0]
             except Exception as e:
                 log.info(e)
+                sample.geometry.fix = fixed_parameters
                 return Isophote(sample, i + 1, False, 3)
 
-            coeffs = coeffs[0]
+            # Mask out coefficients that control fixed ellipse parameters.
+            free_coeffs = ma.masked_array(coeffs[1:], mask=fixed_parameters)
 
-            # largest harmonic in absolute value drives the correction.
-            largest_harmonic_index = np.argmax(np.abs(coeffs[1:]))
-            largest_harmonic = coeffs[1:][largest_harmonic_index]
+            # Largest non-masked harmonic in absolute value drives the correction.
+            largest_harmonic_index = np.argmax(np.abs(free_coeffs))
+            largest_harmonic = free_coeffs[largest_harmonic_index]
 
             # see if the amplitude decreased; if yes, keep the
             # corresponding sample for eventual later use.
@@ -179,15 +186,15 @@ def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
                 # Got a valid solution. But before returning, ensure
                 # that a minimum of iterations has run.
                 if i >= minit - 1:
-                    sample.update()
+                    sample.update(fixed_parameters)
                     return Isophote(sample, i + 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()
+                minimum_amplitude_sample.update(fixed_parameters)
                 return Isophote(minimum_amplitude_sample, i + 1, True, 1)
 
             # pick appropriate corrector code.
@@ -203,21 +210,21 @@ def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
             # (hopefully smaller) price here, by having multiple calls to
             # the EllipseSample constructor.
             sample = corrector.correct(sample, largest_harmonic)
-            sample.update()
+            sample.update(fixed_parameters)
 
             # 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()
+                sample.update(fixed_parameters)
                 return Isophote(sample, i + 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()
+        minimum_amplitude_sample.update(fixed_parameters)
         return Isophote(minimum_amplitude_sample, maxit, True, 2)
 
     @staticmethod
@@ -388,6 +395,8 @@ def fit(self, conver=DEFAULT_CONVERGENCE, minit=DEFAULT_MINIT,
             at the central position.  Thus, most of its attributes are
             hardcoded to `None` or other default value when appropriate.
         """
+        # default values
+        fixed_parameters = np.array([False, False, False, False])
 
-        self._sample.update()
+        self._sample.update(fixed_parameters)
         return CentralPixel(self._sample)

photutils/isophote/geometry.py

@@ -104,11 +104,19 @@ class EllipseGeometry:
         as a relative value (when ``linear_growth=False``).  The default
         is 0.1.
     linear_growth : bool, optional
-        The semimajor axis growing/shrinking mode.  The default is
+        The semimajor axis growing/shrinking mode. The default is
         `False`.
+    fix_center : bool, optional
+        Keep center of ellipse fixed during fit? The default is False.
+    fix_pa : bool, optional
+        Keep position angle of semi-major axis of ellipse fixed during fit?
+        The default is False.
+    fix_eps : bool, optional
+        Keep ellipticity of ellipse fixed during fit? The default is False.
     """
 
-    def __init__(self, x0, y0, sma, eps, pa, astep=0.1, linear_growth=False):
+    def __init__(self, x0, y0, sma, eps, pa, astep=0.1, linear_growth=False,
+                 fix_center=False, fix_pa=False, fix_eps=False):
         self.x0 = x0
         self.y0 = y0
         self.sma = sma
@@ -118,6 +126,11 @@ def __init__(self, x0, y0, sma, eps, pa, astep=0.1, linear_growth=False):
         self.astep = astep
         self.linear_growth = linear_growth
 
+        # Fixed parameters are flagged in here. Note that the
+        # ordering must follow the same ordering used in the
+        # fitter._CORRECTORS list.
+        self.fix = np.array([fix_center, fix_center, fix_pa, fix_eps])
+
         # limits for sector angular width
         self._phi_min = 0.05
         self._phi_max = 0.2

photutils/isophote/sample.py

@@ -278,7 +278,7 @@ def _iter_sigma_clip(self, angles, radii, intensities):
 
         return r_angles, r_radii, r_intensities
 
-    def update(self):
+    def update(self, fixed_parameters):
         """
         Update this `~photutils.isophote.EllipseSample` instance.
 
@@ -288,6 +288,7 @@ def update(self):
         then computes the the mean intensity, local gradient, and other
         associated quantities.
         """
+        self.geometry.fix = fixed_parameters
 
         step = self.geometry.astep
 
@@ -385,11 +386,13 @@ class CentralEllipseSample(EllipseSample):
     the special case of the central pixel in the galaxy image.
     """
 
-    def update(self):
+    def update(self, fixed_parameters):
         """
         Update this `~photutils.isophote.EllipseSample` instance with
         the intensity integrated at the (x0, y0) center position using
         bilinear integration. The local gradient is set to `None`.
+
+        'fixed_parameters' is ignored in this subclass.
         """
 
         s = self.extract()

photutils/isophote/tests/test_ellipse.py

@@ -77,9 +77,12 @@ def test_linear(self):
 
         # difference in sma between successive isohpotes must be constant.
         step = isophote_list[-1].sma - isophote_list[-2].sma
-        assert math.isclose((isophote_list[-2].sma - isophote_list[-3].sma), step, rel_tol=0.01)
-        assert math.isclose((isophote_list[-3].sma - isophote_list[-4].sma), step, rel_tol=0.01)
-        assert math.isclose((isophote_list[2].sma - isophote_list[1].sma), step, rel_tol=0.01)
+        assert math.isclose((isophote_list[-2].sma - isophote_list[-3].sma),
+                            step, rel_tol=0.01)
+        assert math.isclose((isophote_list[-3].sma - isophote_list[-4].sma),
+                            step, rel_tol=0.01)
+        assert math.isclose((isophote_list[2].sma - isophote_list[1].sma),
+                            step, rel_tol=0.01)
 
     def test_fit_one_ellipse(self):
         ellipse = Ellipse(self.data)

photutils/isophote/tests/test_fitter.py

@@ -27,10 +27,12 @@
 DATA = make_test_image(random_state=123)
 DEFAULT_POS = 256
 
+DEFAULT_FIX = np.array([False, False, False, False])
+
 
 def test_gradient():
     sample = EllipseSample(DATA, 40.)
-    sample.update()
+    sample.update(DEFAULT_FIX)
 
     assert_allclose(sample.mean, 200.02, atol=0.01)
     assert_allclose(sample.gradient, -4.222, atol=0.001)
@@ -49,7 +51,7 @@ def test_fitting_raw():
     # pick first guess ellipse that is off in just
     # one of the parameters (eps).
     sample = EllipseSample(DATA, 40., eps=2*0.2)
-    sample.update()
+    sample.update(DEFAULT_FIX)
     s = sample.extract()
 
     harmonics = fit_first_and_second_harmonics(s[0], s[2])

photutils/isophote/tests/test_isophote.py

@@ -14,6 +14,8 @@
 from ..sample import EllipseSample
 from ...datasets import get_path
 
+DEFAULT_FIX = np.array([False, False, False, False])
+
 try:
     import scipy  # noqa
     HAS_SCIPY = True
@@ -126,7 +128,7 @@ def build_list(data, sma0, slen=5):
         iso_list = []
         for k in range(slen):
             sample = EllipseSample(data, float(k + sma0))
-            sample.update()
+            sample.update(DEFAULT_FIX)
             iso_list.append(Isophote(sample, k, True, 0))
         result = IsophoteList(iso_list)
         return result

photutils/isophote/tests/test_sample.py

@@ -11,6 +11,7 @@
 from ..isophote import Isophote
 from ..sample import EllipseSample
 
+DEFAULT_FIX = np.array([False, False, False, False])
 
 DATA = make_test_image(background=100., i0=0., noise=10., random_state=123)
 
@@ -32,7 +33,7 @@ def test_scatter(integrmode, amin, amax):
     """
 
     sample = EllipseSample(DATA, 50., astep=0.2, integrmode=integrmode)
-    sample.update()
+    sample.update(DEFAULT_FIX)
     iso = Isophote(sample, 0, True, 0)
 
     assert iso.pix_stddev < amax
@@ -41,7 +42,7 @@ def test_scatter(integrmode, amin, amax):
 
 def test_coordinates():
     sample = EllipseSample(DATA, 50.)
-    sample.update()
+    sample.update(DEFAULT_FIX)
     x, y = sample.coordinates()
 
     assert isinstance(x, np.ndarray)
@@ -50,7 +51,7 @@ def test_coordinates():
 
 def test_sclip():
     sample = EllipseSample(DATA, 50., nclip=3)
-    sample.update()
+    sample.update(DEFAULT_FIX)
     x, y = sample.coordinates()
 
     assert isinstance(x, np.ndarray)