How is the RMS map calculated and why is it an error?

[mapfiable]

As far as I understand, RMS, or root-mean-square, is just a kind of average, like mean or median e.g. In the wiki article it says it "defined as the square root of the mean square (the arithmetic mean of the squares of a set of numbers)". So how can the RMS serve as an error? Is it possible that instead of the RMS, the RMSE, i.e. the root-mean-square-error is calculated?

In any case, I guess it is calculated from the background, but I couldn't find a reference to how this is done anywhere (neither in the SEP or SE documentation). I even tried looking into the source code of SE, but I didn't manage to figure it out (I don't speak C). I guess the the error in a pixel is somehow related to how far it is off from the avarage local background, but I really don't know. Also the RMS value is different for each pixel, so I guess some kind of smoothing is also applied (supposedly the same natural bicubic spline interpolation that is also applied to the background).

[kboone]

You seem to be on the right track here. To estimate the error, sep first calculates the background level for an image, and then subtracts that off. Once this is done, you should only be left with noise from the uncertainties in each pixel. sep then calculates the RMS of the residual values for each pixel after the background was removed. This RMS is an estimate of the standard deviation since the mean of the image has been removed (see e.g. https://en.wikipedia.org/wiki/Root_mean_square#Relationship_to_other_statistics).

sep estimates the background and RMS in patches across the image. By default, the patches are 64x64 pixels. It then interpolates the background and RMS values between these patches.

[mapfiable]

• Once this is done, you should only be left with noise from the uncertainties in each pixel.

  I guess this is assuming there is no signal in the image? So, if the RMS is high in any pixel, this means there is probably signal there and that's why it can be used for the thresholding during source extraction?

• sep then calculates the RMS of the residual values for each pixel after the background was removed.

  and

  sep estimates the background and RMS in patches across the image. ... It then interpolates the background and RMS values between these patches.

  So in other words, the RMS is really not computed pixelwise, but for a whole patch, while the difference between individual pixel values within a patch is a result of the interpolation?

[kboone]

I guess this is assuming there is no signal in the image? So, if the RMS is high in any pixel, this means there is probably signal there and that's why it can be used for the thresholding during source extraction?

That's right. In many cases astronomical images are sparse, meaning that most of the pixels are background. sep (and source extractor) use a robust estimate of the background level that clips "outliers". Any pixels with signal will look like "outliers", so as long as only a small fraction of pixels have signal this won't be a problem.

If your image isn't very sparse (e.g. crowded fields, larger galaxies), then you can do a two-pass background estimation in sep to clean things up. To do this, use a first pass of sep to identify objects in the image. Then use sep.mask_ellipse to mask out the detected objects. Finally, run sep.Background again with this mask. Here is a code snippet to do this:

# Initial estimate of the sky.
initial_sky = sep.Background(data)

# Do a sky subtraction, detect objects and mask them
initial_skysub_data = data.copy()
initial_sky.subfrom(data)

objects = sep.extract(initial_skysub_data, 3.0*initial_sky.globalrms)

object_mask = np.zeros_like(data)
sep.mask_ellipse(object_mask, objects['x'], objects['y'],
                 cxx=objects['cxx'], cyy=objects['cyy'],
                 cxy=objects['cxy'], r=5.)

# Estimate the sky again with our mask
sky = sep.Background(data, mask=object_mask)

So in other words, the RMS is really not computed pixelwise, but for a whole patch, while the difference between individual pixel values within a patch is a result of the interpolation?

Yes, both the background and error levels are evaluated on a sparse grid (every 64x64 pixels by default) and then interpolated across the image. You can see an example of that here.

[mapfiable]

If your image isn't very sparse (e.g. crowded fields, larger galaxies), then you can do a two-pass background estimation in sep to clean things up.

This is an interesting suggestion, thanks! So the idea behind this is that when you run your background estimation for the second time, you won't overestimate the background signal as much, and thus have more "real" signal left after the backgroud removal, so you will find more of the fainter sources in the image?

[kboone]

The background estimation uses a robust estimator, so the default method will be fine (even at low signal-to-noise) if your image is sparse. You only need to use the masking if there are large objects in your image that will make it so that there are too many "outliers" for the robust estimator.