Added embryonary model-fitting section.

  Position

  PositionWin

  Photom

  Model

.. [#thresh] Pixel values also have to exceed the local detection threshold set with ``DETECT_THRESH``.

.. [#psf_models] PSF models be computed using the |PSFEx|_ package.

.. File Model.rst

.. include:: global.rst

Model fitting

=============

SExtractor can fit models to the images of detected objects since version 2.8. The fit is performed by minimizing the loss function

.. math::

  :label: loss_func

  \lambda(\boldsymbol{q}) = \sum_i \left(g\left(\frac{p_i - \hat{m}_i(\boldsymbol{q})}{\sigma_i}\right)\right)^2 + \sum_j \frac{q_j - \mu_j}{}

with respect to components of the model parameter vector :math:`\boldsymbol{q}`. :math:`\boldsymbol{q}` comprises parameters describing the shape of the model and the model pixel coordinates :math:`\boldsymbol{x}`.

The first term in :eq:`loss_func` is a modified `weighted sum of squares <http://en.wikipedia.org/wiki/Least_squares#Weighted_least_squares>`_ that aims at minimizing the residuals of the fit. :math:`p_i`, :math:`\hat{m}_i(\boldsymbol{q})` and :math:`\sigma_i` are respectively the pixel value above the background, the value of the resampled model, and the pixel value uncertainty at image pixel :math:`i`.

:math:`g(u)` is a derivable monotonous function that reduces the influence of large deviations from the model (e.g., contamination by neighbors):

.. math::

  :label: loss_func

  g(u) = \left\{

    \begin{array}{rl}

       u_0 \log \left(1 + \frac{u}{u_0}\right) & \mbox{if } u \ge 0,\\

      -u_0 \log \left(1 - \frac{u}{u_0}\right) & \mbox{otherwise.}\\

    \end{array}

  \right.

The vector :math:`\hat{\boldsymbol{m}}(\boldsymbol{q})` is obtained by convolving the high resolution model :math:`\boldsymbol{m}(\boldsymbol{q})` with the local PSF model :math:`\boldsymbol{\phi}` and applying a resampling operator :math:`\mathbf{R}(\boldsymbol{x})` to generate the final model raster at position :math:`\boldsymbol{x}` at the nominal image resolution:

.. math::

  :label: model_resampling

  \hat{\boldsymbol{m}}(\boldsymbol{q}) = \mathbf{R}(\boldsymbol{x}) (\boldsymbol{m}(\boldsymbol{q})*\boldsymbol{\phi}).

Levenberg-Marquardt minimization, inside a disk which diameter is scaled to include the isophotal footprint plus a 20 % margin, plus the size of the PSF model image.

The models that can be fit are:

- Exponential disk

  .. math::

      \Sigma_{\tt ExpDisk}(R) = \Sigma(0) \exp \left (- {R\over h}\right ) 

      \label{expdisk}

- Sérsic (:math:`R^{1/n}`) spheroid (bulg)

  .. math::

      \Sigma_{\tt Sersic}(R) = \Sigma(0) \exp \left [- b(n)\,\left({R\over

          R_e}\right)^{1/n}\right ] \ ,

      \label{sersic}

where, for the :raw-latex:`\cite{sersic:1968}` model, :math:`b(n)` is the solution of

  .. math::

      2 \gamma[2\,n,b(n)] = \Gamma(2\,n)

      \label{bofn}

  An accurate approximation for the solution for :math:`b(n)` of equation (bofn) is :raw-latex:`\citep{ciotti:bertin:1999}`

   .. math::

      b(n) = 2\,n - {1\over3} + {4\over 405\,n} + {46\over 25515\,n^2} + {131\over

        1148175\,n^3}

-  :raw-latex:`\cite{devaucouleurs48}` spheroid (bulge, eq. [[sersic]], with :math:`n=4`)

-  Exponential disk + Sérsic (:math:`R^{1/n}`) spheroid (bulge)

-  Point source

-  Background (constant)

For these models, SExtractor can compute fluxes and magnitudes, as well

as sizes (disk scale length for the disks and effective — projected

half-light — radii for the spheroids), characteristic surface

magnitudes, and Sérsic index, as well as their uncertainties.

The models are concentric (they assume the same center) and are all

convolved with the PSF, given by the .psf file, which must be determined

by first running PSFEx (see below).

Unfortunately, the Sérsic profile is very cuspy in the center for

:math:`n>2`. To avoid huge wings in the FFTs when convolving the profile

with the PSF, the profile is split between a 3rd order polynomial,

analytically fit to match, in intensity and its 1st and 2nd spatial

derivatives, the Sérsic profile at :math:`R=4\,\rm pixels`,

:math:`I(r) = I_0 + (r/a)^3`, which has zero first and 2nd derivative at

the center, i.e. a homogeneous core on one hand, and a residual with

finite extent on the other.

For the fit of the spheroid component, the apparent ellipticity allowed

is taken in the range :math:`[0.5, 2]` . This obviously forbids very

flat spheroids to avoid confusion with a flattened disk. By allowing

ellipticities greater than unity, SExtractor avoids dichotomies of

position angle when the ellipticity is very low. The Sérsic index is

allowed values between 1 and 10.

Models are measured according to the following table.

.. math::

   \begin{aligned}

   \hbox{{\tt FLUX\_BACKOFFSET} or {\tt FLUXERR\_BACKOFFSET}} &\to& \hbox{background}

   \nonumber \\

   \hbox{{\tt DISK\_xxx}} &\to& \hbox{exponential disk} \nonumber \\

   \hbox{{\tt SPHEROID\_SERSICN} or {\tt SPHEROID\_SERSICNERR}} &\to&

   \hbox{S\'ersic} \nonumber \\

   \hbox{{\tt SPHEROID\_xxx} without {\tt SPEHEROID\_SERSICN[ERR]}} &\to&

   \hbox{de Vaucouleurs (}n=4 \hbox{ S\'ersic)} \nonumber \\ 

   \hbox{{\tt MODEL\_xxx} only} &\to& \hbox{S\'ersic [???]} \nonumber \\

   \hbox{{\tt SPHEROID\_xxx} and {\tt DISK\_xxx}}&\to& \hbox{S\'ersic spheroid + exponential disk [???]} \nonumber \end{aligned}

Table [modeltriggers] should be interpreted as meaning that if one of

the parameters given in the parameter file (e.g. default.param) includes

the string on the left of the arrow, the model to the right of the arrow

is triggered. For example, when including parameters that contain the

string ‘MODEL’, both galaxies and stars are fit with convolutions of

Sérsic models with the PSF. If no SPHEROID\_xxx or DISK\_xxx parameter

is present, but the model-fitting process is nevertheless triggered by

the presence of other measurement parameters or relevant

CHECKIMAGE\_TYPEs , a single component with Sérsic profile and

adjustable Sérsic index :math:`n` is fitted.

The number of parameters that are fit are 2 for the global center, 4 per

model for the scale, normalization, aspect ratio and position angle,

plus the index for the Sérsic model. For example, fitting a Sérsic +

exponential disk involves a fitting 11 parameters.

Experience shows that the de Vaucouleurs spheroid + exponential disk

combination provides fairly accurate and robust fits for moderately

resolved faint galaxies. An adjustable Sérsic index may offer lower

residuals on spheroids and/or well-resolved galaxies, but makes the fit

less robust and more sensitive to PSF model errors. One might think of

adding some mechanism to lock or unlock the Sérsic index automatically

in future versions of SExtractor.

The measurement parameters related to model-fitting follow the usual

SExtractor rules:

Flux measurements are available in ADUs (FLUX\_xxx parameters) or

magnitudes (MAG\_xxx parameters), Coordinates and radii are available in

pixels or celestial units (provided that the FITS image header contains

the appropriate WCS information).

xxxMODEL\_yyy measurement parameters deal with the global fitted model,

i.e. the sum of all components (e.g. chi-square per d.o.f. CHI2\_MODEL,

PSF-corrected ellipticities E1/2MODEL\_IMAGE, EPS1/2 MODEL\_IMAGE).

:math:`1\,\sigma` error estimates xxxERR\_yyy are provided for most

measurement parameters; they are obtained by marginalizing the full

covariance matrix of the fit.

Since the model fitting involves convolution with the PSF, it is

imperative to launch PSFEx before launching SExtractor. In practice, the

sequence of operations is:

#. Run SExtractor to prepare PSFEx;

#. Run PSFEx to prepare model fits in SExtractor;

#. Run SExtractor with model fit parameters.

   = \overline{x_{\tt WIN}}^{(t)} + 2\,\frac{\sum_{r_i^{(t)} < r_{\rm max}}

   w_i^{(t)} I_i \ (x_i  - \overline{x_{\tt WIN}}^{(t)})}

   {\sum_{r_i^{(t)} < r_{\rm max}} w_i^{(t)} I_i},\\

   \label{eq:ywin}

   {\tt YWIN}^{(t+1)} & = & \overline{y_{\tt WIN}}^{(t+1)}

   = \overline{y_{\tt WIN}}^{(t)} + 2\,\frac{\sum_{r_i^{(t)} < r_{\rm max}}

   w_i^{(t)} I_i\ (y_i - \overline{y_{\tt WIN}}^{(t)})}

# There are two options for replacing |today|: either, you set today to some

# non-false value, then it is used:

#

today = 'Wed Nov 22 2017'

today = 'Wed Dec 13 2017'

#

# Else, today_fmt is used as the format for a strftime call.

#