Repository: astropy/photutils

Issue:
Limit offsets when fitting
As noted in #1799, it is possible to limit the flux to be positive by setting the model flux min to be zero:
```python
    dao_psf_model.flux.min = 0
```

which is useful for rejecting negative sources during fitting.

I'd like a mechanism to do the same for `x_0` and `y_0`, but restricting those parameters doesn't work since they need to be set to the local x,y value.  Instead, the limits need to be set independently for each source.  I think that could be accommodated in `make_psf_model`.

One use case is very crowded fields:  I've seen DAOFind-found peaks close enough that they jump within the 5x5 fit zone onto another star's peak.  

A workaround I'll try (but expect to have limited success with) is grouping.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 52be398e..fcc549f1 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -7,6 +7,12 @@ General
 New Features
 ^^^^^^^^^^^^
 
+- ``photutils.psf``
+
+  - Added new ``xy_bounds`` keyword to ``PSFPhotometry`` and
+    ``IterativePSFPhotometry`` to allow one to bound the x and y
+    model parameters during the fitting. [#1805]
+
 Bug Fixes
 ^^^^^^^^^
 
diff --git a/docs/psf.rst b/docs/psf.rst
index 3f88b124..81e3df97 100644
--- a/docs/psf.rst
+++ b/docs/psf.rst
@@ -502,8 +502,16 @@ The astropy modeling and fitting framework also allows for bounding the
 parameter values during the fitting process. However, not all astropy
 "Fitter" classes support parameter bounds. Please see `Fitting Model to
 Data <https://docs.astropy.org/en/stable/modeling/fitting.html>`_ for
-more details. Currently, parameter bounds apply to all sources in the
-image, thus x and y positions cannot be bounded for individual sources.
+more details.
+
+The model parameter bounds apply to all sources in the image,
+thus this mechanism cannot be used to bound the x and y positions
+of individual sources. However, the x and y positions can be
+bounded for individual sources during the fitting by using the
+``xy_bounds`` keyword in `~photutils.psf.PSFPhotometry` and
+`~photutils.psf.IterativePSFPhotometry`. This keyword accepts a tuple of
+floats representing the maximum distance in pixels that a fitted source
+can be from its initial (x, y) position.
 
 For example, you may want to constrain the flux of a source to be
 between certain values or ensure that it is a non-negative value. This
diff --git a/photutils/psf/photometry.py b/photutils/psf/photometry.py
index 9346a3fa..dc38d816 100644
--- a/photutils/psf/photometry.py
+++ b/photutils/psf/photometry.py
@@ -226,6 +226,14 @@ class PSFPhotometry(ModelImageMixin):
         is not converging for sources (e.g., the output ``flags`` value
         contains 8).
 
+    xy_bounds : `None`, float, or 2-tuple of float, optional
+        The maximum distance in pixels that a fitted source can be from
+        the initial (x, y) position. If a single float, then the same
+        maximum distance is used for both x and y. If a 2-tuple of
+        floats, then the distances are in ``(x, y)`` order. If `None`,
+        then no bounds are applied. Either value can also be `None` to
+        indicate no bound in that direction.
+
     localbkg_estimator : `~photutils.background.LocalBackground` or `None`, optional
         The object used to estimate the local background around each
         source. If `None`, then no local background is subtracted. The
@@ -289,13 +297,13 @@ class PSFPhotometry(ModelImageMixin):
 
     def __init__(self, psf_model, fit_shape, *, finder=None, grouper=None,
                  fitter=LevMarLSQFitter(), fitter_maxiters=100,
-                 localbkg_estimator=None, aperture_radius=None,
+                 xy_bounds=None, localbkg_estimator=None, aperture_radius=None,
                  progress_bar=False):
 
         self._param_maps = self._define_param_maps(psf_model)
         self.psf_model = _validate_psf_model(psf_model)
 
-        self.fit_shape = as_pair('fit_shape', fit_shape, lower_bound=(0, 1),
+        self.fit_shape = as_pair('fit_shape', fit_shape, lower_bound=(1, 0),
                                  check_odd=True)
         self.grouper = self._validate_grouper(grouper, 'grouper')
         self.finder = self._validate_callable(finder, 'finder')
@@ -303,6 +311,7 @@ class PSFPhotometry(ModelImageMixin):
         self.localbkg_estimator = self._validate_localbkg(
             localbkg_estimator, 'localbkg_estimator')
         self.fitter_maxiters = self._validate_maxiters(fitter_maxiters)
+        self.xy_bounds = self._validate_bounds(xy_bounds)
         self.aperture_radius = self._validate_radius(aperture_radius)
         self.progress_bar = progress_bar
 
@@ -421,6 +430,22 @@ class PSFPhotometry(ModelImageMixin):
             maxiters = None
         return maxiters
 
+    def _validate_bounds(self, xy_bounds):
+        if xy_bounds is None:
+            return xy_bounds
+
+        xy_bounds = np.atleast_1d(xy_bounds)
+        if len(xy_bounds) == 1:
+            xy_bounds = np.array((xy_bounds[0], xy_bounds[0]))
+        if len(xy_bounds) != 2:
+            raise ValueError('xy_bounds must have 1 or 2 elements')
+        if xy_bounds.ndim != 1:
+            raise ValueError('xy_bounds must be a 1D array')
+        non_none = [i for i in xy_bounds if i is not None]
+        if np.any(np.array(non_none) <= 0):
+            raise ValueError('xy_bounds must be strictly positive')
+        return xy_bounds
+
     @staticmethod
     def _validate_radius(radius):
         if radius is not None and (not np.isscalar(radius)
@@ -702,6 +727,17 @@ class PSFPhotometry(ModelImageMixin):
                 setattr(model, model_param, value)
                 model.name = source['id']
 
+            if self.xy_bounds is not None:
+                if self.xy_bounds[0] is not None:
+                    x_param = getattr(model, self._param_maps['model']['x'])
+                    x_param.bounds = (x_param.value - self.xy_bounds[0],
+                                      x_param.value + self.xy_bounds[0])
+
+                if self.xy_bounds[1] is not None:
+                    y_param = getattr(model, self._param_maps['model']['y'])
+                    y_param.bounds = (y_param.value - self.xy_bounds[1],
+                                      y_param.value + self.xy_bounds[1])
+
             if index == 0:
                 psf_model = model
             else:
@@ -758,11 +794,14 @@ class PSFPhotometry(ModelImageMixin):
         params = defaultdict(list)
         for model in models:
             for name in param_names:
-                value = getattr(model, name).value
+                param = getattr(model, name)
+                value = param.value
                 if (self.data_unit is not None
                         and name == self._param_maps['model']['flux']):
                     value <<= self.data_unit  # add the flux units
                 params[name].append(value)
+                params[f'{name}_fixed'].append(param.fixed)
+                params[f'{name}_bounds'].append(param.bounds)
 
         table = QTable(params)
         ids = np.arange(len(table)) + 1
@@ -1137,6 +1176,20 @@ class PSFPhotometry(ModelImageMixin):
         except KeyError:
             pass
 
+        # add flag = 32 if x or y fitted value is at the bounds
+        if self.xy_bounds is not None:
+            xcolname = self._param_maps['model']['x']
+            ycolname = self._param_maps['model']['y']
+            for index, row in enumerate(self._fit_model_params):
+                x_bounds = row[f'{xcolname}_bounds']
+                x_bounds = np.array([i for i in x_bounds if i is not None])
+                y_bounds = row[f'{ycolname}_bounds']
+                y_bounds = np.array([i for i in y_bounds if i is not None])
+                dx = x_bounds - row[xcolname]
+                dy = y_bounds - row[ycolname]
+                if np.any(dx == 0) or np.any(dy == 0):
+                    flags[index] += 32
+
         return flags
 
     def _prepare_fit_inputs(self, data, *, mask=None, error=None,
@@ -1302,6 +1355,7 @@ class PSFPhotometry(ModelImageMixin):
                     iterations using the ``fitter_maxiters`` keyword.
                   * 16 : the fitter parameter covariance matrix was not
                     returned
+                  * 32 : the fit x or y position is at the bounded value
         """
         if isinstance(data, NDData):
             data_ = data.data
@@ -1449,6 +1503,14 @@ class IterativePSFPhotometry(ModelImageMixin):
         The maximum number of iterations in which the ``fitter`` is
         called for each source.
 
+    xy_bounds : `None`, float, or 2-tuple of float, optional
+        The maximum distance in pixels that a fitted source can be from
+        the initial (x, y) position. If a single float, then the same
+        maximum distance is used for both x and y. If a 2-tuple of
+        floats, then the distances are in ``(x, y)`` order. If `None`,
+        then no bounds are applied. Either value can also be `None` to
+        indicate no bound in that direction.
+
     maxiters : int, optional
         The maximum number of PSF-fitting/subtraction iterations to
         perform.
@@ -1551,8 +1613,9 @@ class IterativePSFPhotometry(ModelImageMixin):
     """
 
     def __init__(self, psf_model, fit_shape, finder, *, grouper=None,
-                 fitter=LevMarLSQFitter(), fitter_maxiters=100, maxiters=3,
-                 mode='new', localbkg_estimator=None, aperture_radius=None,
+                 fitter=LevMarLSQFitter(), fitter_maxiters=100,
+                 xy_bounds=None, maxiters=3, mode='new',
+                 localbkg_estimator=None, aperture_radius=None,
                  sub_shape=None, progress_bar=False):
 
         if finder is None:
@@ -1566,6 +1629,7 @@ class IterativePSFPhotometry(ModelImageMixin):
         self._psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder,
                                       grouper=grouper, fitter=fitter,
                                       fitter_maxiters=fitter_maxiters,
+                                      xy_bounds=xy_bounds,
                                       localbkg_estimator=localbkg_estimator,
                                       aperture_radius=aperture_radius,
                                       progress_bar=progress_bar)
@@ -1800,6 +1864,7 @@ class IterativePSFPhotometry(ModelImageMixin):
                   * 8 : the fitter may not have converged
                   * 16 : the fitter parameter covariance matrix was not
                     returned
+                  * 32 : the fit x or y position is at the bounded value
         """
         if isinstance(data, NDData):
             data_ = data.data


Test patch:
diff --git a/photutils/psf/tests/test_photometry.py b/photutils/psf/tests/test_photometry.py
index 2214bf92..532a83ed 100644
--- a/photutils/psf/tests/test_photometry.py
+++ b/photutils/psf/tests/test_photometry.py
@@ -74,7 +74,7 @@ def test_invalid_inputs():
         with pytest.raises(ValueError, match=match):
             _ = PSFPhotometry(model, shape)
 
-    match = 'fit_shape must be > 0'
+    match = 'fit_shape must be >= 1'
     with pytest.raises(ValueError, match=match):
         _ = PSFPhotometry(model, (-1, 1))
 
@@ -851,6 +851,67 @@ def test_fitter_no_maxiters_no_metrics(test_data):
         assert np.all(np.isnan(phot['cfit']))
 
 
+@pytest.mark.skipif(not HAS_SCIPY, reason='scipy is required')
+def test_xy_bounds(test_data):
+    data, error, _ = test_data
+
+    psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35)
+    fit_shape = (5, 5)
+    init_params = QTable()
+    init_params['x'] = [65]
+    init_params['y'] = [51]
+    xy_bounds = (1, 1)
+    psfphot = PSFPhotometry(psf_model, fit_shape, finder=None,
+                            aperture_radius=4, xy_bounds=xy_bounds)
+    phot = psfphot(data, error=error, init_params=init_params)
+    assert len(phot) == len(init_params)
+    assert phot['x_fit'] == 64.0  # at lower bound
+    assert phot['y_fit'] == 50.0  # at lower bound
+    assert phot['flags'] == 32
+
+    psfphot2 = PSFPhotometry(psf_model, fit_shape, finder=None,
+                             aperture_radius=4, xy_bounds=1)
+    phot2 = psfphot2(data, error=error, init_params=init_params)
+    cols = ('x_fit', 'y_fit', 'flux_fit')
+    for col in cols:
+        assert np.all(phot[col] == phot2[col])
+
+    xy_bounds = (None, 1)
+    psfphot = PSFPhotometry(psf_model, fit_shape, finder=None,
+                            aperture_radius=4, xy_bounds=xy_bounds)
+    phot = psfphot(data, error=error, init_params=init_params)
+    assert phot['x_fit'] < 64.0
+    assert phot['y_fit'] == 50.0  # at lower bound
+    assert phot['flags'] == 32
+
+    xy_bounds = (1, None)
+    psfphot = PSFPhotometry(psf_model, fit_shape, finder=None,
+                            aperture_radius=4, xy_bounds=xy_bounds)
+    phot = psfphot(data, error=error, init_params=init_params)
+    assert phot['x_fit'] == 64.0  # at lower bound
+    assert phot['y_fit'] < 50.0
+    assert phot['flags'] == 32
+
+    xy_bounds = (None, None)
+    psfphot = PSFPhotometry(psf_model, fit_shape, finder=None,
+                            aperture_radius=4, xy_bounds=xy_bounds)
+    phot = psfphot(data, error=error, init_params=init_params)
+    assert phot['x_fit'] < 64.0
+    assert phot['y_fit'] < 50.0
+    assert phot['flags'] == 0
+
+    # test invalid inputs
+    match = 'xy_bounds must have 1 or 2 elements'
+    with pytest.raises(ValueError, match=match):
+        PSFPhotometry(psf_model, fit_shape, xy_bounds=(1, 2, 3))
+    match = 'xy_bounds must be a 1D array'
+    with pytest.raises(ValueError, match=match):
+        PSFPhotometry(psf_model, fit_shape, xy_bounds=np.ones((1, 1)))
+    match = 'xy_bounds must be strictly positive'
+    with pytest.raises(ValueError, match=match):
+        PSFPhotometry(psf_model, fit_shape, xy_bounds=(-1, 2))
+
+
 @pytest.mark.skipif(not HAS_SCIPY, reason='scipy is required')
 def test_iterative_psf_photometry_mode_new(test_data):
     data, error, sources = test_data