Re: Upper limits for surface brightness levels with noisechisel, mkcatalog

[Mohammad Akhlaghi]

Hi Leslie,

Thank you very much for the very complete details; I think I understand 
the source of the confusion. Let me explain:

If I understood the final verification (4) properly, you measured the 
standard deviation of the "pixel values", or fit with a Gaussian with 
the "pixel values" (within the 1000 apertures). If so, then this is also 
how the surface brightness limit is measured: the standard deviation of 
"pixel values" (of undetected regions). This explains why your two 
verification tests are very similar to the surface brightness limit: 
they are essentially the same as the surface brightness limit.

Ultimately, the surface brightness limit is an extrapolation: 
extrapolating the standard deviation of "pixel values" (over one pixel's 
area) to another arbitrary and shape-agnostic area. Therefore it 
implicitly assumes that there is no correlation between pixels.

But in the upper limit measurement, we measure the standard deviation of 
the distribution of the "sum of pixels within each footprint" (we first 
measure the sum of pixels within each aperture, then measure the 
standard deviation of the sums). Therefore, in the upper limit surface 
brightness, we take into account spatial correlation between pixel 
values within the footprint (which imprints itself in the distribution 
of sums).

To verify this (that the surface brightness limit and upper limit 
surface brightness will be statistically the same with there is no 
correlation between pixel values) I wrote the following script (since 
this exact same discussion has come up in many contexts for me recently!):

https://codeberg.org/Gnuastro/scripts/src/branch/master/sblim-uplimsb-on-mock.sh

You can see the output of one run of this script in the P.S. Everything 
in the script is customizable through the top set of variables, so feel 
free to change them to your particular situation, and run it: it needs 
no inputs ;-).

As you see in the P.S., the upper limit surface brightness and surface 
brightness limit are statistically identical (I used ten thousand random 
apertures for better precision).

As "real world" artifacts are inserted, the nice almost 
symmetric/Gaussian distribution becomes skewed. These real-world 
artifacts can be the following (off the top of my head; other things 
will also be present):

* Varying depth (especially in wide area surveys). For example where we 
have thin regions that are shallower than deeper parts. For example 
compare the first and second panels of the image below from the paper 
describing Gnuastro's 'astscript-pointing-simulate':

https://www.astroexplorer.org/details/rnaasacfc44f1

When the random apertures fall on the shallower regions (with larger 
noise), they will add scatter in the upper-limit surface brightness. The 
extra scatter from this should be symmetric (not causing skewness). 
Depending on the pointing pattern, this type of correlation can be 
more/less significant. One way to remove this type of correlation is to 
use the exposure map of the image and mask (set to NaN) all shallower 
regions in the image that you feed to MakeCatalog for upper-limit 
surface brightness.

* Undetected (or unmasked) pixels: you can mask all your detected pixels 
and visually check the edges of the (larger) detections to see if any 
residual undetected pixels are present or not. If so, you should 
optimize NoiseChisel (and see why they are not detected). This is only 
relevant for very large objects that cover large areas.

* Bad sky subtraction (which can cause over-subtraction in parts of the 
image): you should check the 'SKY' output of NoiseChisel for this.

* Correlated noise due to the warping algorithm (when the single 
exposure images were stacked).

Maybe one nice first that you can do is to plot the distribution of the 
'RANDOM_SUM' column of your '*_upcheck.fits' file (similar to the one in 
the P.S.; see the script). Also, to improve the precision (and see 
problems more clearly in the histogram), you can increase '--upnum' to 
10000 or even 100000. Afterwards, try masking the shallower regions from 
the exposure map and re-run the test to see the effect of the first 
issue above.

For further understanding this large difference, you can again use the 
actual measurements in '*_upcheck.fits' file: find (with 'RANDOM_SUM' 
column values) the apertures that seem to be causing problems in the 
histogram, and using the 'RANDOM_X' and 'RANDOM_Y' values, actually 
check the image: you should be able to see what is causing them.

I hope this helps to better understand those results and improve the 
measurements.

Please keep us updated on the causes as you discover them ;-).

Cheers,
Mohammad





P.S. Output of one run of the script above.

$ ./sblim-uplimsb-on-mock.sh ~/sb-test

Statistics (GNU Astronomy Utilities) 0.22.7-e1f3
-------
Input: ~/sb-test/catalog_upcheck.fits (hdu: 1)
Column: RANDOM_SUM
Unit: input-units
-------
  Number of elements:                      10000
  Minimum:                                 -2.868432e-01
  Maximum:                                 2.889705e-01
  Mode:                                    0.02748960629
  Mode quantile:                           0.5089508951
  Median:                                  2.525770e-02
  Mean:                                    0.02291053996
  Standard deviation:                      0.08211613871
-------
Histogram:
 |                                 **    **
 |                                 ***********
 |                              ****************
 |                             *******************
 |                          * *********************
 |                         **************************
 |                      *******************************
 |                     *********************************
 |                  ***************************************
 |           * ************************************************
 |******************************************************************** *
 |----------------------------------------------------------------------



Results:
Input Surface brightness limit:              29
Measured Surface brightness limit:           29.018350
Measured Upper limit surface brightness:     29.034
Skewness (quantile of mean) of upper-limit:  0.490049