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 7, 2024
1 parent 5f982a1 commit 69a7d17
Showing 1 changed file with 222 additions and 63 deletions.
285 changes: 222 additions & 63 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,87 @@
"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.convolution import convolve\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, make_2dgaussian_kernel, 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 +123,6 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import detect_sources\n",
"from photutils.utils import make_random_cmap"
]
},
Expand All @@ -61,10 +135,12 @@
"source": [
"# Define threshold and minimum object size\n",
"threshold = 5. * std\n",
"npixels = 15\n",
"\n",
"npixels = 10\n",
"fwhm = 5\n",
"kernel = make_2dgaussian_kernel(fwhm, size=5)\n",
"convolved_data = convolve(xdf_image.data, kernel)\n",
"# Create a segmentation image\n",
"segm = detect_sources(xdf_image.data, threshold, npixels)\n",
"segm = detect_sources(convolved_data, threshold, npixels)\n",
"\n",
"print('Found {} sources'.format(segm.max_label))"
]
Expand All @@ -80,8 +156,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 @@ -120,39 +202,56 @@
{
"cell_type": "code",
"execution_count": null,
"id": "fd673ddd-bcf4-4e49-8281-a0bbc75ca080",
"id": "b1d61489-3595-46b9-aaaf-57f11ac2126d",
"metadata": {},
"outputs": [],
"source": [
"fig, axs = plt.subplots(3,3, figsize=(3, 3))\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",
"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",
" 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",
" ax.set_xticks([])\n",
" ax.set_yticks([])"
"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",
"\n",
"# Add an ID column to the find_peaks result\n",
"sources_findpeaks[\"source_id\"] = np.arange(len(sources_findpeaks))"
]
},
{
"cell_type": "markdown",
"id": "40fd8f04-d36e-4d43-b157-60c7b2ee6a3c",
"cell_type": "code",
"execution_count": null,
"id": "1176098d-8a00-402d-9220-a1c54044c0c4",
"metadata": {},
"outputs": [],
"source": [
"<div class=\"alert alert-block alert-info\">\n",
"fig = plt.figure(figsize=(8, 2.5))\n",
"sub_figs = fig.subfigures(nrows=2, ncols=4, wspace=0.01, hspace=0.01)\n",
"\n",
"<h3>Exercises:</h3><br>\n",
"cutout_size = 60\n",
"rng = np.random.default_rng(seed=8675309)\n",
"srcs = rng.permutation(sources_findpeaks)[:len(sub_figs.flatten())]\n",
"for a_fig, src in zip(sub_figs.flatten(), srcs):\n",
" sub_axs = a_fig.subplot_mosaic([[\"orig\", \"segm\"]], sharex=True, sharey=True, )\n",
" left = int(src['x_peak'] - cutout_size)\n",
" bottom = int(src['y_peak'] - cutout_size)\n",
" right = left + 2 * cutout_size\n",
" top = bottom + 2 * cutout_size\n",
" slc = (slice(bottom, top), slice(left, right))\n",
" sources_in_cutout = ((left < sources_findpeaks[\"x_centroid\"]) & (sources_findpeaks[\"x_centroid\"] < right) &\n",
" (bottom < sources_findpeaks[\"y_centroid\"]) & (sources_findpeaks[\"y_centroid\"] < top))\n",
" sub_axs[\"orig\"].imshow(xdf_image.data[slc], cmap=cmap, norm=norm_image, origin=\"lower\")\n",
" sub_axs[\"orig\"].scatter(sources_findpeaks[\"x_centroid\"][sources_in_cutout] - left, \n",
" sources_findpeaks[\"y_centroid\"][sources_in_cutout] - bottom, \n",
" marker=\"x\", color=\"orange\", s=100)\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",
" sub_axs[\"segm\"].imshow(segm.data[slc], cmap=rand_cmap, vmin=1, vmax=len(sources_findpeaks), origin=\"lower\")\n",
" src_id = src[\"source_id\"]\n",
" # ax.text(2, 2, str(src_id), color='w', va='top')\n",
" a_fig.suptitle(f\"Source {src_id}\", fontsize=10)\n",
"\n",
"</div>"
" for ax in sub_axs.values():\n",
" ax.set_xticks([])\n",
" ax.set_yticks([])\n",
"\n",
"fig.suptitle(\"Comparison of find_peaks and segmentation sources\", y=1.1);"
]
},
{
Expand All @@ -166,17 +265,7 @@
"\n",
"Individual objects within the segmentation map can be altered using methods such as `relabel` to change the labels of objects, `remove_labels` to remove objects, or `deblend_sources` to separating overlapping sources that were incorrectly labeled as one source.\n",
"\n",
"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"
"However, perhaps the most powerful aspect of the `SegmentationImage` is the ability to create a catalog using [SourceCatalog](https://photutils.readthedocs.io/en/stable/api/photutils.segmentation.SourceCatalog.html#photutils.segmentation.SourceCatalog) to measure the centroids, photometry, and morphology of the detected objects:"
]
},
{
Expand All @@ -186,7 +275,7 @@
"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 @@ -199,17 +288,7 @@
"source": [
"#### Creating apertures from segmentation data\n",
"\n",
"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"
"We can use this information to create isophotal ellipses for each identified source. These ellipses can also later be used as photometric apertures."
]
},
{
Expand All @@ -225,10 +304,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,24 +322,104 @@
"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",
"cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator())\n",
"\n",
"def format_colorbar(bar):\n",
" # Add minor tickmarks\n",
" bar.ax.yaxis.set_minor_locator(LogLocator(subs=range(1, 10)))\n",
"\n",
" # Force the labels to be displayed as powers of ten and only at exact powers of ten\n",
" bar.ax.set_yticks([1e-4, 1e-3, 1e-2])\n",
" labels = [f'$10^{{{pow:.0f}}}$' for pow in np.log10(bar.ax.get_yticks())]\n",
" bar.ax.set_yticklabels(labels)\n",
"\n",
"format_colorbar(cbar)\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": "19c62f54-3040-4ec8-9c12-f02d271e3ac4",
"metadata": {},
"outputs": [],
"source": [
"daofind = DAOStarFinder(fwhm=fwhm, threshold=threshold)\n",
"sources_dao = daofind(np.ma.masked_where(xdf_image.mask, xdf_image))"
]
},
{
"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",
"apertures[0].plot(color='red', lw=2, alpha=1, axes=ax1, label=\"Image segmentation\")\n",
"for aperture in apertures[1:]:\n",
" aperture.plot(color='red', lw=2, alpha=1, axes=ax1)\n",
"\n",
"# Add find_peaks sources\n",
"ax1.scatter(sources_findpeaks[\"x_centroid\"], sources_findpeaks[\"y_centroid\"], \n",
" marker=\"x\", color=\"orange\", s=200, lw=2, alpha=0.7, label=\"find_peaks\")\n",
"\n",
"# Add find_peaks sources\n",
"ax1.scatter(sources_dao[\"xcentroid\"], sources_dao[\"ycentroid\"], \n",
" marker=\"s\", facecolor=\"none\", edgecolor=\"cyan\", s=200, alpha=0.5, label=\"DAOFind\")\n",
"\n",
"# Define the colorbar\n",
"cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator())\n",
"\n",
"format_colorbar(cbar)\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')"
"\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.legend(ncols=3, loc='lower center', bbox_to_anchor=(0.5, 1))\n",
"\n",
"ax1.set_title('Source detection method comparison', y=1.05);"
]
},
{
"cell_type": "markdown",
"id": "c0b2a1b2-57c0-49a5-8924-5372cd43613c",
"metadata": {},
"source": [
"## Summary\n",
"\n",
"The image above nicely captures the differences between the source detection methods we have discussed. [`DAOStarFinder`](https://photutils.readthedocs.io/en/stable/api/photutils.detection.DAOStarFinder.html#photutils.detection.DAOStarFinder) does an excellent job detecting the star-like sources in the image but not the more extended sources. [`find_peaks`](https://photutils.readthedocs.io/en/stable/api/photutils.detection.find_peaks.html#photutils.detection.find_peaks) detects both extended and point-like sources, but misses several sources. Image segmentation does an excellent job detecting sources, and could detect the smaller sources in this image if we reduced the minimum number of pixels in each source."
]
}
],
Expand All @@ -280,7 +439,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 69a7d17

Please sign in to comment.