GNU Astronomy Utilities

2.2.5 Achieved surface brightness level

In NoiseChisel optimization we customized NoiseChisel for a single-exposure SDSS image of the M51 group and in Image surface brightness limit we measured the surface brightness limit and the upper-limit surface brightness level (which are both measures of the noise level). In this section, let’s do some measurements on the outer-most edges of the M51 group to see how they relate to the noise measurements found in the previous section.

For this measurement, we will need to estimate the average flux on the outer edges of the detection. Fortunately all this can be done with a few simple commands using Arithmetic and MakeCatalog. First, let’s separate each detected region, or give a unique label/counter to all the connected pixels of NoiseChisel’s detection map with the command below. Recall that with the set- operator, the popped operand will be given a name (det in this case) for easy usage later.

$ astarithmetic r_detected.fits -hDETECTIONS set-det \
                det 2 connected-components -olabeled.fits

You can find the label of the main galaxy visually (by opening the image and hovering your mouse over the M51 group’s label). But to have a little more fun, let’s do this automatically (which is necessary in a general scenario). The M51 group detection is by far the largest detection in this image, this allows us to find its ID/label easily. We will first run MakeCatalog to find the area of all the labels, then we will use Table to find the ID of the largest object and keep it as a shell variable (id):

# Run MakeCatalog to find the area of each label.
$ astmkcatalog labeled.fits --ids --geo-area -h1 -ocat.fits

## Sort the table by the area column.
$ asttable cat.fits --sort=AREA_FULL

## The largest object, is the last one, so we will use '--tail'.
$ asttable cat.fits --sort=AREA_FULL --tail=1

## We only want the ID, so let's only ask for that column:
$ asttable cat.fits --sort=AREA_FULL --tail=1 --column=OBJ_ID

## Now, let's put this result in a variable (instead of printing)
$ id=$(asttable cat.fits --sort=AREA_FULL --tail=1 --column=OBJ_ID)

## Just to confirm everything is fine.
$ echo $id

We can now use the id variable to reject all other detections:

$ astarithmetic labeled.fits $id eq -oonly-m51.fits

Open the image and have a look. To separate the outer edges of the detections, we will need to “erode” the M51 group detection. So in the same Arithmetic command as above, we will erode three times (to have more pixels and thus less scatter), using a maximum connectivity of 2 (8-connected neighbors). We will then save the output in eroded.fits.

$ astarithmetic labeled.fits $id eq 2 erode 2 erode 2 erode \

In labeled.fits, we can now set all the 1-valued pixels of eroded.fits to 0 using Arithmetic’s where operator added to the previous command. We will need the pixels of the M51 group in labeled.fits two times: once to do the erosion, another time to find the outer pixel layer. To do this (and be efficient and more readable) we will use the set-i operator (to give this image the name ‘i’). In this way we can use it any number of times afterwards, while only reading it from disk and finding M51’s pixels once.

$ astarithmetic labeled.fits $id eq set-i i \
                i 2 erode 2 erode 2 erode 0 where -oedge.fits

Open the image and have a look. You’ll see that the detected edge of the M51 group is now clearly visible. You can use edge.fits to mark (set to blank) this boundary on the input image and get a visual feeling of how far it extends:

$ astarithmetic r.fits -h0 edge.fits nan where -oedge-masked.fits

To quantify how deep we have detected the low-surface brightness regions (in units of signal to-noise ratio), we will use the command below. In short it just divides all the non-zero pixels of edge.fits in the Sky subtracted input (first extension of NoiseChisel’s output) by the pixel standard deviation of the same pixel. This will give us a signal-to-noise ratio image. The mean value of this image shows the level of surface brightness that we have achieved. You can also break the command below into multiple calls to Arithmetic and create temporary files to understand it better. However, if you have a look at Reverse polish notation and Arithmetic operators, you should be able to easily understand what your computer does when you run this command54.

$ astarithmetic edge.fits -h1                  set-edge \
                r_detected.fits -hSKY_STD      set-skystd \
                r_detected.fits -hINPUT-NO-SKY set-skysub \
                skysub skystd / edge not nan where meanvalue --quiet

We have thus detected the wings of the M51 group down to roughly 1/3rd of the noise level in this image which is a very good achievement! But the per-pixel S/N is a relative measurement. Let’s also measure the depth of our detection in absolute surface brightness units; or magnitudes per square arc-seconds (see Brightness, Flux, Magnitude and Surface brightness). We will also ask for the S/N and magnitude of the full edge we have defined. Fortunately doing this is very easy with Gnuastro’s MakeCatalog:

$ astmkcatalog edge.fits -h1 --valuesfile=r_detected.fits \
               --zeropoint=22.5 --ids --sb --sn --magnitude
$ asttable edge_cat.fits
1      25.6971       55.2406       15.8994

We have thus reached an outer surface brightness of \(25.70\) magnitudes/arcsec\(^2\) (second column in edge_cat.fits) on this single exposure SDSS image! This is very similar to the surface brightness limit measured in Image surface brightness limit (which is a big achievement!). But another point in the result above is very interesting: the total S/N of the edge is \(55.24\) with a total edge magnitude55 of 15.90!!! This is very large for such a faint signal (recall that the mean S/N per pixel was 0.32) and shows a very important point in the study of galaxies: While the per-pixel signal in their outer edges may be very faint (and invisible to the eye in noise), a lot of signal hides deeply buried in the noise.

In interpreting this value, you should just have in mind that NoiseChisel works based on the contiguity of signal in the pixels. Therefore the larger the object, the deeper NoiseChisel can carve it out of the noise (for the same outer surface brightness). In other words, this reported depth, is the depth we have reached for this object in this dataset, processed with this particular NoiseChisel configuration. If the M51 group in this image was larger/smaller than this (the field of view was smaller/larger), or if the image was from a different instrument, or if we had used a different configuration, we would go deeper/shallower.



edge.fits (extension 1) is a binary (0 or 1 valued) image. Applying the not operator on it, just flips all its pixels (from 0 to 1 and vice-versa). Using the where operator, we are then setting all the newly 1-valued pixels (pixels that are not on the edge) to NaN/blank in the sky-subtracted input image (r_detected.fits, extension INPUT-NO-SKY, which we call skysub). We are then dividing all the non-blank pixels (only those on the edge) by the sky standard deviation (r_detected.fits, extension SKY_STD, which we called skystd). This gives the signal-to-noise ratio (S/N) for each of the pixels on the boundary. Finally, with the meanvalue operator, we are taking the mean value of all the non-blank pixels and reporting that as a single number.


You can run MakeCatalog on only-m51.fits instead of edge.fits to see the full magnitude of the M51 group in this image.