2.3.2 Saturated pixels and Segment’s clumps ¶

A constant-depth (flat) image was created in the previous section (Preparing input for extended PSF). As explained in Overview of the PSF scripts, an important step when building the PSF is to mask other sources in the image. Therefore, before going onto selecting stars, let’s detect all significant signal, and identify the clumps of background objects over the wings of the extended PSF.

There is a problem however: the saturated pixels of the bright stars are going to cause problems in the segmentation phase. To see this problem, let’s make a \(1000\times1000\) crop around a bright star to speed up the test (and its solution). Afterwards we will apply the solution to the whole image.

$ astcrop flat/67510.fits --mode=wcs --widthinpix --width=1000 \
          --center=203.3916736,46.7968652 --output=saturated.fits
$ astnoisechisel saturated.fits --output=sat-nc.fits
$ astsegment sat-nc.fits --output=sat-seg.fits
$ astscript-fits-view sat-seg.fits

Have a look at the CLUMPS extension. You will see that instead of a single clump at the center of the bright star, we have many clumps! This has happened because of the saturated pixels! When saturation occurs, the sharp peak of the profile is lost (like cutting off the tip of a mountain to build a telescope!) and all saturated pixels get a noisy value close to the saturation level. To see this saturation noise run the last command again and in SAO DS9, set the “Scale” to “min max” and zoom into the center. You will see the noisy saturation pixels at the center of the star in red.

This noise-at-the-peak disrupts Segment’s assumption to expand clumps from a local maxima: each noisy peak is being treated as a separate local maxima and thus a separate clump. For more on how Segment defines clumps, see Section 3.2.1 and Figure 8 of Akhlaghi and Ichikawa 2015. To have the center identified as a single clump, we should mask these saturated pixels in a way that suites Segment’s non-parametric methodology.

First we need to find the saturation level! The saturation level is usually fixed for any survey or input data that you receive from a certain database, so you will usually have to do this only once (the first time you get data from that database). Let’s make a smaller crop of \(50\times50\) pixels around the star with the first command below. With the next command, please look at the crop with DS9 to visually understand the problem. You will see the saturated pixels as the noisy red pixels in the center of the image. A non-saturated star will have a single pixel as the maximum and will not have such a large area covered by a noisy constant value (find a few stars in the image and see for yourself). Visual and qualitative inspection of the process is very important for understanding the solution.

$ astcrop saturated.fits --mode=wcs --widthinpix --width=50 \
          --center=203.3916736,46.7968652 --output=sat-center.fits
$ astscript-fits-view sat-center.fits --ds9scale=minmax

To quantitatively identify the saturation level in this image, let’s have a look at the distribution of pixels with a value larger than 100 (above the noise level):

$ aststatistics sat-center.fits --greaterequal=100
Histogram:
 |*
 |*
 |*
 |*
 |*                                                              *
 |**                                                             *
 |***                                                           **
 |****                                                          **
 |******                                                        ****
 |********** *    *   *                                        ******
 |************************* ************ * ***  ******* *** ************
 |----------------------------------------------------------------------

The peak you see in the right end (larger values) of the histogram shows the saturated pixels (a constant level, with some scatter due to the large Poisson noise). If there was no saturation, the number of pixels should have decreased at increasing values; until reaching the maximum value of the profile in one pixel. But that is not the case here. Please try this experiment on a non-saturated (fainter) star to see what we mean.

If you still have not experimented on a non-saturated star, please stop reading this tutorial! Please open flat/67510.fits in DS9, select a fainter/smaller star and repeat the last three commands (with a different center). After you have confirmed the point above (visually, and with the histogram), please continue with the rest of this tutorial.

Finding the saturation level is easy with Statistics (by using the --lessthan option until the histogram becomes as expected: only decreasing). First, let’s try --lessthan=3000:

$ aststatistics sat-center.fits --greaterequal=100 --lessthan=3000
-------
Histogram:
 |*
 |*
 |*
 |*
 |*
 |**
 |***                                                                  *
 |****                                                                 *
 |*******                                                             **
 |*********** * *   *   *   *                            *  *       ****
 |************************* *  ***** *******  *****  ** ***** * ********
 |----------------------------------------------------------------------

We still see an increase in the histogram around 3000. Let’s try a threshold of 2500:

$ aststatistics sat-center.fits --greaterequal=100 --lessthan=2500
-------
Histogram:
 |*
 |*
 |**
 |**
 |**
 |**
 |****
 |*****
 |*********
 |*************  *   *  *   *
 |*********************************   ** ** ** *** **  * ****   ** *****
 |----------------------------------------------------------------------

The peak at the large end of the histogram has gone! But let’s have a closer look at the values (the resolution of an ASCII histogram is limited!). To do this, we will ask Statistics to save the histogram into a table with the --histogram option, then look at the last 20 rows:

$ aststatistics sat-center.fits --greaterequal=100 --lessthan=2500 \
                --histogram --output=sat-center-hist.fits
$ asttable sat-center-hist.fits --tail=20
2021.1849112701    1
2045.0495397186    0
2068.9141681671    1
2092.7787966156    1
2116.6434250641    0
2140.5080535126    0
2164.3726819611    0
2188.2373104095    0
2212.101938858     1
2235.9665673065    1
2259.831195755     2
2283.6958242035    0
2307.560452652     0
2331.4250811005    1
2355.289709549     1
2379.1543379974    1
2403.0189664459    2
2426.8835948944    1
2450.7482233429    2
2474.6128517914    2

Since the number of points at the extreme end are increasing (from 1 to 2), We therefore see that a value 2500 is still above the saturation level (the number of pixels has started to increase)! A more reasonable saturation level for this image would be 2200! As an exercise, you can try automating this selection with AWK.

Therefore, we can set the saturation level in this image⁵⁸ to be 2200. Let’s mask all such pixels with the command below:

$ astarithmetic saturated.fits set-i i i 2200 gt nan where \
                --output=sat-masked.fits
$ astscript-fits-view sat-masked.fits --ds9scale=minmax

Please see the peaks of several bright stars, not just the central very bright star. Zoom into each of the peaks you see. Besides the central very bright one that we were looking at closely until now, only one other star is saturated (its center is NaN, or Not-a-Number). Try to find it.

But we are not done yet! Please zoom-in to that central bright star and have another look on the edges of the vertical “bleeding” saturated pixels, there are strong positive/negative values touching it (almost like “waves”). These will also cause problems and have to be masked! So with a small addition to the previous command, let’s dilate the saturated regions (with 2-connectivity, or 8-connected neighbors) four times and have another look:

$ astarithmetic saturated.fits set-i i i 2200 gt \
                2 dilate 2 dilate 2 dilate 2 dilate \
                nan where --output=sat-masked.fits
$ astscript-fits-view sat-masked.fits --ds9scale=minmax

Now that saturated pixels (and their problematic neighbors) have been masked, we can convolve the image (recall that Segment will use the convolved image for identifying clumps) with the command below. However, we will use the Spatial Domain convolution which can account for blank pixels (for more on the pros and cons of spatial and frequency domain convolution, see Spatial vs. Frequency domain). We will also create a Gaussian kernel with \(\rm{FWHM}=2\) pixels, truncated at \(5\times\rm{FWHM}\).

$ astmkprof --kernel=gaussian,2,5 --oversample=1 -okernel.fits
$ astconvolve sat-masked.fits --kernel=kernel.fits --domain=spatial \
              --output=sat-masked-conv.fits
$ astscript-fits-view sat-masked-conv.fits --ds9scale=minmax

Please zoom-in to the star and look closely to see how after spatial-domain convolution, the problematic pixels are still NaN. But Segment requires the profile to start with a maximum value and decrease. So before feeding into Segment, let’s fill the blank values with the maximum value of the neighboring pixels in both the input and convolved images (see Interpolation operators):

$ astarithmetic sat-masked.fits 2 interpolate-maxofregion \
                --output=sat-fill.fits
$ astarithmetic sat-masked-conv.fits 2 interpolate-maxofregion \
                --output=sat-fill-conv.fits
$ astscript-fits-view sat-fill* --ds9scale=minmax

Have a closer look at the opened images. Please zoom-in (you will notice that they are already matched and locked, so they will both zoom-in together). Go to the centers of the saturated stars and confirm how they are filled with the largest non-blank pixel. We can now feed this image to NoiseChisel and Segment as the convolved image:

$ astnoisechisel sat-fill.fits --convolved=sat-fill-conv.fits \
                 --output=sat-nc.fits
$ astsegment sat-nc.fits --convolved=sat-fill-conv.fits \
             --output=sat-seg.fits --rawoutput
$ ds9 -mecube sat-seg.fits -zoom to fit -scale limits -1 1

See the CLUMPS extension. Do you see how the whole center of the star has indeed been identified as a single clump? We thus achieved our aim and did not let the saturated pixels harm the identification of the center!

If the issue was only clumps (like in a normal deep image processing), this was the end of Segment’s special considerations. However, in the scenario here, with the very extended wings of the bright stars, it usually happens that background objects become “clumps” in the outskirts and will rip the bright star outskirts into separate “objects”. In the next section (One object for the whole detection), we will describe how you can modify Segment to avoid this issue.