Skip to content

Commit

Permalink
WIP XDF segmentation
Browse files Browse the repository at this point in the history
  • Loading branch information
@mwcraig
mwcraig committed Apr 6, 2024
1 parent 5f982a1 commit 78b53b8
Showing 1 changed file with 165 additions and 47 deletions.
212 changes: 165 additions & 47 deletions notebooks/photometry/01.04.02-image-segmentation-XDF.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,86 @@
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "4a539fa2-f180-4b88-a620-da47c6af8270",
"metadata": {},
"source": [
"## Data for this notebook \n",
"\n",
"We will be manipulating Hubble eXtreme Deep Field (XDF) data, which was collected using the Advanced Camera for Surveys (ACS) on Hubble between 2002 and 2012. The image we use here is the result of 1.8 million seconds (500 hours!) of exposure time, and includes some of the faintest and most distant galaxies that had ever been observed. \n",
"\n",
"*The methods demonstrated here are also available in the [`photutils.detection` documentation](http://photutils.readthedocs.io/en/stable/detection.html).*\n",
"\n",
"*The original authors of this notebook were Lauren Chambers, Erik Tollerud and Tom Wilson.*\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7e90842d-8f03-408e-8622-532d1fab18c7",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"\n",
"from astropy.nddata import CCDData\n",
"from astropy.stats import sigma_clipped_stats, SigmaClip\n",
"import astropy.units as u\n",
"from astropy.visualization import ImageNormalize, LogStretch\n",
"\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.ticker import LogLocator\n",
"\n",
"import numpy as np\n",
"\n",
"from photutils.aperture import EllipticalAperture\n",
"from photutils.background import Background2D, MeanBackground\n",
"from photutils.centroids import centroid_2dg\n",
"from photutils.detection import find_peaks, DAOStarFinder, IRAFStarFinder\n",
"from photutils.segmentation import detect_sources, SourceCatalog\n",
"\n",
"# Show plots in the notebook\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7fe76ad2-51c6-48fd-b91e-541a1ea9bd90",
"metadata": {},
"outputs": [],
"source": [
"xdf_image = CCDData.read('hlsp_xdf_hst_acswfc-60mas_hudf_f435w_v1_sci.fits')\n",
"# Define the mask\n",
"mask = xdf_image.data == 0\n",
"xdf_image.mask = mask"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "82f387f7-66af-4be6-bd2f-4f54f87ba6ed",
"metadata": {},
"outputs": [],
"source": [
"mean, median, std = sigma_clipped_stats(xdf_image.data, sigma=3.0, maxiters=5, mask=xdf_image.mask)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2b4110f0-a867-42fb-8eac-885a4883c4d6",
"metadata": {},
"outputs": [],
"source": [
"plt.style.use('../photutils_notebook_style.mplstyle')"
]
},
{
"cell_type": "markdown",
"id": "21c15890-d5d2-4f5a-a1dd-7795e0f609a7",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"metadata": {},
"source": [
"Beyond traditional source detection methods, an additional option for identifying sources in image data is a process called **image segmentation**. This method identifies and labels contiguous (connected) objects within an image. \n",
"\n",
Expand Down Expand Up @@ -48,7 +122,6 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import detect_sources\n",
"from photutils.utils import make_random_cmap"
]
},
Expand Down Expand Up @@ -80,8 +153,14 @@
"fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(12, 6))\n",
"plt.tight_layout()\n",
"\n",
"# Plot the data\n",
"# Set up the normalization and colormap\n",
"norm_image = ImageNormalize(vmin=1e-4, vmax=5e-2, stretch=LogStretch(), clip=False)\n",
"cmap = plt.get_cmap('viridis')\n",
"cmap.set_bad('white') # Show masked data as white\n",
"\n",
"# Plot the original data\n",
"fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped),\n",
"fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image),\n",
" norm=norm_image, cmap=cmap)\n",
"ax1.set_ylabel('Y (pixels)')\n",
"ax1.set_title('Original Data')\n",
Expand Down Expand Up @@ -117,23 +196,39 @@
"For a closer look, let's see what the sources we found with `find_peaks` look like in this segmentation map:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b1d61489-3595-46b9-aaaf-57f11ac2126d",
"metadata": {},
"outputs": [],
"source": [
"with warnings.catch_warnings(action=\"ignore\"):\n",
" sources_findpeaks = find_peaks(xdf_image.data - median, mask=xdf_image.mask,\n",
" threshold=20. * std, box_size=29,\n",
" centroid_func=centroid_2dg)\n",
"print(sources_findpeaks)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fd673ddd-bcf4-4e49-8281-a0bbc75ca080",
"metadata": {},
"outputs": [],
"source": [
"fig, axs = plt.subplots(3,3, figsize=(3, 3))\n",
"fig, axs = plt.subplots(3,3, figsize=(8, 8))\n",
"plt.subplots_adjust(wspace=0.1, hspace=0.1)\n",
"\n",
"cutout_size = 20\n",
"\n",
"srcs = np.random.permutation(sources_findpeaks)[:axs.size]\n",
"rng = np.random.default_rng(seed=8675309)\n",
"srcs = rng.permutation(sources_findpeaks)[:axs.size]\n",
"for ax, src in zip(axs.ravel(), srcs):\n",
" slc = (slice(int(src['y_peak'] - cutout_size), int(src['y_peak'] + cutout_size)),\n",
" slice(int(src['x_peak'] - cutout_size), int(src['x_peak'] + cutout_size)))\n",
" ax.imshow(segm.data[slc], cmap=rand_cmap, vmin=1, vmax=len(sources_findpeaks))\n",
" ax.imshow(xdf_image.data[slc], cmap=cmap, norm=norm_image)\n",
" ax.imshow(segm.data[slc], cmap=rand_cmap, vmin=1, vmax=len(sources_findpeaks), alpha=0.5)\n",
" src_id = np.where((sources_findpeaks['x_peak'] == src['x_peak']) &\n",
" (sources_findpeaks['y_peak'] == src['y_peak']))[0][0]\n",
" ax.text(2, 2, str(src_id), color='w', va='top')\n",
Expand All @@ -142,23 +237,19 @@
]
},
{
"cell_type": "markdown",
"id": "40fd8f04-d36e-4d43-b157-60c7b2ee6a3c",
"cell_type": "code",
"execution_count": null,
"id": "1176098d-8a00-402d-9220-a1c54044c0c4",
"metadata": {},
"source": [
"<div class=\"alert alert-block alert-info\">\n",
"\n",
"<h3>Exercises:</h3><br>\n",
"\n",
"Recompute the <code>SegmentationImage</code>, but alter the threshold and the minimum number of pixels in a source. How does changing the threshold affect the results? What about changing the number of pixels?\n",
"\n",
"</div>"
]
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "3e15e7d8-3aa9-4b75-ab03-88d3d6461eb5",
"metadata": {},
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"#### Analyzing `source_properties`\n",
"\n",
Expand All @@ -169,24 +260,14 @@
"However, perhaps the most powerful aspect of the `SegmentationImage` is the ability to create a catalog using `source_properties` to measure the centroids, photometry, and morphology of the detected objects:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6c3b45a9-9bf5-41d1-809b-de6efd8088ae",
"metadata": {},
"outputs": [],
"source": [
"from photutils import source_properties"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "31715b6e-7cc4-47d5-be5b-b49e08740a36",
"metadata": {},
"outputs": [],
"source": [
"catalog = source_properties(xdf_image.data, segm)\n",
"catalog = SourceCatalog(xdf_image.data, segm)\n",
"table = catalog.to_table()\n",
"print(table)\n",
"print(table.colnames)"
Expand All @@ -202,16 +283,6 @@
"We can use this information to create isophotal ellipses for each identified source. These ellipses can also later be used as photometric apertures!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "896f10cb-ebc6-417b-b412-65c726c7c0af",
"metadata": {},
"outputs": [],
"source": [
"from photutils import EllipticalAperture"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -225,10 +296,10 @@
"# Create the apertures\n",
"apertures = []\n",
"for obj in catalog:\n",
" position = (obj.xcentroid.value, obj.ycentroid.value)\n",
" a = obj.semimajor_axis_sigma.value * r\n",
" b = obj.semiminor_axis_sigma.value * r\n",
" theta = obj.orientation.value\n",
" position = (obj.xcentroid, obj.ycentroid)\n",
" a = obj.semimajor_sigma.value * r\n",
" b = obj.semiminor_sigma.value * r\n",
" theta = obj.orientation\n",
" apertures.append(EllipticalAperture(position, a, b, theta=theta))"
]
},
Expand All @@ -243,7 +314,38 @@
"fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
"\n",
"# Plot the data\n",
"fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped),\n",
"fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image),\n",
" norm=norm_image, cmap=cmap)\n",
"\n",
"# Plot the apertures\n",
"for aperture in apertures:\n",
" aperture.plot(color='red', lw=1, alpha=0.7, axes=ax1)\n",
"\n",
"# Define the colorbar\n",
"cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
"labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
"cbar.ax.set_yticklabels(labels)\n",
"\n",
"# Define labels\n",
"cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')),\n",
" rotation=270, labelpad=30)\n",
"ax1.set_xlabel('X (pixels)')\n",
"ax1.set_ylabel('Y (pixels)')\n",
"ax1.set_title('Segmentation Image Apertures')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f3551192-ea06-452e-a5ab-57c322e8c4fc",
"metadata": {},
"outputs": [],
"source": [
"# Set up the figure with subplots\n",
"fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n",
"\n",
"# Plot the data\n",
"fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image),\n",
" norm=norm_image, cmap=cmap)\n",
"\n",
"# Plot the apertures\n",
Expand All @@ -260,8 +362,24 @@
" rotation=270, labelpad=30)\n",
"ax1.set_xlabel('X (pixels)')\n",
"ax1.set_ylabel('Y (pixels)')\n",
"\n",
"top = 1125\n",
"left = 2350\n",
"cutout_size = 500\n",
"\n",
"ax1.set_xlim(left, left + cutout_size)\n",
"ax1.set_ylim(top, top + cutout_size)\n",
"\n",
"ax1.set_title('Segmentation Image Apertures')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b7fb1750-55d3-4eaa-943c-475d5a2db9ad",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand All @@ -280,7 +398,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.3"
"version": "3.11.8"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 78b53b8

Please sign in to comment.