{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global background esimtation\n",
    "\n",
    "`phoptic` uses `photutils` to handle two-dimensional image backgrounds. In this notebook, I will demonstrate how to define backgrounds for use with `phoptic`, as well as explain `phoptic`'s default behaviour when no background is specified.\n",
    "\n",
    "## Test Image\n",
    "\n",
    "First thing's first, let's create and open an image so we can compute its background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:53.299710Z",
     "iopub.status.busy": "2026-01-23T11:18:53.299529Z",
     "iopub.status.idle": "2026-01-23T11:18:56.554364Z",
     "shell.execute_reply": "2026-01-23T11:18:56.554143Z"
    }
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import phoptic\n",
    "\n",
    "out_dir = Path('out')\n",
    "\n",
    "phoptic.generate_observations(\n",
    "    out_directory=out_dir / 'data',\n",
    "    n_images=5,\n",
    "    circular_aperture=False,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:56.555336Z",
     "iopub.status.busy": "2026-01-23T11:18:56.555156Z",
     "iopub.status.idle": "2026-01-23T11:18:56.570361Z",
     "shell.execute_reply": "2026-01-23T11:18:56.570169Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "files = os.listdir(out_dir / 'data')\n",
    "\n",
    "with fits.open(out_dir / 'data' / files[0]) as hdul:\n",
    "    print(repr(hdul[0].header))\n",
    "    image = np.array(hdul[0].data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:56.571049Z",
     "iopub.status.busy": "2026-01-23T11:18:56.570978Z",
     "iopub.status.idle": "2026-01-23T11:18:56.660971Z",
     "shell.execute_reply": "2026-01-23T11:18:56.660719Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.visualization import simple_norm\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "im = ax.imshow(\n",
    "    image,\n",
    "    norm=simple_norm(image, stretch=\"sqrt\"),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"Greys\",\n",
    "    )\n",
    "\n",
    "ax.set_xlabel(\"X\")\n",
    "ax.set_ylabel(\"Y\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Default Background\n",
    "\n",
    "`phoptic`'s default background estimator is the default `Background2D()` estimator from `photutils` with some added convenience tailored to OPTICAM. Let's look at the background image produced by the default estimator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:56.662074Z",
     "iopub.status.busy": "2026-01-23T11:18:56.661966Z",
     "iopub.status.idle": "2026-01-23T11:18:56.800481Z",
     "shell.execute_reply": "2026-01-23T11:18:56.800249Z"
    }
   },
   "outputs": [],
   "source": [
    "default_background = phoptic.DefaultBackground(\n",
    "    box_size=image.shape[0] // 32,  # phoptic sets this parameter to image.shape[0] // 32 by default\n",
    "    )\n",
    "\n",
    "bkg = default_background(image)  # compute the background\n",
    "bkg_image = bkg.background  # get the background image\n",
    "\n",
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "im = ax.imshow(\n",
    "    bkg_image,\n",
    "    norm=simple_norm(bkg_image, stretch=\"sqrt\"),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"Greys\",\n",
    "    )\n",
    "\n",
    "fig.colorbar(im)\n",
    "\n",
    "ax.set_xlabel(\"X\")\n",
    "ax.set_ylabel(\"Y\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the background is around 100 counts, which is the expected value. We can also look at the background mesh to see if any regions of the image have been excluded by the background estimator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:56.801326Z",
     "iopub.status.busy": "2026-01-23T11:18:56.801232Z",
     "iopub.status.idle": "2026-01-23T11:18:57.146300Z",
     "shell.execute_reply": "2026-01-23T11:18:57.145981Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "im = ax.imshow(\n",
    "    image,\n",
    "    norm=simple_norm(image, stretch=\"sqrt\"),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"Greys\",\n",
    "    )\n",
    "bkg.plot_meshes(\n",
    "    outlines=True,\n",
    "    marker='.',\n",
    "    color='r',\n",
    "    alpha=0.3,\n",
    "    ax=ax,\n",
    "    )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default size of the background \"pixels\" for `phoptic.DefaultBackground` is the width of the image divided by 16. This value is generally good across a range of observing conditions, but it can, of course, be changed on a case-by-case basis. The `photutils` documentation suggests setting `box_size` to a value which is small, but larger than the typical source size. For these simulated data, we could probably get away with a smaller `box_size`, but such fine-tunings are left to the user.\n",
    "\n",
    "We can see that brighter sources have been clipped, which fainter sources haven't. A custom background estimator may therefore be useful if more aggressive clipping is required.\n",
    "\n",
    "Let's subtract the background and compare the histograms of pixel values before and after removing the background from the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.147681Z",
     "iopub.status.busy": "2026-01-23T11:18:57.147581Z",
     "iopub.status.idle": "2026-01-23T11:18:57.206631Z",
     "shell.execute_reply": "2026-01-23T11:18:57.206411Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "ax.hist(\n",
    "    image.flatten(),\n",
    "    bins=100,\n",
    "    histtype=\"step\",\n",
    "    label=\"Raw\",\n",
    "    )\n",
    "\n",
    "ax.hist(\n",
    "    (image - bkg_image).flatten(),\n",
    "    bins=100,\n",
    "    histtype=\"step\",\n",
    "    label=\"Background subtracted\",\n",
    "    )\n",
    "\n",
    "ax.set_xlabel(\"Pixel value [ADU]\")\n",
    "ax.set_ylabel(\"Number of pixels\")\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the background has been massively reduced. After subtracting the background, the pixel values appear approximately Gaussian distributed about zero.\n",
    "\n",
    "Using `phoptic`'s default background should be \"good enough\" most of the time (at least for OPTICAM-MX data). However, in some cases, we may see better results if we implement a custom background estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom Backgrounds\n",
    "\n",
    "Let's now define a custom background estimator. Custom background estimators must be `callable`s that take an image (`numpy` `NDArray`) as input and return a `photutils.background.Background2D` instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.207438Z",
     "iopub.status.busy": "2026-01-23T11:18:57.207343Z",
     "iopub.status.idle": "2026-01-23T11:18:57.209422Z",
     "shell.execute_reply": "2026-01-23T11:18:57.209240Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy.stats import SigmaClip\n",
    "from photutils.background import Background2D, BiweightLocationBackground, BiweightScaleBackgroundRMS\n",
    "\n",
    "class CustomBackground:\n",
    "    \n",
    "    def __init__(\n",
    "        self,\n",
    "        box_size,\n",
    "        bkg_estimator=BiweightLocationBackground(),\n",
    "        bkg_rms_estimator=BiweightScaleBackgroundRMS(),\n",
    "        ):\n",
    "        \n",
    "        self.box_size = box_size\n",
    "        self.bkg_estimator=bkg_estimator\n",
    "        self.bkg_rms_estimator=bkg_rms_estimator\n",
    "    \n",
    "    def __call__(\n",
    "        self,\n",
    "        image,\n",
    "        ):\n",
    "        \n",
    "        return Background2D(\n",
    "            image,\n",
    "            self.box_size,\n",
    "            bkg_estimator=self.bkg_estimator,\n",
    "            bkg_rms_estimator=self.bkg_rms_estimator,\n",
    "            sigma_clip=SigmaClip(\n",
    "                sigma=3,\n",
    "                maxiters=None,\n",
    "                ),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, I have defined `CustomBackground` as a class that requires a `box_size` parameter, just like `phoptic.DefaultBackground`. Unlike `phoptic.DefaultBackground`, however, `CustomBackground` takes the background and background RMS estimators as additional arguments. By default, `CustomBackground` uses the `photutils.background.BiweightLocationBackground` background estimator and `photutils.background.BiweightScaleBackgroundRMS` background RMS estimator, while `phoptic.DefaultBackground` is limited to the `photutils.background.SExtractorBackground` and `photutils.background.StdBackgroundRMS` estimators. Of course, custom background and background RMS estimators can also be defined, though this is more advanced (see https://photutils.readthedocs.io/en/stable/user_guide/background.html#d-background-and-noise-estimation for more details). A final crucial difference of our custom background estimator is the use of a custom `astropy.stats.SigmaClip()` instance with `maxiters=None`, meaning sigma clipping will be run until convergence is reached. In contrast, `phoptic.DefaultBackground` performs no more than 10 sigma clipping iterations.\n",
    "\n",
    "We could have equally defined a custom background estimator using a function instead of a class, which is much simpler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.210022Z",
     "iopub.status.busy": "2026-01-23T11:18:57.209950Z",
     "iopub.status.idle": "2026-01-23T11:18:57.211447Z",
     "shell.execute_reply": "2026-01-23T11:18:57.211250Z"
    }
   },
   "outputs": [],
   "source": [
    "%%writefile out/custom_routines.py\n",
    "from astropy.stats import SigmaClip\n",
    "from photutils.background import Background2D, BiweightLocationBackground, BiweightScaleBackgroundRMS\n",
    "\n",
    "\n",
    "def custom_background(\n",
    "    image,\n",
    "    box_size,\n",
    "    ):\n",
    "    \n",
    "    return Background2D(\n",
    "            image,\n",
    "            box_size,\n",
    "            bkg_estimator=BiweightLocationBackground(),\n",
    "            bkg_rms_estimator=BiweightScaleBackgroundRMS(),\n",
    "            sigma_clip=SigmaClip(\n",
    "                sigma=3,\n",
    "                maxiters=None,\n",
    "                ),\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to compatibility issues, routines defined inside of IPython notebooks cannot be used by `multiprocessing` (used by `phoptic`). To overcome this issue, we have used the `%%writefile` magic command to write our `custom_background()` function to a temporary module called `custom_routines.py` in our `out` directory. This is required in Python $\\geq$ 3.14.\n",
    "\n",
    "Let's now compare `CustomBackground` to `phoptic.DefaultBackground`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.212113Z",
     "iopub.status.busy": "2026-01-23T11:18:57.212043Z",
     "iopub.status.idle": "2026-01-23T11:18:57.318993Z",
     "shell.execute_reply": "2026-01-23T11:18:57.318743Z"
    }
   },
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "from out.custom_routines import custom_background\n",
    "\n",
    "# class:\n",
    "custom_background_class = CustomBackground(box_size=image.shape[0] // 32)\n",
    "\n",
    "# or, for our function, we want to fix the box_size parameter using functools.partial:\n",
    "custom_background_func = partial(\n",
    "    custom_background,\n",
    "    box_size=image.shape[0] // 32,  # set box size to image width / 32, similarly to phoptic's default behaviour\n",
    "    )\n",
    "\n",
    "custom_bkg_class_image = custom_background_class(image).background\n",
    "custom_bkg_func_image = custom_background_func(image).background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.320095Z",
     "iopub.status.busy": "2026-01-23T11:18:57.320012Z",
     "iopub.status.idle": "2026-01-23T11:18:57.533822Z",
     "shell.execute_reply": "2026-01-23T11:18:57.533570Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(\n",
    "    ncols=3,\n",
    "    tight_layout=True,\n",
    "    figsize=(15, 5),\n",
    "    )\n",
    "\n",
    "im = ax[0].imshow(\n",
    "    bkg_image,\n",
    "    norm=simple_norm(bkg_image, stretch=\"sqrt\"),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"Greys\",\n",
    "    )\n",
    "ax[0].set_title(\"Default background\")\n",
    "\n",
    "im = ax[1].imshow(\n",
    "    custom_bkg_class_image,\n",
    "    norm=simple_norm(bkg_image, stretch=\"sqrt\"),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"Greys\",\n",
    "    )\n",
    "ax[1].set_title(\"Custom background (class)\")\n",
    "\n",
    "im = ax[2].imshow(\n",
    "    custom_bkg_func_image,\n",
    "    norm=simple_norm(bkg_image, stretch=\"sqrt\"),\n",
    "    origin=\"lower\",\n",
    "    cmap=\"Greys\",\n",
    "    )\n",
    "ax[2].set_title(\"Custom background (func)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the default background and our custom backgrounds are subtly different (note that they also share the same normalisation). Let's compare the histograms between these background estimators to better quantify the differences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.534883Z",
     "iopub.status.busy": "2026-01-23T11:18:57.534787Z",
     "iopub.status.idle": "2026-01-23T11:18:57.610447Z",
     "shell.execute_reply": "2026-01-23T11:18:57.610252Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(tight_layout=True)\n",
    "\n",
    "ax.hist(\n",
    "    bkg_image.flatten(),\n",
    "    bins='auto',\n",
    "    histtype=\"step\",\n",
    "    label=\"Default background\",\n",
    "    )\n",
    "ax.hist(\n",
    "    custom_bkg_class_image.flatten(),\n",
    "    bins='auto',\n",
    "    histtype=\"step\",\n",
    "    label=\"Custom background (class)\",\n",
    "    )\n",
    "ax.hist(\n",
    "    custom_bkg_func_image.flatten(),\n",
    "    bins='auto',\n",
    "    histtype=\"step\",\n",
    "    label=\"Custom background (func)\",\n",
    "    ls='--',\n",
    "    )\n",
    "\n",
    "ax.set_xlabel(\"Pixel value [ADU]\")\n",
    "ax.set_ylabel(\"Number of pixels\")\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see that the default estimator is very different to our custom estimators, while both our class-based and function-based background estimators are, of course, identical.\n",
    "\n",
    "Let's see how to pass our custom background estimators to `phoptic.Reducer`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:18:57.611263Z",
     "iopub.status.busy": "2026-01-23T11:18:57.611176Z",
     "iopub.status.idle": "2026-01-23T11:19:02.204731Z",
     "shell.execute_reply": "2026-01-23T11:19:02.204426Z"
    }
   },
   "outputs": [],
   "source": [
    "reducer = phoptic.Reducer(\n",
    "    out_directory=out_dir / 'reduced' / 'custom_background_function',\n",
    "    data_directory=out_dir / 'data',\n",
    "    background=custom_background_func,  # pass custom background estimator\n",
    "    remove_cosmic_rays=True,\n",
    ")\n",
    "\n",
    "reducer.create_catalogs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, our custom background estimator infers an average background of $\\sim 100$ counts with an RMS of $\\sim 10$, both of which are the expected values. We could have equally used our class-based background estimator, similar to the implementation of the default background estimator, instead of the function-based one.\n",
    "\n",
    "Let's check how this compares to `phoptic`'s default background estimator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-01-23T11:19:02.205675Z",
     "iopub.status.busy": "2026-01-23T11:19:02.205587Z",
     "iopub.status.idle": "2026-01-23T11:19:06.504435Z",
     "shell.execute_reply": "2026-01-23T11:19:06.504163Z"
    }
   },
   "outputs": [],
   "source": [
    "reducer = phoptic.Reducer(\n",
    "    out_directory=out_dir / 'reduced' / 'default_background',\n",
    "    data_directory=out_dir / 'data',\n",
    "    remove_cosmic_rays=True,\n",
    ")\n",
    "\n",
    "reducer.create_catalogs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsurprisingly, we again find that the default background estimator infers an average background of $\\sim 100$ counts with an RMS of $\\sim 10$.\n",
    "\n",
    "That concludes the backgrounds tutorial for `phoptic`! In most cases, `phoptic.DefaultBackground` should be \"good enough\", but we have seen how custom background estimators could be implemented if necessary. For more details on implementing custom background estimators, I refer to the excellent `photutils` documentation: https://photutils.readthedocs.io/en/stable/user_guide/background.html#d-background-and-noise-estimation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opticam7",
   "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.14.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}