{ "cells": [ { "cell_type": "markdown", "id": "3e52c34d", "metadata": {}, "source": [ "# Image Fitting\n", "\n", "\n", "Most galaxy light profiles can be well described by PSF-convolved models like the Sérsic profile. PetroFit uses the `astropy` `modeling` sub-module to provide tools to perform two-dimensional fits of galaxy light profiles. To this end, we use the PetroFit ` PSFConvolvedModel2D` class, which applies PSF convolution to and handles oversampling for `astropy` based models. \n", "\n", "In this section, we demonstrate the basics of light profile modeling on a galaxy using a single component Sérsic profile. \n", "\n", "To start with PetroFit, simply import it as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "c6987e5f", "metadata": {}, "outputs": [], "source": [ "import petrofit as pf " ] }, { "cell_type": "markdown", "id": "a68cbd8e", "metadata": {}, "source": [ "## Loading Example Data\n", "\n", "The dataset we're using is a synthetic image of a galaxy, created using astropy's `Sersic2D` model.\n", "This galaxy representation is convolved with a PSF for the F105W filter using petrofit's PSFConvolvedModel2D to simulate observational data.\n", "We also added noise to the data and provide a corresponding RMS map. \n", "\n", "Key features of the synthetic galaxy:\n", "\n", "- Sérsic index of 1 (exponential profile).\n", "- Effective radius of 15 pixels.\n", "- Positioned at (100, 75) pixels.\n", "- Rotated by $\\frac{\\pi}{4}$.\n", "- With ellip=0.1\n", "\n", "\n", "### Loading Data and RMS Images\n", "We first use `astropy`'s ``CCDData`` to load the example data and visualize it through `matplotlib`. The RMS image is loaded using `astropy`'s ``fits`` sub-module.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "91575418", "metadata": {}, "outputs": [], "source": [ "from astropy.nddata import CCDData\n", "from astropy.io import fits \n", "\n", "image = CCDData.read('data/example_sersic.fits.gz', unit='electron s-1')\n", "rms = fits.getdata('data/example_rms.fits.gz')" ] }, { "cell_type": "code", "execution_count": null, "id": "7c75fd46", "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "# Hidden cell\n", "\n", "%matplotlib inline\n", "\n", "# Stop Fit Model to Data section warnings\n", "import warnings\n", "warnings.filterwarnings('ignore', append=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "ffb27660", "metadata": { "scrolled": false }, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "plt.rcParams['figure.figsize'] = [6, 6]\n", "plt.rcParams['image.origin'] = 'lower'\n", "plt.rcParams['font.size'] = 12\n", "\n", "vmax = 0.005 # Use the image std as max and min of all plots \n", "vmin = - vmax \n", "\n", "fig, axs = plt.subplots(1,2, figsize=[12, 6])\n", "plt.sca(axs[0])\n", "plt.imshow(image.data, vmin=vmin, vmax=vmax)\n", "plt.title(\"Mock Galaxy\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "\n", "plt.sca(axs[1])\n", "plt.imshow(rms)\n", "plt.title(\"RMS Image\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "84dec3c3", "metadata": {}, "source": [ "## PSF \n", "\n", "A Point Spread Function (PSF) describes how light from a point source is distributed on detector due to optical effects such as diffraction. Images or cutouts of stars are good approximations of PSFs because stars are single-point sources and their images describe how their light is distributed on the detector. To make cutouts of stars in an image, use the ` astropy.nddata.Cutout2D` function.\n", "\n", "The following PSF is a cutout of a star in the Hubble Frontier Fields image of Abell 2744 (same dataset as the example image). Since we will be using the PSF image as a convolution kernel, it is **very important** that the following requirements are satisfied:\n", "\n", "- The image of the PSF should be at the same resolution as the data.\n", "- The star or PSF is centered in the image.\n", "- The PSF image does not contain other sources. \n", "- The image is normalized so that the sum of the PSF image is near or equal to 1.0. \n", "- The PSF image should have odd dimensions on each side (because it is a convolution kernel). \n" ] }, { "cell_type": "code", "execution_count": null, "id": "8dc9f08a", "metadata": {}, "outputs": [], "source": [ "from astropy.io import fits\n", "\n", "# Load PSF image (2D array)\n", "PSF = fits.getdata('data/f105w_psf.fits.gz')\n", "\n", "# Normalize PSF \n", "PSF = PSF / PSF.sum()\n", "\n", "# Note that the PSF shape is odd on all sides\n", "print(\"PSF Shape = {}\".format(PSF.shape))\n", "\n", "# Plot PSF and use vmax and vmin to show difraction spikes\n", "plt.imshow(PSF, vmin=0, vmax=PSF.std()/10)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ab39545f", "metadata": {}, "source": [ "## Sérsic Model\n", "\n", "### Sérsic Parameters \n", "\n", "The `amplitude`, `r_eff`, `n`, `x_0`, `y_0`, `ellip`, and `theta` represent the galaxy's brightness, \n", "effective radius, Sérsic index, position, ellipticity, and orientation, respectively. Here we make rough estimates of the parameters:" ] }, { "cell_type": "code", "execution_count": null, "id": "63efe398", "metadata": {}, "outputs": [], "source": [ "amplitude=0.2\n", "r_eff=20\n", "n=1\n", "x_0=107\n", "y_0=70\n", "ellip=0.1\n", "theta=0.1" ] }, { "cell_type": "markdown", "id": "33f984f4", "metadata": {}, "source": [ "### AstroPy Sérsic Model\n", "\n", "Here, we are setting up a 2D galaxy light profile model using astropy's Sersic2D model. \n", "The Sersic2D model is a widely-used representation of the light distribution of elliptical galaxies.\n", "We also define a set of `bounds`, a dictionary of lower and upper bounds of parameters. \n", "Keys are parameter names. The values are a list or a tuple of length 2 giving the desired range for \n", "the parameter and a value of `None` means no bounds. The default bounds can be provided using the \n", "PetroFit `get_default_sersic_bounds` function. For example, we restrain the fitter from exploring \n", "half-light radii that are negative by adding `'r_eff': (0, None)`. \n", "We also apply a custom restriction for the center of the model to be within a range (`center_slack`) from the initial guess." ] }, { "cell_type": "code", "execution_count": null, "id": "3cb0021c", "metadata": {}, "outputs": [], "source": [ "from astropy.modeling import models \n", "\n", "center_slack = 20\n", "\n", "sersic_model = models.Sersic2D(\n", " \n", " amplitude=amplitude,\n", " r_eff=r_eff,\n", " n=n,\n", " x_0=x_0,\n", " y_0=y_0,\n", " ellip=ellip, \n", " theta=theta,\n", " bounds = pf.get_default_sersic_bounds({\n", " 'x_0': (x_0 - center_slack/2, x_0 + center_slack/2),\n", " 'y_0': (y_0 - center_slack/2, y_0 + center_slack/2),\n", " }),\n", ") \n", "\n", "sersic_model" ] }, { "cell_type": "markdown", "id": "efbbdd0e", "metadata": {}, "source": [ "## PSFConvolvedModel2D\n", "\n", "The `petrofit` `PSFConvolvedModel2D` is a `Fittable2DModel` that adds PSF convolution and model to image sampling to `astropy` core models. `PSFConvolvedModel2D` makes an image of the underlying model and samples it onto a grid. The model image is then convolved with a PSF if one is provided. Since `PSFConvolvedModel2D` is a `Fittable2DModel`, it can be used to fit model images to data images. For example, we wrap an `astropy` `Sersic2D` model in this doc with `PSFConvolvedModel2D`, which produces an oversampled and PSF convolved version of the Sérsic profile at each iteration of the fitting algorithm. **Note that `PSFModel` is deprecated and replaced by `PSFConvolvedModel2D`.**\n" ] }, { "cell_type": "raw", "id": "3653d689", "metadata": {}, "source": [ "
\n", "

Note

\n", "

PSFConvolvedModel2D is agnostic to the models it wraps and can handle complex multi-component astropy models.

\n", "
" ] }, { "cell_type": "markdown", "id": "f826cd76", "metadata": {}, "source": [ "### Pixel Centering in PSFConvolvedModel2D\n", "PSFConvolvedModel2D adopts the DS9 coordinate system, where the pixel index corresponds to its center. Thus, an index of 0 designates the center of the first pixel. This is distinct from the GALFIT convention, and users should note this difference when comparing results between tools.\n" ] }, { "cell_type": "markdown", "id": "817344e0", "metadata": {}, "source": [ "### Oversampling \n", "\n", "One of the advantages of using `PSFConvolvedModel2D` is its ability to sample models onto model images. Sometimes the models have regions that have to be oversampled to produce better estimates of the data. `PSFConvolvedModel2D` can oversample the entire model image or a specific pixel region of the image. The oversampling factor and region can be specified in the `oversample` keyword argument when wrapping an `astropy` model or during run time by setting the `PSFConvolvedModel2D.oversample` attribute. \n" ] }, { "cell_type": "markdown", "id": "09315692", "metadata": {}, "source": [ "**Disable Oversampling (Defailt)**\n", "\n", "To disable oversampling, set the `oversampling` argument or attribute to `None`" ] }, { "cell_type": "code", "execution_count": null, "id": "5b1918ae", "metadata": {}, "outputs": [], "source": [ "# Disable Oversampling\n", "oversample = None " ] }, { "cell_type": "markdown", "id": "70d6d88a", "metadata": {}, "source": [ "**Oversample Entire Model Image**\n", "\n", "To oversample the image by a factor, you can pass a single integer value. For example:" ] }, { "cell_type": "code", "execution_count": null, "id": "f7d7d660", "metadata": {}, "outputs": [], "source": [ "# Oversample the entire image by a factor of 4\n", "oversample = 4" ] }, { "cell_type": "markdown", "id": "a3565e79", "metadata": {}, "source": [ "**Oversample a Fixed Region**\n", "\n", "To oversample a fixed region of finite size, specify the center pixel, the length of the square region and thee oversampling factor. This means passing a tuple of `(center_x, center_y, box_length, oversample_factor)`. For example:" ] }, { "cell_type": "code", "execution_count": null, "id": "1980fae8", "metadata": {}, "outputs": [], "source": [ "# Replace the pixel values in a box of \n", "# length 20 cented at (x=50, y=60) with a box of \n", "# the same size that has been oversampled by a factor of 5 \n", "# i.e (x=50 y=60, box_length=20, oversample_factor=5)\n", "\n", "oversample = (50, 60, 20, 5)" ] }, { "cell_type": "markdown", "id": "ceac6bad", "metadata": {}, "source": [ "**Oversample a Moving Region**\n", "\n", "If the model is being fit, the center of the model is likely to move around. To account for this, we can specify the names of the model parameters that define the center of the box that we are interested in oversampling as strings. This means passing a tuple of `(model_param_x, model_param_y, box_length, oversample_factor)`. For example:" ] }, { "cell_type": "code", "execution_count": null, "id": "b7a45dc2", "metadata": {}, "outputs": [], "source": [ "# Replace the pixel values in a box of \n", "# length 20 cented at (x=model.x_0, y=model.y_0) with a box of \n", "# the same size that has been oversampled by a factor of 5 \n", "# i.e (model.x_0, model.y_0, box_length=20, oversample_factor=5)\n", "\n", "oversample = ('x_0', 'y_0', 20, 5)" ] }, { "cell_type": "markdown", "id": "cdd76eb0", "metadata": {}, "source": [ "### Oversampled PSF \n", "\n", "The PSF can have intricate details and variations that are not well-captured if we simply sample at the same rate as the data image. \n", "This is where the concept of an oversampled PSF comes into play.\n", "An oversampled PSF is essentially a higher-resolution representation of the PSF, capturing its subtle variations with more detail. \n", "This is beneficial because, during convolution, these details interact with the underlying data, ensuring a more accurate representation of the light distribution.\n", "`PSFConvolvedModel2D` facilitates this by allowing users to specify an oversampled PSF alongside the model. \n", "The `psf_oversample` keyword argument, or attribute, controls the oversampling factor of the PSF. \n", "It's essential to remember that when working with both oversampled models and PSFs, compatibility is key. \n", "The `PSFConvolvedModel2D` class ensures that the model's oversampling rate (oversample) is always an integer multiple of the PSF's oversampling rate (`psf_oversample`). " ] }, { "cell_type": "code", "execution_count": null, "id": "c47b65b1", "metadata": {}, "outputs": [], "source": [ "# The star image PSF is at the \n", "# same resolution as the data\n", "psf_oversample = 1" ] }, { "cell_type": "markdown", "id": "f8a178f5", "metadata": {}, "source": [ "### Create PetroFit Model\n", "\n", "Now that we have an `astropy` model, PSF and oversampling rule, we can create a `PSFConvolvedModel2D` model as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "ef24ad08", "metadata": {}, "outputs": [], "source": [ "psf_sersic_model = pf.PSFConvolvedModel2D(sersic_model, psf=PSF, oversample=4, psf_oversample=1)" ] }, { "cell_type": "markdown", "id": "1fe26506", "metadata": {}, "source": [ " The `PSFConvolvedModel2D` etherates all of the parameters, fixed-parameter rules and parameter bounds from the input `astropy` model. Notice that a new parameter, `psf_pa` is added to enable PSF rotation. " ] }, { "cell_type": "code", "execution_count": null, "id": "31e610be", "metadata": {}, "outputs": [], "source": [ "print(psf_sersic_model.param_names)" ] }, { "cell_type": "code", "execution_count": null, "id": "6e4bd2a9", "metadata": {}, "outputs": [], "source": [ "print(psf_sersic_model.bounds)" ] }, { "cell_type": "markdown", "id": "a98e6f6e", "metadata": {}, "source": [ "### PSF Rotation \n", "\n", "`PSFConvolvedModel2D` can to rotate the PSF image until an optimal rotation angle is found. This is useful for when the PSF comes from a dataset where the orientation of the diffraction spikes are not the same as the image being fit. `psf_pa` is in degrees.\n", "\n", "To restrict the bounds of the rotation or disable the PSF rotation, you can set the psf_pa to fixed:" ] }, { "cell_type": "code", "execution_count": null, "id": "8014d5a1", "metadata": {}, "outputs": [], "source": [ "# Limit PSF rotation to -5 to 5 degrees\n", "psf_sersic_model.bounds['psf_pa'] = (-5, 5)\n", "\n", "# To disable the PSF rotation, \n", "# you can set the psf_pa to fixed.\n", "psf_sersic_model.fixed['psf_pa'] = True" ] }, { "cell_type": "markdown", "id": "63e5eaf2", "metadata": {}, "source": [ "### Accessing the Underlying Model \n", "\n", "At any point, a copy of the input model with the same parameter values as the corresponding `PSFConvolvedModel2D` can be accessed using the `PSFConvolvedModel2D.model` attribute:" ] }, { "cell_type": "code", "execution_count": null, "id": "68426b46", "metadata": {}, "outputs": [], "source": [ "psf_sersic_model.model" ] }, { "cell_type": "markdown", "id": "fc16830d", "metadata": {}, "source": [ "### Visualize Inital Guess Model\n", "\n", "Here we visualize the inital guess model using the `plot_fit` function:" ] }, { "cell_type": "code", "execution_count": null, "id": "bdfaf554", "metadata": {}, "outputs": [], "source": [ "pf.plot_fit(psf_sersic_model, image.data, vmax=vmax, vmin=vmin, figsize=[3*6, 6])\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "9abeba8b", "metadata": {}, "source": [ "Looks like we'd better fit this model to optimize its paramters..." ] }, { "cell_type": "markdown", "id": "f53b0ee0", "metadata": {}, "source": [ "## Fitting Models\n", "\n", "PetroFit uses the Levenberg-Marquardt, Trust Region Reflective algorithm, and linear least-squares algorithms to fit parametrized models. To achieve this, it uses `astropy` fitting and provides wrappers to fit models to images. One such function is `fit_model`, which takes any `Fittable2DModel` model and an image to fit, and returns a fitted copy of the model and the `fit_info` dictionary. If the image to be fit contains pixels that are set to `np.nan`, those pixels are ignored by the fitter. The `fit_model` function also allows us to define parameters, such as ` maxiter`, for the `astropy` fitter.\n", "\n", "Before we fit the image, we compute the weights of each pixel using rms data as follows (please note that weights are optional and set to `None` by defualt):" ] }, { "cell_type": "code", "execution_count": null, "id": "646afb23", "metadata": {}, "outputs": [], "source": [ "fitting_weights = 1 / rms" ] }, { "cell_type": "markdown", "id": "aff6a5c9", "metadata": {}, "source": [ "To fit the galaxy we prepared with the `PSFConvolvedModel2D` we constructed, we call the `fit_model` as follows: " ] }, { "cell_type": "code", "execution_count": null, "id": "990fe9a4", "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "fitted_model, fit_info = pf.fit_model(\n", " image.data, psf_sersic_model,\n", " weights=fitting_weights,\n", " calc_uncertainties=True,\n", " maxiter=10000,\n", " epsilon=1.4901161193847656e-08,\n", " acc=1e-09,\n", ")" ] }, { "cell_type": "markdown", "id": "58e7d19c", "metadata": {}, "source": [ "That’s it! The retuned `fitted_model` is a copy of the input model (`psf_sersic_model`) but with the optimized parameter values. We can inspect the parameters of any `astropy` model using the `print_model_params`:" ] }, { "cell_type": "code", "execution_count": null, "id": "544b714d", "metadata": {}, "outputs": [], "source": [ "pf.print_model_params(fitted_model)" ] }, { "cell_type": "markdown", "id": "274ac604", "metadata": {}, "source": [ "### Paramter Errors\n", "When `calc_uncertainties` is enabled in the `fit_model` function, Astropy's fitter calculates the parameter uncertainties using the covariance matrix. \n", "To extract the standard deviation of the parameters, given that the covariance matrix is available:" ] }, { "cell_type": "code", "execution_count": null, "id": "c308b609", "metadata": {}, "outputs": [], "source": [ "# covariance matrix dict:\n", "fitted_model.cov_matrix" ] }, { "cell_type": "code", "execution_count": null, "id": "cf0b9cd4", "metadata": {}, "outputs": [], "source": [ "param_stds = fitted_model.stds\n", "for param, std in zip(param_stds.param_names, param_stds.stds):\n", " print(\"{:<10} {}\".format(param, std))" ] }, { "cell_type": "markdown", "id": "92341da0", "metadata": {}, "source": [ "## Generate Model Image\n", "\n", "To generate a model image we use the `plot_fit` function. The function, given a 2D model and fitted image, converts the model into a model-image we can visualize and manipulate. " ] }, { "cell_type": "code", "execution_count": null, "id": "73fd5f44", "metadata": { "scrolled": false }, "outputs": [], "source": [ "pf.plot_fit(fitted_model, image.data, vmax=vmax, vmin=vmin, figsize=[3*6, 6])\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.5" } }, "nbformat": 4, "nbformat_minor": 5 }