This is gnuastro.info, produced by makeinfo version 7.2 from
gnuastro.texi.

This book documents version 0.23.119-586f of the GNU Astronomy Utilities
(Gnuastro).  Gnuastro provides various programs and libraries for
astronomical data manipulation and analysis.

   Copyright © 2015-2025 Free Software Foundation, Inc.

     Permission is granted to copy, distribute and/or modify this
     document under the terms of the GNU Free Documentation License,
     Version 1.3 or any later version published by the Free Software
     Foundation; with no Invariant Sections, no Front-Cover Texts, and
     no Back-Cover Texts.  A copy of the license is included in the
     section entitled "GNU Free Documentation License".
INFO-DIR-SECTION Astronomy
START-INFO-DIR-ENTRY
* Gnuastro: (gnuastro).       GNU Astronomy Utilities.
* libgnuastro: (gnuastro)Gnuastro library. Full Gnuastro library doc.

* help-gnuastro: (gnuastro)help-gnuastro mailing list. Getting help.

* bug-gnuastro: (gnuastro)Report a bug. How to report bugs

* Arithmetic: (gnuastro)Arithmetic. Arithmetic operations on pixels.
* astarithmetic: (gnuastro)Invoking astarithmetic. Options to Arithmetic.

* BuildProgram: (gnuastro)BuildProgram. Compile and run programs using Gnuastro's library.
* astbuildprog: (gnuastro)Invoking astbuildprog. Options to BuildProgram.

* ConvertType: (gnuastro)ConvertType. Convert different file types.
* astconvertt: (gnuastro)Invoking astconvertt. Options to ConvertType.

* Convolve: (gnuastro)Convolve. Convolve an input file with kernel.
* astconvolve: (gnuastro)Invoking astconvolve. Options to Convolve.

* CosmicCalculator: (gnuastro)CosmicCalculator. For cosmological params.
* astcosmiccal: (gnuastro)Invoking astcosmiccal. Options to CosmicCalculator.

* Crop: (gnuastro)Crop. Crop region(s) from image(s).
* astcrop: (gnuastro)Invoking astcrop. Options to Crop.

* Fits: (gnuastro)Fits. View and manipulate FITS extensions and keywords.
* astfits: (gnuastro)Invoking astfits. Options to Fits.

* MakeCatalog: (gnuastro)MakeCatalog. Make a catalog from labeled image.
* astmkcatalog: (gnuastro)Invoking astmkcatalog. Options to MakeCatalog.

* MakeProfiles: (gnuastro)MakeProfiles. Make mock profiles.
* astmkprof: (gnuastro)Invoking astmkprof. Options to MakeProfiles.

* Match: (gnuastro)Match. Match two separate catalogs.
* astmatch: (gnuastro)Invoking astmatch. Options to Match.

* NoiseChisel: (gnuastro)NoiseChisel. Detect signal in noise.
* astnoisechisel: (gnuastro)Invoking astnoisechisel. Options to NoiseChisel.

* Segment: (gnuastro)Segment. Segment detections based on signal structure.
* astsegment: (gnuastro)Invoking astsegment. Options to Segment.

* Query: (gnuastro)Query. Access remote databases for downloading data.
* astquery: (gnuastro)Invoking astquery. Options to Query.

* Statistics: (gnuastro)Statistics. Get image Statistics.
* aststatistics: (gnuastro)Invoking aststatistics. Options to Statistics.

* Table: (gnuastro)Table. Read and write FITS binary or ASCII tables.
* asttable: (gnuastro)Invoking asttable. Options to Table.

* Warp: (gnuastro)Warp. Warp a dataset to a new grid.
* astwarp: (gnuastro)Invoking astwarp. Options to Warp.

* astscript: (gnuastro)Installed scripts. Gnuastro's installed scripts.
* astscript-ds9-region: (gnuastro)Invoking astscript-ds9-region. Options to this script
* astscript-fits-view: (gnuastro)Invoking astscript-fits-view. Options to this script
* astscript-pointing-simulate: (gnuastro)Invoking astscript-pointing-simulate. Options to this script
* astscript-psf-scale-factor: (gnuastro)Invoking astscript-psf-scale-factor. Options to this script
* astscript-psf-select-stars: (gnuastro)Invoking astscript-psf-select-stars. Options to this script
* astscript-psf-stamp: (gnuastro)Invoking astscript-psf-stamp. Options to this script
* astscript-psf-subtract: (gnuastro)Invoking astscript-psf-subtract. Options to this script
* astscript-psf-unite: (gnuastro)Invoking astscript-psf-unite. Options to this script
* astscript-radial-profile: (gnuastro)Invoking astscript-radial-profile. Options to this script
* astscript-sort-by-night: (gnuastro)Invoking astscript-sort-by-night. Options to this script
* astscript-zeropoint: (gnuastro)Invoking astscript-zeropoint. Options to this script
END-INFO-DIR-ENTRY


File: gnuastro.info,  Node: NoiseChisel,  Next: Segment,  Prev: Statistics,  Up: Data analysis

7.2 NoiseChisel
===============

Once instrumental signatures are removed from the raw data (image) in
the initial reduction process (see *note Data manipulation::).  You are
naturally eager to start answering the scientific questions that
motivated the data collection in the first place.  However, the raw
dataset/image is just an array of values/pixels, that is all!  These raw
values cannot directly be used to answer your scientific questions; for
example, "how many galaxies are there in the image?"  and "What is their
magnitude?".

   The first high-level step your analysis will therefore be to
classify, or label, the dataset elements (pixels) into two classes: 1)
Noise, where random effects are the major contributor to the value, and
2) Signal, where non-random factors (for example, light from a distant
galaxy) are present.  This classification of the elements in a dataset
is formally known as _detection_.

   In an observational/experimental dataset, signal is always buried in
noise: only mock/simulated datasets are free of noise.  Therefore
detection, or the process of separating signal from noise, determines
the number of objects you study and the accuracy of any higher-level
measurement you do on them.  Detection is thus the most important step
of any analysis and is not trivial.  In particular, the most
scientifically interesting astronomical targets are faint, can have a
large variety of morphologies, along with a large distribution in
magnitude and size.  Therefore when noise is significant, proper
detection of your targets is a uniquely decisive step in your final
scientific analysis/result.

   NoiseChisel is Gnuastro's program for detection of targets that do
not have a sharp border (almost all astronomical objects).  When the
targets have sharp edges/borders (for example, cells in biological
imaging), a simple threshold is enough to separate them from noise and
each other (if they are not touching).  To detect such sharp-edged
targets, you can use Gnuastro's Arithmetic program in a command like
below (assuming the threshold is ‘100’, see *note Arithmetic::):

     $ astarithmetic in.fits 100 gt 2 connected-components

   Since almost no astronomical target has such sharp edges, we need a
more advanced detection methodology.  NoiseChisel uses a new noise-based
paradigm for detection of very extended and diffuse targets that are
drowned deeply in the ocean of noise.  It was initially introduced in
Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664) and
improvements after the first four were published in Akhlaghi 2019
(https://arxiv.org/abs/1909.11230).  Please take the time to go through
these papers to most effectively understand the need of NoiseChisel and
how best to use it.

   The name of NoiseChisel is derived from the first thing it does after
thresholding the dataset: to erode it.  In mathematical morphology,
erosion on pixels can be pictured as carving-off boundary pixels.
Hence, what NoiseChisel does is similar to what a wood chisel or stone
chisel do.  It is just not a hardware, but a software.  In fact, looking
at it as a chisel and your dataset as a solid cube of rock will greatly
help in effectively understanding and optimally using it: with
NoiseChisel you literally carve your targets out of the noise.  Try
running it with the ‘--checkdetection’ option, and open the temporary
output as a multi-extension cube, to see each step of the carving
process on your input dataset (see *note Viewing FITS file contents with
DS9 or TOPCAT::).

   NoiseChisel's primary output is a binary detection map with the same
size as the input but its pixels only have two values: 0 (background)
and 1 (foreground).  Pixels that do not harbor any detected signal
(noise) are given a label (or value) of zero and those with a value of 1
have been identified as hosting signal.

   Segmentation is the process of classifying the signal into
higher-level constructs.  For example, if you have two separate galaxies
in one image, NoiseChisel will give a value of 1 to the pixels of both
(each forming an "island" of touching foreground pixels).  After
segmentation, the connected foreground pixels will get separate labels,
enabling you to study them individually.  NoiseChisel is only focused on
detection (separating signal from noise), to _segment_ the signal (into
separate galaxies for example), Gnuastro has a separate specialized
program *note Segment::.  NoiseChisel's output can be directly/readily
fed into Segment.

   For more on NoiseChisel's output format and its benefits (especially
in conjunction with *note Segment:: and later *note MakeCatalog::),
please see Akhlaghi 2016 (https://arxiv.org/abs/1611.06387).  Just note
that when that paper was published, Segment was not yet spun-off into a
separate program, and NoiseChisel done both detection and segmentation.

   NoiseChisel's output is designed to be generic enough to be easily
used in any higher-level analysis.  If your targets are not touching
after running NoiseChisel and you are not interested in their
sub-structure, you do not need the Segment program at all.  You can ask
NoiseChisel to find the connected pixels in the output with the
‘--label’ option.  In this case, the output will not be a binary image
any more, the signal will have counters/labels starting from 1 for each
connected group of pixels.  You can then directly feed NoiseChisel's
output into MakeCatalog for measurements over the detections and the
production of a catalog (see *note MakeCatalog::).

   Thanks to the published papers mentioned above, there is no need to
provide a more complete introduction to NoiseChisel in this book.
However, published papers cannot be updated any more, but the software
has evolved/changed.  The changes since publication are documented in
*note NoiseChisel changes after publication::.  In *note Invoking
astnoisechisel::, the details of running NoiseChisel and its options are
discussed.

   As discussed above, detection is one of the most important steps for
your scientific result.  It is therefore very important to obtain a good
understanding of NoiseChisel (and afterwards *note Segment:: and *note
MakeCatalog::).  We strongly recommend reviewing two tutorials of *note
General program usage tutorial:: and *note Detecting large extended
targets::.  They are designed to show how to most effectively use
NoiseChisel for the detection of small faint objects and large extended
objects.  In the meantime, they also show the modular principle behind
Gnuastro's programs and how they are built to complement, and build
upon, each other.

   *note General program usage tutorial:: culminates in using
NoiseChisel to detect galaxies and use its outputs to find the galaxy
colors.  Defining colors is a very common process in most science-cases.
Therefore it is also recommended to (patiently) complete that tutorial
for optimal usage of NoiseChisel in conjunction with all the other
Gnuastro programs.  *note Detecting large extended targets:: shows you
can optimize NoiseChisel's settings for very extended objects to
successfully carve out to signal-to-noise ratio levels of below 1/10.
After going through those tutorials, play a little with the settings (in
the order presented in the paper and *note Invoking astnoisechisel::) on
a dataset you are familiar with and inspect all the check images
(options starting with ‘--check’) to see the effect of each parameter.

   Below, in *note Invoking astnoisechisel::, we will review
NoiseChisel's input, detection, and output options in *note NoiseChisel
input::, *note Detection options::, and *note NoiseChisel output::.  If
you have used NoiseChisel within your research, please run it with
‘--cite’ to list the papers you should cite and how to acknowledge its
funding sources.

* Menu:

* NoiseChisel changes after publication::  Updates since published papers.
* Invoking astnoisechisel::     Options and arguments for NoiseChisel.


File: gnuastro.info,  Node: NoiseChisel changes after publication,  Next: Invoking astnoisechisel,  Prev: NoiseChisel,  Up: NoiseChisel

7.2.1 NoiseChisel changes after publication
-------------------------------------------

NoiseChisel was initially introduced in Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664) and updates after the first four
years were published in Akhlaghi 2019
(https://arxiv.org/abs/1909.11230).  To help in understanding how it
works, those papers have many figures showing every step on multiple
mock and real examples.  We recommended to read these papers for a good
understanding of what it does and how each parameter influences the
output.

   However, the papers cannot be updated anymore, but NoiseChisel has
evolved (and will continue to do so): better algorithms or steps have
been found and implemented and some options have been added, removed or
changed behavior.  This book is thus the final and definitive guide to
NoiseChisel.  The aim of this section is to make the transition from the
papers above to the installed version on your system, as smooth as
possible with the list below.  For a more detailed list of changes in
each Gnuastro version, please see the ‘NEWS’ file(1).

   • An improved outlier rejection for identifying tiles without any
     signal has been implemented in the quantile-threshold phase: Prior
     to version 0.14, outliers were defined globally: the distribution
     of all tiles with an acceptable ‘--meanmedqdiff’ was inspected and
     outliers were found and rejected.  However, this caused problems
     when there are strong gradients over the image (for example, an
     image prior to flat-fielding, or in the presence of a large
     foreground galaxy).  In these cases, the faint wings of
     galaxies/stars could be mistakenly identified as Sky (leaving a
     footprint of the object on the Sky output) and wrongly subtracted.

     It was possible to play with the parameters to correct this for
     that particular dataset, but that was frustrating.  Therefore from
     version 0.14, instead of finding outliers from the full tile
     distribution, we now measure the _slope_ of the tile's nearby tiles
     and find outliers locally.  Three options have been added to
     configure this part of NoiseChisel: ‘--outliernumngb’,
     ‘--outliersclip’ and ‘--outliersigma’.  For more on the local
     outlier-by-distance algorithm and the definition of _slope_
     mentioned above, see *note Quantifying signal in a tile::.  In our
     tests, this gave a much improved estimate of the quantile
     thresholds and final Sky values with default values.

   ---------- Footnotes ----------

   (1) The ‘NEWS’ file is present in the released Gnuastro tarball, see
*note Release tarball::.


File: gnuastro.info,  Node: Invoking astnoisechisel,  Prev: NoiseChisel changes after publication,  Up: NoiseChisel

7.2.2 Invoking NoiseChisel
--------------------------

NoiseChisel will detect signal in noise producing a multi-extension
dataset containing a binary detection map which is the same size as the
input.  Its output can be readily used for input into *note Segment::,
for higher-level segmentation, or *note MakeCatalog:: to do measurements
and generate a catalog.  The executable name is ‘astnoisechisel’ with
the following general template

     $ astnoisechisel [OPTION ...] InputImage.fits

One line examples:

     ## Detect signal in input.fits.
     $ astnoisechisel input.fits

     ## Inspect all the detection steps after changing a parameter.
     $ astnoisechisel input.fits --qthresh=0.4 --checkdetection

     ## Detect signal assuming input has 4 amplifier channels along first
     ## dimension and 1 along the second. Also set the regular tile size
     ## to 100 along both dimensions:
     $ astnoisechisel --numchannels=4,1 --tilesize=100,100 input.fits

If NoiseChisel is to do processing (for example, you do not want to get
help, or see the values to each input parameter), an input image should
be provided with the recognized extensions (see *note Arguments::).
NoiseChisel shares a large set of common operations with other Gnuastro
programs, mainly regarding input/output, general processing steps, and
general operating modes.  To help in a unified experience between all of
Gnuastro's programs, these operations have the same command-line
options, see *note Common options:: for a full list/description (they
are not repeated here).

   As in all Gnuastro programs, options can also be given to NoiseChisel
in configuration files.  For a thorough description on Gnuastro's
configuration file parsing, please see *note Configuration files::.  All
of NoiseChisel's options with a short description are also always
available on the command-line with the ‘--help’ option, see *note
Getting help::.  To inspect the option values without actually running
NoiseChisel, append your command with ‘--printparams’ (or ‘-P’).

   NoiseChisel's input image may contain blank elements (see *note Blank
pixels::).  Blank elements will be ignored in all steps of NoiseChisel.
Hence if your dataset has bad pixels which should be masked with a mask
image, please use Gnuastro's *note Arithmetic:: program (in particular
its ‘where’ operator) to convert those pixels to blank pixels before
running NoiseChisel.  Gnuastro's Arithmetic program has bitwise
operators helping you select specific kinds of bad-pixels when
necessary.

   A convolution kernel can also be optionally given.  If a value (file
name) is given to ‘--kernel’ on the command-line or in a configuration
file (see *note Configuration files::), then that file will be used to
convolve the image prior to thresholding.  Otherwise a default kernel
will be used.  For a 2D image, the default kernel is a 2D Gaussian with
a FWHM of 2 pixels truncated at 5 times the FWHM. This choice of the
default kernel is discussed in Section 3.1.1 of Akhlaghi and Ichikawa
2015 (https://arxiv.org/abs/1505.01664).  For a 3D cube, it is a
Gaussian with FWHM of 1.5 pixels in the first two dimensions and 0.75
pixels in the third dimension.  See *note Convolution kernel:: for
kernel related options.  Passing ‘none’ to ‘--kernel’ will disable
convolution.  On the other hand, through the ‘--convolved’ option, you
may provide an already convolved image, see descriptions below for more.

   NoiseChisel defines two tessellations over the input (see *note
Tessellation::).  This enables it to deal with possible gradients in the
input dataset and also significantly improve speed by processing each
tile on different threads simultaneously.  Tessellation related options
are discussed in *note Processing options::.  In particular, NoiseChisel
uses two tessellations (with everything between them identical except
the tile sizes): a fine-grained one with smaller tiles (used in
thresholding and Sky value estimations) and another with larger tiles
which is used for pseudo-detections over non-detected regions of the
image.  The common Tessellation options described in *note Processing
options:: define all parameters of both tessellations.  The large tile
size for the latter tessellation is set through the ‘--largetilesize’
option.  To inspect the tessellations on your input dataset, run
NoiseChisel with ‘--checktiles’.

*Usage TIP:* Frequently use the options starting with ‘--check’.  Since
the noise properties differ between different datasets, you can often
play with the parameters/options for a better result than the default
parameters.  You can start with ‘--checkdetection’ for the main steps.
For the full list of NoiseChisel's checking options please run:
     $ astnoisechisel --help | grep check

*Not detecting wings of bright galaxies:* In such cases, probably the
best solution is to increase ‘--outliernumngb’ (to reject tiles that are
affected by very flat diffuse signal).  For more, see *note Quantifying
signal in a tile::.

   When working on 3D data cubes, the tessellation options need three
values and updating them every time can be annoying/buggy.  To simplify
the job, NoiseChisel also installs a ‘astnoisechisel-3d.conf’
configuration file (see *note Configuration files::).  You can use this
for default values on data cubes.  For example, if you installed
Gnuastro with the prefix ‘/usr/local’ (the default location, see *note
Installation directory::), you can benefit from this configuration file
by running NoiseChisel like the example below.

     $ astnoisechisel cube.fits \
                  --config=/usr/local/etc/gnuastro/astnoisechisel-3d.conf

   To further simplify the process, you can define a shell alias in any
startup file (for example, ‘~/.bashrc’, see *note Installation
directory::).  Assuming that you installed Gnuastro in ‘/usr/local’, you
can add this line to the startup file (you may put it all in one line,
it is broken into two lines here for fitting within page limits).

     alias astnoisechisel-3d="astnoisechisel \
                    --config=/usr/local/etc/gnuastro/astnoisechisel-3d.conf"

Using this alias, you can call NoiseChisel with the name
‘astnoisechisel-3d’ (instead of ‘astnoisechisel’).  It will
automatically load the 3D specific configuration file first, and then
parse any other arguments, options or configuration files.  You can
change the default values in this 3D configuration file by calling them
on the command-line as you do with ‘astnoisechisel’(1).  For example:

     $ astnoisechisel-3d --numchannels=3,3,1 cube.fits

   Below, we will discuss NoiseChisel's options, classified into
separate sub-sections to help in easy navigation.  *note NoiseChisel
input:: discusses the basic options relating to input file(s) and data;
these have no effect on the the detection process.  Afterwards, *note
Detection options:: fully describes every configuration parameter
(option) related to detection and how they affect the final result.  The
order of options in this section follow the logical order within
NoiseChisel.  On first reading (while you are still new to NoiseChisel),
it is therefore strongly recommended to read the options in the given
order below.  The output of ‘--printparams’ (or ‘-P’) also has this
order.  However, the output of ‘--help’ is sorted alphabetically.
Finally, in *note NoiseChisel output:: the format of NoiseChisel's
output is discussed.

* Menu:

* NoiseChisel input::           NoiseChisel's input options.
* Detection options::           Configure detection in NoiseChisel.
* NoiseChisel output::          NoiseChisel's output options and format.

   ---------- Footnotes ----------

   (1) Recall that for single-invocation options, the last command-line
invocation takes precedence over all previous invocations (including
those in the 3D configuration file).  See the description of ‘--config’
in *note Operating mode options::.


File: gnuastro.info,  Node: NoiseChisel input,  Next: Detection options,  Prev: Invoking astnoisechisel,  Up: Invoking astnoisechisel

7.2.2.1 NoiseChisel input
.........................

The options here can be used to configure the inputs and output of
NoiseChisel, along with some general processing options.  Recall that
you can always see the full list of Gnuastro's options with the ‘--help’
(see *note Getting help::), or ‘--printparams’ (or ‘-P’) to see their
values (see *note Operating mode options::).

‘-k FITS’
‘--kernel=FITS’
     File name of kernel to smooth the image before applying the
     threshold, see *note Convolution kernel::.  If no convolution is
     needed, give this option a value of ‘none’.

     The first step of NoiseChisel is to convolve/smooth the image and
     use the convolved image in multiple steps including the finding and
     applying of the quantile threshold (see ‘--qthresh’).  The
     ‘--kernel’ option is not mandatory.  If not called, for a 2D, image
     a 2D Gaussian profile with a FWHM of 2 pixels truncated at 5 times
     the FWHM is used.  This choice of the default kernel is discussed
     in Section 3.1.1 of Akhlaghi and Ichikawa [2015].

     For a 3D cube, when no file name is given to ‘--kernel’, a Gaussian
     with FWHM of 1.5 pixels in the first two dimensions and 0.75 pixels
     in the third dimension will be used.  The reason for this
     particular configuration is that commonly in astronomical
     applications, 3D datasets do not have the same nature in all three
     dimensions, commonly the first two dimensions are spatial (RA and
     Dec) while the third is spectral (for example, wavelength).  The
     samplings are also different, in the default case, the spatial
     sampling is assumed to be larger than the spectral sampling, hence
     a wider FWHM in the spatial directions, see *note Sampling
     theorem::.

     You can use MakeProfiles to build a kernel with any of its
     recognized profile types and parameters.  For more details, please
     see *note MakeProfiles output dataset::.  For example, the command
     below will make a Moffat kernel (with $\beta=2.8$) with FWHM of 2
     pixels truncated at 10 times the FWHM.

          $ astmkprof --oversample=1 --kernel=moffat,2,2.8,10

     Since convolution can be the slowest step of NoiseChisel, for large
     datasets, you can convolve the image once with Gnuastro's Convolve
     (see *note Convolve::), and use the ‘--convolved’ option to feed it
     directly to NoiseChisel.  This can help getting faster results when
     you are playing/testing the higher-level options.

‘--khdu=STR’
     HDU containing the kernel in the file given to the ‘--kernel’
     option.

‘--convolved=FITS’
     Use this file as the convolved image and do not do convolution
     (ignore ‘--kernel’).  NoiseChisel will just check the size of the
     given dataset is the same as the input's size.  If a wrong image
     (with the same size) is given to this option, the results (errors,
     bugs, etc.)  are unpredictable.  So please use this option with
     care and in a highly controlled environment, for example, in the
     scenario discussed below.

     In almost all situations, as the input gets larger, the single most
     CPU (and time) consuming step in NoiseChisel (and other programs
     that need a convolved image) is convolution.  Therefore minimizing
     the number of convolutions can save a significant amount of time in
     some scenarios.  One such scenario is when you want to segment
     NoiseChisel's detections using the same kernel (with *note
     Segment::, which also supports this ‘--convolved’ option).  This
     scenario would require two convolutions of the same dataset: once
     by NoiseChisel and once by Segment.  Using this option in both
     programs, only one convolution (prior to running NoiseChisel) is
     enough.

     Another common scenario where this option can be convenient is when
     you are testing NoiseChisel (or Segment) for the best parameters.
     You have to run NoiseChisel multiple times and see the effect of
     each change.  However, once you are happy with the kernel,
     re-convolving the input on every change of higher-level parameters
     will greatly hinder, or discourage, further testing.  With this
     option, you can convolve the input image with your chosen kernel
     once before running NoiseChisel, then feed it to NoiseChisel on
     each test run and thus save valuable time for better/more tests.

     To build your desired convolution kernel, you can use *note
     MakeProfiles::.  To convolve the image with a given kernel you can
     use *note Convolve::.  Spatial domain convolution is mandatory: in
     the frequency domain, blank pixels (if present) will cover the
     whole image and gradients will appear on the edges, see *note
     Spatial vs. Frequency domain::.

     Below you can see an example of the second scenario: you want to
     see how variation of the growth level (through the ‘--detgrowquant’
     option) will affect the final result.  Recall that you can ignore
     all the extra spaces, new lines, and backslash's ('‘\’') if you are
     typing in the terminal.  In a shell script, remove the ‘$’ signs at
     the start of the lines.

          ## Make the kernel to convolve with.
          $ astmkprof --oversample=1 --kernel=gaussian,2,5

          ## Convolve the input with the given kernel.
          $ astconvolve input.fits --kernel=kernel.fits                \
                        --domain=spatial --output=convolved.fits

          ## Run NoiseChisel with seven growth quantile values.
          $ for g in 60 65 70 75 80 85 90; do                          \
              astnoisechisel input.fits --convolved=convolved.fits     \
                             --detgrowquant=0.$g --output=$g.fits;     \
            done

‘--chdu=STR’
     The HDU/extension containing the convolved image in the file given
     to ‘--convolved’.

‘-w FITS’
‘--widekernel=FITS’
     File name of a wider kernel to use in estimating the difference of
     the mode and median in a tile (this difference is used to identify
     the significance of signal in that tile, see *note Quantifying
     signal in a tile::).  As displayed in Figure 4 of Akhlaghi and
     Ichikawa 2015 (https://arxiv.org/abs/1505.01664), a wider kernel
     will help in identifying the skewness caused by data in noise.  The
     image that is convolved with this kernel is _only_ used for this
     purpose.  Once the mode is found to be sufficiently close to the
     median, the quantile threshold is found on the image convolved with
     the sharper kernel (‘--kernel’), see ‘--qthresh’).

     Since convolution will significantly slow down the processing, this
     feature is optional.  When it is not given, the image that is
     convolved with ‘--kernel’ will be used to identify good tiles _and_
     apply the quantile threshold.  This option is mainly useful in
     conditions were you have a very large, extended, diffuse signal
     that is still present in the usable tiles when using ‘--kernel’.
     See *note Detecting large extended targets:: for a practical
     demonstration on how to inspect the tiles used in identifying the
     quantile threshold.

‘--whdu=STR’
     HDU containing the kernel file given to the ‘--widekernel’ option.

‘-L INT[,INT]’
‘--largetilesize=INT[,INT]’
     The size of each tile for the tessellation with the larger tile
     sizes.  Except for the tile size, all the other parameters for this
     tessellation are taken from the common options described in *note
     Processing options::.  The format is identical to that of the
     ‘--tilesize’ option that is discussed in that section.


File: gnuastro.info,  Node: Detection options,  Next: NoiseChisel output,  Prev: NoiseChisel input,  Up: Invoking astnoisechisel

7.2.2.2 Detection options
.........................

Detection is the process of separating the pixels in the image into two
groups: 1) Signal, and 2) Noise.  Through the parameters below, you can
customize the detection process in NoiseChisel.  Recall that you can
always see the full list of NoiseChisel's options with the ‘--help’ (see
*note Getting help::), or ‘--printparams’ (or ‘-P’) to see their values
(see *note Operating mode options::).

‘-Q FLT’
‘--meanmedqdiff=FLT’
     The maximum acceptable distance between the quantiles of the mean
     and median in each tile, see *note Quantifying signal in a tile::.
     The quantile threshold estimates are measured on tiles where the
     quantiles of their mean and median are less distant than the value
     given to this option.  For example, ‘--meanmedqdiff=0.01’ means
     that only tiles where the mean's quantile is between 0.49 and 0.51
     (recall that the median's quantile is 0.5) will be used.

‘-a INT’
‘--outliernumngb=INT’
     Number of neighboring tiles to use for outlier rejection (mostly
     the wings of bright stars or galaxies).  For optimal detection of
     the wings of bright stars or galaxies, this is *the most important*
     option in NoiseChisel.  This is because the extended wings of
     bright galaxies or stars (the PSF) can become flat over the tile.
     In this case, they will satisfy the ‘--meanmedqdiff’ condition and
     pass that step.  Therefore, to correctly identify such bad tiles,
     we need to look at the neighboring nearby tiles.  A tile that is on
     the wing of a bright galaxy/star will clearly be an outlier when
     looking at the neighbors.  For more on the details of the outlier
     rejection algorithm, see the latter half of *note Quantifying
     signal in a tile::.  If this option is given a value of zero, no
     outlier rejection will take place.

‘--outliersclip=FLT,FLT’
     $\sigma$-clipping parameters for the outlier rejection of the
     quantile threshold.  The format of the given values is similar to
     ‘--sigmaclip’ below.  In NoiseChisel, outlier rejection on tiles is
     used when identifying the quantile thresholds (‘--qthresh’,
     ‘--noerodequant’, and ‘detgrowquant’).

     Outlier rejection is useful when the dataset contains a large and
     diffuse (almost flat within each tile) signal.  The flatness of the
     profile will cause it to successfully pass the mean-median quantile
     difference test, so we will need to use the distribution of
     successful tiles for removing these false positives.  For more, see
     the latter half of *note Quantifying signal in a tile::.

‘--outliersigma=FLT’
     Multiple of sigma to define an outlier.  If this option is given a
     value of zero, no outlier rejection will take place.  For more see
     ‘--outliersclip’ and the latter half of *note Quantifying signal in
     a tile::.

‘-t FLT’
‘--qthresh=FLT’
     The quantile threshold to apply to the convolved image.  The
     detection process begins with applying a quantile threshold to each
     of the tiles in the small tessellation.  The quantile is only
     calculated for tiles that do not have any significant signal within
     them, see *note Quantifying signal in a tile::.  Interpolation is
     then used to give a value to the unsuccessful tiles and it is
     finally smoothed.

     The quantile value is a floating point value between 0 and 1.
     Assume that we have sorted the $N$ data elements of a distribution
     (the pixels in each mesh on the convolved image).  The quantile
     ($q$) of this distribution is the value of the element with an
     index of (the nearest integer to) $q\times{N}$ in the sorted data
     set.  After thresholding is complete, we will have a binary (two
     valued) image.  The pixels above the threshold are known as
     foreground pixels (have a value of 1) while those which lie below
     the threshold are known as background (have a value of 0).

‘--smoothwidth=INT’
     Width of flat kernel used to smooth the interpolated quantile
     thresholds, see ‘--qthresh’ for more.

‘--checkqthresh’
     Check the quantile threshold values on the mesh grid.  A
     multi-extension FITS file, suffixed with ‘_qthresh.fits’ will be
     created showing each step of how the final quantile threshold is
     found.  With this option, NoiseChisel will abort as soon as
     quantile estimation has been completed, allowing you to inspect the
     steps leading to the final quantile threshold, this can be disabled
     with ‘--continueaftercheck’.  By default the output will have the
     same pixel size as the input, but with the ‘--oneelempertile’
     option, only one pixel will be used for each tile (see *note
     Processing options::).

     The key things to remember are:
        • The measurements to find the thresholds are done on tiles that
          cover the whole image in a tessellation.  Recall that you can
          set the size of tiles with ‘--tilesize’ and check them with
          ‘--checktiles’.  Therefore except for the first and last
          extensions, the rest only show tiles.
        • NoiseChisel ultimately has three thresholds: the quantile
          threshold (that you set with ‘--qthresh’), the no-erode
          quantile (set with ‘--noerodequant’) and the growth quantile
          (set with ‘--detgrowquant’).  Therefore for each step, we have
          three extensions.

     The output file will have the following extensions.  Below, the
     extensions are put in the same order as you see in the file, with
     their name.

     ‘CONVOLVED’
          This is the input image after convolution with the kernel
          (which is a FWHM=2 Gaussian by default, but you can change
          with ‘--kernel’).  Recall that the thresholds are defined on
          the convolved image.

     ‘QTHRESH_ERODE’
     ‘QTHRESH_NOERODE’
     ‘QTHRESH_EXPAND’
          In these three extensions, the tiles that have a
          quantile-of-mean more/less than 0.5 (quantile of median) $\pm
          d$ are set to NaN ($d$ is the value given to ‘--meanmedqdiff’,
          see *note Quantifying signal in a tile::).  Therefore the
          non-NaN tiles that you see here are the tiles where there is
          no significant skewness (changing signal) within that tile.
          The only differing thing between the three extensions is the
          values of the non-NaN tiles.  These values will be used to
          construct the final threshold map over the whole image.

     ‘VALUE1_NO_OUTLIER’
     ‘VALUE2_NO_OUTLIER’
     ‘VALUE3_NO_OUTLIER’
          All outlier tiles have been masked.  The reason for removing
          outliers is that the quantile-of-mean is only sensitive to
          signal that varies on a scale that is smaller than the tile
          size.  Therefore the extended wings of large galaxies or
          bright stars (which vary on scales much larger than the tile
          size) will pass that test.  As described in *note Quantifying
          signal in a tile:: outlier rejection is customized through
          ‘--outliernumngb’, ‘--outliersclip’ and ‘--outliersigma’.

     ‘THRESH1_INTERP’
     ‘THRESH2_INTERP’
     ‘THRESH3_INTERP’
          Using the successful values that remain after the previous
          step, give values to all (interpolate) the tiles in the image.
          The interpolation is done using the nearest-neighbor method:
          for each tile, the N nearest neighbors are found and the
          median of their values is used to fill it.  You can set the
          value of N through the ‘--interpnumngb’ option.

     ‘THRESH1_SMOOTH’
     ‘THRESH2_SMOOTH’
     ‘THRESH3_SMOOTH’
          Smooth the interpolated image to remove the strong differences
          between touching tiles.  Because we used the median value of
          the N nearest neighbors in the previous step, there can be
          strong discontinuities on the edges of tiles (which can
          directly show in the image after applying the threshold).  The
          scale of the smoothing (number of nearby tiles to smooth with)
          is set with the ‘--smoothwidth’ option.

     ‘QTHRESH-APPLIED’
          The pixels in this image can only have three values:

          ‘0’
               These pixels had a value below the quantile threshold.
          ‘1’
               These pixels had a value above the quantile threshold,
               but below the threshold for no erosion.  Therefore in the
               next step, NoiseChisel will erode (set them to 0) these
               pixels if they are touching a 0-valued pixel.
          ‘2’
               These pixels had a value above the no-erosion threshold.
               So NoiseChisel will not erode these pixels, it will only
               apply Opening to them afterwards.  Recall that this was
               done to avoid loosing sharp point-sources (like stars in
               space-based imaging).

‘--blankasforeground’
     In the erosion and opening steps below, treat blank elements as
     foreground (regions above the threshold).  By default, blank
     elements in the dataset are considered to be background, so if a
     foreground pixel is touching it, it will be eroded.  This option is
     irrelevant if the datasets contains no blank elements.

     When there are many blank elements in the dataset, treating them as
     foreground will systematically erode their regions less, therefore
     systematically creating more false positives.  So use this option
     (when blank values are present) with care.

‘-e INT’
‘--erode=INT’
     The number of erosions to apply to the binary thresholded image.
     Erosion is simply the process of flipping (from 1 to 0) any of the
     foreground pixels that neighbor a background pixel.  In a 2D image,
     there are two kinds of neighbors, 4-connected and 8-connected
     neighbors.  In a 3D dataset, there are three: 6-connected,
     18-connected, and 26-connected.  You can specify which class of
     neighbors should be used for erosion with the ‘--erodengb’ option,
     see below.

     Erosion has the effect of shrinking the foreground pixels.  To put
     it another way, it expands the holes.  This is a founding principle
     in NoiseChisel: it exploits the fact that with very low thresholds,
     the holes in the very low surface brightness regions of an image
     will be smaller than regions that have no signal.  Therefore by
     expanding those holes, we are able to separate the regions
     harboring signal.

‘--erodengb=INT’
     The type of neighborhood (structuring element) used in erosion, see
     ‘--erode’ for an explanation on erosion.  If the input is a 2D
     image, only two integer values are acceptable: 4 or 8.  For a 3D
     input data cube, the acceptable values are: 6, 18 and 26.

     In 2D 4-connectivity, the neighbors of a pixel are defined as the
     four pixels on the top, bottom, right and left of a pixel that
     share an edge with it.  The 8-connected neighbors on the other hand
     include the 4-connected neighbors along with the other 4 pixels
     that share a corner with this pixel.  See Figure 6 (a) and (b) in
     Akhlaghi and Ichikawa (2015) for a demonstration.  A similar
     argument applies to 3D data cubes.

‘--noerodequant’
     Pure erosion is going to carve off sharp and small objects
     completely out of the detected regions.  This option can be used to
     avoid missing such sharp and small objects (which have significant
     pixels, but not over a large area).  All pixels with a value larger
     than the significance level specified by this option will not be
     eroded during the erosion step above.  However, they will undergo
     the erosion and dilation of the opening step below.

     Like the ‘--qthresh’ option, the significance level is determined
     using the quantile (a value between 0 and 1).  Just as a reminder,
     in the normal distribution, $1\sigma$, $1.5\sigma$, and $2\sigma$
     are approximately on the 0.84, 0.93, and 0.98 quantiles.

‘-p INT’
‘--opening=INT’
     Depth of opening to be applied to the eroded binary image.  Opening
     is a composite operation.  When opening a binary image with a depth
     of $n$, $n$ erosions (explained in ‘--erode’) are followed by $n$
     dilations.  Simply put, dilation is the inverse of erosion.  When
     dilating an image any background pixel is flipped (from 0 to 1) to
     become a foreground pixel.  Dilation has the effect of fattening
     the foreground.  Note that in NoiseChisel, the erosion which is
     part of opening is independent of the initial erosion that is done
     on the thresholded image (explained in ‘--erode’).  The structuring
     element for the opening can be specified with the ‘--openingngb’
     option.  Opening has the effect of removing the thin foreground
     connections (mostly noise) between separate foreground 'islands'
     (detections) thereby completely isolating them.  Once opening is
     complete, we have _initial_ detections.

‘--openingngb=INT’
     The structuring element used for opening, see ‘--erodengb’ for more
     information about a structuring element.

‘--skyfracnoblank’
     Ignore blank pixels when estimating the fraction of undetected
     pixels for Sky estimation.  NoiseChisel only measures the Sky over
     the tiles that have a sufficiently large fraction of undetected
     pixels (value given to ‘--minskyfrac’).  By default this fraction
     is found by dividing number of undetected pixels in a tile by the
     tile's area.  But this default behavior ignores the possibility of
     blank pixels.  In situations that blank/masked pixels are scattered
     across the image and if they are large enough, all the tiles can
     fail the ‘--minskyfrac’ test, thus not allowing NoiseChisel to
     proceed.  With this option, such scenarios can be fixed: the
     denominator of the fraction will be the number of non-blank
     elements in the tile, not the total tile area.

‘-B FLT’
‘--minskyfrac=FLT’
     Minimum fraction (value between 0 and 1) of Sky (undetected) areas
     in a tile.  Only tiles with a fraction of undetected pixels (Sky)
     larger than this value will be used to estimate the Sky value.
     NoiseChisel uses this option value twice to estimate the Sky value:
     after initial detections and in the end when false detections have
     been removed.

     Because of the PSF and their intrinsic amorphous properties,
     astronomical objects (except cosmic rays) never have a clear cutoff
     and commonly sink into the noise very slowly.  Even below the very
     low thresholds used by NoiseChisel.  So when a large fraction of
     the area of one mesh is covered by detections, it is very plausible
     that their faint wings are present in the undetected regions (hence
     causing a bias in any measurement).  To get an accurate measurement
     of the above parameters over the tessellation, tiles that harbor
     too many detected regions should be excluded.  The used tiles are
     visible in the respective ‘--check’ option of the given step.

‘--checkdetsky’
     Check the initial approximation of the sky value and its standard
     deviation in a FITS file ending with ‘_detsky.fits’.  With this
     option, NoiseChisel will abort as soon as the sky value used for
     defining pseudo-detections is complete.  This allows you to inspect
     the steps leading to the final quantile threshold, this behavior
     can be disabled with ‘--continueaftercheck’.  By default the output
     will have the same pixel size as the input, but with the
     ‘--oneelempertile’ option, only one pixel will be used for each
     tile (see *note Processing options::).

‘-s FLT,FLT’
‘--sigmaclip=FLT,FLT’
     The $\sigma$-clipping parameters for measuring the initial and
     final Sky values from the undetected pixels, see *note Sigma
     clipping::.

     This option takes two values which are separated by a comma (<,>).
     Each value can either be written as a single number or as a
     fraction of two numbers (for example, ‘3,1/10’).  The first value
     to this option is the multiple of $\sigma$ that will be clipped
     ($\alpha$ in that section).  The second value is the exit criteria.
     If it is less than 1, then it is interpreted as tolerance and if it
     is larger than one it is assumed to be the fixed number of
     iterations.  Hence, in the latter case the value must be an
     integer.

‘-R FLT’
‘--dthresh=FLT’
     The detection threshold: a multiple of the initial Sky standard
     deviation added with the initial Sky approximation (which you can
     inspect with ‘--checkdetsky’).  This flux threshold is applied to
     the initially undetected regions on the unconvolved image.  The
     background pixels that are completely engulfed in a 4-connected
     foreground region are converted to background (holes are filled)
     and one opening (depth of 1) is applied over both the initially
     detected and undetected regions.  The Signal to noise ratio of the
     resulting 'pseudo-detections' are used to identify true vs.  false
     detections.  See Section 3.1.5 and Figure 7 in Akhlaghi and
     Ichikawa (2015) for a very complete explanation.

‘--dopening=INT’
     The number of openings to do after applying ‘--dthresh’.

‘--dopeningngb=INT’
     The connectivity used in the opening of ‘--dopening’.  In a 2D
     image this must be either 4 or 8.  The stronger the connectivity,
     the more smaller regions will be discarded.

‘--holengb=INT’
     The connectivity (defined by the number of neighbors) to fill holes
     after applying ‘--dthresh’ (above) to find pseudo-detections.  For
     example, in a 2D image it must be 4 (the neighbors that are most
     strongly connected) or 8 (all neighbors).  The stronger the
     connectivity, the stronger the hole will be enclosed.  So setting a
     value of 8 in a 2D image means that the walls of the hole are
     4-connected.  If standard (near Sky level) values are given to
     ‘--dthresh’, setting ‘--holengb=4’, might fill the complete dataset
     and thus not create enough pseudo-detections.

‘--pseudoconcomp=INT’
     The connectivity (defined by the number of neighbors) to find
     individual pseudo-detections.  If it is a weaker connectivity (4 in
     a 2D image), then pseudo-detections that are connected on the
     corners will be treated as separate.

‘-m INT’
‘--snminarea=INT’
     The minimum area to calculate the Signal to noise ratio on the
     pseudo-detections of both the initially detected and undetected
     regions.  When the area in a pseudo-detection is too small, the
     Signal to noise ratio measurements will not be accurate and their
     distribution will be heavily skewed to the positive.  So it is best
     to ignore any pseudo-detection that is smaller than this area.  Use
     ‘--detsnhistnbins’ to check if this value is reasonable or not.

‘--checksn’
     Save the S/N values of the pseudo-detections (and possibly grown
     detections if ‘--cleangrowndet’ is called) into separate tables.
     If ‘--tableformat’ is a FITS table, each table will be written into
     a separate extension of one file suffixed with ‘_detsn.fits’.  If
     it is plain text, a separate file will be made for each table
     (ending in ‘_detsn_sky.txt’, ‘_detsn_det.txt’ and
     ‘_detsn_grown.txt’).  For more on ‘--tableformat’ see *note Input
     output options::.

     You can use these to inspect the S/N values and their distribution
     (in combination with the ‘--checkdetection’ option to see where the
     pseudo-detections are).  You can use Gnuastro's *note Statistics::
     to make a histogram of the distribution or any other analysis you
     would like for better understanding of the distribution (for
     example, through a histogram).

‘--minnumfalse=INT’
     The minimum number of 'pseudo-detections' over the undetected
     regions to identify a Signal-to-Noise ratio threshold.  The Signal
     to noise ratio (S/N) of false pseudo-detections in each tile is
     found using the quantile of the S/N distribution of the
     pseudo-detections over the undetected pixels in each mesh.  If the
     number of S/N measurements is not large enough, the quantile will
     not be accurate (can have large scatter).  For example, if you set
     ‘--snquant=0.99’ (or the top 1 percent), then it is best to have at
     least 100 S/N measurements.

‘-c FLT’
‘--snquant=FLT’
     The quantile of the Signal to noise ratio distribution of the
     pseudo-detections in each mesh to use for filling the large mesh
     grid.  Note that this is only calculated for the large mesh grids
     that satisfy the minimum fraction of undetected pixels (value of
     ‘--minbfrac’) and minimum number of pseudo-detections (value of
     ‘--minnumfalse’).

‘--snthresh=FLT’
     Manually set the signal-to-noise ratio of true pseudo-detections.
     With this option, NoiseChisel will not attempt to find
     pseudo-detections over the noisy regions of the dataset, but will
     directly go onto applying the manually input value.

     This option is useful in crowded images where there is no blank sky
     to find the sky pseudo-detections.  You can get this value on a
     similarly reduced dataset (from another region of the Sky with more
     undetected regions spaces).

‘-d FLT’
‘--detgrowquant=FLT’
     Quantile limit to "grow" the final detections.  As discussed in the
     previous options, after applying the initial quantile threshold,
     layers of pixels are carved off the objects to identify true
     signal.  With this step you can return those low surface brightness
     layers that were carved off back to the detections.  To disable
     growth, set the value of this option to ‘1’.

     The process is as follows: after the true detections are found, all
     the non-detected pixels above this quantile will be put in a list
     and used to "grow" the true detections (seeds of the growth).  Like
     all quantile thresholds, this threshold is defined and applied to
     the convolved dataset.  Afterwards, the dataset is dilated once
     (with minimum connectivity) to connect very thin regions on the
     boundary: imagine building a dam at the point rivers spill into an
     open sea/ocean.  Finally, all holes are filled.  In the geography
     metaphor, holes can be seen as the closed (by the dams) rivers and
     lakes, so this process is like turning the water in all such rivers
     and lakes into soil.  See ‘--detgrowmaxholesize’ for configuring
     the hole filling.

     Note that since the growth occurs on all neighbors of a data
     element, the quantile for 3D detection must be must larger than
     that of 2D detection.  Recall that in 2D each element has 8
     neighbors while in 3D there are 27 neighbors.

‘--detgrowmaxholesize=INT’
     The maximum hole size to fill during the final expansion of the
     true detections as described in ‘--detgrowquant’.  This is
     necessary when the input contains many smaller objects and can be
     used to avoid marking blank sky regions as detections.

     For example, multiple galaxies can be positioned such that they
     surround an empty region of sky.  If all the holes are filled, the
     Sky region in between them will be taken as a detection which is
     not desired.  To avoid such cases, the integer given to this option
     must be smaller than the hole between such objects.  However, we
     should caution that unless the "hole" is very large, the combined
     faint wings of the galaxies might actually be present in between
     them, so be very careful in not filling such holes.

     On the other hand, if you have a very large (and extended) galaxy,
     the diffuse wings of the galaxy may create very large holes over
     the detections.  In such cases, a large enough value to this option
     will cause all such holes to be detected as part of the large
     galaxy and thus help in detecting it to extremely low surface
     brightness limits.  Therefore, especially when large and extended
     objects are present in the image, it is recommended to give this
     option (very) large values.  For one real-world example, see *note
     Detecting large extended targets::.

‘--cleangrowndet’
     After dilation, if the signal-to-noise ratio of a detection is less
     than the derived pseudo-detection S/N limit, that detection will be
     discarded.  In an ideal/clean noise, a true detection's S/N should
     be larger than its constituent pseudo-detections because its area
     is larger and it also covers more signal.  However, on a false
     detections (especially at lower ‘--snquant’ values), the increase
     in size can cause a decrease in S/N below that threshold.

     This will improve purity and not change completeness (a true
     detection will not be discarded).  Because a true detection has
     flux in its vicinity and dilation will catch more of that flux and
     increase the S/N. So on a true detection, the final S/N cannot be
     less than pseudo-detections.

     However, in many real images bad processing creates artifacts that
     cannot be accurately removed by the Sky subtraction.  In such
     cases, this option will decrease the completeness (will
     artificially discard true detections).  So this feature is not
     default and should to be explicitly called when you know the noise
     is clean.

‘--checkdetection’
     Every step of the detection process will be added as an extension
     to a file with the suffix ‘_det.fits’.  Going through each would
     just be a repeat of the explanations above and also of those in
     Akhlaghi and Ichikawa (2015).  The extension label should be
     sufficient to recognize which step you are observing.  Viewing all
     the steps can be the best guide in choosing the best set of
     parameters.  With this option, NoiseChisel will abort as soon as a
     snapshot of all the detection process is saved.  This behavior can
     be disabled with ‘--continueaftercheck’.

‘--checksky’
     Check the derivation of the final sky and its standard deviation
     values on the mesh grid.  With this option, NoiseChisel will abort
     as soon as the sky value is estimated over the image (on each
     tile).  This behavior can be disabled with ‘--continueaftercheck’.
     By default the output will have the same pixel size as the input,
     but with the ‘--oneelempertile’ option, only one pixel will be used
     for each tile (see *note Processing options::).


File: gnuastro.info,  Node: NoiseChisel output,  Prev: Detection options,  Up: Invoking astnoisechisel

7.2.2.3 NoiseChisel output
..........................

NoiseChisel's output is a multi-extension FITS file.  The main
extension/dataset is a (binary) detection map.  It has the same size as
the input but with only two possible values for all pixels: 0 (for
pixels identified as noise) and 1 (for those identified as
signal/detections).  The detection map is followed by a Sky and Sky
standard deviation dataset (which are calculated from the binary image).
By default (when ‘--rawoutput’ is not called), NoiseChisel will also
subtract the Sky value from the input and save the sky-subtracted input
as the first extension in the output with data.  The zero-th extension
(that contains no data), contains NoiseChisel's configuration as FITS
keywords, see *note Output FITS files::.

   The name of the output file can be set by giving a value to
‘--output’ (this is a common option between all programs and is
therefore discussed in *note Input output options::).  If ‘--output’ is
not used, the input name will be suffixed with ‘_detected.fits’ and used
as output, see *note Automatic output::.  If any of the options starting
with ‘--check*’ are given, NoiseChisel will not complete and will abort
as soon as the respective check images are created.  For more
information on the different check images, see the description for the
‘--check*’ options in *note Detection options:: (this can be disabled
with ‘--continueaftercheck’).

   The last two extensions of the output are the Sky and its Standard
deviation, see *note Sky value:: for a complete explanation.  They are
calculated on the tile grid that you defined for NoiseChisel.  By
default these datasets will have the same size as the input, but with
all the pixels in one tile given one value.  To be more space-efficient
(keep only one pixel per tile), you can use the ‘--oneelempertile’
option, see *note Tessellation::.

   To inspect any of NoiseChisel's output files, assuming you use SAO
DS9, you can configure your Graphic User Interface (GUI) to open
NoiseChisel's output as a multi-extension data cube.  This will allow
you to flip through the different extensions and visually inspect the
results.  This process has been described for the GNOME GUI (most common
GUI in GNU/Linux operating systems) in *note Viewing FITS file contents
with DS9 or TOPCAT::.

   NoiseChisel's output configuration options are described in detail
below.

‘--continueaftercheck’
     Continue NoiseChisel after any of the options starting with
     ‘--check’ (see *note Detection options::.  NoiseChisel involves
     many steps and as a result, there are many checks, allowing you to
     inspect the status of the processing.  The results of each step
     affect the next steps of processing.  Therefore, when you want to
     check the status of the processing at one step, the time spent to
     complete NoiseChisel is just wasted/distracting time.

     To encourage easier experimentation with the option values, when
     you use any of the NoiseChisel options that start with ‘--check’,
     NoiseChisel will abort once its desired extensions have been
     written.  With ‘--continueaftercheck’ option, you can disable this
     behavior and ask NoiseChisel to continue with the rest of the
     processing, even after the requested check files are complete.

‘--ignoreblankintiles’
     Do not set the input's blank pixels to blank in the tiled outputs
     (for example, Sky and Sky standard deviation extensions of the
     output).  This is only applicable when the tiled output has the
     same size as the input, in other words, when ‘--oneelempertile’ is
     not called.

     By default, blank values in the input (commonly on the edges which
     are outside the survey/field area) will be set to blank in the
     tiled outputs also.  But in other scenarios this default behavior
     is not desired; for example, if you have masked something in the
     input, but want the tiled output under that also.

‘-l’
‘--label’
     Run a connected-components algorithm on the finally detected pixels
     to identify which pixels are connected to which.  By default the
     main output is a binary dataset with only two values: 0 (for noise)
     and 1 (for signal/detections).  See *note NoiseChisel output:: for
     more.

     The purpose of NoiseChisel is to detect targets that are extended
     and diffuse, with outer parts that sink into the noise very
     gradually (galaxies and stars for example).  Since NoiseChisel digs
     down to extremely low surface brightness values, many such targets
     will commonly be detected together as a single large body of
     connected pixels.

     To properly separate connected objects, sophisticated segmentation
     methods are commonly necessary on NoiseChisel's output.  Gnuastro
     has the dedicated *note Segment:: program for this job.  Since
     input images are commonly large and can take a significant volume,
     the extra volume necessary to store the labels of the connected
     components in the detection map (which will be created with this
     ‘--label’ option, in 32-bit signed integer type) can thus be a
     major waste of space.  Since the default output is just a binary
     dataset, an 8-bit unsigned dataset is enough.

     The binary output will also encourage users to segment the result
     separately prior to doing higher-level analysis.  As an alternative
     to ‘--label’, if you have the binary detection image, you can use
     the ‘connected-components’ operator in Gnuastro's Arithmetic
     program to identify regions that are connected with each other.
     For example, with this command (assuming NoiseChisel's output is
     called ‘nc.fits’):

          $ astarithmetic nc.fits 2 connected-components -hDETECTIONS

‘--rawoutput’
     Do not include the Sky-subtracted input image as the first
     extension of the output.  By default, the Sky-subtracted input is
     put in the first extension of the output.  The next extensions are
     NoiseChisel's main outputs described above.

     The extra Sky-subtracted input can be convenient in checking
     NoiseChisel's output and comparing the detection map with the
     input: visually see if everything you expected is detected
     (reasonable completeness) and that you do not have too many false
     detections (reasonable purity).  This visual inspection is
     simplified if you use SAO DS9 to view NoiseChisel's output as a
     multi-extension data-cube, see *note Viewing FITS file contents
     with DS9 or TOPCAT::.

     When you are satisfied with your NoiseChisel configuration
     (therefore you do not need to check on every run), or you want to
     archive/transfer the outputs, or the datasets become large, or you
     are running NoiseChisel as part of a pipeline, this Sky-subtracted
     input image can be a significant burden (take up a large volume).
     The fact that the input is also noisy, makes it hard to compress it
     efficiently.

     In such cases, this ‘--rawoutput’ can be used to avoid the extra
     sky-subtracted input in the output.  It is always possible to
     easily produce the Sky-subtracted dataset from the input (assuming
     it is in extension ‘1’ of ‘in.fits’) and the ‘SKY’ extension of
     NoiseChisel's output (let's call it ‘nc.fits’) with a command like
     below (assuming NoiseChisel was not run with ‘--oneelempertile’,
     see *note Tessellation::):

          $ astarithmetic in.fits nc.fits - -h1 -hSKY

*Save space:* with the ‘--rawoutput’ and ‘--oneelempertile’,
NoiseChisel's output will only be one binary detection map and two much
smaller arrays with one value per tile.  Since none of these have noise
they can be compressed very effectively (without any loss of data) with
exceptionally high compression ratios.  This makes it easy to archive,
or transfer, NoiseChisel's output even on huge datasets.  To compress it
with the most efficient method (take up less volume), run the following
command:

     $ gzip --best noisechisel_output.fits

The resulting ‘.fits.gz’ file can then be fed into any of Gnuastro's
programs directly, or viewed in viewers like SAO DS9, without having to
decompress it separately (they will just take a little longer, because
they have to internally decompress it before starting).  See *note
NoiseChisel optimization for storage:: for an example on a real dataset.


File: gnuastro.info,  Node: Segment,  Next: MakeCatalog,  Prev: NoiseChisel,  Up: Data analysis

7.3 Segment
===========

Once signal is separated from noise (for example, with *note
NoiseChisel::), you have a binary dataset: each pixel is either signal
(1) or noise (0).  Signal (for example, every galaxy in your image) has
been "detected", but all detections have a label of 1.  Therefore while
we know which pixels contain signal, we still cannot find out how many
galaxies they contain or which detected pixels correspond to which
galaxy.  At the lowest (most generic) level, detection is a kind of
segmentation (segmenting the whole dataset into signal and noise, see
*note NoiseChisel::).  Here, we will define segmentation only on signal:
to separate sub-structure within the detections.

   If the targets are clearly separated, or their detected regions are
not touching, a simple connected components(1) algorithm (very basic
segmentation) is enough to separate the regions that are
touching/connected.  This is such a basic and simple form of
segmentation that Gnuastro's Arithmetic program has an operator for it:
see ‘connected-components’ in *note Arithmetic operators::.  Assuming
the binary dataset is called ‘binary.fits’, you can use it with a
command like this:

     $ astarithmetic binary.fits 2 connected-components

You can even do a very basic detection (a threshold, say at value ‘100’)
_and_ segmentation in Arithmetic with a single command like below:

     $ astarithmetic in.fits 100 gt 2 connected-components

   However, in most astronomical situations our targets are not nicely
separated or have a sharp boundary/edge (for a threshold to suffice):
they touch (for example, merging galaxies), or are simply in the same
line-of-sight (which is much more common).  This causes their images to
overlap.

   In particular, when you do your detection with NoiseChisel, you will
detect signal to very low surface brightness limits: deep into the faint
wings of galaxies or bright stars (which can extend very far and
irregularly from their center).  Therefore, it often happens that
several galaxies are detected as one large detection.  Since they are
touching, a simple connected components algorithm will not suffice.  It
is therefore necessary to do a more sophisticated segmentation and break
up the detected pixels (even those that are touching) into multiple
target objects as accurately as possible.

   Segment will use a detection map and its corresponding dataset to
find sub-structure over the detected areas and use them for its
segmentation.  Until Gnuastro version 0.6 (released in 2018), Segment
was part of *note NoiseChisel::.  Therefore, similar to NoiseChisel, the
best place to start reading about Segment and understanding what it does
(with many illustrative figures) is Section 3.2 of Akhlaghi and Ichikawa
2015 (https://arxiv.org/abs/1505.01664), and continue with Akhlaghi 2019
(https://arxiv.org/abs/1909.11230).

   As a summary, Segment first finds true _clump_s over the detections.
Clumps are associated with local maxima/minima(2) and extend over the
neighboring pixels until they reach a local minimum/maximum
(_river_/_watershed_).  By default, Segment will use the distribution of
clump signal-to-noise ratios over the undetected regions as reference to
find "true" clumps over the detections.  Using the undetected regions
can be disabled by directly giving a signal-to-noise ratio to
‘--clumpsnthresh’.

   The true clumps are then grown to a certain threshold over the
detections.  Based on the strength of the connections
(rivers/watersheds) between the grown clumps, they are considered parts
of one _object_ or as separate _object_s.  See Section 3.2 of Akhlaghi
and Ichikawa 2015 (https://arxiv.org/abs/1505.01664) for more.
Segment's main output are thus two labeled datasets: 1) clumps, and 2)
objects.  See *note Segment output:: for more.

   To start learning about Segment, especially in relation to detection
(*note NoiseChisel::) and measurement (*note MakeCatalog::), the
recommended references are Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664), Akhlaghi 2016
(https://arxiv.org/abs/1611.06387) and Akhlaghi 2019
(https://arxiv.org/abs/1909.11230).  If you have used Segment within
your research, please run it with ‘--cite’ to list the papers you should
cite and how to acknowledge its funding sources.

   Those papers cannot be updated any more but the software will evolve.
For example, Segment became a separate program (from NoiseChisel) in
2018 (after those papers were published).  Therefore this book is the
definitive reference.  Finally, in *note Invoking astsegment::, we will
discuss Segment's inputs, outputs and configuration options.

* Menu:

* Invoking astsegment::         Inputs, outputs and options to Segment

   ---------- Footnotes ----------

   (1) <https://en.wikipedia.org/wiki/Connected-component_labeling>

   (2) By default the maximum is used as the first clump pixel, to
define clumps based on local minima, use the ‘--minima’ option.


File: gnuastro.info,  Node: Invoking astsegment,  Prev: Segment,  Up: Segment

7.3.1 Invoking Segment
----------------------

Segment will identify substructure within the detected regions of an
input image.  Segment's output labels can be directly used for
measurements (for example, with *note MakeCatalog::).  The executable
name is ‘astsegment’ with the following general template

     $ astsegment [OPTION ...] InputImage.fits

One line examples:

     ## Segment NoiseChisel's detected regions.
     $ astsegment default-noisechisel-output.fits

     ## Use a hand-input S/N value for keeping true clumps
     ## (avoid finding the S/N using the undetected regions).
     $ astsegment nc-out.fits --clumpsnthresh=10

     ## Inspect all the segmentation steps after changing a parameter.
     $ astsegment input.fits --snquant=0.9 --checksegmentaion

     ## Use the fixed value of 0.01 for the input's Sky standard deviation
     ## (in the units of the input), and assume all the pixels are a
     ## detection (for example, a large structure extending over the whole
     ## image), and only keep clumps with S/N>10 as true clumps.
     $ astsegment in.fits --std=0.01 --detection=all --clumpsnthresh=10

If Segment is to do processing (for example, you do not want to get
help, or see the values of each option), at least one input dataset is
necessary along with detection and error information, either as separate
datasets (per-pixel) or fixed values, see *note Segment input::.
Segment shares a large set of common operations with other Gnuastro
programs, mainly regarding input/output, general processing steps, and
general operating modes.  To help in a unified experience between all of
Gnuastro's programs, these common operations have the same names and
defined in *note Common options::.

   As in all Gnuastro programs, options can also be given to Segment in
configuration files.  For a thorough description of Gnuastro's
configuration file parsing, please see *note Configuration files::.  All
of Segment's options with a short description are also always available
on the command-line with the ‘--help’ option, see *note Getting help::.
To inspect the option values without actually running Segment, append
your command with ‘--printparams’ (or ‘-P’).

   To help in easy navigation between Segment's options, they are
separately discussed in the three sub-sections below: *note Segment
input:: discusses how you can customize the inputs to Segment.  *note
Segmentation options:: is devoted to options specific to the high-level
segmentation process.  Finally, in *note Segment output::, we will
discuss options that affect Segment's output.

* Menu:

* Segment input::               Input files and options.
* Segmentation options::        Parameters of the segmentation process.
* Segment output::              Outputs of Segment


File: gnuastro.info,  Node: Segment input,  Next: Segmentation options,  Prev: Invoking astsegment,  Up: Invoking astsegment

7.3.1.1 Segment input
.....................

Besides the input dataset (for example, astronomical image), Segment
also needs to know the Sky standard deviation and the regions of the
dataset that it should segment.  The values dataset is assumed to be Sky
subtracted by default.  If it is not, you can ask Segment to subtract
the Sky internally by calling ‘--sky’.  For the rest of this discussion,
we will assume it is already sky subtracted.

   The Sky and its standard deviation can be a single value (to be used
for the whole dataset) or a separate dataset (for a separate value per
pixel).  If a dataset is used for the Sky and its standard deviation,
they must either be the size of the input image, or have a single value
per tile (generated with ‘--oneelempertile’, see *note Processing
options:: and *note Tessellation::).

   The detected regions/pixels can be specified as a detection map (for
example, see *note NoiseChisel output::).  If ‘--detection=all’, Segment
will not read any detection map and assume the whole input is a single
detection.  For example, when the dataset is fully covered by a large
nearby galaxy/globular cluster.

   When dataset are to be used for any of the inputs, Segment will
assume they are multiple extensions of a single file by default (when
‘--std’ or ‘--detection’ are not called).  For example, NoiseChisel's
default output *note NoiseChisel output::.  When the Sky-subtracted
values are in one file, and the detection and Sky standard deviation are
in another, you just need to use ‘--detection’: in the absence of
‘--std’, Segment will look for both the detection labels and Sky
standard deviation in the file given to ‘--detection’.  Ultimately, if
all three are in separate files, you need to call both ‘--detection’ and
‘--std’.

   The extensions of the three mandatory inputs can be specified with
‘--hdu’, ‘--dhdu’, and ‘--stdhdu’.  For a full discussion on what to
give to these options, see the description of ‘--hdu’ in *note Input
output options::.  To see their default values (along with all the other
options), run Segment with the ‘--printparams’ (or ‘-P’) option.  Just
recall that in the absence of ‘--detection’ and ‘--std’, all three are
assumed to be in the same file.  If you only want to see Segment's
default values for HDUs on your system, run this command:

     $ astsegment -P | grep hdu

   By default Segment will convolve the input with a kernel to improve
the signal-to-noise ratio of true peaks.  If you already have the
convolved input dataset, you can pass it directly to Segment for faster
processing (using the ‘--convolved’ and ‘--chdu’ options).  Just do not
forget that the convolved image must also be Sky-subtracted before
calling Segment.  If a value/file is given to ‘--sky’, the convolved
values will also be Sky subtracted internally.  Alternatively, if you
prefer to give a kernel (with ‘--kernel’ and ‘--khdu’), Segment can do
the convolution internally.  To disable convolution, use
‘--kernel=none’.

‘--sky=STR/FLT’
     The Sky value(s) to subtract from the input.  This option can
     either be given a constant number or a file name containing a
     dataset (multiple values, per pixel or per tile).  By default,
     Segment will assume the input dataset is Sky subtracted, so this
     option is not mandatory.

     If the value cannot be read as a number, it is assumed to be a file
     name.  When the value is a file, the extension can be specified
     with ‘--skyhdu’.  When it is not a single number, the given dataset
     must either have the same size as the output or the same size as
     the tessellation (so there is one pixel per tile, see *note
     Tessellation::).

     When this option is given, its value(s) will be subtracted from the
     input and the (optional) convolved dataset (given to ‘--convolved’)
     prior to starting the segmentation process.

‘--skyhdu=STR/INT’
     The HDU/extension containing the Sky values.  This is mandatory
     when the value given to ‘--sky’ is not a number.  Please see the
     description of ‘--hdu’ in *note Input output options:: for the
     different ways you can identify a special extension.

‘--std=STR/FLT’
     The Sky standard deviation value(s) corresponding to the input.
     The value can either be a constant number or a file name containing
     a dataset (multiple values, per pixel or per tile).  The Sky
     standard deviation is mandatory for Segment to operate.

     If the value cannot be read as a number, it is assumed to be a file
     name.  When the value is a file, the extension can be specified
     with ‘--skyhdu’.  When it is not a single number, the given dataset
     must either have the same size as the output or the same size as
     the tessellation (so there is one pixel per tile, see *note
     Tessellation::).

     When this option is not called, Segment will assume the standard
     deviation is a dataset and in a HDU/extension (‘--stdhdu’) of
     another one of the input file(s).  If a file is given to
     ‘--detection’, it will assume that file contains the standard
     deviation dataset, otherwise, it will look into input filename (the
     main argument, without any option).

‘--stdhdu=INT/STR’
     The HDU/extension containing the Sky standard deviation values,
     when the value given to ‘--std’ is a file name.  Please see the
     description of ‘--hdu’ in *note Input output options:: for the
     different ways you can identify a special extension.

‘--variance’
     The input Sky standard deviation value/dataset is actually
     variance.  When this option is called, the square root of input Sky
     standard deviation (see ‘--std’) is used internally, not its raw
     value(s).

‘-d FITS’
‘--detection=FITS’
     Detection map to use for segmentation.  If given a value of ‘all’,
     Segment will assume the whole dataset must be segmented, see below.
     If a detection map is given, the extension can be specified with
     ‘--dhdu’.  If not given, Segment will assume the desired
     HDU/extension is in the main input argument (input file specified
     with no option).

     The final segmentation (clumps or objects) will only be over the
     non-zero pixels of this detection map.  The dataset must have the
     same size as the input image.  Only datasets with an integer type
     are acceptable for the labeled image, see *note Numeric data
     types::.  If your detection map only has integer values, but it is
     stored in a floating point container, you can use Gnuastro's
     Arithmetic program (see *note Arithmetic::) to convert it to an
     integer container, like the example below:

          $ astarithmetic float.fits int32 --output=int.fits

     It may happen that the whole input dataset is covered by signal,
     for example, when working on parts of the Andromeda galaxy, or
     nearby globular clusters (that cover the whole field of view).  In
     such cases, segmentation is necessary over the complete dataset,
     not just specific regions (detections).  By default Segment will
     first use the undetected regions as a reference to find the proper
     signal-to-noise ratio of "true" clumps (give a purity level
     specified with ‘--snquant’).  Therefore, in such scenarios you also
     need to manually give a "true" clump signal-to-noise ratio with the
     ‘--clumpsnthresh’ option to disable looking into the undetected
     regions, see *note Segmentation options::.  In such cases, is
     possible to make a detection map that only has the value ‘1’ for
     all pixels (for example, using *note Arithmetic::), but for
     convenience, you can also use ‘--detection=all’.

‘--dhdu’
     The HDU/extension containing the detection map given to
     ‘--detection’.  Please see the description of ‘--hdu’ in *note
     Input output options:: for the different ways you can identify a
     special extension.

‘-k FITS’
‘--kernel=FITS’
     The name of file containing kernel that will be used to convolve
     the input image.  The usage of this option is identical to
     NoiseChisel's ‘--kernel’ option (*note NoiseChisel input::).
     Please see the descriptions there for more.  To disable
     convolution, you can give it a value of ‘none’.

‘--khdu’
     The HDU/extension containing the kernel used for convolution.  For
     acceptable values, please see the description of ‘--hdu’ in *note
     Input output options::.

‘--convolved=FITS’
     The convolved image's file name to avoid internal convolution by
     Segment.  The usage of this option is identical to NoiseChisel's
     ‘--convolved’ option.  Please see *note NoiseChisel input:: for a
     thorough discussion of the usefulness and best practices of using
     this option.

     If you want to use the same convolution kernel for detection (with
     *note NoiseChisel::) and segmentation, with this option, you can
     use the same convolved image (that is also available in
     NoiseChisel) and avoid two convolutions.  However, just be careful
     to use the input to NoiseChisel as the input to Segment also, then
     use the ‘--sky’ and ‘--std’ to specify the Sky and its standard
     deviation (from NoiseChisel's output).  Recall that when
     NoiseChisel is not called with ‘--rawoutput’, the first extension
     of NoiseChisel's output is the _Sky-subtracted_ input (see *note
     NoiseChisel output::).  So if you use the same convolved image that
     you fed to NoiseChisel, but use NoiseChisel's output with Segment's
     ‘--convolved’, then the convolved image will not be Sky subtracted.

‘--chdu’
     The HDU/extension containing the convolved image (given to
     ‘--convolved’).  For acceptable values, please see the description
     of ‘--hdu’ in *note Input output options::.

‘-L INT[,INT]’
‘--largetilesize=INT[,INT]’
     The size of the large tiles to use for identifying the clump S/N
     threshold over the undetected regions.  The usage of this option is
     identical to NoiseChisel's ‘--largetilesize’ option (*note
     NoiseChisel input::).  Please see the descriptions there for more.

     The undetected regions can be a significant fraction of the dataset
     and finding clumps requires sorting of the desired regions, which
     can be slow.  To speed up the processing, Segment finds clumps in
     the undetected regions over separate large tiles.  This allows it
     to have to sort a much smaller set of pixels and also to treat them
     independently and in parallel.  Both these issues greatly speed it
     up.  Just be sure to not decrease the large tile sizes too much
     (less than 100 pixels in each dimension).  It is important for them
     to be much larger than the clumps.


File: gnuastro.info,  Node: Segmentation options,  Next: Segment output,  Prev: Segment input,  Up: Invoking astsegment

7.3.1.2 Segmentation options
............................

The options below can be used to configure every step of the
segmentation process in the Segment program.  For a more complete
explanation (with figures to demonstrate each step), please see Section
3.2 of Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664),
and also *note Segment::.  By default, Segment will follow the procedure
described in the paper to find the S/N threshold based on the noise
properties.  This can be disabled by directly giving a trustable
signal-to-noise ratio to the ‘--clumpsnthresh’ option.

   Recall that you can always see the full list of Gnuastro's options
with the ‘--help’ (see *note Getting help::), or ‘--printparams’ (or
‘-P’) to see their values (see *note Operating mode options::).

‘-B FLT’
‘--minskyfrac=FLT’
     Minimum fraction (value between 0 and 1) of Sky (undetected) areas
     in a large tile.  Only (large) tiles with a fraction of undetected
     pixels (Sky) greater than this value will be used for finding
     clumps.  The clumps found in the undetected areas will be used to
     estimate a S/N threshold for true clumps.  Therefore this is an
     important option (to decrease) in crowded fields.  Operationally,
     this is almost identical to NoiseChisel's ‘--minskyfrac’ option
     (*note Detection options::).  Please see the descriptions there for
     more.

‘--minima’
     Build the clumps based on the local minima, not maxima.  By
     default, clumps are built starting from local maxima (see Figure 8
     of Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664)).
     Therefore, this option can be useful when you are searching for
     true local minima (for example, absorption features).

‘-m INT’
‘--snminarea=INT’
     The minimum area which a clump in the undetected regions should
     have in order to be considered in the clump Signal to noise ratio
     measurement.  If this size is set to a small value, the Signal to
     noise ratio of false clumps will not be accurately found.  It is
     recommended that this value be larger than the value to
     NoiseChisel's ‘--snminarea’.  Because the clumps are found on the
     convolved (smoothed) image while the pseudo-detections are found on
     the input image.  You can use ‘--checksn’ and ‘--checksegmentation’
     to see if your chosen value is reasonable or not.

‘--checksn’
     Save the S/N values of the clumps over the sky and detected regions
     into separate tables.  If ‘--tableformat’ is a FITS format, each
     table will be written into a separate extension of one file
     suffixed with ‘_clumpsn.fits’.  If it is plain text, a separate
     file will be made for each table (ending in ‘_clumpsn_sky.txt’ and
     ‘_clumpsn_det.txt’).  For more on ‘--tableformat’ see *note Input
     output options::.

     You can use these tables to inspect the S/N values and their
     distribution (in combination with the ‘--checksegmentation’ option
     to see where the clumps are).  You can use Gnuastro's *note
     Statistics:: to make a histogram of the distribution (ready for
     plotting in a text file, or a crude ASCII-art demonstration on the
     command-line).

     With this option, Segment will abort as soon as the two tables are
     created.  This allows you to inspect the steps leading to the final
     S/N quantile threshold, this behavior can be disabled with
     ‘--continueaftercheck’.

‘--minnumfalse=INT’
     The minimum number of clumps over undetected (Sky) regions to
     identify the requested Signal-to-Noise ratio threshold.
     Operationally, this is almost identical to NoiseChisel's
     ‘--minnumfalse’ option (*note Detection options::).  Please see the
     descriptions there for more.

‘-c FLT’
‘--snquant=FLT’
     The quantile of the signal-to-noise ratio distribution of clumps in
     undetected regions, used to define true clumps.  After identifying
     all the usable clumps in the undetected regions of the dataset, the
     given quantile of their signal-to-noise ratios is used to define
     the signal-to-noise ratio of a "true" clump.  Effectively, this can
     be seen as an inverse p-value measure.  See Figure 9 and Section
     3.2.1 of Akhlaghi and Ichikawa 2015
     (https://arxiv.org/abs/1505.01664) for a complete explanation.  The
     full distribution of clump signal-to-noise ratios over the
     undetected areas can be saved into a table with ‘--checksn’ option
     and visually inspected with ‘--checksegmentation’.

‘-v’
‘--keepmaxnearriver’
     Keep a clump whose maximum (minimum if ‘--minima’ is called) flux
     is 8-connected to a river pixel.  By default such clumps over
     detections are considered to be noise and are removed irrespective
     of their significance measure; see Akhlaghi 2019
     (https://arxiv.org/abs/1909.11230).  Over large profiles, that sink
     into the noise very slowly, noise can cause part of the profile
     (which was flat without noise) to become a very large and with a
     very high Signal to noise ratio.  In such cases, the pixel with the
     maximum flux in the clump will be immediately touching a river
     pixel.

‘-s FLT’
‘--clumpsnthresh=FLT’
     The signal-to-noise threshold for true clumps.  If this option is
     given, then the segmentation options above will be ignored and the
     given value will be directly used to identify true clumps over the
     detections.  This can be useful if you have a large dataset with
     similar noise properties.  You can find a robust signal-to-noise
     ratio based on a (sufficiently large) smaller portion of the
     dataset.  Afterwards, with this option, you can speed up the
     processing on the whole dataset.  Other scenarios where this option
     may be useful is when, the image might not contain enough/any Sky
     regions.

‘-G FLT’
‘--gthresh=FLT’
     Threshold (multiple of the sky standard deviation added with the
     sky) to stop growing true clumps.  Once true clumps are found, they
     are set as the basis to segment the detected region.  They are
     grown until the threshold specified by this option.

‘-y INT’
‘--minriverlength=INT’
     The minimum length of a river between two grown clumps for it to be
     considered in signal-to-noise ratio estimations.  Similar to
     ‘--snminarea’, if the length of the river is too short, the
     signal-to-noise ratio can be noisy and unreliable.  Any existing
     rivers shorter than this length will be considered as non-existent,
     independent of their Signal to noise ratio.  The clumps are grown
     on the input image, therefore this value can be smaller than the
     value given to ‘--snminarea’.  Recall that the clumps were defined
     on the convolved image so ‘--snminarea’ should be larger.

‘-O FLT’
‘--objbordersn=FLT’
     The maximum Signal to noise ratio of the rivers between two grown
     clumps in order to consider them as separate 'objects'.  If the
     Signal to noise ratio of the river between two grown clumps is
     larger than this value, they are defined to be part of one
     'object'.  Note that the physical reality of these 'objects' can
     never be established with one image, or even multiple images from
     one broad-band filter.  Any method we devise to define 'object's
     over a detected region is ultimately subjective.

     Two very distant galaxies or satellites in one halo might lie in
     the same line of sight and be detected as clumps on one detection.
     On the other hand, the connection (through a spiral arm or tidal
     tail for example) between two parts of one galaxy might have such a
     low surface brightness that they are broken up into multiple
     detections or objects.  In fact if you have noticed, exactly for
     this purpose, this is the only Signal to noise ratio that the user
     gives into NoiseChisel.  The 'true' detections and clumps can be
     objectively identified from the noise characteristics of the image,
     so you do not have to give any hand input Signal to noise ratio.

‘--checksegmentation’
     A file with the suffix ‘_seg.fits’ will be created.  This file
     keeps all the relevant steps in finding true clumps and segmenting
     the detections into multiple objects in various extensions.  Having
     read the paper or the steps above.  Examining this file can be an
     excellent guide in choosing the best set of parameters.  Note that
     calling this function will significantly slow NoiseChisel.  In
     verbose mode (without the ‘--quiet’ option, see *note Operating
     mode options::) the important steps (along with their extension
     names) will also be reported.

     With this option, NoiseChisel will abort as soon as the two tables
     are created.  This behavior can be disabled with
     ‘--continueaftercheck’.


File: gnuastro.info,  Node: Segment output,  Prev: Segmentation options,  Up: Invoking astsegment

7.3.1.3 Segment output
......................

The main output of Segment are two label datasets (with integer types,
separating the dataset's elements into different classes).  They have
HDU/extension names of ‘CLUMPS’ and ‘OBJECTS’.

   Similar to all Gnuastro's FITS outputs, the zero-th extension/HDU of
the main output file only contains header keywords and image or table.
It contains the Segment input files and parameters (option names and
values) as FITS keywords.  Note that if an option name is longer than 8
characters, the keyword name is the second word.  The first word is
‘HIERARCH’.  Also note that according to the FITS standard, the keyword
names must be in capital letters, therefore, if you want to use Grep to
inspect these keywords, use the ‘-i’ option, like the example below.

     $ astfits image_segmented.fits -h0 | grep -i snquant

   By default, besides the ‘CLUMPS’ and ‘OBJECTS’ extensions, Segment's
output will also contain the (technically redundant) sky-subtracted
input dataset (‘INPUT-NO-SKY’) and the sky standard deviation dataset
(‘SKY_STD’, if it was not a constant number).  This can help in visually
inspecting the result when viewing the images as a "Multi-extension data
cube" in SAO DS9 for example, (see *note Viewing FITS file contents with
DS9 or TOPCAT::).  You can simply flip through the extensions and see
the same region of the image and its corresponding clumps/object labels.
It also makes it easy to feed the output (as one file) into MakeCatalog
when you intend to make a catalog afterwards (see *note MakeCatalog::.
To remove these redundant extensions from the output (for example, when
designing a pipeline), you can use ‘--rawoutput’.

   The ‘OBJECTS’ and ‘CLUMPS’ extensions can be used as input into *note
MakeCatalog:: to generate a catalog for higher-level analysis.  If you
want to treat each clump separately, you can give a very large value (or
even a NaN, which will always fail) to the ‘--gthresh’ option (for
example, ‘--gthresh=1e10’ or ‘--gthresh=nan’), see *note Segmentation
options::.

   For a complete definition of clumps and objects, please see Section
3.2 of Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664) and
*note Segmentation options::.  The clumps are "true" local maxima
(minima if ‘--minima’ is called) and their surrounding pixels until a
local minimum/maximum (caused by noise fluctuations, or another "true"
clump).  Therefore it may happen that some of the input detections are
not covered by clumps at all (very diffuse objects without any strong
peak), while some objects may contain many clumps.  Even in those that
have clumps, there will be regions that are too diffuse.  The diffuse
regions (within the input detected regions) are given a negative label
(-1) to help you separate them from the undetected regions (with a value
of zero).

   Each clump is labeled with respect to its host object.  Therefore, if
an object has three clumps for example, the clumps within it have labels
1, 2 and 3.  As a result, if an initial detected region has multiple
objects, each with a single clump, all the clumps will have a label of
1.  The total number of clumps in the dataset is stored in the ‘NCLUMPS’
keyword of the ‘CLUMPS’ extension and printed in the verbose output of
Segment (when ‘--quiet’ is not called).

   The ‘OBJECTS’ extension of the output will give a positive
counter/label to every detected pixel in the input.  As described in
Akhlaghi and Ichikawa [2015], the true clumps are grown until a certain
threshold.  If the grown clumps touch other clumps and the connection is
strong enough, they are considered part of the same _object_.  Once
objects (grown clumps) are identified, they are grown to cover the whole
detected area.

   The options to configure the output of Segment are listed below:

‘--continueaftercheck’
     Do not abort Segment after producing the check image(s).  The usage
     of this option is identical to NoiseChisel's ‘--continueaftercheck’
     option (*note NoiseChisel input::).  Please see the descriptions
     there for more.

‘--noobjects’
     Abort Segment after finding true clumps and do not continue with
     finding options.  Therefore, no ‘OBJECTS’ extension will be present
     in the output.  Each true clump in ‘CLUMPS’ will get a unique
     label, but diffuse regions will still have a negative value.

     To make a catalog of the clumps, the input detection map (where all
     the labels are one) can be fed into *note MakeCatalog:: along with
     the input detection map to Segment (that only had a value of ‘1’
     for all detected pixels) with ‘--clumpscat’.  In this way,
     MakeCatalog will assume all the clumps belong to a single "object".

‘--grownclumps’
     In the output ‘CLUMPS’ extension, store the grown clumps.  If a
     detected region contains no clumps or only one clump, then it will
     be fully given a label of ‘1’ (no negative valued pixels).

‘--rawoutput’
     Only write the ‘CLUMPS’ and ‘OBJECTS’ datasets in the output file.
     Without this option (by default), the first and last extensions of
     the output will the Sky-subtracted input dataset and the Sky
     standard deviation dataset (if it was not a number).  When the
     datasets are small, these redundant extensions can make it
     convenient to inspect the results visually or feed the output to
     *note MakeCatalog:: for measurements.  Ultimately both the input
     and Sky standard deviation datasets are redundant (you had them
     before running Segment).  When the inputs are large/numerous, these
     extra dataset can be a burden.

*Save space:* with the ‘--rawoutput’, Segment's output will only be two
labeled datasets (only containing integers).  Since they have no noise,
such datasets can be compressed very effectively (without any loss of
data) with exceptionally high compression ratios.  You can use the
following command to compress it with the best ratio:

     $ gzip --best segment_output.fits

The resulting ‘.fits.gz’ file can then be fed into any of Gnuastro's
programs directly, without having to decompress it separately (it will
just take them a little longer, because they have to decompress it
internally before use).

   When the input is a 2D image, to inspect NoiseChisel's output you can
configure SAO DS9 in your Graphic User Interface (GUI) to open
NoiseChisel's output as a multi-extension data cube.  This will allow
you to flip through the different extensions and visually inspect the
results.  This process has been described for the GNOME GUI (most common
GUI in GNU/Linux operating systems) in *note Viewing FITS file contents
with DS9 or TOPCAT::.


File: gnuastro.info,  Node: MakeCatalog,  Next: Match,  Prev: Segment,  Up: Data analysis

7.4 MakeCatalog
===============

At the lowest level, a dataset (for example, an image) is just a
collection of values, placed after each other in any number of
dimensions (for example, an image is a 2D dataset).  Each data-element
(pixel) just has two properties: its position (relative to the rest) and
its value.  In higher-level analysis, an entire dataset (an image for
example) is rarely treated as a singular entity(1).  You usually want to
know/measure the properties of the (separate) scientifically interesting
targets that are embedded in it.  For example, the magnitudes, positions
and elliptical properties of the galaxies that are in the image.

   MakeCatalog is Gnuastro's program for localized measurements over a
dataset.  In other words, MakeCatalog is Gnuastro's program to convert
low-level datasets (like images), to high level catalogs.  The role of
MakeCatalog in a scientific analysis and the benefits of its model
(where detection/segmentation is separated from measurement) is
discussed in Akhlaghi 2016 (https://arxiv.org/abs/1611.06387v1)(2) and
summarized in *note Detection and catalog production::.  We strongly
recommend reading this short paper for a better understanding of this
methodology.  Understanding the effective usage of MakeCatalog, will
thus also help effective use of other (lower-level) Gnuastro's programs
like *note NoiseChisel:: or *note Segment::.

   Following the Unix philosophy(3), MakeCatalog is specialized in doing
measurements accurately and efficiently, not to do detection,
segmentation, or defining apertures on requested positions in your
dataset.  It is therefore important to define the pixels (of each
"object" in the image for example) _before_ running MakeCatalog.  There
are separate, highly specialized and customizable programs in Gnuastro
for these other jobs as discussed below, their outputs can be fed into
MakeCatalog (for a usage example in a real-world analysis, see *note
General program usage tutorial:: and *note Detecting large extended
targets::).

   • *note Arithmetic::: Detection with a simple threshold.

   • *note NoiseChisel::: Advanced detection.

   • *note Segment::: Segmentation (substructure over detections).

   • *note MakeProfiles::: Aperture creation for known positions.

   To identify which pixels should be measured together, a "label"
dataset is necessary (which is the main input to MakeCatalog).  The
number of labels in the labels dataset will be the number of rows in the
output.  Having identified which pixels to use for each row of the
output catalog, the "values" dataset will be used to read the value of
each pixel to use in the measurement.  The labeled dataset should have
the same size/dimensions as the values image, but with integer valued
pixels.  Any pixel with the same label/counter for each sub-set of
pixels that must be measured together.  For example, all the pixels
covering one galaxy in an image, get the same label.

   The final result is a catalog where each row corresponds to the
measurements on pixels with a specific label.  For example, the flux
weighted average position of all the pixels with a label of 42 will be
written into the 42nd row of the output catalog/table's central position
column(4).  Similarly, the sum of all these pixels will be the 42nd row
in the sum column, etc.  Pixels with labels equal to, or smaller than,
zero will be ignored by MakeCatalog.  In other words, the number of rows
in MakeCatalog's output is already known before running it (the maximum
value of the labeled dataset).

   Before getting into the details of running MakeCatalog (in *note
Invoking astmkcatalog:: (last sub-section), we will start with some
basics that are good to review if you are new to optical/astronomical
data analysis.  The basics of MakeCatalog's approach to separating
detection/segmentation from measurements in is first discussed in *note
Detection and catalog production::.  We then introduce the core
units/concepts of brightness measurements in optical astronomy (which
inherits a lot from its multi-millennial history) in *note Brightness
flux magnitude::.  Having reviewed the core concepts, *note MakeCatalog
measurements on each label:: provides a categorized list of all
measurements that are currently available in MakeCatalog.

   A very important factor in any measurement is understanding its
validity range, or limits.  Therefore in *note Metameasurements on full
input::, we will discuss how to estimate the reliability of the
detection and basic measurements.  This section will continue with a
derivation of elliptical parameters from the labeled datasets in *note
Morphology measurements elliptical::.  For those who feel MakeCatalog's
existing measurements/columns are not enough and would like to add
further measurements, in *note Adding new columns to MakeCatalog::, a
checklist of steps is provided for readily adding your own new
measurements/columns.

* Menu:

* Detection and catalog production::  Discussing why/how to treat these separately.
* Brightness flux magnitude::   More on Magnitudes, surface brightness, etc.
* Standard deviation vs Standard error::  To avoid confusions with error measurements.
* MakeCatalog measurements on each label::  The columns that you can request in output catalog.
* Metameasurements on full input::  Measurements on/about measurements.
* Manual metameasurements::     These need custom runs of MakeCatalog and/or other programs.
* Adding new columns to MakeCatalog::  How to add new columns.
* Invoking astmkcatalog::       Options and arguments to MakeCatalog.

   ---------- Footnotes ----------

   (1) You can derive the over-all properties of a complete dataset (1D
table column, 2D image, or 3D data-cube) treated as a single entity with
Gnuastro's Statistics program (see *note Statistics::).

   (2) A published paper cannot undergo any more change, so this manual
is the definitive guide.

   (3) Make each program do one thing well.  To do a new job, build
afresh rather than complicate old programs by adding new "features".

   (4) See *note Morphology measurements elliptical:: for a discussion
on this and the derivation of positional parameters, which includes the
center.


File: gnuastro.info,  Node: Detection and catalog production,  Next: Brightness flux magnitude,  Prev: MakeCatalog,  Up: MakeCatalog

7.4.1 Detection and catalog production
--------------------------------------

Most existing common tools in low-level astronomical data-analysis (for
example, SExtractor(1)) merge the two processes of detection and
measurement (catalog production) in one program.  However, in light of
Gnuastro's modularized approach (modeled on the Unix system) detection
is separated from measurements and catalog production.  This modularity
is therefore new to many experienced astronomers and deserves a short
review here.  Further discussion on the benefits of this methodology can
be seen in Akhlaghi 2016 (https://arxiv.org/abs/1611.06387v1).

   As discussed in the introduction of *note MakeCatalog::, detection
(identifying which pixels to do measurements on) can be done with
different programs.  Their outputs (a labeled dataset) can be directly
fed into MakeCatalog to do the measurements and write the result as a
catalog/table.  Beyond that, Gnuastro's modular approach has many
benefits that will become clear as you get more experienced in
astronomical data analysis and want to be more creative in using your
valuable data for the exciting scientific project you are working on.
In short the reasons for this modularity can be classified as below:

   • Simplicity/robustness of independent, modular tools: making a
     catalog is a logically separate process from labeling (detection,
     segmentation, or aperture production).  A user might want to do
     certain operations on the labeled regions before creating a catalog
     for them.  Another user might want the properties of the same
     pixels/objects in another image (another filter for example) to
     measure the colors or SED fittings.

     Here is an example of doing both: suppose you have images in
     various broad band filters at various resolutions and orientations.
     The image of one color will thus not lie exactly on another or even
     be in the same scale.  However, it is imperative that the same
     pixels be used in measuring the colors of galaxies.

     To solve the problem, NoiseChisel can be run on the reference image
     to generate the labeled detection image.  Afterwards, the labeled
     image can be warped into the grid of the other color (using *note
     Warp::).  MakeCatalog will then generate the same catalog for both
     colors (with the different labeled images).  It is currently
     customary to warp the images to the same pixel grid, however,
     modification of the scientific dataset is very harmful for the data
     and creates correlated noise.  It is much more accurate to do the
     transformations on the labeled image.

   • Complexity of a monolith: Adding in a catalog functionality to the
     detector program will add several more steps (and many more
     options) to its processing that can equally well be done outside of
     it.  This makes following what the program does harder for the
     users and developers, it can also potentially add many bugs.

     As an example, if the parameter you want to measure over one
     profile is not provided by the developers of MakeCatalog.  You can
     simply open this tiny little program and add your desired
     calculation easily.  This process is discussed in *note Adding new
     columns to MakeCatalog::.  However, if making a catalog was part of
     NoiseChisel for example, adding a new column/measurement would
     require a lot of energy to understand all the steps and internal
     structures of that huge program.  It might even be so intertwined
     with its processing, that adding new columns might cause
     problems/bugs in its primary job (detection).

   ---------- Footnotes ----------

   (1) <https://www.astromatic.net/software/sextractor>


File: gnuastro.info,  Node: Brightness flux magnitude,  Next: Standard deviation vs Standard error,  Prev: Detection and catalog production,  Up: MakeCatalog

7.4.2 Brightness, Flux, Magnitude and Surface brightness
--------------------------------------------------------

After taking an image with your camera (or a large telescope), the value
in each pixel of the output file is a proxy for the amount of _energy_
that accumulated in it (while it was exposed to light).  In
astrophysics, the centimeter–gram–second (cgs) units are commonly used
so energy is commonly reported in units of $erg$.  In an image, the
energy of a galaxy for example will be distributed over many pixels.
Therefore, the collected energy of an object in the image is the total
sum of the values in the pixels associated to the object.

   To be able to compare our scientific data with other data, optical
astronomers have a unique terminology based on the concept of
"Magnitude"s.  But before getting to those, let's review the following
basic physical concepts first:

Energy ($erg$; also in _counts_ or _ADU_s)
     Within the electromagnetic regime, we measure the received energy
     of astronomical source by counting the number of photons that have
     been converted to electrons (electric potential) in our detectors.

     When counting/measuring the electric potential changes, the optical
     (but also near ultra-violet and near infra-red) detectors do not
     actually count individual electrons but bundles/packages of
     electrons known as the analog-to-digital unit (ADU). The number of
     electrons in each ADU is known as the _gain_ of the instrument and
     is measured as part of its calibration.

Power ($erg/s$)
     The amount of energy in a fixed interval of time (1 second) is
     known as power.  Power is used in two contexts within astronomy
     which are listed below.  Both have the same units of energy per
     time, but their difference is very important to understand in
     physical interpretation:

     Brightness
          The _received_ power from a source (the thing we measure).  To
          be able to compare data taken with different exposure times,
          we define the received power of the source in a detector as
          the _brightness_.

     Luminosity
          The total _emitted_ power of a source in _all directions_.

     Unlike brightness (a measured property), luminosity is an inherent
     property of the object that is calculated from the combination of
     multiple measurements (flux and distance; see below).

Flux ($erg/s/cm^2$)
     To be able to compare with data from different telescopes (with
     different collecting areas), we define the _flux_ which is defined
     by dividing the brightness by the exposed aperture of our
     telescope.  Because we are using the cgs units, the collecting area
     is reported in $cm^2$.  Knowing the flux ($f$) and distance to the
     object ($r$), we can derive its _luminosity_: $L=4{\pi}r^2f$.

Spectral flux density ($erg/s/cm^2/Hz$ or $erg/s/cm^2/\AA$)
     To take into account the different spectral coverage of filters and
     detectors, we define the _spectral flux density_, which is defined
     in either of these units (based on context): $erg/s/cm^2/Hz$
     (frequency-based) $erg/s/cm^2/\AA$ (wavelength-based).

     As in other objects in nature, astronomical objects do not emit or
     reflect the same flux at all wavelengths.  On the other hand, our
     detector techologies are different for different wavelength ranges.
     Therefore, even if we wanted to, there is no way to measure the
     "total" (at all wavelengths; also known as "bolometric") luminosity
     of an object with a single tool.  To be able to analyze objects
     with different spectral features (compare measurements of the same
     object taken in different spectral regimes), it is therefore
     important to account for the wavelength (or frequency) range of the
     photons that we measured through the spectral flux density.

     Jansky ($10^{-23}erg/s/cm^2/Hz$)
          A "Jansky" is a predefined/nominal level of frequency flux
          density that is commonly used in radio astronomy.  The AB
          magnitude system (see below; used in optical astronomy) is
          also in frequency-based so there is a simple conversion
          between the two.

          Janskys can be converted to wavelength flux density using the
          ‘jy-to-wavelength-flux-density’ operator of Gnuastro's
          Arithmetic program, see the derivation under this operator's
          description in *note Unit conversion operators::.

   Having summarized the relevant basic physical concepts above, let's
review the terminology that is used in optical astronomy.  The reason
optical astronomers don't use modern physical terminology is that
optical astronomy precedes modern physical concepts by thousands of
years!

   Once the modern physical concepts where mature enough, optical
astronomers found the correct conversion factors to better define their
own terminology (and easily use previous results) instead of abandoning
them.  Other fields of astronomy (for example X-ray or radio) were
discovered in the last century when modern physical concepts had already
matured and were being extensively used, so for those fields, the
concepts above are enough.

Magnitude
     The observed spectral flux density of astronomical objects span
     over a very large range: the Sun (as the brightest object) is
     roughly $10^{24}$ times brighter than the fainter galaxies we can
     currently detect in our deepest images.  Therefore the scale that
     was originally used from the ancient times to measure the incoming
     light (used by Hipparchus of Nicaea; 190-120 BC) can be
     parametrized as a logarithmic function of the spectral flux
     density.

     But the logarithm can only be usable with a value which is always
     positive and has no units.  Fortunately brightness is always
     positive.  To remove the units, we divide the spectral flux density
     of the object ($F$) by a reference spectral flux density ($F_r$).
     We then define a logarithmic scale through the relation below and
     call it the _magnitude_.  The $-2.5$ factor is also a legacy of our
     ancient origins: was necessary to approximately match the used
     magnitude system of Hipparchus.

          $$m-m_r=-2.5\log_{10} \left( F \over F_r \right)$$

     In the equation above, $m$ is the magnitude of the object and $m_r$
     is the pre-defined magnitude of the reference spectral flux
     density.  For estimating the error in measuring a magnitude, see
     *note Metameasurements on full input::.  The equation above is
     ultimately a relative relation.  To tie it to physical units,
     astronomers use the concept of a zero point which is discussed in
     the next item.

     See the ‘mag-to-luminosity’ operator of Arithmetic in *note Unit
     conversion operators:: for more on the conversion of the observed
     magnitudes (described below) of an object to luminosity.  The
     received brightness of two object with the same luminosity but at
     different distances is going to be different (the closer one will
     be brighter).  Therefore, astronomers have defined the following
     terminology to help avoid confusing distant-dependent and
     distance-independent magnitudes.

     Apparent magnitude
          The apparent magnitude is directly related (through the
          equation above) to the received spectral flux density (that we
          measure in our detectors).  Therefore, it depends on the
          distance to the object (and any absorption that may occur in
          the light path).
     Absolute magnitude
          Knowing the distance of an object and absorptions in the light
          path, we can obtain the luminosity.  From the luminosity, we
          can measure the apparent magnitude if the object was a point
          at a nominal (or fixed, or standard, or absolute) distance.
          The magnitude at this absolute distance is known as the
          absolute magnitude.  The standard (or abstract: just to help
          in comparisons) distance is historically defined as 10
          parsecs.  By reporting the magnitude at a fixed distance for
          all objects, the absolute magnitude therefore helps to compare
          the intrinsic (independent of distance) magnitude of
          astronomical objects.

Zero point
     A unique situation in the magnitude equation above occurs when the
     reference spectral flux density is unity ($F_r=1$).  In other
     words, the increase in spectral flux density that produces an
     increment in the detector's native measurement units (ADUs).

     The word "increment" above is used intentionally: because ADUs are
     discrete and measured as integers.  In other words, an increase in
     spectral flux density that is below $F_r$ will not be measured by
     the device.  The reference magnitde ($m_r$) that corresponds to
     $F_r$ is known as the _Zero point_ magnitude of that detector +
     filter + atmosphere (for on-ground observations).

     Therefore, the increase in spectral flux density (from an
     astrophysical source) that produces an increment in ADUs depends on
     all hardware and observational parameters that the image was taken
     in.  These include the quantum efficiency of the dector, the
     detector's coating, the filter transmission curve, the transmission
     of the optical path and the atmospheric absorption (for
     ground-based images; for example observations at different
     altitudes from the horizon where the thickness of the atmosphere is
     different).

     The rest of the absoptions (for example due to the interstellar
     medium, or ISM) are not considered in the zero point definition
     because for most purposes, they are not related to our observing
     conditions, but position on the sky.  In other words, while ISM
     absorption should be taken into account when measuring the
     luminosity of the source for example, ISM absorption is not in the
     zero point.  If we can later observe the universe from outside the
     MilkyWay, the ISM absorption should also be included (it would
     become like the atmosphere).  But the farthest we have got so far
     for scientific observations beyond the Solar system is the L2 orbit
     of Earth (for instruments like Euclid, Gaia or JWST).

     The zero point therefore allows us to summarize all these
     "observational" (non-astrophysical) factors into a single number
     and compare different observations from different instruments at
     different observing conditions (which are critical to do science).
     Defining the zero point magnitude as $m_r=Z$ in the magnitude
     equation, we can write it in simpler format (recall that $F_r=1$):

                     $$m = -2.5\log_{10}(F) + Z$$

     The zero point is found through comparison of measurements with
     pre-defined standards (in other words, it is a calibration of the
     pixel values).  Gnuastro has an installed script with a complete
     tutorial to estimate the zero point of any image, see *note Zero
     point estimation::.

     Historically, the reference was defined to be measurements of the
     star Vega, producing the _vega magnitude_ system.  In this system,
     the star Vega had a magnitude of zero (similar to the catalog of
     Hipparchus of Nicaea).  However, this caused many problems because
     Vega itself has its unique spectral features which are not in other
     stars and it is a variable star when measured precisely.

     Therefore, based on previous efforts, in 1983 Oke & Gunn proposed
     (https://ui.adsabs.harvard.edu/abs/1983ApJ...266..713O) the AB
     (absolute) magnitude system from accurate spectroscopy of Vega.  To
     avoid confusion with the "absolute magnitude" of a source (at a
     fixed distance), this magnitude system is always written as AB
     magnitude.  The AB magnitude zero point (when the input is
     frequency flux density; $F_\nu$ with units of $erg/s/cm^2/Hz$) was
     defined such that a star with a flat spectra around $5480\AA$ have
     a similar magnitude in the AB and Vega-based systems:

               $$m_{AB} = -2.5\log_{10}(F_\nu) + 48.60$$

     Reversing this equation and using Janskys, an object with a
     magnitude of zero ($m_{AB}=0$) has a spectral flux density of
     $3631Jy$.  Once the AB magnitude zero point of an image is found,
     you can directly convert any measurement on it from instrument ADUs
     to Janskys.  In Gnuastro, the Arithmetic program has an operator
     called ‘counts-to-jy’ which will do this though a given AB
     Magnitude-based zero point like below (SDSS data have a fixed zero
     point of 22.5 in the AB magnitude system):

          $ astarithmetic sdss.fits 22.5 counts-to-jy

     *Verify the zero point definition on new databases:* observational
     factors like the exposure time, the gain, telescope aperture,
     filter transmission curve and other factors are usually taken into
     account in the reduction pipeline that produces high-level science
     products.  But some reduction pipelines may not account for some of
     these for special reasons: for example not accounting for the gain
     or exposure time.  To avoid annoyingly strange results, when using
     a new database, verify (in the documentation of the database) that
     the zero points they provide directly converts pixel values to
     Janskys (is an AB magnitude zero point), or not.  If they not, you
     need to apply corrections your self.

     Let's look at one example where the given zero point has not
     accounted for the exposure time (in other words it is only for a
     fixed exposure time: $Z_E$), but the pixel values ($p$) have been
     corrected for the exposure time.  One solution would be to first
     multiply the pixels by the exposure time, use that zero point to
     get your desired measurement and delete the temporary file.  But a
     more optimal way (in terms of storage, execution and clean code)
     would be to correct the zero point.  Let's take $t$ to show time in
     units of seconds and $p_E$ to be the pixel value that would be
     measured after the fixed exposure time (in other words $p_E=p\times
     t$).  We then have the following:

   $$m = -2.5\log_{10}(p_E) + Z_E = -2.5\log_{10}(p\times t) + Z_E$$

     From the properties of logarithms, we can then derive the correct
     zero point ($Z$) to use directly (without touching the original
     pixels):

$$m = -2.5\log_{10}(p) + Z \quad\rm{where}\quad Z = Z_E - 2.5\log_{10}(t)$$

Surface brightness
     The definition of magnitude above was for the total spectral flux
     density coming from an object (recall how we mentioned at the start
     of this section that the total energy of an object is calculated by
     summing all its pixels).  The total flux is (mostly!)  independent
     of the angular size of your pixels, so we didn't need to account
     for the pixel area.  But when you want to study extended structures
     where the total magnitude is not desired (for example the
     sub-structure of a galaxy, or the brightness of the background
     sky), you need to report values that are independent of the area
     that total spectral flux density was measured on.

     For this, we define the _surface brightness_ to be the magnitude of
     an object's brightness divided by its solid angle over the
     celestial sphere (or coverage in the sky, commonly in units of
     arcsec$^2$).  The solid angle is expressed in units of arcsec$^2$
     because astronomical targets are usually much smaller than one
     steradian.  Recall that the steradian is the dimension-less SI unit
     of a solid angle and 1 steradian covers $1/4\pi$ (almost $8\%$) of
     the full celestial sphere.

     Surface brightness is therefore most commonly expressed in units of
     mag/arcsec$^2$.  For example, when the spectral flux density is
     measured over an area of A arcsec$^2$, the surface brightness is
     calculated by:

$$S = -2.5\log_{10}(F/A) + Z = -2.5\log_{10}(F) + 2.5\log_{10}(A) + Z$$

     In other words, the surface brightness (in units of mag/arcsec$^2$)
     is related to the object's magnitude ($m$) and area ($A$, in units
     of arcsec$^2$) through this equation:

                      $$S = m + 2.5\log_{10}(A)$$

     A common mistake is to follow the mag/arcsec$^2$ unit literally,
     and divide the object's magnitude by its area.  But this is wrong
     because magnitude is a logarithmic scale while area is linear.  It
     is the spectral flux density that should be divided by the solid
     angle because both have linear scales.  The magnitude of that ratio
     is then defined to be the surface brightness.

     Besides applications in catalogs and the resulting scientific
     analysis, converting pixels to surface brightness is usually a good
     way to display a FITS file in a publication!  See *note FITS images
     in a publication:: for a fully working tutorial on how to do this.

*Do not warp or convolve magnitude or surface brightness images:*
Warping an image involves calculating new pixel values (of the new pixel
grid) from the input grid's pixel values.  Convolution is also a process
of finding the weighted mean of pixel values.  During these processes,
many arithmetic operations are done on the original pixel values, for
example, addition or multiplication.  However, $log_{10}(a+b)\ne
log_{10}(a)+log_{10}(b)$.  Therefore if you generate color, magnitude or
surface brightness images (where pixels are in units of mangnitudes), do
not apply any such operations on them!  If you need to warp or convolve
the image, do it _before_ the conversion to magnitude-based units.


File: gnuastro.info,  Node: Standard deviation vs Standard error,  Next: MakeCatalog measurements on each label,  Prev: Brightness flux magnitude,  Up: MakeCatalog

7.4.3 Standard deviation vs Standard error
------------------------------------------

The standard deviation and standard error are different concepts to
convey different aspects of a measurement, but they can be easily
confused with each other: for example, the standard deviation is also
called as the "error".  A nice description of this difference is given
in the following quote from Wikipedia.  In this section, we'll show the
concept described in this quote with a hands-on and practical set of
commands to clarify this important distinction.

     The standard deviation of the sample data is a description of the
     variation in measurements, while the standard error of the mean is
     a probabilistic statement about how the sample size will provide a
     better bound on estimates of the population mean, in light of the
     central limit theorem.  Put simply, the _standard error_ of the
     sample mean is an estimate of how far the sample mean is likely to
     be from the population mean, whereas the _standard deviation_ of
     the sample is the degree to which individuals within the sample
     differ from the sample mean.
   -- _Wikipedia page on "Standard error" (retrieved on 2024-08-09)_

   Let's simulate an observation of the sky, but without any
astronomical sources to simplify the logic.  In other words, we only
have a background flux level (assume you want to measure the brightness
of the twilight sky that is not yet faint enough to let stars be
visible).  With the first command below, let's make an image called
‘1.fits’ that contains $200\times200$ pixels that are filled with random
noise from a Poisson distribution with a mean of 100 counts (the flux
from the background sky).  With the second command, we'll have a look at
the image.  Recall that the Poisson distribution is equal to a normal
distribution for large mean values (as in this case).

     $ astarithmetic 200 200 2 makenew 100 mknoise-poisson \
                     --output=1.fits

     $ astscript-fits-view 1.fits

   The standard deviation ($\sigma$) of the Poisson distribution is the
square root of the mean, see *note Photon counting noise::.  Note that
due to the random nature of the noise, the values reported in the next
steps on your computer will be slightly different (which is one of the
points of this section!).  To reproducible exactly the same values in
different runs, see *note Generating random numbers::, and for more on
the first command, see *note Arithmetic::.

   Each pixel shows the result of one sampling from the Poisson
distribution.  In other words, assuming the sky emission in our
simulation is constant over our field of view, each pixel's value shows
one measurement of the same sky emission.  Statistically speaking, a
"measurement" is a sampling from an underlying distribution of values.
Through our measurements, we aim to identify that underlying
distribution (the "truth")!  With the command below, let's look at the
pixel statistics of ‘1.fits’ (output is shown immediately under it).

     $ aststatistics 1.fits
     Statistics (GNU Astronomy Utilities) 0.23.119-586f
     -------
     Input: 1.fits (hdu: 1)
     -------
       Number of elements:                      40000
       Minimum:                                 61
       Maximum:                                 155
       Median:                                  100
       Mean:                                    100.044925
       Standard deviation:                      10.00066032
     -------
     Histogram:
      |                          *  *
      |                          *  *  *
      |                       *  *  *  *
      |                       *  *  *  *  *
      |                       *  *  *  *  *
      |                    *  * ********  * *
      |                    *  ************* *
      |                 *  ******************  *
      |                 ************************  *
      |              *********************************
      |* ********************************************************** **      *
      |----------------------------------------------------------------------

   As expected, you see that the ASCII histogram nicely resembles a
normal distribution.  The measured mean and standard deviation are also
very similar to the input (mean of 100, standard deviation of 10).  But
the measured mean (and standard deviation) aren't exactly equal to the
input!

   Every time we make a different simulated image from the same
distribution, the measured mean and standard deviation will slightly
differ.  Run the commands above one more time (this time calling the
output ‘2.fits’) and check for your self (actually open the two FITS
images and check visually, don't just rely on the statistics).

   Now that you have a good feeling of the change, let's automate this
and scale it up for some nice statistics.  With the commands below, you
will build 500 images like above and measure their mean and standard
deviation and save each measurment into a file (‘mean-stds.txt’.  In the
first command we are deleting it to make sure we write into an empty
file within the first run of the loop.  With the third command, let's
view the top 10 rows:

     $ rm -f mean-stds.txt
     $ for i in $(seq 500); do \
           astarithmetic 200 200 2 makenew 100 mknoise-poisson \
                         --output=$i.fits --quiet; \
           aststatistics $i.fits --mean --std >> mean-stds.txt; \
           echo "$i: complete"; \
       done

     $ asttable mean-stds.txt -Y --head=10
     99.989381               9.936407
     100.036622              10.059997
     100.006054              9.985470
     99.944535               9.960069
     100.050318              9.970116
     100.002718              9.905395
     100.067555              9.964038
     100.027167              10.018562
     100.051951              9.995859
     100.000212              9.970293

   From this table, you see that each simulation has produced a slightly
different measured mean and measured standard deviation.  They are just
fluctuating around the input mean (which was 100) and input standard
deviation (which was 10).  Let's have a look at the distribution of mean
measurements:

     $ aststatistics mean-stds.txt -c1
     Statistics (GNU Astronomy Utilities) 0.23.119-586f
     -------
     Input: mean-stds.txt
     Column: 1
     -------
       Number of elements:                      500
       Minimum:                                 9.98183528700191e+01
       Maximum:                                 1.00146490891332e+02
       Mode:                                    99.99709739
       Mode quantile:                           0.49498998
       Median:                                  9.99977393190436e+01
       Mean:                                    99.99891826
       Standard deviation:                      0.04901635275
     -------
     Histogram:
      |                                       *
      |                                   *   **
      |                               ****** **** * *
      |                               ****** **** * *    *
      |                          *  * ************* *    *
      |                          *  ******************   **
      |                   *      *********************  ***   *
      |                   *   ***************************** ***
      |                   *** **********************************      *
      |             ***  *******************************************  **
      |           * ************************************************* **    *
      |----------------------------------------------------------------------

   The standard deviation you see above (approximately $0.05$) shows the
scatter in measuring the mean with an image of this size and is
different from the standard deviation of the Poisson distribution that
the values were drawn from.  This is therefore defined as the _standard
error of the mean_, or "standard error" for short (since most
measurements are actually the mean of a population).

   From the example above, you see that the error is smaller than the
standard deviation (smaller when you have a larger sample).  In fact, it
can be shown (https://en.wikipedia.org/wiki/Standard_error#Derivation)
that this "error of the mean" ($\sigma_{\bar{x}}$; recall that $\bar{x}$
represents the mean) is related to the distribution standard deviation
($\sigma$) through the following equation.  Where $N$ is the number of
points used to measure the mean in one sample ($N=200\times200=40000$ in
this case).  Note that the $10.0$ below was reported as "standard
deviation" in the first run of ‘aststatistics’ on ‘1.fits’ above):

$$\sigma_{\bar{x}}\approx\frac{\sigma}{\sqrt{N}} = \frac{10.0}{200} = 0.05$$

   Therefore the standard error of the mean is directly related to the
number of pixels you used to measure the mean.  You can test this by
changing the ‘200’s in the commands above to smaller or larger values.
As you make larger and larger images, you will be able to measure the
mean much more precisely (the standard error of the mean will go to
zero).  But no matter how many pixels you use, the standard deviation
will always be the same.

   Within MakeCatalog, the options related to dispersion/error in the
measurements have the following conventions:
‘--*std’
     For example ‘--std’ or ‘--sigclip-std’.  These return the standard
     deviation of the values within a label.  If the underlying object
     (in the noise) is flat, then this will be the ‘\sigma’ that is
     mentioned above.  However, no object in astronomy is flat!  So this
     option should be used with extreme care!  It only makes sense in
     special contexts like measuring the radial profile where we assume
     that the values at a certain radius have the same flux (see *note
     Generate radial profile::).

‘--*error’
     For example ‘--mag-error’, ‘--mean-error’ or ‘--sum-error’.  These
     options should be used when when pixel values are different.  When
     the pixels do not have the same value (for example different parts
     of one galaxy), their standard deviation is meaningless.  To
     measure the total error in such cases, we need to know the standard
     deviation of each pixel separately.  Therefore, for these columns,
     MakeCatalog needs a separate dataset that contains the underlying
     sky standard deviation for those pixels.  That dataset should have
     the same size (number and dimension of pixels) as the values
     dataset.  You can use *note NoiseChisel:: to generate such an
     image.

     If the underlying profile and sky standard deviations is flat, then
     ‘--sum-error’ will be the standard deviation that we discussed in
     this section and ‘--mean-error’ will be the standard error.  When
     the values are different, the combined error is calculated by
     adding the variances (second power of the standard deviation) of
     each pixel, added with its value.  When the values are smaller than
     one, a correction is applied (that is defined in Section 3.3 of
     Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664)).


File: gnuastro.info,  Node: MakeCatalog measurements on each label,  Next: Metameasurements on full input,  Prev: Standard deviation vs Standard error,  Up: MakeCatalog

7.4.4 MakeCatalog measurements on each label
--------------------------------------------

MakeCatalog's output measurements can be classified into two types:
columns (independent measurements on each label) or metadata (of the
whole image, see *note Metameasurements on full input::).  The precise
measurements that it will actually do is specified using command-line
options that are described in the subsections below (*note Options::).

   The majority of the first group are those which only produce a
single-valued measurement for each label (for example, its magnitude: a
single number for every label/object).  There are also multi-valued
measurements per label (like spectra on 3D cubes) which keep their
values in a single vector column of the same table (see *note Vector
columns::).  Both types of measurements will be written as one column in
a final table/catalog that also contains all the other requested
columns/measurements.  Measurements that produce one value for the whole
input (not a specific label/object) are stored as metadata in the 0-th
(first) HDU of the output with predefined keywords described in *note
Metameasurements on full input::.

   The majority of this section is devoted to MakeCatalog's
single-valued measurements.  However, MakeCatalog can also do
measurements that produce more than one value for each label.  Currently
the only such measurement is generation of spectra from 3D cubes with
the ‘--spectrum’ option and it is discussed in the end of this section.

   Command-line options are used to identify which measurements you want
in the final catalog(s) and in what order.  If any of the options below
is called on the command-line or in any of the configuration files, it
will be included as a column in the output catalog.  The order of the
columns is in the same order as the options were seen by MakeCatalog
(see *note Configuration file precedence::).  Some of the columns apply
to both "objects" and "clumps" and some are particular to only one of
them (for the definition of "objects" and "clumps", see *note
Segment::).  Columns/options that are unique to one catalog (only
objects, or only clumps), are explicitly marked with [Objects] or
[Clumps] to specify the catalog they will be placed in.

* Menu:

* Identifier columns::          Identifying labels of each row (object/clumps).
* Position measurements in pixels::  Containing image/pixel (X/Y) measurements.
* Position measurements in WCS::  Containing WCS (for example RA/Dec) measurements.
* Brightness measurements::     Using pixel values of each label.
* Surface brightness measurements::  Various ways to measure surface brightness.
* Upper limit measurements::
* Morphology measurements nonparametric::  Non-parametric morphology.
* Morphology measurements elliptical::  Elliptical morphology measurements.
* Measurements per slice spectra::  Measurements on each slice (like spectra).


File: gnuastro.info,  Node: Identifier columns,  Next: Position measurements in pixels,  Prev: MakeCatalog measurements on each label,  Up: MakeCatalog measurements on each label

7.4.4.1 Identifier columns
..........................

The identifier of each row (group of measurements) is usually the first
thing you will be requesting from MakeCatalog.  Without the identifier,
it is not clear which measurement corresponds to which label for the
input.

   Since MakeCatalog can also optionally take sub-structure label
(clumps; see *note Segment::), there are various identifiers in general
that are listed below.  The most generic (and shortest and easiest to
type!)  is the ‘--ids’ option which can be used in object-only or
object-clump catalogs.

‘--i’
‘--ids’
     This is a unique option which can add multiple columns to the final
     catalog(s).  Calling this option will put the object IDs
     (‘--obj-id’) in the objects catalog and host-object-ID
     (‘--host-obj-id’) and ID-in-host-object (‘--id-in-host-obj’) into
     the clumps catalog.  Hence if only object catalogs are required, it
     has the same effect as ‘--obj-id’.

‘--obj-id’
     [Objects] ID of this object.

‘-j’
‘--host-obj-id’
     [Clumps] The ID of the object which hosts this clump.

‘--id-in-host-obj’
     [Clumps] The ID of this clump in its host object.


File: gnuastro.info,  Node: Position measurements in pixels,  Next: Position measurements in WCS,  Prev: Identifier columns,  Up: MakeCatalog measurements on each label

7.4.4.2 Position measurements in pixels
.......................................

The position of a labeled region within your input dataset (in its own
units) can be measured with the options in this section.  By "in its own
units" we mean pixels in a 2D image or voxels in a 3D cube.  For example
if the flux-weighted center of a label lies 123 pixels on the horizontal
and 456 pixels on the vertical, the ‘--x’ and ‘--y’ options will put a
value of 123 and 456 in their respective columns.  As you see below,
there are various ways to define the "position" of an object, so read
the differences carefully to choose the one that corresponds best to
your usage.

‘-x’
‘--x’
     The flux weighted center of all objects and clumps along the first
     FITS axis (horizontal when viewed in SAO DS9), see $\overline{x}$
     in *note Morphology measurements elliptical::.  The weight has to
     have a positive value (pixel value larger than the Sky value) to be
     meaningful!  Specially when doing matched photometry, this might
     not happen: no pixel value might be above the Sky value.  For such
     detections, the geometric center will be reported in this column
     (see ‘--geo-x’).  You can use ‘--weight-area’ to see which was
     used.

‘-y’
‘--y’
     The flux weighted center of all objects and clumps along the second
     FITS axis (vertical when viewed in SAO DS9).  See ‘--x’.

‘-z’
‘--z’
     The flux weighted center of all objects and clumps along the third
     FITS axis.  See ‘--x’.

‘--geo-x’
     The geometric center of all objects and clumps along the first FITS
     axis axis.  The geometric center is the average pixel positions
     irrespective of their pixel values.

‘--geo-y’
     The geometric center of all objects and clumps along the second
     FITS axis axis, see ‘--geo-x’.

‘--geo-z’
     The geometric center of all objects and clumps along the third FITS
     axis axis, see ‘--geo-x’.

‘--min-val-x’
     Position of pixel with minimum value in objects and clumps, along
     the first FITS axis.

‘--max-val-x’
     Position of pixel with maximum value in objects and clumps, along
     the first FITS axis.

‘--min-val-y’
     Position of pixel with minimum value in objects and clumps, along
     the first FITS axis.

‘--max-val-y’
     Position of pixel with maximum value in objects and clumps, along
     the first FITS axis.

‘--min-val-z’
     Position of pixel with minimum value in objects and clumps, along
     the first FITS axis.

‘--max-val-z’
     Position of pixel with maximum value in objects and clumps, along
     the first FITS axis.

‘--min-x’
     The minimum position of all objects and clumps along the first FITS
     axis.

‘--max-x’
     The maximum position of all objects and clumps along the first FITS
     axis.

‘--min-y’
     The minimum position of all objects and clumps along the second
     FITS axis.

‘--max-y’
     The maximum position of all objects and clumps along the second
     FITS axis.

‘--min-z’
     The minimum position of all objects and clumps along the third FITS
     axis.

‘--max-z’
     The maximum position of all objects and clumps along the third FITS
     axis.

‘--clumps-x’
     [Objects] The flux weighted center of all the clumps in this object
     along the first FITS axis.  See ‘--x’.

‘--clumps-y’
     [Objects] The flux weighted center of all the clumps in this object
     along the second FITS axis.  See ‘--x’.

‘--clumps-z’
     [Objects] The flux weighted center of all the clumps in this object
     along the third FITS axis.  See ‘--x’.

‘--clumps-geo-x’
     [Objects] The geometric center of all the clumps in this object
     along the first FITS axis.  See ‘--geo-x’.

‘--clumps-geo-y’
     [Objects] The geometric center of all the clumps in this object
     along the second FITS axis.  See ‘--geo-x’.

‘--clumps-geo-z’
     [Objects] The geometric center of all the clumps in this object
     along the third FITS axis.  See ‘--geo-z’.


File: gnuastro.info,  Node: Position measurements in WCS,  Next: Brightness measurements,  Prev: Position measurements in pixels,  Up: MakeCatalog measurements on each label

7.4.4.3 Position measurements in WCS
....................................

The position of a labeled region within your input dataset (in the World
Coordinate System, or WCS) can be measured with the options in this
section.  As you see below, there are various ways to define the
"position" of an object, so read the differences carefully to choose the
one that corresponds best to your usage.

   The most common WCS coordinates are Right Ascension (RA) and
Declination in an equatorial system.  Therefore, to simplify their
usage, we have special ‘--ra’ and ‘--dec’ options.  However, the WCS of
datasets are in Galactic coordinates, so to be generic, you can use the
‘--w1’, ‘--w2’ or ‘--w3’ (if you have a 3D cube) options.  In case your
dataset's WCS is not in your desired system (for example it is Galactic,
but you want equatorial 2000), you can use the ‘--wcscoordsys’ option of
Gnuastro's Fits program on the labeled image before running MakeCatalog
(see *note Keyword inspection and manipulation::).

‘-r’
‘--ra’
     Flux weighted right ascension of all objects or clumps, see ‘--x’.
     This is just an alias for one of the lower-level ‘--w1’ or ‘--w2’
     options.  Using the FITS WCS keywords (‘CTYPE’), MakeCatalog will
     determine which axis corresponds to the right ascension.  If no
     ‘CTYPE’ keywords start with ‘RA’, an error will be printed when
     requesting this column and MakeCatalog will abort.

‘-d’
‘--dec’
     Flux weighted declination of all objects or clumps, see ‘--x’.
     This is just an alias for one of the lower-level ‘--w1’ or ‘--w2’
     options.  Using the FITS WCS keywords (‘CTYPE’), MakeCatalog will
     determine which axis corresponds to the declination.  If no ‘CTYPE’
     keywords start with ‘DEC’, an error will be printed when requesting
     this column and MakeCatalog will abort.

‘--w1’
     Flux weighted first WCS axis of all objects or clumps, see ‘--x’.
     The first WCS axis is commonly used as right ascension in images.

‘--w2’
     Flux weighted second WCS axis of all objects or clumps, see ‘--x’.
     The second WCS axis is commonly used as declination in images.

‘--w3’
     Flux weighted third WCS axis of all objects or clumps, see ‘--x’.
     The third WCS axis is commonly used as wavelength in integral field
     unit data cubes.

‘--geo-w1’
     Geometric center in first WCS axis of all objects or clumps, see
     ‘--geo-x’.  The first WCS axis is commonly used as right ascension
     in images.

‘--geo-w2’
     Geometric center in second WCS axis of all objects or clumps, see
     ‘--geo-x’.  The second WCS axis is commonly used as declination in
     images.

‘--geo-w3’
     Geometric center in third WCS axis of all objects or clumps, see
     ‘--geo-x’.  The third WCS axis is commonly used as wavelength in
     integral field unit data cubes.

‘--clumps-w1’
     [Objects] Flux weighted center in first WCS axis of all clumps in
     this object, see ‘--x’.  The first WCS axis is commonly used as
     right ascension in images.

‘--clumps-w2’
     [Objects] Flux weighted declination of all clumps in this object,
     see ‘--x’.  The second WCS axis is commonly used as declination in
     images.

‘--clumps-w3’
     [Objects] Flux weighted center in third WCS axis of all clumps in
     this object, see ‘--x’.  The third WCS axis is commonly used as
     wavelength in integral field unit data cubes.

‘--clumps-geo-w1’
     [Objects] Geometric center right ascension of all clumps in this
     object, see ‘--geo-x’.  The first WCS axis is commonly used as
     right ascension in images.

‘--clumps-geo-w2’
     [Objects] Geometric center declination of all clumps in this
     object, see ‘--geo-x’.  The second WCS axis is commonly used as
     declination in images.

‘--clumps-geo-w3’
     [Objects] Geometric center in third WCS axis of all clumps in this
     object, see ‘--geo-x’.  The third WCS axis is commonly used as
     wavelength in integral field unit data cubes.


File: gnuastro.info,  Node: Brightness measurements,  Next: Surface brightness measurements,  Prev: Position measurements in WCS,  Up: MakeCatalog measurements on each label

7.4.4.4 Brightness measurements
...............................

Within an image, pixels have both a position and a value.  In the
sections above all the measurements involved position (see *note
Position measurements in pixels:: or *note Position measurements in
WCS::).  The measurements in this section only deal with pixel values
and ignore the pixel positions completely.  In other words, for the
options of this section each labeled region within the input is just a
group of values (and their associated error values if necessary), and
they let you do various types of measurements on the resulting
distribution of values.  For more on the difference between the
‘--*error’ or ‘--*std’ columns see *note Standard deviation vs Standard
error::.

‘--sum’
     The sum of all pixel values associated to this label (object or
     clump).  Note that if a sky value or image has been given, it will
     be subtracted before any column measurement.  For clumps, the
     ambient values (average of river pixels around the clump,
     multiplied by the area of the clump) is subtracted, see
     ‘--river-mean’.  So the sum of all the clump-sums in the clump
     catalog of one object will be smaller than the ‘--clumps-sum’
     column of the objects catalog.

     If no usable pixels are present over the clump or object (for
     example, they are all blank), the returned value will be NaN (note
     that zero is meaningful).

‘--sum-error’
     The standard deviation of the sum of values of a label (objects or
     clumps).  The value is calculated by using the values image (for
     signal above the sky level) and the sky standard deviation image
     (extension ‘--stdhdu’ of file given to ‘--instd’); which you can
     derive for any image using *note NoiseChisel::.  This column is
     internally used to measure the signal-to-noise (‘--sn’).

     For objects this is calculated by adding the sky variance (second
     power of the sky standard deviation image) of each pixel in the
     label, with the value of the pixel if the value is not negative
     (error only increases).  This is done to account for brighter
     pixels which have higher noise in the Poisson distribution (its
     side effect is that the error will always be slightly
     over-estimated due to the positive values close to the noise).  A
     correction may be applied if the sky standard deviation is
     negative; see Section 3.3 of Akhlaghi & Ichikawa 2015
     (https://arxiv.org/abs/1505.01664).  For clumps, the variance of
     the rivers (which are subtracted from the value of pixels in
     calclulating the sum) are also added to generate the final standard
     deviation.

     The returned value will be NaN when the label covers only NaN
     pixels in the values image, or a pixel is NaN in the ‘--instd’
     image, but non-NaN in the values image.  The latter situation
     usually happens when there is a bug in the previous steps of your
     analysis.  This is because the sky standard deviation should have a
     value in all pixels.  In such cases, it is important to find the
     cause and fix it because those pixels with a NaN in the ‘--instd’
     image may contribute significantly to the final error.  If you want
     to ignore those pixels in the error measurement, set them to zero
     (which is a meaningful number in such scenarios).

‘--clumps-sum’
     [Objects] The total sum of the pixels covered by clumps (before
     subtracting the river) within each object.  This is simply the sum
     of ‘--sum-no-river’ in the clumps catalog (see below).  If no
     usable pixels are present over the clump or object (for example,
     they are all blank), the stored value will be NaN (note that zero
     is meaningful).

‘--sum-no-river’
     [Clumps] The sum of Sky (not river) subtracted clump pixel values.
     By definition, for the clumps, the average value of the rivers
     surrounding it are subtracted from it for a first order accounting
     for contamination by neighbors.

     If no usable pixels are present over the clump or object (for
     example, they are all blank), the stored value will be NaN (note
     that zero is meaningful).

‘--mean’
     The mean sky subtracted value of pixels within the object or clump.
     For clumps, the average river flux is subtracted from the sky
     subtracted mean.

‘--mean-error’
     The error in measuring the mean; using both the values file and the
     sky standard deviation image.  In case the given standard deviation
     or variance image already contains the contributions from the pixel
     values (it is not just the sky standard deviation), use
     ‘--novalinerror’).

‘--std’
     The standard deviation of the pixels within the object or clump.
     For clumps, the river pixels are not subtracted because that is a
     constant (per pixel) value and should not affect the standard
     deviation.

‘--median’
     The median sky subtracted value of pixels within the object or
     clump.  For clumps, the average river flux is subtracted from the
     sky subtracted median.

‘--maximum’
     The maximum value of pixels within the object or clump.  When the
     label (object or clump) is larger than three pixels, the maximum is
     actually derived by the mean of the brightest three pixels, not the
     largest pixel value of the same label.  This is because noise
     fluctuations can be very strong in the extreme values of the
     objects/clumps due to Poisson noise (which gets stronger as the
     mean gets higher).  Simply using the maximum pixel value will
     create a strong scatter in results that depend on the maximum (for
     example, the ‘--fwhm’ option also uses this value internally).

‘--sigclip-number’
     The number of elements/pixels in the dataset after sigma-clipping
     the object or clump.  The sigma-clipping parameters can be set with
     the ‘--sigmaclip’ option described in *note MakeCatalog inputs and
     basic settings::.  For more on Sigma-clipping, see *note Sigma
     clipping::.

‘--sigclip-median’
     The sigma-clipped median value of the object of clump's pixel
     distribution.  For more on sigma-clipping and how to define it, see
     ‘--sigclip-number’.

‘--sigclip-mean’
     The sigma-clipped mean value of the object of clump's pixel
     distribution.  For more on sigma-clipping and how to define it, see
     ‘--sigclip-number’.

‘--sigclip-std’
     The sigma-clipped standard deviation of the object of clump's pixel
     distribution.  For more on sigma-clipping and how to define it, see
     ‘--sigclip-number’.

‘-m’
‘--magnitude’
     The magnitude of clumps or objects.  It is derived through the
     pixel counts over the label (which you can see in the ‘--sum’
     column) and the value given to the ‘--zeropoint’ through the
     definition of the magnitude described in *note Brightness flux
     magnitude::.

‘--magnitude-error’
     The magnitude error of clumps or objects.  The method used is
     described below.  As we see there, this error assumes uncorrelated
     pixel values and also does not include the error in estimating the
     aperture (or error in generating the labeled image).  See the
     status of implementation of such factors in Task 14124
     (https://savannah.gnu.org/task/index.php?14124).

     The returned value will be NaN when the label covers only NaN
     pixels in the values image, or a pixel is NaN in the ‘--instd’
     image, but non-NaN in the values image.  The latter situation
     usually happens when there is a bug in the previous steps of your
     analysis, and is important because those pixels with a NaN in the
     ‘--instd’ image may contribute significantly to the final error.
     If you want to ignore those pixels in the error measurement, set
     them to zero (which is a meaningful number in such scenarios).

     The raw error in measuring the magnitude is only meaningful when
     the object's magnitude is brighter than the upper-limit magnitude
     (see below).  As discussed in *note Brightness flux magnitude::,
     the magnitude ($M$) of an object with brightness $B$ and zero point
     magnitude $z$ can be written as:

                       $$M=-2.5\log_{10}(B)+z$$

     Calculating the derivative with respect to $B$, we get:

            $${dM\over dB} = {-2.5\over {B\times ln(10)}}$$

     From the Tailor series ($\Delta{M}=dM/dB\times\Delta{B}$), we can
     write:

$$\Delta{M} = \left|{-2.5\over ln(10)}\right|\times{\Delta{B}\over{B}}$$

     But, $\Delta{B}/B$ is just the inverse of the Signal-to-noise ratio
     ($S/N$), so we can write the error in magnitude in terms of the
     signal-to-noise ratio:

             $$\Delta{M} = {2.5\over{S/N\times ln(10)}} $$

     MakeCatalog uses this relation to estimate the magnitude errors.
     The signal-to-noise ratio is calculated in different ways for
     clumps and objects, see Akhlaghi and Ichikawa 2015
     (https://arxiv.org/abs/1505.01664)), but this single equation can
     be used to estimate the measured magnitude error afterwards for any
     type of target.

‘--clumps-magnitude’
     [Objects] The magnitude of all clumps in this object, see
     ‘--clumps-sum’.

‘--river-mean’
     [Clumps] The average of the river pixel values around this clump.
     River pixels were defined in Akhlaghi and Ichikawa 2015
     (https://arxiv.org/abs/1505.01664).  In short they are the pixels
     immediately outside of the clumps.  This value is used internally
     to find the sum (or magnitude) and signal to noise ratio of the
     clumps.  It can generally also be used as a scale to gauge the base
     (ambient) flux surrounding the clump.  In case there was no river
     pixels, then this column will have the value of the Sky under the
     clump.  So note that this value is _not_ sky subtracted.

‘--river-num’
     [Clumps] The number of river pixels around this clump, see
     ‘--river-mean’.

‘--river-min’
     [Clumps] Minimum river value around this clump, see ‘--river-mean’.

‘--river-max’
     [Clumps] Maximum river value around this clump, see ‘--river-mean’.

‘--sn’
     The Signal to noise ratio (S/N) of all clumps or objects.  See
     Akhlaghi and Ichikawa (2015) for the exact equations used.

     The returned value will be NaN when the label covers only NaN
     pixels in the values image, or a pixel is NaN in the ‘--instd’
     image, but non-NaN in the values image.  The latter situation
     usually happens when there is a bug in the previous steps of your
     analysis, and is important because those pixels with a NaN in the
     ‘--instd’ image may contribute significantly to the final error.
     If you want to ignore those pixels in the error measurement, set
     them to zero (which is a meaningful number).

‘--sky’
     The sky flux (per pixel) value under this object or clump.  This is
     actually the mean value of all the pixels in the sky image that lie
     on the same position as the object or clump.

‘--sky-std’
     The sky value standard deviation (per pixel) for this clump or
     object.  This is the square root of the mean variance under the
     object, or the root mean square.


File: gnuastro.info,  Node: Surface brightness measurements,  Next: Upper limit measurements,  Prev: Brightness measurements,  Up: MakeCatalog measurements on each label

7.4.4.5 Surface brightness measurements
.......................................

In astronomy, Surface brightness is most commonly measured in units of
magnitudes per arcsec$^2$ (for the formal definition, see *note
Brightness flux magnitude::).  Therefore it involves both the values of
the pixels within each input label (or output row) and their position.

‘--sb’
     The surface brightness (in units of mag/arcsec$^2$) of the labeled
     region (objects or clumps).  For more on the definition of the
     surface brightness, see *note Brightness flux magnitude::.

‘--sb-error’
     Error in measuring the surface brightness (the ‘--sb’ column) as
     derived below.

     We can derive the error in measuring the surface brightness based
     on the surface brightness (SB) equation of *note Brightness flux
     magnitude:: and the generic magnitude error ($\Delta{M}$) that is
     described under ‘--magnitude-error’ in *note Brightness
     measurements::.  Let's set $A$ to represent the area and
     $\Delta{A}$ to represent the error in measuring the area.  For more
     on $\Delta{A}$, see the description of ‘--spatialresolution’ in
     *note MakeCatalog inputs and basic settings::.

$$\Delta{(SB)} = \Delta{M} + \left|{-2.5\over ln(10)}\right|\times{\Delta{A}\over{A}}$$

     In the surface brightness equation mentioned above, $A$ is in units
     of arcsecond squared and the conversion between arcseconds to
     pixels is a multiplication factor.  Therefore as long as $A$ and
     $\Delta{A}$ have the same units, it does not matter if they are in
     arcseconds or pixels.  Since the measure of spatial resolution (or
     area error) is the FWHM of the PSF which is usually defined in
     terms of pixels, its more intuitive to use pixels for $A$ and
     $\Delta{A}$.

‘--upperlimit-sb’
     The upper-limit surface brightness (in units of mag/arcsec$^2$) of
     this labeled region (object or clump).  In other words, this option
     measures the surface brightness of noise within the footprint of
     each input label.

     This is just a simple wrapper over lower-level columns: setting B
     and A as the value in the columns ‘--upperlimit’ and
     ‘--area-arcsec2’, we fill this column by simply use the surface
     brightness equation of *note Brightness flux magnitude::.

‘--half-sum-sb’
     Surface brightness (in units of mag/arcsec$^2$) within the area
     that contains half the total sum of the label's pixels (object or
     clump).  This is useful as a measure of the sharpness of an
     astronomical object: for example a star will have very few pixels
     at half the maximum, so its ‘--half-sum-sb’ will be much brighter
     than a galaxy at the same magnitude.  Also consider ‘--half-max-sb’
     below.

     This column just plugs in the values of half the value of the
     ‘--sum’ column and the ‘--half-sum-area’ column, into the surface
     brightness equation.  Therefore please see the description in
     ‘--half-sum-area’ to understand the systematics of this column and
     potential biases (see *note Morphology measurements
     nonparametric::).

‘--half-max-sb’
     The surface brightness (in units of mag/arcsec$^2$) within the
     region that contains half the maximum value of the labeled region.
     Like ‘--half-sum-sb’ this option this is a good way to identify the
     "central" surface brightness of an object.  To know when this
     measurement is reasonable, see the description of ‘--fwhm’ in *note
     Morphology measurements nonparametric::.

‘--sigclip-mean-sb’
     Surface brightness (over 1 pixel's area in arcsec$^2$) of the
     sigma-clipped mean value of the pixel values distribution
     associated to each label (object or clump).  This is useful in
     scenarios where your labels have approximately _constant_ surface
     brightness values _after_ after removing outliers: for example in a
     radial profile, see *note Invoking astscript-radial-profile::).

     In other scenarios it should be used with extreme care.  For
     example over the full area of a galaxy/star the pixel distribution
     is not constant (or symmetric after adding noise), their pixel
     distributions are inherently skewed (with fewer pixels in the
     center, having a very large value and many pixels in the outer
     parts having lower values).  Therefore, sigma-clipping is not
     meaningful for such objects!  For more on the definition of the
     surface brightness, see *note Brightness flux magnitude::, for more
     on sigma-clipping, see *note Sigma clipping::.

     The error in this magnitude can be retrieved from the
     ‘--sigclip-mean-sb-delta’ column described below, and you can use
     the ‘--sigclip-std-sb’ column to find when the magnitude has become
     noise-dominated (signal-to-noise ratio is roughly 1).  See the
     description of these two options for more.

‘--sigclip-mean-sb-delta’
     Scatter in the ‘--sigclip-mean-sb’ without using the standard
     deviation of each pixel (that is given by ‘--instd’ in *note
     MakeCatalog inputs and basic settings::).  The scatter here is
     measured from the values of the label themselves.  This measurement
     is therefore most meaningful when you expect the flux across one
     label to be constant (as in a radial profile for example).

     This is calculated using surface brightness error derived under
     ‘--sb-error’ in *note Surface brightness measurements::, but with
     $\Delta{A}=0$ (since sigma-clip is calculated per pixel and there
     is no error in a single pixel).  Within the equation to derive
     $\Delta{M}$ (the error in magnitude, derived in the description of
     ‘--magnitude-error’ in *note Brightness measurements::), the
     signal-to-noise ratio is defined by dividing the sigma-clipped mean
     by the sigma-clipped standard deviation.

‘--sigclip-std-sb’
     The surface brightness of the sigma-clipped standard deviation of
     all the pixels with the same label.  For labels that are expected
     to have the same value in all their pixels (for example each
     annulus of a radial profile) this can be used to find the reliable
     ($1\sigma$) surface brightness for that label.  In other words, if
     ‘--sigclip-mean-sb’ is fainter than the value of this column, you
     know that noise is becoming significant.  However, as described in
     ‘--sigclip-mean-sb’, the sigma-clipped measurements of MakeCatalog
     should only be used in certain situations like radial profiles, see
     the description there for more.


File: gnuastro.info,  Node: Upper limit measurements,  Next: Morphology measurements nonparametric,  Prev: Surface brightness measurements,  Up: MakeCatalog measurements on each label

7.4.4.6 Upper limit measurements
................................

Due to the noisy nature of data, it is possible to get arbitrarily faint
magnitudes, especially when you use labels from another image (for
example see *note Working with catalogs estimating colors::).  Given the
scatter caused by the dataset's noise, values fainter than a certain
level are meaningless: another similar depth observation will give a
radically different value.  In such cases, measurements like the image
magnitude limit are not useful because it is estimated for a certain
morphology and is given for the whole image (it is a crude
generalization; see *note Metameasurements on full input::.  You want a
quality measure that is specific to each object.

   For example, assume that you have done your detection and
segmentation on one filter and now you do measurements over the same
labeled regions, but on other filters to measure colors (as we did in
the tutorial *note Segmentation and making a catalog::).  Some objects
are not going to have any significant signal in the other filters, but
for example, you measure magnitude of 36 for one of them!  This is
clearly unreliable (no dataset in current astronomy is able to detect
such a faint signal).  In another image with the same depth, using the
same filter, you might measure a magnitude of 30 for it, and yet another
might give you 33.  Furthermore, the total sum of pixel values might
actually be negative in some images of the same depth (due to noise).
In these cases, no magnitude can be defined and MakeCatalog will place a
NaN there (recall that a magnitude is a base-10 logarithm).

   Using such unreliable measurements will directly affect our analysis,
so we must not use the raw measurements.  When approaching the limits of
your detection method, it is therefore important to be able to identify
such cases.  But how can we know how reliable a measurement of one
object on a given dataset is?

   When we confront such unreasonably faint magnitudes, there is one
thing we can deduce: that if something actually exists under our labeled
pixels (possibly buried deep under the noise), it's inherent magnitude
is fainter than an _upper limit magnitude_.  To find this upper limit
magnitude, we place the object's footprint (segmentation map) over a
random part of the image where there are no detections, and measure the
sum of pixel values within the footprint.  Doing this a large number of
times will give us a distribution of measurements of the sum.  The
standard deviation ($\sigma$) of that distribution can be used to
quantify the upper limit magnitude for that particular object (given its
particular shape and area):

$$M_{up,n\sigma}=-2.5\times\log_{10}{(n\sigma_m)}+z \quad\quad [mag/target]$$

   Traditionally, faint/small object photometry was done using fixed
circular apertures (for example, with a diameter of $N$ arc-seconds) and
there was not much processing involved (to make a deep coadd).  Hence,
the upper limit was synonymous with the surface brightness limit
discussed above: one value for the whole image.  The problem with this
simplified approach is that the number of pixels in the aperture
directly affects the final distribution and thus magnitude.  Also the
image correlated noise might actually create certain patterns, so the
shape of the object can also affect the final result.  Fortunately, with
the much more advanced hardware and software of today, we can make
customized segmentation maps (footprint) for each object and have enough
computing power to actually place that footprint over many random
places.  As a result, the per-target upper-limit magnitude and general
surface brightness limit have diverged.

   When any of the upper-limit-related columns requested, MakeCatalog
will randomly place each target's footprint over the undetected parts of
the dataset as described above, and estimate the required properties.
The procedure is fully configurable with the options in *note
Upper-limit settings::.  You can get the full list of upper-limit
related columns of MakeCatalog with this command (the extra ‘--’ before
‘--upperlimit’ is necessary(1)):

     $ astmkcatalog --help | grep -- --upperlimit

‘--upperlimit’
     The upper limit value (in units of the input image) for this object
     or clump.  This is the sigma-clipped standard deviation of the
     random distribution, multiplied by the value of ‘--upnsigma’).
     This is very important for the fainter and smaller objects in the
     image where the measured magnitudes are not reliable.

‘--upperlimit-mag’
     The upper limit magnitude for this object or clump.  This is very
     important for the fainter and smaller objects in the image where
     the measured magnitudes are not reliable.

‘--upperlimit-onesigma’
     The $1\sigma$ upper limit value (in units of the input image) for
     this object or clump.  When ‘--upnsigma=1’, this column's values
     will be the same as ‘--upperlimit’.

‘--upperlimit-sigma’
     The position of the label's sum measured within the distribution of
     randomly placed upperlimit measurements in units of the
     distribution's $\sigma$ or standard deviation.

‘--upperlimit-quantile’
     The position of the label's sum within the distribution of randomly
     placed upperlimit measurements as a quantile (value between 0 or
     1).  If the object is brighter than the brightest randomly placed
     profile, a value of ‘inf’ is returned.  If it is less than the
     minimum, a value of ‘-inf’ is reported.

‘--upperlimit-skew’
     This column contains the non-parametric skew of the
     $\sigma$-clipped random distribution that was used to estimate the
     upper-limit magnitude.  Taking $\mu$ as the mean, $\nu$ as the
     median and $\sigma$ as the standard deviation, the traditional
     definition of skewness is defined as: $(\mu-\nu)/\sigma$.

     This can be a good measure to see how much you can trust the random
     measurements, or in other words, how accurately the regions with
     signal have been masked/detected.  If the skewness is strong (and
     to the positive), then you can tell that you have a lot of
     undetected signal in the dataset, and therefore that the
     upper-limit measurement (and other measurements) are not reliable.

   ---------- Footnotes ----------

   (1) Without the extra ‘--’, grep will assume that ‘--upperlimit’ is
one of its own options, and will thus abort, complaining that it has no
option with this name.


File: gnuastro.info,  Node: Morphology measurements nonparametric,  Next: Morphology measurements elliptical,  Prev: Upper limit measurements,  Up: MakeCatalog measurements on each label

7.4.4.7 Morphology measurements (non-parametric)
................................................

Morphology defined as a way to quantify the "shape" of an object in your
input image.  This includes both the position and value of the pixels
within your input labels.  There are many ways to define the morphology
of an object.  In this section, we will review the available
non-parametric measures of morphology.  By non-parametric, we mean that
no functional shape is assumed for the measurement.

   In *note Morphology measurements elliptical:: you can see some
parametric elliptical measurements (which are only valid when the object
is actually an ellipse).

‘--num-clumps’
     [Objects] The number of clumps in this object.

‘--area’
     The raw area (number of pixels/voxels) in any clump or object
     independent of what pixel it lies over (if it is NaN/blank or
     unused for example).

‘--arcsec2-area’
     The used (non-blank in values image) area of the labeled region in
     units of arc-seconds squared.  This column is just the value of the
     ‘--area’ column, multiplied by the area of each pixel in the input
     image (in units of arcsec^2).  Similar to the ‘--ra’ or ‘--dec’
     columns, for this option to work, the objects extension used has to
     have a WCS structure.

‘--area-min-val’
     The number of pixels that are equal to the minimum value of the
     labeled region (clump or object).

‘--area-max-val’
     The number of pixels that are equal to the maximum value of the
     labeled region (clump or object).

‘--area-xy’
     Similar to ‘--area’, when the clump or object is projected onto the
     first two dimensions.  This is only available for 3-dimensional
     datasets.  When working with Integral Field Unit (IFU) datasets,
     this projection onto the first two dimensions would be a
     narrow-band image.

‘--fwhm’
     The full width at half maximum (in units of pixels, along the
     semi-major axis) of the labeled region (object or clump).  The
     maximum value is estimated from the mean of the top-three pixels
     with the highest values, see the description under ‘--maximum’.
     The number of pixels that have half the value of that maximum are
     then found (value in the ‘--half-max-area’ column) and a radius is
     estimated from the area.  See the description under
     ‘--half-sum-radius’ for more on converting area to radius along
     major axis.

     Because of its non-parametric nature, this column is most reliable
     on clumps and should only be used in objects with great caution.
     This is because objects can have more than one clump (peak with
     true signal) and multiple peaks are not treated separately in
     objects, so the result of this column will be biased.

     Also, because of its non-parametric nature, this FWHM it does not
     account for the PSF, and it will be strongly affected by noise if
     the object is faint/diffuse So when half the maximum value (which
     can be requested using the ‘--maximum’ column) is too close to the
     local noise level (which can be requested using the ‘--sky-std’
     column), the value returned in this column is meaningless (its just
     noise peaks which are randomly distributed over the area).  You can
     therefore use the ‘--maximum’ and ‘--sky-std’ columns to remove, or
     flag, unreliable FWHMs.  For example, if a labeled region's maximum
     is less than 2 times the sky standard deviation, the value will
     certainly be unreliable (half of that is $1\sigma$!).  For a more
     reliable value, this fraction should be around 4 (so half the
     maximum is 2$\sigma$).

‘--half-max-area’
     The number of pixels with values larger than half the maximum flux
     within the labeled region.  This option is used to estimate
     ‘--fwhm’, so please read the notes there for the caveats and
     necessary precautions.

‘--half-max-radius’
     The radius of region containing half the maximum flux within the
     labeled region.  This is just half the value reported by ‘--fwhm’.

‘--half-max-sum’
     The sum of the pixel values containing half the maximum flux within
     the labeled region (or those that are counted in ‘--halfmaxarea’).
     This option uses the pixels within ‘--fwhm’, so please read the
     notes there for the caveats and necessary precautions.

‘--half-sum-area’
     The number of pixels that contain half the object or clump's total
     sum of pixels (half the value in the ‘--sum’ column).  To count
     this area, all the non-blank values associated with the given label
     (object or clump) will be sorted and summed in order (starting from
     the maximum), until the sum becomes larger than half the total sum
     of the label's pixels.

     This option is thus good for clumps (which are defined to have a
     single peak in their morphology), but for objects you should be
     careful: if the object includes multiple peaks/clumps at roughly
     the same level, then the area reported by this option will be
     distributed over all the peaks.

‘--half-sum-radius’
     Radius (in units of pixels) derived from the area that contains
     half the total sum of the label's pixels (value reported by
     ‘--halfsumarea’).  If the area is $A_h$ and the axis ratio is $q$,
     then the value returned in this column is $\sqrt{A_h/({\pi}q)}$.
     This option is a good measure of the concentration of the
     _observed_ (after PSF convolution and noisy) object or clump, But
     as described below it underestimates the effective radius.  Also,
     it should be used in caution with objects that may have multiple
     clumps.  It is most reliable with clumps or objects that have one
     or zero clumps, see the note under ‘--halfsumarea’.

     Recall that in general, for an ellipse with semi-major axis $a$,
     semi-minor axis $b$, and axis ratio $q=b/a$ the area ($A$) is
     $A={\pi}ab={\pi}qa^2$.  For a circle (where $q=1$), this simplifies
     to the familiar $A={\pi}a^2$.

     This option should not be confused with the _effective radius_ for
     Sérsic profiles, commonly written as $r_e$.  For more on the Sérsic
     profile and $r_e$, please see *note Galaxies::.  Therefore, when
     $r_e$ is meaningful for the target (the target is elliptically
     symmetric and can be parameterized as a Sérsic profile), $r_e$
     should be derived from fitting the profile with a Sérsic function
     which has been convolved with the PSF. But from the equation above,
     you see that this radius is derived from the raw image's labeled
     values (after convolution, with no parametric profile), so this
     column's value will generally be (much) smaller than $r_e$,
     depending on the PSF, depth of the dataset, the morphology, or if a
     fraction of the profile falls on the edge of the image.

     In other words, this option can only be interpreted as an effective
     radius if there is no noise and no PSF and the profile within the
     image extends to infinity (or a very large multiple of the
     effective radius) and it not near the edge of the image.

‘--frac-max1-area’
‘--frac-max2-area’
     Number of pixels brighter than the given fraction(s) of the maximum
     pixel value.  For the maximum value, see the description of
     ‘--maximum’ column.  The fraction(s) are given through the
     ‘--frac-max’ option (that can take two values) and is described in
     *note MakeCatalog inputs and basic settings::.  Recall that in
     ‘--halfmaxarea’, the fraction is fixed to 0.5.  Hence, added with
     these two columns, you can sample three parts of the profile area.

‘--frac-max1-sum’
‘--frac-max2-sum’
     Sum of pixels brighter than the given fraction(s) of the maximum
     pixel value.  For the maximum value, see the description of
     ‘--maximum’ column below.  The fraction(s) are given through the
     ‘--frac-max’ option (that can take two values) and is described in
     *note MakeCatalog inputs and basic settings::.  Recall that in
     ‘--halfmaxsum’, the fraction is fixed to 0.5.  Hence, added with
     these two columns, you can sample three parts of the profile's sum
     of pixels.

‘--frac-max1-radius’
‘--frac-max2-radius’
     Radius (in units of pixels) derived from the area that contains the
     given fractions of the maximum valued pixel(s) of the label's
     pixels (value reported by ‘--frac-max1-area’ or
     ‘--frac-max2-area’).  For the maximum value, see the description of
     ‘--maximum’ column below.  The fractions are given through the
     ‘--frac-max’ option (that can take two values) and is described in
     *note MakeCatalog inputs and basic settings::.  Recall that in
     ‘--fwhm’, the fraction is fixed to 0.5.  Hence, added with these
     two columns, you can sample three parts of the profile's radius.

‘--clumps-area’
     [Objects] The total area of all the clumps in this object.

‘--weight-area’
     The area (number of pixels) used in the flux weighted position
     calculations.

‘--geo-area’
     The area of all the pixels labeled with an object or clump.  Note
     that unlike ‘--area’, pixel values are completely ignored in this
     column.  For example, if a pixel value is blank, it will not be
     counted in ‘--area’, but will be counted here.

‘--geo-area-xy’
     Similar to ‘--geo-area’, when the clump or object is projected onto
     the first two dimensions.  This is only available for 3-dimensional
     datasets.  When working with Integral Field Unit (IFU) datasets,
     this projection onto the first two dimensions would be a
     narrow-band image.


File: gnuastro.info,  Node: Morphology measurements elliptical,  Next: Measurements per slice spectra,  Prev: Morphology measurements nonparametric,  Up: MakeCatalog measurements on each label

7.4.4.8 Morphology measurements (elliptical)
............................................

When your target objects are sufficiently ellipse-like, you can use the
measurements below to quantify the various parameters of the ellipse.
The derivation of the elliptical parameters are described in detail
below and followed by the available MakeCatalog column options.

   The shape or morphology of a target is one of the most commonly
desired parameters of a target.  Here, we will review the derivation of
the most basic/simple morphological parameters: the elliptical
parameters for a set of labeled pixels.  The elliptical parameters are:
the (semi-)major axis, the (semi-)minor axis and the position angle
along with the central position of the profile.  The derivations below
follow the SExtractor manual derivations with some added explanations
for easier reading.

   Let's begin with one dimension for simplicity: Assume we have a set
of $N$ values $B_i$ (for example, showing the spatial distribution of a
target's brightness), each at position $x_i$.  The simplest parameter we
can define is the geometric center of the object ($x_g$) (ignoring the
brightness values): $x_g=(\sum_ix_i)/N$.  _Moments_ are defined to
incorporate both the value (brightness) and position of the data.  The
first moment can be written as:

            $$\overline{x}={\sum_iB_ix_i \over \sum_iB_i}$$

This is essentially the weighted (by $B_i$) mean position.  The
geometric center ($x_g$, defined above) is a special case of this with
all $B_i=1$.  The second moment is essentially the variance of the
distribution:

      $$\overline{x^2}\equiv{\sum_iB_i(x_i-\overline{x})^2 \over
            \sum_iB_i} = {\sum_iB_ix_i^2 \over \sum_iB_i} -
      2\overline{x}{\sum_iB_ix_i\over\sum_iB_i} + \overline{x}^2
         ={\sum_iB_ix_i^2 \over \sum_iB_i} - \overline{x}^2$$

The last step was done from the definition of $\overline{x}$.  Hence,
the square root of $\overline{x^2}$ is the spatial standard deviation
(along the one-dimension) of this particular brightness distribution
($B_i$).  Crudely (or qualitatively), you can think of its square root
as the distance (from $\overline{x}$) which contains a specific amount
of the flux (depending on the $B_i$ distribution).  Similar to the first
moment, the geometric second moment can be found by setting all $B_i=1$.
So while the first moment quantified the position of the brightness
distribution, the second moment quantifies how that brightness is
dispersed about the first moment.  In other words, it quantifies how
"sharp" the object's image is.

   Before continuing to two dimensions and the derivation of the
elliptical parameters, let's pause for an important implementation
technicality.  You can ignore this paragraph and the next two if you do
not want to implement these concepts.  The basic definition (first
definition of $\overline{x^2}$ above) can be used without any major
problem.  However, using this fraction requires two runs over the data:
one run to find $\overline{x}$ and another run to find $\overline{x^2}$
from $\overline{x}$, this can be slow.  The advantage of the last
fraction above, is that we can estimate both the first and second
moments in one run (since the $-\overline{x}^2$ term can easily be added
later).

   The logarithmic nature of floating point number digitization creates
a complication however: suppose the object is located between pixels
10000 and 10020.  Hence the target's pixels are only distributed over 20
pixels (with a standard deviation $<20$), while the mean has a value of
$\sim10000$.  The $\sum_iB_i^2x_i^2$ will go to very very large values
while the individual pixel differences will be orders of magnitude
smaller.  This will lower the accuracy of our calculation due to the
limited accuracy of floating point operations.  The variance only
depends on the distance of each point from the mean, so we can shift all
position by a constant/arbitrary $K$ which is much closer to the mean:
$\overline{x-K}=\overline{x}-K$.  Hence we can calculate the second
order moment using:

        $$\overline{x^2}={\sum_iB_i(x_i-K)^2 \over \sum_iB_i} -
                         (\overline{x}-K)^2 $$

The closer $K$ is to $\overline{x}$, the better (the sums of squares
will involve smaller numbers), as long as $K$ is within the object
limits (in the example above: $10000\leq{K}\leq10020$), the floating
point error induced in our calculation will be negligible.  For the most
simplest implementation, MakeCatalog takes $K$ to be the smallest
position of the object in each dimension.  Since $K$ is arbitrary and an
implementation/technical detail, we will ignore it for the remainder of
this discussion.

   In two dimensions, the mean and variances can be written as:

             $$\overline{x}={\sum_iB_ix_i\over B_i}, \quad
           \overline{x^2}={\sum_iB_ix_i^2 \over \sum_iB_i} -
                           \overline{x}^2$$
             $$\overline{y}={\sum_iB_iy_i\over B_i}, \quad
           \overline{y^2}={\sum_iB_iy_i^2 \over \sum_iB_i} -
                           \overline{y}^2$$
            $$\quad\quad\quad\quad\quad\quad\quad\quad\quad
           \overline{xy}={\sum_iB_ix_iy_i \over \sum_iB_i} -
                   \overline{x}\times\overline{y}$$

   If an elliptical profile's major axis exactly lies along the $x$
axis, then $\overline{x^2}$ will be directly proportional with the
profile's major axis, $\overline{y^2}$ with its minor axis and
$\overline{xy}=0$.  However, in reality we are not that lucky and
(assuming galaxies can be parameterized as an ellipse) the major axis of
galaxies can be in any direction on the image (in fact this is one of
the core principles behind weak-lensing by shear estimation).  So the
purpose of the remainder of this section is to define a strategy to
measure the position angle and axis ratio of some randomly positioned
ellipses in an image, using the raw second moments that we have
calculated above in our image coordinates.

   Let's assume we have rotated the galaxy by $\theta$, the new second
order moments are:

        $$\overline{x_\theta^2} = \overline{x^2}\cos^2\theta +
                     \overline{y^2}\sin^2\theta -
                 2\overline{xy}\cos\theta\sin\theta $$
        $$\overline{y_\theta^2} = \overline{x^2}\sin^2\theta +
                     \overline{y^2}\cos^2\theta +
                 2\overline{xy}\cos\theta\sin\theta$$
     $$\overline{xy_\theta} = \overline{x^2}\cos\theta\sin\theta -
                 \overline{y^2}\cos\theta\sin\theta +
              \overline{xy}(\cos^2\theta-\sin^2\theta)$$

The best $\theta$ ($\theta_0$, where major axis lies along the
$x_\theta$ axis) can be found by:

$$\left.{\partial \overline{x_\theta^2} \over \partial \theta}\right|_{\theta_0}=0$$
   Taking the derivative, we get:
     $$2\cos\theta_0\sin\theta_0(\overline{y^2}-\overline{x^2}) +
        2(\cos^2\theta_0-\sin^2\theta_0)\overline{xy}=0$$ When
   $\overline{x^2}\neq\overline{y^2}$, we can write:
                           $$\tan2\theta_0 =
        2{\overline{xy} \over \overline{x^2}-\overline{y^2}}.$$

MakeCatalog uses the standard C math library's ‘atan2’ function to
estimate $\theta_0$, which we define as the position angle of the
ellipse.  To recall, this is the angle of the major axis of the ellipse
with the $x$ axis.  By definition, when the elliptical profile is
rotated by $\theta_0$, then $\overline{xy_{\theta_0}}=0$,
$\overline{x_{\theta_0}^2}$ will be the extent of the maximum variance
and $\overline{y_{\theta_0}^2}$ the extent of the minimum variance
(which are perpendicular for an ellipse).  Replacing $\theta_0$ in the
equations above for $\overline{x_\theta}$ and $\overline{y_\theta}$, we
can get the semi-major ($A$) and semi-minor ($B$) lengths:

        $$A^2\equiv\overline{x_{\theta_0}^2}= {\overline{x^2} +
\overline{y^2} \over 2} + \sqrt{\left({\overline{x^2}-\overline{y^2} \over 2}\right)^2 + \overline{xy}^2}$$

        $$B^2\equiv\overline{y_{\theta_0}^2}= {\overline{x^2} +
\overline{y^2} \over 2} - \sqrt{\left({\overline{x^2}-\overline{y^2} \over 2}\right)^2 + \overline{xy}^2}$$

   As a summary, it is important to remember that the units of $A$ and
$B$ are in pixels (the standard deviation of a positional distribution)
and that they represent the spatial light distribution of the object in
both image dimensions (rotated by $\theta_0$).  When the object cannot
be represented as an ellipse, this interpretation breaks down:
$\overline{xy_{\theta_0}}\neq0$ and $\overline{y_{\theta_0}^2}$ will not
be the direction of minimum variance.

   The list of options/columns to measure the elliptical properties of
an object's light distribution is given below.  Those that start with
‘--geo-*’ ignore the pixel values and just do the measurements on the
label's "geometric" shape.

‘--semi-major’
     The pixel-value weighted root mean square (RMS) along the
     semi-major axis of the profile (assuming it is an ellipse) in units
     of pixels.

‘--semi-minor’
     The pixel-value weighted root mean square (RMS) along the
     semi-minor axis of the profile (assuming it is an ellipse) in units
     of pixels.

‘--axis-ratio’
     The pixel-value weighted axis ratio (semi-minor/semi-major) of the
     object or clump.

‘--position-angle’
     The pixel-value weighted angle of the semi-major axis with the
     first FITS axis in degrees.

‘--geo-semi-major’
     The geometric (ignoring pixel values) root mean square (RMS) along
     the semi-major axis of the profile, assuming it is an ellipse, in
     units of pixels.

‘--geo-semi-minor’
     The geometric (ignoring pixel values) root mean square (RMS) along
     the semi-minor axis of the profile, assuming it is an ellipse, in
     units of pixels.

‘--geo-axis-ratio’
     The geometric (ignoring pixel values) axis ratio of the profile,
     assuming it is an ellipse.

‘--geo-position-angle’
     The geometric (ignoring pixel values) angle of the semi-major axis
     with the first FITS axis in degrees.


File: gnuastro.info,  Node: Measurements per slice spectra,  Prev: Morphology measurements elliptical,  Up: MakeCatalog measurements on each label

7.4.4.9 Measurements per slice (spectra)
........................................

When the input is a 3D data cube, MakeCatalog also has the following
multi-valued measurements per label (as well as the single-valued
measurements).  These will be stored as vector columns (*note Vector
columns::) in the same output table as the single-valued measurements.
For a tutorial on how to use these options and interpret their values,
see *note Detecting lines and extracting spectra in 3D data::.

   These options will do measurements on each 2D slice of the input 3D
cube; hence the common the format of ‘--*-in-slice’.  Each slice usually
corresponds to a certain wavelength, you can also think of these
measurements as spectra.

*No in-slice measurement for clumps:* Unfortunately due to
time-constrains, the columns described in this section are not yet
implemented for clumps.  For the status of this issue, see bug 66286
(https://savannah.gnu.org/bugs/?66286).

   For each row (input label), each of the columns described here will
contain multiple values as a vector column.  The number of measurements
in each column is the number of slices in the cube, or the size of the
cube along the third dimension.  To learn more about vector columns and
how to manipulate them, see *note Vector columns::.  For example usage
of these columns in the tutorial above, see *note 3D measurements and
spectra:: and *note Extracting a single spectrum and plotting it::.

There are two ways to do each measurement on a slice for each label:
Only label
     The measurement will only be done on the voxels in the slice that
     are associated to that label.  These types of per-slice measurement
     therefore have the following properties:
        • This will only be a measurement of that label and will not be
          affected by any other label.
        • The number of voxels used in each slice can be different
          (usually only one or two voxels at the two extremes of the
          label (along the third dimension), and many in the middle.
        • Since most labels are localized along the third dimension
          (maybe only covering 20 slices out of thousands!), many of the
          measurements (on slices where the label doesn't exist) will be
          NaN (for the sum measurements for example) or 0 (for the area
          measurements).
Projected label
     MakeCatalog will first project the 3D label into a 2D surface
     (along the third dimension) to get its 2D footprint.  Afterwards,
     all the voxels in that 2D footprint will be measured all slices.
     All these measurements will have a ‘-proj-’ component in their
     name.  These types of per-slice measurement therefore has the
     following properties:

        • A measurement will be done on each slice of the cube.
        • All measurements will be done on the same surface area.
        • Labels can overlap when they are projected onto the first two
          FITS dimensions (the spatial coordinates, not spectral).  As a
          result, other emission lines or objects may contaminate the
          resulting spectrum for each label.

     To help separate other labels, MakeCatalog can do a third type of
     measurement on each slice: measurements on the voxels that belong
     to other labels but overlap with the 2D projection.  This can be
     used to see how much your projected measurement is affected by
     other emission sources (on the projected spectra) and also if
     multiple lines (labeled regions) belong to the same physical
     object.  These measurements contain ‘-other-’ in their name.

‘--sum-in-slice’
     [Only label] Sum of values in each slice.

‘--sum-err-in-slice’
     [Only label] Error in '-sum-in-slice'.

‘--area-in-slice’
     [Only label] Number of labeled in each slice.

‘--sum-proj-in-slice’
     [Projected label] Sum of projected area in each slice.

‘--area-proj-in-slice:’
     [Projected label] Number of voxels that are used in
     ‘--sum-proj-in-slice’.

‘--sum-proj-err-in-slice’
     [Projected label] Error of ‘--sum-proj-in-slice’.

‘--area-other-in-slice’
     [Projected label] Area of other label in projected area on each
     slice.

‘--sum-other-in-slice’
     [Projected label] Sum of other label in projected area on each
     slice.

‘--sum-other-err-in-slice:’
     [Projected label] Area in ‘--sum-other-in-slice’.


File: gnuastro.info,  Node: Metameasurements on full input,  Next: Manual metameasurements,  Prev: MakeCatalog measurements on each label,  Up: MakeCatalog

7.4.5 Metameasurements on full input
------------------------------------

The data (rows and columns) produced by MakeCatalog are independent
measurements of each label and were the focus of *note MakeCatalog
measurements on each label::.  However, to correctly interpret those
measurements especially in light of other studies, we also need
statistical analysis of the full input or catalog.  These
metameasurements (measurements on/about independent/individual
measurements) are a single value for the whole input and are therefore
stored as meta-data (keywords) in the 0th (first) HDU of the
MakeCatalog's output.  The most common type of such metameasurements are
commonly known as "limits" or "depth" of a dataset.  These are also the
only type that is currently implemented in MakeCatalog; so let's
continue with them.

   No measurement on a real dataset can be perfect: you can only reach a
certain level/limit of accuracy and a meaningful (scientific) analysis
on the output catalog requires an understanding of these limits.
Different datasets have different noise properties and different
detection methods(1) will have different abilities to detect or measure
certain kinds of signal (astronomical objects) and their properties in
the dataset.  Hence, quantifying the detection and measurement
limitations with a particular dataset and analysis tool is the most
crucial/critical aspect of any high-level analysis.  Due to their
importance, we have already touched upon some of these in two tutorials
on real data: see *note Measuring the dataset limits:: and *note Image
surface brightness limit::.

   The sections below describe each of the metameasurements that
MakeCatalog will calculate and report when given the ‘--meta-measures’
option.  Here, we just focus on the concept behind each measurement and
its derviation.  The actual keywords that are written for each is
described in *note MakeCatalog output keywords::.

* Menu:

* Surface brightness limit of image::  A portable measure of the noise level.
* Noise based magnitude limit of image::  Magnitude of objects above a certain threshold of the noise.
* Confusion limit of image::    A measure of density in the image.

   ---------- Footnotes ----------

   (1) A method/algorithm/software that is run with a different set of
parameters is considered as a different detection method


File: gnuastro.info,  Node: Surface brightness limit of image,  Next: Noise based magnitude limit of image,  Prev: Metameasurements on full input,  Up: Metameasurements on full input

7.4.5.1 Surface brightness limit of image
.........................................

As we make more observations on one region of the sky and add/combine
the observations into one dataset, both the signal and the noise
increase.  However, the signal increases much faster than the noise:
Assuming you add $N$ datasets with equal exposure times, the signal will
increases as a multiple of $N$, while noise increases as $\sqrt{N}$.
Therefore the signal-to-noise ratio increases by a factor of $\sqrt{N}$.
Visually, fainter (per pixel) parts of the objects/signal in the image
will become more visible/detectable.  The noise-level is known as the
dataset's surface brightness limit.

   You can think of the noise as muddy water in a pond (the "sky" will
be the bottom of the pond).  The signal (coming from astronomical
objects in real data) will be the rocks and their peaks can sometimes
reach above the muddy water.  Let's assume that in your first
observation the muddy water has just been stirred and except a few small
peaks, you cannot see anything through the mud.  As you wait (and make
more observations/exposures), the mud settles down and the _depth_ of
the transparent water increases.  As a result, more and more summits
become visible and the lower parts of the rocks (parts with lower
surface brightness) can be seen more clearly.  In this analogy(1),
height (from the ground) is the _surface brightness_ and the height of
the muddy water at the moment you do your measurements, is your _surface
brightness limit_.

   On different instruments, pixels cover different spatial angles over
the sky.  For example, the width of each pixel on the ACS camera on the
Hubble Space Telescope (HST) is roughly 0.05 seconds of arc, while the
pixels of SDSS are each 0.396 seconds of arc (almost eight times
wider(2)).  Nevertheless, irrespective of its sky coverage, a pixel is
our unit of data collection.

   To start with, we define the low-level Surface brightness limit or
_depth_, in units of magnitude/pixel with the equation below (assuming
the image has zero point magnitude $z$ and we want the $n$th multiple of
$\sigma_p$).  Where $\sigma_p$ is the per-pixel standard deviation (for
example median of the measured sky standard deviation in the image).

$$SB_{n\sigma,\rm pixel}=-2.5\times\log_{10}{(n\sigma_p)}+z \quad\quad [mag/pixel]$$

   As an example, the XDF survey covers part of the sky that the HST has
observed the most (for 85 orbits) and is consequently very small
($\sim4$ minutes of arc, squared).  On the other hand, the CANDELS
survey, is one of the widest multi-color surveys done by the HST
covering several fields (about 720 arcmin$^2$) but its deepest fields
have only 9 orbits observation.  The $1\sigma$ depth of the XDF and
CANDELS-deep surveys in the near infrared WFC3/F160W filter are
respectively 34.40 and 32.45 magnitudes/pixel.  In a single orbit image,
this same field has a $1\sigma$ depth of 31.32 magnitudes/pixel.  Recall
that a larger magnitude corresponds to fainter objects, see *note
Brightness flux magnitude::.

   The low-level magnitude/pixel measurement above is only useful when
all the datasets you want to use, or compare, have the same pixel size.
However, you will often find yourself using, or comparing, datasets from
various instruments with different pixel scales (projected pixel width,
in arc-seconds).  If we know the pixel scale, we can obtain a more
easily comparable surface brightness limit in units of:
magnitude/arcsec$^2$.  But another complication is that astronomical
objects are usually larger than 1 arcsec$^2$.  As a result, it is common
to measure the surface brightness limit over a larger (but fixed,
depending on context) area.

   Let's assume that every pixel is $p$ arcsec$^2$ and we want the
surface brightness limit for an object covering A arcsec$^2$ (so $A/p$
is the number of pixels that cover an area of $A$ arcsec$^2$).  On the
other hand, noise is added in RMS(3), hence the noise level in $A$
arcsec$^2$ is $n\sigma_p\sqrt{A/p}$.  But we want the result in units of
arcsec$^2$, so we should divide this by $A$ arcsec$^2$:
$n\sigma_p\sqrt{A/p}/A=n\sigma_p\sqrt{A/(pA^2)}=n\sigma_p/\sqrt{pA}$.
Plugging this into the magnitude equation, we get the $n\sigma$ surface
brightness limit, over an area of A arcsec$^2$, in units of
magnitudes/arcsec$^2$ (see the end of this section on how to derive this
using Gnuastro):

$$SB_{{n\sigma,\rm A arcsec}^2}=-2.5\times\log_{10}{\left(n\sigma_p\over \sqrt{pA}\right)+z} \quad\quad [mag/arcsec^2]$$

   As you saw in its derivation, the calculation above extrapolates the
noise in one pixel over all the input's pixels!  In other words, all
pixels are treated independently in the measurement of the standard
deviation.  It therefore implicitly assumes that the noise is the same
in all of the pixels.  But this only happens in individual exposures:
reduced data will have correlated noise because they are a coadd of many
individual exposures that have been warped (thus mixing the pixel
values).  A more accurate measure which will provide a realistic value
for every labeled region is known as the _upper-limit magnitude_, which
is discussed in the next section (*note Upper limit surface brightness
of image::).

   Within Gnuastro, measuring the surface brighntess limit of the image
is very easy: the outputs of NoiseChisel include the Sky standard
deviation ($\sigma$) on every group of pixels (a tile) that were
calculated from the undetected pixels in each tile, see *note
Tessellation:: and *note NoiseChisel output::.  The $\sigma_p$ used
above is recorded in the ‘MEDSTD’ keyword of the ‘SKY_STD’ extension of
NoiseChisel's output.  Therefore, the first thing you need to do is to
run NoiseChisel on the image.

   When MakeCatalog is called with ‘--meta-measurements’, it will
calculate the input dataset's $SB_{{n\sigma,\rm A arcsec}^2}$ and will
write it ‘SBL’ keywords the 0th HDU of the output, see *note MakeCatalog
output keywords::.  You can set your desired $n$-th multiple of $\sigma$
and the $A$ arcsec$^2$ area using the following two options
respectively: ‘--sbl-sigma’ and ‘--sbl-area’ (see *note MakeCatalog
output keywords::).  Just note that $SB_{{n\sigma,\rm A arcsec}^2}$ is
only calculated if the input has World Coordinate System (WCS). Without
WCS, the pixel scale cannot be derived.

   In case you just want the surface brightness limit of an image
without separating the objects and clumps within the image or doing any
measurements, you can use the series of commands below.  They assume
your image is called ‘image.fits’ and its zero point magnitude is 25.
For a more detailed hands-on understanding of the surface brightness
limit, see *note Image surface brightness limit:: (which is part of a
larger tutorial).

     $ zp=25.0
     $ astnoisechisel image.fits -o nc.fits
     $ astmkcatalog nc.fits -hDETECTIONS --sn --zeropoint=$zp \
                    --output=cat.fits --meta-measures
     $ astfits cat.fits --keyvalue=SBL --quiet | asttable -Y

   Because the surface brightness limit is just an extrapolation of the
noise level, we can extend the concept to use a reference surface
brightness limit to find the expected surface brightness limit of
another telescope (in a similar location), or a different exposure time
for the same telescope.  This is discussed more fully in *note Expected
surface brightness limit::.

   ---------- Footnotes ----------

   (1) Note that this muddy water analogy is not perfect, because while
the water-level remains the same all over a peak, in data analysis, the
Poisson noise increases with the level of data.

   (2) Ground-based instruments like the SDSS suffer from strong
smoothing due to the atmosphere.  Therefore, increasing the pixel
resolution (or decreasing the width of a pixel) will not increase the
received information).

   (3) If you add three datasets with noise $\sigma_1$, $\sigma_2$ and
$\sigma_3$, the resulting noise level is
$\sigma_t=\sqrt{\sigma_1^2+\sigma_2^2+\sigma_3^2}$, so when
$\sigma_1=\sigma_2=\sigma_3\equiv\sigma$, then
$\sigma_t=\sigma\sqrt{3}$.  In this case, the area $A$ is covered by
$A/p$ pixels, so the noise level is $\sigma_t=\sigma\sqrt{A/p}$.


File: gnuastro.info,  Node: Noise based magnitude limit of image,  Next: Confusion limit of image,  Prev: Surface brightness limit of image,  Up: Metameasurements on full input

7.4.5.2 Noise based magnitude limit of image
............................................

Assuming $\sigma_p$ to be per-pixel noise level, and we apply a
threshold of $n\sigma_p$ on the image, what would be the faintest
magnitude of an object that covers more than $A$ square arcseconds?
Answering this question is the goal of the noise-based magnitude limit
(nML or NML) is defined for (for more on $\sigma_p$, see *note Surface
brightness limit of image::).  The commonly used parameters are: $n=5$
and $A=7$ square arcseconds (approximate area of a circle with radius
1.5 arcseconds).

   Let's show the area of a pixel (in square arcseconds) as $p$.  The
arbitrary area $A$ covers $A/p$ pixels (which is just a counter; without
any units).  On the other hand, $n\sigma_p$ (the threshold we want to
apply) is the minimum flux that each pixel should have.  Assuming all
the pixels above the threshold have the same $n\sigma$ flux, then the
total flux above the threshold would be $n\sigma\times A/p$.  Plugging
this into the definition of magnitudes (see *note Brightness flux
magnitude::), we get the noise based magnitude limit:

  $$nML = -2.5 \log_{10}\left(n\sigma_p\times\frac{A}{p}\right) + z$$

There are several caveats to this measure:

   • Astronomical objects are not flat!  Their flux gradually sinks
     deeply into the noise; see Figure 1 of Akhlaghi & Ichikawa 2015
     (https://arxiv.org/abs/1505.01664).  Therefore, the centers of
     stars and galaxies will be (much) brighter than the threshold and
     major fraction of the object's flux will not be counted in this
     measure (more as $n$ increases).
   • The area used in the equation is just a number (without any
     shape!).  In other words, it is an extrapolated/abstract area that
     is shapeless.
   • Since this depends on a per-pixel measurement of $\sigma_p$, it
     does not account for correlated noise.  For more on this, see *note
     Surface brightness limit of image::.

   Given the caveats above, it may be strage to see that the noise based
magnitude limit is still reported in many reports (it is sometimes
called simply as "magnitude limit" or even "detection limit").  So let's
review the history of this metameasure for a better understanding of its
usage.  When computers started to enter astronomical data analysis
(1970s) their processing powers were very low (compared to today!).
Also, images were taken on photographic plates and digitized/scanned
later (not too clean noise).  Therefore, they needed to apply high
thresholds ($+2\sigma_p$), find groups of pixels that were contiguously
above that threshold, and define the "real" objects based on their
contiguous area above that threshold.  Therefore, this limit was also
known as "detection limit".

   Because this lost flux was a problem, authors like Petrosian 1976
(https://ui.adsabs.harvard.edu/abs/1976ApJ...209L...1P) or Kron
<https://ui.adsabs.harvard.edu/abs/1980ApJS...43..305K> suggested ways
of finding the object's "total" flux (accounting for the flux below the
threshold also).  But with new paradigms to detection coming in the
current century (for example Gnuastro's *note NoiseChisel::), the
"detection limit" name will cause confusion and mis-interpretation.
Therefore in Gnuastro we call this limit the "noise-based magnitude
limit" and include it only for legacy/historical reasons (it is just an
extrapolatoin).

   The reason for the "noise-based" prefix is that the equation above is
purely based on the noise (does not depend on the actual signal!).
However, different objects have different profiles: some are
sharp/concentrated (like stars or elliptical galaxies), others are
diffuse (like ultra diffuse galaxies) and others are clumpy (like a
globular cluster) and etc.  Therefore, we classify the more robust
"magnitude limit" (that actually takes morphology into account) as a
_manual_ metameasurement (that you should do after calling MakeCatalog
based on your desired targets/goals), it is described in *note Magnitude
limit for certain objects::.

   When MakeCatalog is called with ‘--meta-measurements’ and
‘--clumpscat’, it will calculate the input dataset's noise-based
magnitude limit and will write it ‘NML’ keyword the 0th HDU of the
output.  You can set your desired $n$-th multiple of $\sigma$ and the
$A$ arcsec$^2$ area using the following two options respectively:
‘--nml-sigma’ and ‘--nml-area’.  See *note MakeCatalog output keywords::
for more.


File: gnuastro.info,  Node: Confusion limit of image,  Prev: Noise based magnitude limit of image,  Up: Metameasurements on full input

7.4.5.3 Confusion limit of image
................................

The confusion limit is expressed in units of distance (pixels or
arcseconds), derived from the distribution of the distances of all
resolved sources with their nearest neighbor.  The confusion limit
therefore depends on several factors, including the pointing on the sky.
For example, a shallow image taken on the disk of M31, or the center of
a large globular cluster, or close to the disk of the Milky Way will be
heavily contaminated by confusion limit (resolved sources will be very
close to each other) even if it has a low exposure time on a small
telescope.  But at exactly the same imaging conditions (exposure time,
telescope aperture, PSF and etc), a random pointing away from the disk
of the Milky Way will have a good confusion limit (resolved sources are
more distant).  But as we go very deep (more exposure time or larger
telescope) without a similar improvement in the pixel scale, the
confusion limit will become mostly uniform (assuming a relatively
uniform background galaxy population).  In wide field imaging, it is
also important to account for the projection effects when calculating
the projected distances (so centers have to be in RA and Dec; and
distances have to be measured on a sphere, not on a plane).

   Therefore, unlike the previously discussed automatic
meta-measurements like *note Surface brightness limit of image:: and
*note Noise based magnitude limit of image:: the confusion limit depends
on the pointing on sky.  This is because those are defined purely based
on the noise of the image, but the confusion limit comes from the
resolved signal's distribution in the image.  But as we go very deep in
extra-galactic fields (high galactic latitudes: away from the disk of
the milky wayw) the density of sources as you go deeper should become
mostly uniform (assuming a relatively uniform background galaxy
population).

   The confusion limit is found by finding the nearest resolved source
to each source (Segment's clumps by default, see *note Segment::).
Technically (when calculating the distance between the clump and its
nearest neighbor), it is also important to account for the spherical
nature of RA and Dec and to use spherical distance measures (see
description of ‘gal_wcs_angular_distance_deg’ in *note World Coordinate
System::).  Otherwise (if using pixel coordinates/distances) the
confusion limit will be affected by projection/distortion issues.

   The distribution of distances to the nearest neighbor (in arcsec)
will naturally be heavily skewed (since distance is always positive).
The distribution relies heavily on the depth, pixel scale and PSF of the
image (which can be very different from one image to another).
Therefore, the best measure we have found so far is to use the "width"
of the distribution (difference between the 0.75 and 0.25 quantiles).
However, this distance is most relevant in pixel units, not arcseconds.
This is because the confusion limit depends on detection and
segmentation which are done on the pixel grid not in the abstract,
continuous or instrument-agnostic units of arcseconds.

   When MakeCatalog is called with ‘--meta-measurements’ and
‘--clumpscat’, it will measure the confusion limit and store it in the
‘CNL’ output keyword (see *note MakeCatalog output keywords::).  Further
information (five percentiles) about the distribution of distances to
the nearest neighbor are provided in the ‘CNLP*’.  The percentile
information is necessary because the distribution to the nearest label
is by nature a very non-symmetric distribution (skewed to the positive).
As a result, it cannot be simply parametrized as in the other
metameasures here.

   Therefore as described in , MakeCatalog doesn't just return the
median (at quantile of 0.5 or percentile of 50), but several other
quantiles.  Ultimately, you can see each clump's nearest neighbor as
well as the spherical distance to its nearest clump by activating the
‘--cnl-check’ option (see *note MakeCatalog output HDUs::).


File: gnuastro.info,  Node: Manual metameasurements,  Next: Adding new columns to MakeCatalog,  Prev: Metameasurements on full input,  Up: MakeCatalog

7.4.6 Manual metameasurements
-----------------------------

The concept of metameasurements, and the automatically calculated
metameasurements of MakeCatalog were introduced in *note
Metameasurements on full input::.  There are other metameasures that
cannot (as of this writing) be done automatically in the same run of
MakeCatalog and need a customized execution of MakeCatalog.  In this
section these types of metameasures are introduced.

* Menu:

* Expected surface brightness limit::
* Upper limit surface brightness of image::  Necessary for diffuse regions of image.
* Magnitude limit for certain objects::
* Completeness limit for certain objects::


File: gnuastro.info,  Node: Expected surface brightness limit,  Next: Upper limit surface brightness of image,  Prev: Manual metameasurements,  Up: Manual metameasurements

7.4.6.1 Expected surface brightness limit
.........................................

The Surface brighntess limit was defined in *note Surface brightness
limit of image::.  As we saw there, it is ultimately a theoretical
extrapolation of one pixel's noise level.  Therefore, we can use a
reference image's surface brightness limit ($SB_r$) to derive an
expected surface brightness limit of another image (on another telescope
with another exposure time, but with _the same_ filter; let's call it
$SB_i$).  Assuming no correlated noise (which is only valid on a single
exposure!), the $\sigma_p$ described above is purely due to the Poisson
noise of the background signal (for example Zodiacal light, light
pollution or etc).  Therefore, taking $S$ to be the exposed surface of
the telescope's primary mirror (or lens) and $t$ to be the exposure
time, the signal of the background will increase with $S\times t$ and
therefore $\sigma_p\propto \sqrt{St}$.  Plugging this into two instances
of the equation above, allows us to derive $SB_i$ from $SB_r$

$$SB_i = SB_r + 2.5\log_{10}\left( {n_i \over n_r} \sqrt{{S_i \over S_r}{t_i \over t_r}{A_i \over A_r}} \right)$$

   Since almost all mirrors are circular, we can simplify the relation
above by replacing the exposed surface with the exposed radius
(accounting for the non-exposed area in most reflective telescopes due
to the secondary mirror) as shown in the quation below.  Note that we
didn't simply say "radius" (and the equation does not have $r$), but
"exposed radius" ($r_e$ in the equation).  This is a very important
factor to have in mind.  For example in the Vera C. Rubin Observatory
(that will be conducting the LSST survey) the primary mirror has a
diameter of 8.4 meters, however, the inner cicular area of radius 5
meters is not used (due to the second and third mirrors)(1).  Therefore,
the useful surface of the Vera C. Rubin telescope is
$\pi(8.4^2-5.0^2)=\pi45.56=\pi6.75^2$, giving it an exposed radius of
$r_e=6.75m$ (for an easy implementation of this equation in Gnuastro,
see the ‘sblim-diff’ operator of *note Unit conversion operators::).

$$SB_i = SB_r + 2.5\log_{10}\left( {r_{ei} \over r_{er}}{n_i \over n_r} \sqrt{{t_i \over t_r}{A_i \over A_r}} \right)$$

   ---------- Footnotes ----------

   (1) <https://commons.wikimedia.org/wiki/File:LSSToptics.jpg>


File: gnuastro.info,  Node: Upper limit surface brightness of image,  Next: Magnitude limit for certain objects,  Prev: Expected surface brightness limit,  Up: Manual metameasurements

7.4.6.2 Upper limit surface brightness of image
...............................................

As mentioned in *note Surface brightness limit of image::, the surface
brightness limit assumes independent pixels when deriving the standard
deviation (the main input in the equation).  It just extrapolates the
standard deviation derived from one pixel to the requested area.  But as
mentioned at the end of that section, we have correlated noise in our
science-ready (deep) images and the noise of the pixels are not
independent.

   Because of this, the surface brightness limit will always
under-estimate the surface brightness (give fainter values than what is
statistically possible in the data for the requested area).  To account
for the correlated noise in the images, we need to derive the standard
deviation over a group of pixels that fall within a certain
footprint/shape.  For example over a circular aperture of radius 5.6419
arcsec, or a square with a side length of $10$ arcsec.  Depending on the
correlated noise systematics, the limit can be (very) different for
different shapes, even if they have the same area (as in the circle and
square mentioned in the previous sentence: both have an area of 100
arcsec$^2$).

   Therefore we need to derive the standard deviation that goes into the
surface brightness limit equation over a certain footprint/shape.  To do
that, we should:
  1. Place the desired footprint many times randomly over all the
     undetected pixels in an image.  In MakeCatalog, the number of these
     random positions can be configured with ‘--upnum’ and you can check
     their positions with ‘--checkuplim’.
  2. Calculate the sum of pixel values in each randomly placed
     footprint.
  3. Calculate the sigma-clipped standard deviation of the resulting
     distribution (of sum of pixel values in the randomly placed
     apertures).  Therefore, each footprint's measurement is independent
     of the other.
  4. Calculate the surface brightness of that standrad deviation (after
     multiplying it with your desired multiple of sigma).  For the
     definition of surface brightness, see *note Brightness flux
     magnitude::.

   The measurements over randomly placed apertures may remind you of
*note Upper limit measurements::.  Generally, the "upper limit" prefix
is given to all measurements with this way of measurement.  Therefore
this limit is called "Upper limit surface brightness" of an image (for a
multiple of sigma, over a certain shape).

   Traditionally a circular aperture of a fixed radius (in arcseconds)
has been used.  In Gnuastro, a labeled image containing the desired
shape/aperture can be generated with MakeProfiles.  The position of the
label is irrelevant because the upper limit measurements are done on the
many randomly placed footprints in undetected regions (independent of
where the label is positioned).  That labeled image should then be given
to MakeCatalog, while requesting ‘--upperlimit-sb’.  Of course, all
detected signal in the image needs to be masked (set to blank/NaN) so
MakeCatalog doesn't use randomly placed apertures that overlap with
detected signal in the image.

   Going into the implementation details can get pretty hard to follow
in English, so a full hands-on tutorial is available in the second half
of *note Image surface brightness limit::.  Read that tutorial with the
same input images and run the commands, and see each output image to get
a good understanding of how to properly measure the upper limit surface
brightness of your images.


File: gnuastro.info,  Node: Magnitude limit for certain objects,  Next: Completeness limit for certain objects,  Prev: Upper limit surface brightness of image,  Up: Manual metameasurements

7.4.6.3 Magnitude limit for certain objects
...........................................

Suppose we have taken two images of the same field of view with the same
CCD, once with a smaller telescope, and once with a larger one.  Because
we used the same CCD, the noise will be very similar.  However, the
larger telescope has gathered more light, therefore the same star or
galaxy will have a higher signal-to-noise ratio (S/N) in the image taken
with the larger one.  The same applies for a coadded image of the field
compared to a single-exposure image of the same telescope.

   This concept is used by some researchers to define the "magnitude
limit" or "detection limit" at a certain S/N (sometimes 10, 5 or 3 for
example, also written as $10\sigma$, $5\sigma$ or $3\sigma$).  To do
this, they measure the magnitude and signal-to-noise ratio of all the
objects that have similar within an image and measure the mean (or
median) magnitude of objects at the desired S/N. A fully working example
of deriving the magnitude limit is available in the tutorials section:
*note Measuring the dataset limits::.

   However, this method should be used with extreme care!  This is
because the shape of the object becomes important in this method: a
sharper object will have a higher _measured_ S/N compared to a more
diffuse object at the same original magnitude.  Besides the inherent
shape/sharpness of the object, issues like the PSF also become important
in this method (because the finally observed shapes of objects are
important here): two surveys with the same surface brightness limit (see
*note Surface brightness limit of image::) will have different magnitude
limits if one is taken from space and the other from the ground.


File: gnuastro.info,  Node: Completeness limit for certain objects,  Prev: Magnitude limit for certain objects,  Up: Manual metameasurements

7.4.6.4 Completeness limit for certain objects
..............................................

As the surface brightness of the objects decreases, the ability to
detect them will also decrease.  An important statistic is thus the
fraction of objects of similar morphology and magnitude that will be
detected with our detection algorithm/parameters in a given image.  This
fraction is known as _completeness_.  -For brighter objects,
completeness is 1: all bright objects that might exist over the image
will be detected.  However, as we go to objects of lower overall surface
brightness, we will fail to detect a fraction of them, and fainter than
a certain surface brightness level (for each morphology), nothing will
be detectable in the image: you will need more data to construct a
"deeper" image.  For a given profile and dataset, the magnitude where
the completeness drops below a certain level (usually above $90\%$) is
known as the completeness limit.

   Another important parameter in measuring completeness is purity: the
fraction of true detections to all detections.  In effect purity is the
measure of contamination by false detections: the higher the purity, the
lower the contamination.  Completeness and purity are anti-correlated:
if we can allow a large number of false detections (that we might be
able to remove by other means), we can significantly increase the
completeness limit.

   One traditional way to measure the completeness and purity of a given
sample is by embedding mock profiles in regions of the image with no
detection.  However in such a study we must be really careful to choose
model profiles as similar to the target of interest as possible.


File: gnuastro.info,  Node: Adding new columns to MakeCatalog,  Next: Invoking astmkcatalog,  Prev: Manual metameasurements,  Up: MakeCatalog

7.4.7 Adding new columns to MakeCatalog
---------------------------------------

MakeCatalog is designed to allow easy addition of different measurements
over a labeled image; see Akhlaghi 2016
(https://arxiv.org/abs/1611.06387v1).  A check-list style description of
necessary steps to do that is described in this section.  The common
development characteristics of MakeCatalog and other Gnuastro programs
is explained in *note Developing::.  We strongly encourage you to have a
look at that chapter to greatly simplify your navigation in the code.
After adding and testing your column, you are most welcome (and
encouraged) to share it with us so we can add to the next release of
Gnuastro for everyone else to also benefit from your efforts.

   MakeCatalog will first pass over each label's pixels two times and do
necessary raw/internal calculations.  Once the passes are done, it will
use the raw information for filling the final catalog's columns.  In the
first pass it will gather mainly object information and in the second
run, it will mainly focus on the clumps, or any other measurement that
needs an output from the first pass.  These two passes are designed to
be raw summations: no extra processing.  This will allow parallel
processing and simplicity/clarity.  So if your new calculation, needs
new raw information from the pixels, then you will need to also modify
the respective ‘mkcatalog_first_pass’ and ‘mkcatalog_second_pass’
functions (both in ‘bin/mkcatalog/mkcatalog.c’) and define new raw table
columns in ‘main.h’ (hopefully the comments in the code are clear
enough).

   In all these different places, the final columns are sorted in the
same order (same order as *note Invoking astmkcatalog::).  This allows a
particular column/option to be easily found in all steps.  Therefore in
adding your new option, be sure to keep it in the same relative place in
the list in all the separate places (it does not necessarily have to be
in the end), and near conceptually similar options.

‘main.h’
     The ‘objectcols’ and ‘clumpcols’ enumerated variables (‘enum’)
     define the raw/internal calculation columns.  If your new column
     requires new raw calculations, add a row to the respective list.
     If your calculation requires any other settings parameters, you
     should add a variable to the ‘mkcatalogparams’ structure.

‘ui.c’
     If the new column needs raw calculations (an entry was added in
     ‘objectcols’ and ‘clumpcols’), specify which inputs it needs in
     ‘ui_necessary_inputs’, similar to the other options.  Afterwards,
     if your column includes any particular settings (you needed to add
     a variable to the ‘mkcatalogparams’ structure in ‘main.h’), you
     should do the sanity checks and preparations for it here.

‘ui.h’
     The ‘option_keys_enum’ associates a unique value for each option to
     MakeCatalog.  The options that have a short option version, the
     single character short comment is used for the value.  Those that
     do not have a short option version, get a large integer
     automatically.  You should add a variable here to identify your
     desired column.

‘args.h’
     This file specifies all the parameters for the GNU C library, Argp
     structure that is in charge of reading the user's options.  To
     define your new column, just copy an existing set of parameters and
     change the first, second and 5th values (the only ones that differ
     between all the columns), you should use the macro you defined in
     ‘ui.h’ here.

‘columns.c’
     This file contains the main definition and high-level calculation
     of your new column through the ‘columns_define_alloc’ and
     ‘columns_fill’ functions.  In the first, you specify the basic
     information about the column: its name, units, comments, type (see
     *note Numeric data types::) and how it should be printed if the
     output is a text file.  You should also specify the raw/internal
     columns that are necessary for this column here as the many
     existing examples show.  Through the types for objects and rows,
     you can specify if this column is only for clumps, objects or both.

     The second main function (‘columns_fill’) writes the final value
     into the appropriate column for each object and clump.  As you can
     see in the many existing examples, you can define your processing
     on the raw/internal calculations here and save them in the output.

‘mkcatalog.c’
     This file contains the low-level parsing functions.  To be
     optimized, the parsing is done in parallel through the
     ‘mkcatalog_single_object’ function.  This function initializes the
     necessary arrays and calls the lower-level ‘parse_objects’ and
     ‘parse_clumps’ for actually going over the pixels.  They are all
     heavily commented, so you should be able to follow where to add
     your necessary low-level calculations.

‘doc/gnuastro.texi’
     Update this manual and add a description for the new column.


File: gnuastro.info,  Node: Invoking astmkcatalog,  Prev: Adding new columns to MakeCatalog,  Up: MakeCatalog

7.4.8 Invoking MakeCatalog
--------------------------

MakeCatalog will do measurements and produce a catalog from a labeled
dataset and optional values dataset(s).  The executable name is
‘astmkcatalog’ with the following general template

     $ astmkcatalog [OPTION ...] InputImage.fits

One line examples:

     ## Create catalog with RA, Dec, Magnitude and Magnitude error,
     ## from Segment's output:
     $ astmkcatalog --ra --dec --magnitude seg-out.fits

     ## Same catalog as above (using short options):
     $ astmkcatalog -rdm seg-out.fits

     ## Write the catalog to a text table:
     $ astmkcatalog -rdm seg-out.fits --output=cat.txt

     ## Output columns specified in `columns.conf':
     $ astmkcatalog --config=columns.conf seg-out.fits

     ## Use object and clump labels from a K-band image, but pixel values
     ## from an i-band image.
     $ astmkcatalog K_segmented.fits --hdu=DETECTIONS --clumpscat     \
                    --clumpsfile=K_segmented.fits --clumpshdu=CLUMPS  \
                    --valuesfile=i_band.fits

If MakeCatalog is to do processing (not printing help or option values),
an input labeled image should be provided.  The options described in
this section are those that are particular to MakeProfiles.  For
operations that MakeProfiles shares with other programs (mainly
involving input/output or general processing steps), see *note Common
options::.  Also see *note Common program behavior:: for some general
characteristics of all Gnuastro programs including MakeCatalog.

   The various measurements/columns of MakeCatalog are requested as
options, either on the command-line or in configuration files, see *note
Configuration files::.  The full list of available columns is available
in *note MakeCatalog measurements on each label::.  Depending on the
requested columns, MakeCatalog needs more than one input dataset, for
more details, please see *note MakeCatalog inputs and basic settings::.
The upper-limit measurements in particular need several configuration
options which are thoroughly discussed in *note Upper-limit settings::.
Finally, in *note MakeCatalog output HDUs:: and *note MakeCatalog output
keywords:: the output HDUs and keywords that are written by MakeCatalog
are discussed.

* Menu:

* MakeCatalog inputs and basic settings::  Input files and basic settings.
* Upper-limit settings::        Settings for upper-limit measurements.
* MakeCatalog output HDUs::     Options that specify the output HDUs.
* MakeCatalog output keywords::  Meta data and metameasures and related options.


File: gnuastro.info,  Node: MakeCatalog inputs and basic settings,  Next: Upper-limit settings,  Prev: Invoking astmkcatalog,  Up: Invoking astmkcatalog

7.4.8.1 MakeCatalog inputs and basic settings
.............................................

MakeCatalog works by using a localized/labeled dataset (see *note
MakeCatalog::).  This dataset maps/labels pixels to a specific target
(row number in the final catalog) and is thus the only necessary input
dataset to produce a minimal catalog in any situation.  Because it only
has labels/counters, it must have an integer type (see *note Numeric
data types::), see below if your labels are in a floating point
container.  When the requested measurements only need this dataset (for
example, ‘--geo-x’, ‘--geo-y’, or ‘--geo-area’), MakeCatalog will not
read any more datasets.

   Low-level measurements that only use the labeled image are rarely
sufficient for any high-level science case.  Therefore necessary input
datasets depend on the requested columns in each run.  For example,
let's assume you want the brightness/magnitude and signal-to-noise ratio
of your labeled regions.  For these columns, you will also need to
provide an extra dataset containing values for every pixel of the
labeled input (to measure magnitude) and another for the Sky standard
deviation (to measure error).  All such auxiliary input files have to
have the same size (number of pixels in each dimension) as the input
labeled image.  Their numeric data type is irrelevant (they will be
converted to 32-bit floating point internally).  For the full list of
available measurements, see *note MakeCatalog measurements on each
label::.

   The "values" dataset is used for measurements like
brightness/magnitude, or flux-weighted positions.  If it is a real
image, by default it is assumed to be already Sky-subtracted prior to
running MakeCatalog.  If it is not, you should use the ‘--subtractsky’
option so MakeCatalog reads and subtracts the Sky dataset before any
processing.  To obtain the Sky value, you can use the ‘--sky’ option of
*note Statistics::, but the best recommended method is *note
NoiseChisel::, see *note Sky value::.

   MakeCatalog can also do measurements on sub-structures of detections.
In other words, it can produce two catalogs.  Following the nomenclature
of Segment (see *note Segment::), the main labeled input dataset is
known as "object" labels and the (optional) sub-structure input dataset
is known as "clumps".  If MakeCatalog is run with the ‘--clumpscat’
option, it will also need a labeled image containing clumps, similar to
what Segment produces (see *note Segment output::).  Since clumps are
defined within detected regions (they exist over signal, not noise),
MakeCatalog uses their boundaries to subtract the level of signal under
them.

   There are separate options to explicitly request a file name and
HDU/extension for each of the required input datasets as fully described
below (with the ‘--*file’ format).  When each dataset is in a separate
file, these options are necessary.  However, one great advantage of the
FITS file format (that is heavily used in astronomy) is that it allows
the storage of multiple datasets in one file.  So in most situations
(for example, if you are using the outputs of *note NoiseChisel:: or
*note Segment::), all the necessary input datasets can be in one file.

   When none of the ‘--*file’ options are given (for example
‘--clumpsfile’ or ‘--valuesfile’), MakeCatalog will assume the necessary
input datasets are available as HDUs in the file given as its argument
(without any option).  When the Sky or Sky standard deviation datasets
are necessary and the only ‘--*file’ option called is ‘--valuesfile’,
MakeCatalog will search for these datasets (with the default/given HDUs)
in the file given to ‘--valuesfile’ (before looking into the main
argument file).

   It may happen that your labeled objects image was created with a
program that only outputs floating point files.  However, you know it
only has integer valued pixels that are stored in a floating point
container.  In such cases, you can use Gnuastro's Arithmetic program
(see *note Arithmetic::) to change the numerical data type of the image
(‘float.fits’) to an integer type image (‘int.fits’) with a command like
below:

     $ astarithmetic float.fits int32 --output=int.fits

   To summarize: if the input file to MakeCatalog is the default/full
output of Segment (see *note Segment output::) you do not have to worry
about any of the ‘--*file’ options below.  You can just give Segment's
output file to MakeCatalog as described in *note Invoking
astmkcatalog::.  To feed NoiseChisel's output into MakeCatalog, just
change the labeled dataset's header (with ‘--hdu=DETECTIONS’).  The full
list of input dataset options and general setting options are described
below.

‘-l FITS’
‘--clumpsfile=FITS’
     The FITS file containing the labeled clumps dataset when
     ‘--clumpscat’ is called (see *note MakeCatalog output HDUs::).
     When ‘--clumpscat’ is called, but this option is not, MakeCatalog
     will look into the main input file (given as an argument) for the
     required extension/HDU (value to ‘--clumpshdu’).

‘--clumpshdu=STR’
     The HDU/extension of the clump labels dataset.  Only pixels with
     values above zero will be considered.  The clump labels dataset has
     to be an integer data type (see *note Numeric data types::) and
     only pixels with a value larger than zero will be used.  See *note
     Segment output:: for a description of the expected format.

‘-v FITS’
‘--valuesfile=FITS’
     The file name of the (sky-subtracted) values dataset.  When any of
     the columns need values to associate with the input labels (for
     example, to measure the sum of pixel values or magnitude of a
     galaxy, see *note Brightness flux magnitude::), MakeCatalog will
     look into a "values" for the respective pixel values.  In most
     common processing, this is the actual astronomical image that the
     labels were defined, or detected, over.  The HDU/extension of this
     dataset in the given file can be specified with ‘--valueshdu’.  If
     this option is not called, MakeCatalog will look for the given
     extension in the main input file.

‘--valueshdu=STR/INT’
     The name or number (counting from zero) of the extension containing
     the "values" dataset, see the descriptions above and those in
     ‘--valuesfile’ for more.

‘-s FITS/FLT’
‘--insky=FITS/FLT’
     Sky value as a single number, or the file name containing a dataset
     (different values per pixel or tile).  The Sky dataset is only
     necessary when ‘--subtractsky’ is called or when a column directly
     related to the Sky value is requested (currently ‘--sky’).  This
     dataset may be a tessellation, with one element per tile (see
     ‘--oneelempertile’ of NoiseChisel's *note Processing options::).

     When the Sky dataset is necessary but this option is not called,
     MakeCatalog will assume it is an HDU/extension (specified by
     ‘--skyhdu’) in one of the already given files.  First it will look
     for it in the ‘--valuesfile’ (if it is given) and then the main
     input file (given as an argument).

     By default the values dataset is assumed to be already Sky
     subtracted, so this dataset is not necessary for many of the
     columns.

‘--skyhdu=STR’
     HDU/extension of the Sky dataset, see ‘--skyfile’.

‘--subtractsky’
     Subtract the sky value or dataset from the values file prior to any
     processing.

‘-t STR/FLT’
‘--instd=STR/FLT’
     Sky standard deviation value as a single number, or the file name
     containing a dataset (different values per pixel or tile).  With
     the ‘--variance’ option you can tell MakeCatalog to interpret this
     value/dataset as a variance image, not standard deviation.

     *Important note:* This must only be the SKY standard deviation or
     variance (not including the signal's contribution to the error).
     In other words, the final standard deviation of a pixel depends on
     how much signal there is in it.  MakeCatalog will find the amount
     of signal within each pixel (while subtracting the Sky, if
     ‘--subtractsky’ is called) and account for the extra error due to
     it's value (signal).  Therefore if the input standard deviation (or
     variance) image also contains the contribution of signal to the
     error, then the final error measurements will be over-estimated.

‘--stdhdu=STR’
     The HDU of the Sky value standard deviation image.

‘--variance’
     The value/file given to ‘--instd’ (and ‘--stdhdu’ has the Sky
     variance of every pixel, not the Sky standard deviation.

‘--novalinerror’
     The value/file given to ‘--instd’ is not just due to the sky
     (noise), but also contains the contribution of the signal to each
     pixel's standard deviation or variance.  If this option is given,
     the pixel values will be ignored when measuring the ‘--*-error’
     columns.

‘--forcereadstd’
     Read the input STD image even if it is not required by any of the
     requested columns.  This is because some of the output catalog's
     metadata may need it, for example, to calculate the dataset's
     surface brightness limit (see *note Metameasurements on full
     input::, configured with ‘--sbl-area’ and ‘--sbl-sigma’ in *note
     MakeCatalog output keywords::).

     Furthermore, if the input STD image does not have the ‘MEDSTD’
     keyword (that is meant to contain the representative standard
     deviation of the full image), with this option, the median will be
     calculated and used for the surface brightness limit.

‘-z FLT’
‘--zeropoint=FLT’
     The zero point magnitude for the input image, see *note Brightness
     flux magnitude::.

‘--sigmaclip FLT,FLT’
     The sigma-clipping parameters when any of the sigma-clipping
     related columns are requested (for example, ‘--sigclip-median’ or
     ‘--sigclip-number’).

     This option takes two values: the first is the multiple of
     $\sigma$, and the second is the termination criteria.  If the
     latter is larger than 1, it is read as an integer number and will
     be the number of times to clip.  If it is smaller than 1, it is
     interpreted as the tolerance level to stop clipping.  See *note
     Sigma clipping:: for a complete explanation.

‘--frac-max=FLT[,FLT]’
     The fractions (one or two) of maximum value in objects or clumps to
     be used in the related columns, for example, ‘--frac-max1-area’,
     ‘--frac-max1-sum’ or ‘--frac-max1-radius’.  For the maximum value,
     see the description of ‘--maximum’ column below.  The value(s) of
     this option must be larger than 0 and smaller than 1 (they are a
     fraction).  When only ‘--frac-max1-area’ or ‘--frac-max1-sum’ is
     requested, one value must be given to this option, but if
     ‘--frac-max2-area’ or ‘--frac-max2-sum’ are also requested, two
     values must be given to this option.  The values can be written as
     simple floating point numbers, or as fractions, for example,
     ‘0.25,0.75’ and ‘0.25,3/4’ are the same.

‘--spatialresolution=FLT’
     The error in measuring spatial properties (for example, the area)
     in units of pixels.  You can think of this as the FWHM of the
     dataset's PSF and is used in measurements like the error in surface
     brightness (‘--sb-error’ described in *note Surface brightness
     measurements::).

     But because most objects in astronomy are either concentrated (like
     stars or galaxies) or very large (like interstellar dust filaments,
     also known as Galactic cirrus), our tests
     (https://savannah.gnu.org/task/?16586) showed that this factor
     should be set to zero to give a reasonable error measurement.

     In summary, this is because in both situations the error in the
     number of pixels used in the measurement only effects the few
     pixels in the circumference of the label and they usually have the
     lowest values (contribute least to the total sum or average).
     Therefore, the default value to this option is 0.  But in case your
     targets are very small (a handful of pixels wide) and their light
     profiles are flat (for example a cosmic ray!), it may be necessary
     to increase the value to this option.

‘--inbetweenints’
     Output will contain one row for all integers between 1 and the
     largest label in the input (irrespective of their existance in the
     input image).  By default, MakeCatalog's output will only contain
     rows with integers that actually corresponded to at least one pixel
     in the input dataset.

     For example, if the input's only labeled pixel values are 11 and
     13, MakeCatalog's default output will only have two rows.  If you
     use this option, it will have 13 rows and all the columns
     corresponding to integer identifiers that did not correspond to any
     pixel will be 0 or NaN (depending on context).


File: gnuastro.info,  Node: Upper-limit settings,  Next: MakeCatalog output HDUs,  Prev: MakeCatalog inputs and basic settings,  Up: Invoking astmkcatalog

7.4.8.2 Upper-limit settings
............................

The upper-limit magnitude was discussed in *note Metameasurements on
full input::.  Unlike other measured values/columns in MakeCatalog, the
upper limit magnitude needs several extra parameters which are discussed
here.  All the options specific to the upper-limit measurements start
with ‘up’ for "upper-limit".  The only exception is ‘--envseed’ that is
also present in other programs and is general for any job requiring
random number generation in Gnuastro (see *note Generating random
numbers::).

   One very important consideration in Gnuastro is reproducibility.
Therefore, the values to all of these parameters along with others (like
the random number generator type and seed) are also reported in the
comments of the final catalog when the upper limit magnitude column is
desired.  The random seed that is used to define the random positions
for each object or clump is unique and set based on the (optionally)
given seed, the total number of objects and clumps and also the labels
of the clumps and objects.  So with identical inputs, an identical
upper-limit magnitude will be found.  However, even if the seed is
identical, when the ordering of the object/clump labels differs between
different runs, the result of upper-limit measurements will not be
identical.

   MakeCatalog will randomly place the object/clump footprint over the
dataset.  When the randomly placed footprint does not fall on any object
or masked region (see ‘--upmaskfile’) it will be used in the final
distribution.  Otherwise that particular random position will be ignored
and another random position will be generated.  Finally, when the
distribution has the desired number of successfully measured random
samples (‘--upnum’) the distribution's properties will be measured and
placed in the catalog.

   When the profile is very large or the image is significantly covered
by detections, it might not be possible to find the desired number of
samplings in a reasonable time.  MakeProfiles will continue searching
until it is unable to find a successful position (since the last
successful measurement(1)), for a large multiple of ‘--upnum’
(currently(2) this is 10).  If ‘--upnum’ successful samples cannot be
found until this limit is reached, MakeCatalog will set the upper-limit
magnitude for that object to NaN (blank).

   MakeCatalog will also print a warning if the range of positions
available for the labeled region is smaller than double the size of the
region.  In such cases, the limited range of random positions can
artificially decrease the standard deviation of the final distribution.
If your dataset can allow it (it is large enough), it is recommended to
use a larger range if you see such warnings.

‘--upmaskfile=FITS’
     File name of mask image to use for upper-limit calculation.  In
     some cases (especially when doing matched photometry), the object
     labels specified in the main input and mask image might not be
     adequate.  In other words they do not necessarily have to cover
     _all_ detected objects: the user might have selected only a few of
     the objects in their labeled image.  This option can be used to
     ignore regions in the image in these situations when estimating the
     upper-limit magnitude.  All the non-zero pixels of the image
     specified by this option (in the ‘--upmaskhdu’ extension) will be
     ignored in the upper-limit magnitude measurements.

     For example, when you are using labels from another image, you can
     give NoiseChisel's objects image output for this image as the value
     to this option.  In this way, you can be sure that regions with
     data do not harm your distribution.  See *note Metameasurements on
     full input:: for more on the upper limit magnitude.

‘--upmaskhdu=STR’
     The extension in the file specified by ‘--upmask’.

‘--upnum=INT’
     The number of random samples to take for all the objects.  A larger
     value to this option will give a more accurate result
     (asymptotically), but it will also slow down the process.  When a
     randomly positioned sample overlaps with a detected/masked pixel it
     is not counted and another random position is found until the
     object completely lies over an undetected region.  So you can be
     sure that for each object, this many samples over undetected
     objects are made.  See the upper limit magnitude discussion in
     *note Metameasurements on full input:: for more.

‘--uprange=INT,INT’
     The range/width of the region (in pixels) to do random sampling
     along each dimension of the input image around each object's
     position.  This is not a mandatory option and if not given (or
     given a value of zero in a dimension), the full possible range of
     the dataset along that dimension will be used.  This is useful when
     the noise properties of the dataset vary gradually.  In such cases,
     using the full range of the input dataset is going to bias the
     result.  However, note that decreasing the range of available
     positions too much will also artificially decrease the standard
     deviation of the final distribution (and thus bias the upper-limit
     measurement).

‘--envseed’
     Read the random number generator type and seed value from the
     environment (see *note Generating random numbers::).  Random
     numbers are used in calculating the random positions of different
     samples of each object.

‘--upsigmaclip=FLT,FLT’
     The raw distribution of random values will not be used to find the
     upper-limit magnitude, it will first be $\sigma$-clipped (see *note
     Sigma clipping::) to avoid outliers in the distribution (mainly the
     faint undetected wings of bright/large objects in the image).  This
     option takes two values: the first is the multiple of $\sigma$, and
     the second is the termination criteria.  If the latter is larger
     than 1, it is read as an integer number and will be the number of
     times to clip.  If it is smaller than 1, it is interpreted as the
     tolerance level to stop clipping.  See *note Sigma clipping:: for a
     complete explanation.

‘--upnsigma=FLT’
     The multiple of the final ($\sigma$-clipped) standard deviation (or
     $\sigma$) used to measure the upper-limit sum or magnitude.

‘--checkuplim=INT[,INT]’
     Write a table of positions and measured values for the full random
     distribution used in the upper-limit positions for one particular
     object or clump.  The table will be placed as a HDU in the output
     file with a name of ‘CHECK-UPPERLIMIT’.  If only one integer is
     given to this option, it is interpreted to be an object's label.
     If two values are given, the first is the object label and the
     second is the ID of requested clump within it.

     The output is a table with three columns on a 2D image input (and
     four columns on a 3D cube input).  The first two columns are the
     pixel X,Y positions of the center of each label's tile (see next
     paragraph), in each random sampling of this particular object/clump
     (on a 3D cube, the third column will be the position along the
     third dimension).  The final column is the measured flux over that
     region.

     If a randomly placed region overlapped with a detection or masked
     pixel, its measured value will be a NaN (not-a-number).  The total
     number of rows is thus unknown before running.  However, if an
     upper-limit measurement was made in the main output of MakeCatalog,
     you can be sure that the number of rows with non-NaN measurements
     is the number given to the ‘--upnum’ option.

     The "tile" of each label is defined by the minimum and maximum
     positions of each label: values of the ‘--min-x’, ‘--max-x’,
     ‘--min-y’ and ‘--max-y’ columns in the main output table for each
     label.  Therefore, the tile center position that is recorded in the
     output of this column ignores the distribution of labeled pixels
     within the tile.  Note that this is only about the center position,
     not the measurement.

     Precise interpretation of the position is only relevant when the
     footprint of your label is highly un-symmetrical and you want to
     use this catalog to insert your object into the image.  In such a
     case, you can also ask for ‘--min-x’ and ‘--min-y’ and manually
     calculate their difference with the following two positional
     measurements of your desired label: ‘--geo-x’ and ‘--geo-y’ (which
     report the label's "geometric" center; only using the label
     positions ignoring any "values") or ‘--x’ and ‘--y’ (which report
     the value-weighted center of the label).  Adding the difference
     with the position reported by this column, will let you define
     alternative "center"s for your label in particular situations (this
     will usually not be necessary!).  For more on these positional
     columns, see *note Position measurements in pixels::.

   ---------- Footnotes ----------

   (1) The counting of failed positions restarts on every successful
measurement.

   (2) In Gnuastro's source, this constant number is defined as the
‘MKCATALOG_UPPERLIMIT_MAXFAILS_MULTIP’ macro in ‘bin/mkcatalog/main.h’,
see *note Downloading the source::.


File: gnuastro.info,  Node: MakeCatalog output HDUs,  Next: MakeCatalog output keywords,  Prev: Upper-limit settings,  Up: Invoking astmkcatalog

7.4.8.3 MakeCatalog output HDUs
...............................

After it has completed all the requested measurements (see *note
MakeCatalog measurements on each label::), MakeCatalog writes them in
table(s) that are stored into a FITS file.  This is necessary because
the number of output tables and the metadata that MakeCatalogs provides
can be numerous.  If any of the output tables are necessary in another
format (for example plain-text), you can use Gnuastro's Table program
(with executable name ‘asttable’, see *note Table::).

   This section focuses on the HDUs of the output, for the keywords, see
*note MakeCatalog output keywords::.  The name of the output FITS table
can be given to the ‘--output’ option, with a recognized FITS suffix (as
defined in *note Arguments::).  When it is not given, the input name
will be appended with a ‘-cat.fits’ suffix (see *note Automatic
output::) and its format (ASCII or Binary FITS table) will be determined
from the ‘--tableformat’ option, which is also discussed in *note Input
output options::.  By default (when ‘--spectrum’, ‘--checkuplim’ or
‘--clumpscat’ are not called) only a single catalog/table will be
created for the labeled objects.

   • When ‘--clumpscat’ is called, a secondary catalog/table HDU will
     also be created for "clumps" (one of the outputs of the Segment
     program, for more on "objects" and "clumps", see *note Segment::).
     In short, if you only have one labeled image, you do not have to
     worry about clumps and just ignore this.
   • When ‘--checkuplim’ is called, a HDU is added to the output FITS
     file and is fully described in the description of this option in
     *note Upper-limit settings::.

   The full list of MakeCatalog's options relating to the output file
format and keywords are listed below.  See *note MakeCatalog
measurements on each label:: for specifying which columns you want in
the final catalog.

‘-C’
‘--clumpscat’
     Do measurements on clumps and produce a second catalog (only
     devoted to clumps).  When this option is given, MakeCatalog will
     also look for a secondary labeled dataset (identifying
     substructure) and produce a catalog from that.  For more on the
     definition on "clumps", see *note Segment::.

     When the output is a FITS file, the objects and clumps
     catalogs/tables will be stored as multiple extensions of one FITS
     file.  You can use *note Table:: to inspect the column meta-data
     and contents in this case.  However, in plain text format (see
     *note Gnuastro text table format::), it is only possible to keep
     one table per file.  Therefore, if the output is a text file, two
     output files will be created, ending in ‘_o.txt’ (for objects) and
     ‘_c.txt’ (for clumps).

‘--noclumpsort’
     Do not sort the clumps catalog based on object ID (only relevant
     with ‘--clumpscat’).  This option will benefit the performance(1)
     of MakeCatalog when it is run on multiple threads _and_ the
     position of the rows in the clumps catalog is irrelevant (for
     example, you just want the number-counts).

     MakeCatalog does all its measurements on each _object_
     independently and in parallel.  As a result, while it is writing
     the measurements on each object's clumps, it does not know how many
     clumps there were in previous objects.  Each thread will just fetch
     the first available row and write the information of clumps (in
     order) starting from that row.  After all the measurements are
     done, by default (when this option is not called), MakeCatalog will
     reorder/permute the clumps catalog to have both the object and
     clump ID in an ascending order.

     If you would like to order the catalog later (when it is a plain
     text file), you can run the following command to sort the rows by
     object ID (and clump ID within each object), assuming they are
     respectively the first and second columns:

          $ awk '!/^#/' out_c.txt | sort -g -k1,1 -k2,2

‘--cnl-check’
     Add an extra HDU to the output of MakeCatalog that contains the
     positions of every clump as well as the positions of that clump's
     nearest clump and the distance between the two (in both units of
     arcseconds and pixels).

   ---------- Footnotes ----------

   (1) The performance boost due to ‘--noclumpsort’ can only be felt
when there are a huge number of objects.  Therefore, by default the
output is sorted to avoid miss-understandings or bugs in the user's
scripts when the user forgets to sort the outputs.


File: gnuastro.info,  Node: MakeCatalog output keywords,  Prev: MakeCatalog output HDUs,  Up: Invoking astmkcatalog

7.4.8.4 MakeCatalog output keywords
...................................

The columns and rows that include the various measurements for the
variuos labels of the input in its HDUs were described separately in
*note MakeCatalog output HDUs::.  But those raw numbers are not the only
thing that MakeCatalog writes in its output!  MakeCatalog will also
write metadata (header keywords) in the 0th (first) HDU of the output
FITS file which add a lot of value and help to interpret the raw
numbers.  The only keywords in the other HDUs are the column names,
units and comments; generic metadata are written in the 0th HDU.

   You can see the full list of keywords written by MakeCatalog in its
output with Gnuastro's Fits program, for example ‘astfits out.fits -h0’.
If you only want the value of certain keywords (in a script/pipeline for
example), its ‘--keyvalue’ option is pretty convenient see *note Keyword
inspection and manipulation::.

   The keywords are grouped in the output based on context with a title
above each group.  We'll follow the same structure here, skipping the
first three groups (that are generic to all Gnuastro's programs, see
*note Output FITS files::).  Some of the values reported below will be
repeated (are the same as the respective option keyword in the "Option
values" group of keywords).  The "Option value" group contains the raw
option names and values: the option (keyword) names are just written in
full-caps according to the FITS standard.  This is done for the
following reason: Human readability is important for option names, so
they tend to be long and force a ‘HIERARCH’ string at the start of the
line.  On the contrary, for automatic extraction and human readability,
it helps a lot to have all keywords related to a certain metameasure
start with the same characters and because the keyword description is
written just after it, there is no problem if the name is cryptic.

Input file(s) and HDUs
     ‘INLAB’ and ‘INLABHDU’: main (object) label dataset.
     ‘INCLU’ and ‘INCLUHDU’: clump label dataset.
     ‘INVAL’ and ‘INVALHDU’: value dataset.
     ‘INSTD’ and ‘INSTDHDU’: standard deviation dataset.
     ‘INVAR’ and ‘INVARHDU’: variance dataset.
     ‘INUPM’ and ‘INUPMHDU’: upper-limit mask dataset.
          The file name and HDU of all the possible inputs to
          MakeCatalog.  Only the keywords that correspond to inputs
          which are actually used in a given run of MakeCatalog will be
          written in the output.  This is based on the columns or
          metameasurements that you request, not if you gave them on the
          command-line (described in *note MakeCatalog inputs and basic
          settings::), or the HDUs that exist in the main argument.

          For the standard deviation or variance, in case a single
          number was given (instead of a dataset), that number will be
          written for ‘INSTD’ or ‘INVAR’ (without any ‘INSTDHDU’ and
          ‘INVALHDU’).

Input pixel grid and value properties
     Basic information about the pixel values and grid properties of the
     input:

     ‘PIXWIDTH’
          The width of one pixel on the sky (in units of arcseconds).

     ‘PIXAREA’
          The area of one pixel on the sky (in units of arcseconds
          squared) on the reference point.  The difference in pixel area
          across the image will be negligible in most science images.
          In case you would like to check this for your input images,
          use the ‘--pixelareaonwcs’ option of the Fits program (see
          *note Pixel information images::).

     ‘ZEROPNT’
          The zero point of the values image (used to convert This is
          the same value you gave to the ‘--zeropoint’ option of *note
          MakeCatalog inputs and basic settings::.

     ‘STDUSED’
          Per-pixel standard deviation (used in noise-based
          metameasurements like *note Surface brightness limit of
          image:: and *note Noise based magnitude limit of image::).
          This keyword will only be present when a standard deviation
          image has been loaded (done automatically for any column
          measurement that involves noise, for example, ‘--sn’,
          ‘--magnitude-error’ or ‘--sky-std’).  In case your catalog
          does not include any such columns and you want this keyword,
          you can use the ‘--meta-measures’ option (see *note
          Metameasurements on full input::).

          If the ‘MEDSTD’ keyword is present in the standard deviation
          dataset (see *note NoiseChisel output::), it will be used.
          Otherwise, the median of the standard deviation input is
          calculated, used for the metameasures and written in this
          keyword.

Upper-limit parameters
     When any of the upper-limit measurements are requested, the input
     parameters for the upper-limit measurement are stored in the
     following keywords (see *note Upper limit measurements::).
     ‘UPSIGMA’
          The multiple of sigma to measure the upper-limit.  This is the
          same value given to the ‘--upnsigma’ option of *note
          Upper-limit settings::.

     ‘UPNUMBER’
          The number of random positions with a successful reading.
          This is the same value given to the ‘--upnum’ option of *note
          Upper-limit settings::.

     ‘UPRNGNAM’
          Name of the random number generator used for finding the
          random positions; see *note Generating random numbers::.

     ‘UPRNGSEE’
          Seed used for the random number generator.  This will be
          differnet on every run, unless ‘--envseed’ is called.  For
          more details, see *note Generating random numbers::

     ‘UPSCMLTP’
          $\sigma$-clipping paramter: multiple of sigma.  Clipping is
          necessary to reject strong outliers that can affect the
          statistics.  This is the first value given to the
          ‘--upsigmaclip’ option of *note Upper-limit settings::.

     ‘UPSCTOL’
          $\sigma$-clipping paramter: tolerance level.  Clipping is
          necessary to reject strong outliers that can affect the
          statistics.  This is the second value given to the
          ‘--upsigmaclip’ option of *note Upper-limit settings::.

Noise-based metameasures
     The following metamesurement are calculated purely based on the
     measured noise level.  But they are not written by default, if you
     want them run with ‘--meta-measure’.

     ‘SBL’
          Measured surface brighntess limit (in units of
          mag/arcsec$^2$); as described in *note Surface brightness
          limit of image::.
     ‘NML’
          Measured surface brighntess limit (in units of
          mag/arcsec$^2$); as described in *note Noise based magnitude
          limit of image::.

Confusion limit
     The confusion limit (CNL) is a measure of density of resolved
     sources in the input image.  For a complete review on its goals and
     how to interpret the values to these keywords, see *note Confusion
     limit of image::.  In particular, if you would like to see the full
     distribution of nearest neighbors and their distances, you can use
     ‘--cnl-check’ as described in *note MakeCatalog output HDUs::.  By
     default, the clumps catalog is used for the distribution, however
     if you would like to use the objects for any reason, you can use
     ‘--cnl-with-objects’ as described below in this section (not
     recommended unless you understand the risks).

     ‘CNL’
          The confusion limit: difference between ‘CNLP75’ and ‘CNLP25’.

     ‘CNLP05’
     ‘CNLP25’
     ‘CNLP50’
     ‘CNLP75’
     ‘CNLP95’
          Various percentiles of the distribution of distances to the
          nearest neighbor in units of pixels.  The ‘P’ before each
          number is for percentile, so ‘CNLP50’ shows the median value
          of the distribution to nearest neighbors.  Note that a
          quantile is just the percentile after division by 100.
          Percentiles are used in the keyword names because they are
          simple integers and do not need a floating point.

The following MakeCatalog options are specifically related to the
various keywords above.

‘--sbl-sigma=FLT’
     Value to multiply with the median standard deviation (from a
     ‘MEDSTD’ keyword in the Sky standard deviation image) for
     estimating the surface brightness limit.  Note that the surface
     brightness limit is only reported when a standard deviation image
     is read, in other words a column using it is requested (for
     example, ‘--sn’) or ‘--forcereadstd’ is called.

     This value is a per-pixel value, not per object/clump and is not
     found over an area or aperture, like the common $5\sigma$ values
     that are commonly reported as a measure of depth or the upper-limit
     measurements (see *note Metameasurements on full input::).

‘--sbl-area=FLT’
     Area (in arc-seconds squared) to convert the per-pixel estimation
     of ‘--sbl-sigma’ in the comments section of the output tables.
     Note that the surface brightness limit is only reported when a
     standard deviation image is read, in other words a column using it
     is requested (for example, ‘--sn’) or ‘--forcereadstd’ is called.

     Note that this is just a unit conversion using the World Coordinate
     System (WCS) information in the input's header.  It does not
     actually do any measurements on this area.  For random measurements
     on any area, please use the upper-limit columns of MakeCatalog (see
     the discussion on upper-limit measurements in *note
     Metameasurements on full input::).

‘--cnl-with-objects’
     Use the object positions instead of clumps for measuing the
     distance to the nearest label.  This is only useful if you have
     generated your labels image with something other than Gnuastro's
     *note Segment::, and that program doesn't have the capacity to
     identify individual peaks and extended signal at the same time.
     Therefore, in case you have generated the input labels with
     Segment, we do not recommend using this option.


File: gnuastro.info,  Node: Match,  Prev: MakeCatalog,  Up: Data analysis

7.5 Match
=========

High-level data (catalogs) of a single astronomical source can come come
from different telescopes, filters, software (detection, segmentation
and cataloging tools on the same image) and even different
configurations for a single software.  As a result, one of the first
things we usually do after generating or querying catalog data (for
example, with *note MakeCatalog:: or *note Query::), is to find which
sources in one catalog correspond to which in the other(s).  In other
words, to 'match' the two catalogs with each other.  Within Gnuastro,
the Match program is in charge of such operations.  The matching rows
are defined within an aperture, which can be a circle or an ellipse with
any orientation.

   Before digging into the usage and command-line execution details (in
*note Invoking astmatch::), it is important to discuss several concepts
in detail.  If this is the first time you are using Match, please take
the time to read these three sections (before *note Invoking astmatch::)
in detail to optimize your usage of Match and understanding its outputs.
The first is *note Arranging match output:: which introduces how the
matched rows should be arranged in the outputs and is determined by the
purpose/context of your analysis (there are examples in that section).
When you want an "inner" or "full" match, some matches can become
ambiguous; in *note Unambiguous matching:: Gnuastro's Match strategy for
those is discussed.  We will then review the *note Matching
algorithms::.

* Menu:

* Arranging match output::      Various meanings of matching.
* Unambiguous matching::        How we avoid multiple matches.
* Matching algorithms::         Different ways to find the match.
* Invoking astmatch::           Inputs, outputs and options of Match.


File: gnuastro.info,  Node: Arranging match output,  Next: Unambiguous matching,  Prev: Match,  Up: Match

7.5.1 Arranging match output
----------------------------

Once the unambiguous match (see *note Unambiguous matching::) between
two catalogs is found the thing you want to do with the matched rows can
be different based on context.  For example, sometimes you exclusively
want the matched rows of the two input catalogs.  In other times you
want to find the nearest row in a reference catalog for every row of a
query catalog.  Yet other times, you want to merge two catalogs (having
all the matched and non-matched rows in one table).  Within Match, these
are defined as different "arrangements" of the output.

   The core matching algorithm (to find the nearest points between two
catalogs *note Matching algorithms::) is the same in all the
arrangements.  The thing that is different is how you want to use that
information in constructing/arranging the output.  The various
arrangements (that should be given to the ‘--arrange’ option) are
formally defined below.

‘inner’
     The output of an inner match will contain umambiguously/exclusively
     matched rows in both catalogs within a given aperture.  In other
     words, the number of rows in the output will be fewer (or equal) to
     the input with fewest rows, _and_ that no two rows of any input can
     be repeated.  Another way to say this is that the order of the
     inputs does not matter for the final number of matched rows in an
     inner match.  An important aspect of an inner match is the
     necessity of an aperture to define the maximum/largest acceptable
     distance to have a match.  Without an aperture to define a
     reasonable match, there is always bound to be a "nearest" item in
     both catalogs!

     A real-world example (in astrophysics) of inner matching is when
     you have two catalogs of galaxies from different filters/telescopes
     and want to find the same galaxy in both catalogs.  The simplified
     example below shows inner matching in action.  The values in each
     row are intentionally abstract/ideal to help understanding the
     concept (for example, the first two columns can be the RA and Dec
     in a real-world catalog).  First, let's look at the contents of the
     two tables we want to match:

          $ cat a.txt
          # Column 1: X    [pix,  f32]   X axis position
          # Column 2: Y    [pix,  f32]   Y axis position
          # Column 3: MAGa [unit, f32] Electrong counts
          1.2   1.2   18.2
          3.2   3.2   20.8
          8.2   8.2   15.2

          $ cat b.txt
          # Column 1: X    [pix, f32] X axis position
          # Column 2: Y    [pix, f32] Y axis position
          # Column 3: MAGb [pix, f32] Radius of object in pixels
          1   1   18.5
          2   2   23.1
          3   3   16.2
          4   4   22.7
          5   5   23.4

     In this hypothetical example, the first catalog has the magnitude
     of three objects in filter a, and the second catalog has the
     magnitude of 5 objects in filter b.  But the coordinates are not
     identical (which is natural when the catalogs come from different
     source) and not all of the galaxies in one catalog match the other!

     With the first command below, we will find the matching rows with
     an inner arrangement, and put the radius and magnitude of the
     matched rows in the same row (using the ‘--outcols’, see *note
     Invoking astmatch::).  The second command simply prints the
     contents of the table on the command-line in easy-to-read format
     (simplified here by removing trailing zeros to help visual
     comparison with above).

          $ astmatch a.txt --ccol1=X,Y b.txt --ccol2=X,Y \
                     --aperture=0.5 --outcols=a1,b1,a2,b2,a3,b3 \
                     --arrange=inner --output=inner.fits

          $ asttable inner.fits -Y
          1.2    1    1.2    1    18.2    18.5
          3.2    3    3.2    3    20.8    16.2

     As you see, while both inputs had more than two rows, the output of
     an inner arrangement only contains the rows that matched in both
     catalogs (two in this case).  For inner matching, we rarely need
     both coordinate columns; instead, we can simply only select the
     coordinates of the catalog that had a better precision (in other
     words, a better point spread function, or PSF). For example,
     assuming that ‘a.fits’ had better precision, we could simplify the
     output above with ‘--outcols=a1,a2,a3,b3’ (which has removed ‘b1’
     and ‘b2’ compared to the one in the example).  Try this out for
     yourself as an exercise to visually understand this important
     feature of inner arrangements.

‘full’
     The output of a full match will contain matching and non-matching
     rows in both catalogs.  In other words, when there is a match
     between the two catalogs, it is unambiguous and exclusive (as in
     the "Inner" mode above).  However, if a row from any of the two
     inputs does not match, it is still present in the output catalog,
     but rows that belong to the other input are given blank values for
     it.  When there is no matching rows between the two tables, the
     output will effectively be a concatenation of the two inputs.  In
     SQL jargon (https://www.w3schools.com/sql/sql_join.asp) this is
     known as "full-outer joining".

     A "full" output arrangement is therefore effectively an inner
     arrangement, but with the non-matching rows from both catalogs
     appended to it.  Therefore, similar to an inner match, an aperture
     is mandatory.

     Let's use the same two example demo inputs of the inner arrangement
     above, but call the Match program with a ‘--arrange=full’ this
     time:

          $ astmatch a.txt --ccol1=X,Y b.txt --ccol2=X,Y \
                     --aperture=0.5 --outcols=a1,b1,a2,b2,a3,b3 \
                     --arrange=full --output=full-raw.fits

          $ asttable full-raw.fits -Y
          1.2       1       1.2     1       18.2      18.5
          3.2       3       3.2     3       20.8      16.2
          8.2       nan     8.2     nan     15.2      nan
          nan       2       nan     2       nan       23.1
          nan       4       nan     4       nan       22.7
          nan       5       nan     5       nan       23.4

     You see that first two rows are the same as the inner example
     above.  But now we have all the non-matching rows of both catalogs
     also.

     In a real world example, the problem with these non-matching rows
     is that the coordinates are NaN/blank!  To fix that, you can use
     the ‘where’ operator and *note Column arithmetic:: like the command
     below (where the coordinates of ‘b.fits’ are put in the first and
     third columns and the coordinate columns of ‘b.fits’ are removed).
     With the last two ‘--colmetadata’ options, we are giving names,
     units and comments to the newly synthesized columns.  This is
     because after column arithmetic, the metadata is lost by default
     and it is up to you to set it: always do so (data without metadata
     is not too useful!):

          $ asttable full-raw.fits --output=full.fits \
                     -c'arith $1 set-i i i isblank $2 where' \
                     -c'arith $3 set-i i i isblank $4 where' \
                     -cMAGa,MAGb \
                     --colmetadata=1,X,pix,"X axis position (merged)" \
                     --colmetadata=2,Y,pix,"Y axis position (merged)"

          $ asttable full.fits -Y
          1.2     1.2     18.2     18.5
          3.2     3.2     20.8     16.2
          8.2     8.2     15.2     nan
          2.0     2.0     nan      23.1
          4.0     4.0     nan      22.7
          5.0     5.0     nan      23.4

     The ‘full.fits’ table above now has the matching rows; but it also
     has entries for the non-matching objects of both catalogs and only
     one X and Y position.  To find objects that were only detected in
     one or the other filter, the user can simply use the ‘--noblank’
     option like the first example below.  To get the inner match
     output, the user of the catalog can simply use the same option with
     ‘_all’, like the second example below.

          $ asttable full.fits --noblank=MAGa
          $ asttable full.fits --noblank=_all

‘outer’
     An outer row arrangement will produce a table with the same number
     of rows as the second input to the Match program: every row of the
     second input is matched with the nearest row of the first.  An
     outer match is useful when you want to find the nearest row of a
     _reference_ catalog (first input) for each row of a _query_
     catalog's entries (second input).  In this scenario, if there are
     multiple entries of the reference catalog in the output, there is
     no problem (in other words, the match is not exclusive).  Also,
     since you simply want the nearest entry, no aperture is necessary.

     In SQL jargon (https://www.w3schools.com/sql/sql_join.asp), there
     is two types of outer arrangements (or "joining" as SQL calls it):
     left-outer and right-outer.  But in Gnuastro's Match program, we
     current only have the right-outer join.  This is because unlike SQL
     (that does many more things than matching/joining) Gnuastro's Match
     program is not a programming language!  So the user can always
     change the order of their inputs to achieve left-outer matching!
     In other words, within the scope of Gnuastro's Match, there is no
     need for the extra complexity.

     For example, let's assume you want to find the Galactic extinction
     for the galaxies in your catalog.  In this scenario, the reference
     catalog is a table of Galactic extinction values for a grid of
     positions in the sky and the check catalog contains one row for
     each of your target galaxies.  If two or more galaxies are near the
     same entry in the reference catalog, you expect them to have the
     same extinction column value.

     To see it in practice, let's use the ‘b.fits’ of the inner
     arrangement example above as our query catalog.  For the reference
     catalog, let's make a new ‘c.fits’; assuming that it has galactic
     extinction for some points you have already measured:

          $ cat c.txt
          # Column 1: X    [pix, f32] X axis position
          # Column 2: Y    [pix, f32] Y axis position
          # Column 3: EXT  [mag, f32] Galactic extinction at (X,Y)
          1  1  0.04
          4  4  0.05

          $ astmatch c.txt --ccol1=X,Y b.txt --ccol2=X,Y \
                     --outcols=bX,bY,bMAGb,aEXT --arrange=outer \
                     --kdtree=internal --output=outer.fits

          $ asttable outer.fits -Y
          1   1   18.5   0.04
          2   2   23.1   0.04
          3   3   16.2   0.05
          4   4   22.7   0.05
          5   5   23.4   0.05

     The output has the same number of rows as the query catalog, and
     rows from the reference catalog are repeated when they are the
     nearest to each query item.

‘outer-within-aperture’
     Similar to the outer match described above, with the difference
     that if the nearest point is farther than ‘--aperture’, all
     reference table (first input) columns will be NaN in the output.
     For more on ‘--aperture’, see *note Invoking astmatch::.  This is
     useful in scenarios where the distance matters and you do not want
     reference points that are too distant from the query catalog.

     As an example, let's repeat the command from the example for the
     outer arrangement, but use this arrangement instead:

          $ astmatch c.txt --ccol1=X,Y b.txt --ccol2=X,Y \
                     --aperture=0.5 --outcols=bX,bY,bMAGb,aEXT \
                     --arrange=outer-within-aperture --kdtree=internal \
                     --output=outer-wa.fits

          $ asttable outer-wa.fits -Y
          1   1   18.5   0.040
          2   2   23.1   nan
          3   3   16.2   nan
          4   4   22.7   0.050
          5   5   23.4   nan

     You can confirm that only the reference rows that were within the
     given aperture, are not NaN/blank.


File: gnuastro.info,  Node: Unambiguous matching,  Next: Matching algorithms,  Prev: Arranging match output,  Up: Match

7.5.2 Unambiguous matching
--------------------------

When you want an "inner" or "full" match (discussed previously in *note
Arranging match output::), it is important that the matches are
unambiguous.  Within Gnuastro's Match program, an ambiguous match is
defined when there is more than one "match" within the given aperture or
the output depends on the order of the inputs.  In other words, to be
unambiguous, there should only a single row/entry in both catalogs
within the given aperture and the match should not depend on the order
of the inputs (which one is given first on the command-line).  If there
is more than one match in either catalog to the other, Gnuastro's Match
program will remove the match and put all the entries of the second
catalog that have multiple matches in a HDU of the output FITS file
called ‘FLAGGED-IN-2ND’.  These are necessary for the match output to be
reliable: simply reporting the nearest match when there is more than one
match within the given aperture will lead to an ambiguous result.

   The issue is best displayed in practice.  Let's assume you have two
catalogs called ‘x.fits’ and ‘o.fits’ and the points are distributed as
follows (in four "clusters"):

     $ asttable x.fits
     -2    -1
      1     4
      7     0
      8     4
      8    -1
      2     1
      9     0

     $ asttable o.fits
     -2    -2
      2     3
      4     1
      7     4
      8     0
      9     4

     ## ASCII visualization of the points in the two tables
     |
     |            x                o  x  o
     |              o
     |
     |              x   o
     |                             x  o  x
     |   x                            x
     |   o
     --------------------------------------->

   When you run Gnuastro's Match program on these two tables (and your
aperture is large enough to fit the various "clusters" of points) only
the first two points of the two catalogs (on the bottom-left side of the
chart) will match.  The rest will be flagged and the row-counter of the
second input for those matches will be saved in the ‘FLAGGED-IN-2ND’ HDU
of the output.  This is very important, especially when you consider the
matching aperture to be analogous to your error in the two positions.
In other words, while the other points do indeed have a closest match,
they have other points that are also within the error/aperture of the
match and therefore they will contaminate/pollute any higher-level
analysis you want to do with the match.

   You can check if your match command produced flagged/ambiguous
matches:
   • On the lines that Match prints on the terminal, the number of
     flagged sources will be printed just under the number of
     unambiguous matches.  But this will require human checking (which
     is not easy in large pipelines!)  and is primarily good when you
     are developing/debugging your higher-level scripts or just
     understanding your data.
   • Read the number of rows of the ‘FLAGGED-IN-2ND’ HDU of the output
     and put that in a conditional, like the minimal working example
     below (you can replace the ‘echo GOOD’ and ‘echo BAD’ parts with
     any series of complex operations):

          $ nflag=$(asttable red-first.fits -hFLAGGED-IN-2ND \
                             --info-num-rows)
          $ if [ $nflag = 0 ]; then echo GOOD; else echo BAD; fi

   As mentioned above, the ‘FLAGGED-IN-2ND’ HDU will only contain the
row counter (starting from 1) of the ambiguous matches of the second
input.  To get the full row of the second input from that, you can use
the commands below.  It uses Table's *note Column arithmetic:: and a
second call on Match with the column counters (a one-dimensional match):

     $ asttable second.fits --output=2nd-with-counter.fits
                -c'arith $1 counter swap free',_all
     $ astmatch original-output.fits  --hdu=FLAGGED-IN-2ND --ccol1=1 \
                2nd-with-counter.fits --hdu2=1             --ccol2=1 \
                --aperture=0.1 --outcols=b_all \
                --output=flagged-in-2nd.fits

   The solution to having flagged matches depends on the context, so
there is not one easy solution unfortunately.  For example one situation
where this can occur is the following: You have two catalogs of galaxies
and one is from a much deeper image (with a higher spatial resolution)
than the other.  Due to its depth and higher resolution, the deeper
catalog's points will be much more dense in the same region of the sky
and the probability of such ambiguous matches will increase.

   To fix the problem above, before doing the match, you can remove all
the points in the deep catalog that would not be detectable in the
shallow one.  To find the points of the deep catalog that should be
removed, you can generate the number-counts of your two catalogs
(histogram of their magnitudes) and plot them together.  This will show
your shallower catalog's completeness limit (the magnitude where the
shallower catalog start to deviate (drop) compared to the deeper
catalog.  You can then remove the sources from the deeper catalog that
are fainter than the completeness limit of the shallower catalog (for
example, when the limit is on the 24th magnitude, you can use ‘asttable
deep.fits --range=MAG,-inf,24’).  Finally, you can do the match on the
the resulting sub-set of the deeper catalog.

