gnuastro (0.22)

Browse

Build Log

Usage

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

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

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: Operand storage in memory or a file, Prev: Building new dataset and stack management, Up: Arithmetic operators

6.2.4.21 Operand storage in memory or a file
............................................

In your early days of using Gnuastro, to do multiple operations, it is
likely that you will simply call Arithmetic (or Table, with column
arithmetic) multiple times: feed the output file of the first call to
the second call. But as you get more proficient in the reverse polish
notation, you will find yourself combining many operations into one
call. This greatly speeds up your operation, because instead of writing
the dataset to a file in one command, and reading it in the next
command, it will just keep the intermediate dataset in memory!

But adding more complexity to your operations, can make them much
harder to debug, or extend even further. Therefore in this section we
have some special operators that behave differently from the rest: they
do not touch the contents of the data, only where/how they are stored.
They are designed to do complex operations, without necessarily having a
complex command.

‘set-AAA’
Set the characters after the dash (‘AAA’ in the case shown here) as
a name for the first popped operand on the stack. The named
dataset will be freed from memory as soon as it is no longer
needed, or if the name is reset to refer to another dataset later
in the command. This operator thus enables reusability of a
dataset without having to reread it from a file every time it is
necessary during a process. When a dataset is necessary more than
once, this operator can thus help simplify reading/writing on the
command-line (thus avoiding potential bugs), while also speeding up
the processing.

Like all operators, this operator pops the top operand off of the
main processing stack, but unlike other operands, it will not add
anything back to the stack immediately. It will keep the popped
dataset in memory through a separate list of named datasets (not on
the main stack). That list will be used to add/copy any requested
dataset to the main processing stack when the name is called.

The name to give the popped dataset is part of the operator's name.
For example, the ‘set-a’ operator of the command below, gives the
name "‘a’" to the contents of ‘image.fits’. This name is then used
instead of the actual filename to multiply the dataset by two.

$ astarithmetic image.fits set-a a 2 x

The name can be any string, but avoid strings ending with standard
filename suffixes (for example, ‘.fits’)(1).

One example of the usefulness of this operator is in the ‘where’
operator. For example, let's assume you want to mask all pixels
larger than ‘5’ in ‘image.fits’ (extension number 1) with a NaN
value. Without setting a name for the dataset, you have to read
the file two times from memory in a command like this:

$ astarithmetic image.fits image.fits 5 gt nan where -g1

But with this operator you can simply give ‘image.fits’ the name
‘i’ and simplify the command above to the more readable one below
(which greatly helps when the filename is long):

$ astarithmetic image.fits set-i i i 5 gt nan where

‘repeat’
Add N copies of the second popped operand to the stack of operands.
N is the first popped operand. For example, let's assume
‘image.fits’ is a $100\times100$ image. The output of the command
below will be a 3D datacube of size $100\times100\times20$ voxels
(volume-pixels):

$ astarithmetic image.fits 20 repeat 20 add-dimension-slow

‘tofile-AAA’
Write the top operand on the operands stack into a file called
‘AAA’ (can be any FITS file name) without changing the operands
stack. If you do not need the dataset any more and would like to
free it, see the ‘tofilefree’ operator below.

By default, any file that is given to this operator is deleted
before Arithmetic actually starts working on the input datasets.
The deletion can be deactivated with the ‘--dontdelete’ option (as
in all Gnuastro programs, see *note Input output options::). If
the same FITS file is given to this operator multiple times, it
will contain multiple extensions (in the same order that it was
called.

For example, the operator ‘tofile-check.fits’ will write the top
operand to ‘check.fits’. Since it does not modify the operands
stack, this operator is very convenient when you want to debug, or
understanding, a string of operators and operands given to
Arithmetic: simply put ‘tofile-AAA’ anywhere in the process to see
what is happening behind the scenes without modifying the overall
process.

‘tofilefree-AAA’
Similar to the ‘tofile’ operator, with the only difference that the
dataset that is written to a file is popped from the operand stack
and freed from memory (cannot be used any more).

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

(1) A dataset name like ‘a.fits’ (which can be set with ‘set-a.fits’)
will cause confusion in the initial parser of Arithmetic. It will
assume this name is a FITS file, and if it is used multiple times,
Arithmetic will abort, complaining that you have not provided enough
HDUs.

File: gnuastro.info, Node: Invoking astarithmetic, Prev: Arithmetic operators, Up: Arithmetic

6.2.5 Invoking Arithmetic
-------------------------

Arithmetic will do pixel to pixel arithmetic operations on the
individual pixels of input data and/or numbers. For the full list of
operators with explanations, please see *note Arithmetic operators::.
Any operand that only has a single element (number, or single pixel FITS
image) will be read as a number, the rest of the inputs must have the
same dimensions. The general template is:

$ astarithmetic [OPTION...] ASTRdata1 [ASTRdata2] OPERATOR ...

One line examples:

## Calculate (10.32-3.84)^2.7 quietly (will just print 155.329):
$ astarithmetic -q 10.32 3.84 - 2.7 pow

## Inverse the input image (1/pixel):
$ astarithmetic 1 image.fits / --out=inverse.fits

## Multiply each pixel in image by -1:
$ astarithmetic image.fits -1 x --out=negative.fits

## Subtract extension 4 from extension 1 (counting from zero):
$ astarithmetic image.fits image.fits - --out=skysub.fits \
--hdu=1 --hdu=4

## Add two images, then divide them by 2 (2 is read as floating point):
## Note that without the '.0', the '2' will be read/used as an integer.
$ astarithmetic image1.fits image2.fits + 2.0 / --out=average.fits

## Use Arithmetic's average operator:
$ astarithmetic image1.fits image2.fits average --out=average.fits

## Calculate the median of three images in three separate extensions:
$ astarithmetic img1.fits img2.fits img3.fits median \
-h0 -h1 -h2 --out=median.fits

Arithmetic's notation for giving operands to operators is fully
described in *note Reverse polish notation::. The output dataset is
last remaining operand on the stack. When the output dataset a single
number, and ‘--output’ is not called, it will be printed on the standard
output (command-line). When the output is an array, it will be stored
as a file.

The name of the final file can be specified with the ‘--output’
option, but if it is not given (and the output dataset has more than one
element), Arithmetic will use "automatic output" on the name of the
first FITS image encountered to generate an output file name, see *note
Automatic output::. By default, if the output file already exists, it
will be deleted before Arithmetic starts operation. However, this can
be disabled with the ‘--dontdelete’ option (see below). At any point
during Arithmetic's operation, you can also write the top operand on the
stack to a file, using the ‘tofile’ or ‘tofilefree’ operators, see *note
Arithmetic operators::.

By default, the world coordinate system (WCS) information of the
output dataset will be taken from the first input image (that contains a
WCS) on the command-line. This can be modified with the ‘--wcsfile’ and
‘--wcshdu’ options described below. When the ‘--quiet’ option is not
given, the name and extension of the dataset used for the output's WCS
is printed on the command-line.

Through operators like those starting with ‘collapse-’, the
dimensionality of the inputs may not be the same as the outputs. By
default, when the output is 1D, Arithmetic will write it as a table, not
an image/array. The format of the output table (plain text or FITS
ASCII or binary) can be set with the ‘--tableformat’ option, see *note
Input output options::). You can disable this feature (write 1D arrays
as FITS images/arrays, or to the standard output) with the
‘--onedasimage’ or ‘--onedonstdout’ options.

See *note Common options:: for a review of the options in all
Gnuastro programs. Arithmetic just redefines the ‘--hdu’ and
‘--dontdelete’ options as explained below.

‘--arguments=STR’
A plain-text file containing the command-line arguments that will
be used by Arithmetic. This option is only relevant when no
arguments are given on the command-line: if any arguments are
given, this option is ignored.

This is necessary when the set of of input files and operators
(arguments; see *note Arguments and options::) are very long
(thousands of long file names for example; usually generated within
large pipelines). Such long arguments will cause the shell to
abort with an ‘Argument list too long’ error. In such cases, you
can put the list into a plain-text file and use this option like
below (assuming you want to stack all the files in a certain
directory with the ‘mean’ operator; see *note Stacking
operators::):

$ counter=0;
$ for f in $(pwd)/*.fits; do \
echo $f; counter=$((counter+1)); \
done > arguments.txt; \
echo "$counter mean" >> arguments.txt
$ astarithmetic --arguments=list.txt -g1

‘-h INT/STR’
‘--hdu INT/STR’
The header data unit of the input FITS images, see *note Input
output options::. Unlike most options in Gnuastro (which will
ultimately only have one value for this option), Arithmetic allows
‘--hdu’ to be called multiple times and the value of each
invocation will be stored separately (for the unlimited number of
input images you would like to use). Recall that for other
programs this (common) option only takes a single value. So in
other programs, if you specify it multiple times on the
command-line, only the last value will be used and in the
configuration files, it will be ignored if it already has a value.

The order of the values to ‘--hdu’ has to be in the same order as
input FITS images. Options are first read from the command-line
(from left to right), then top-down in each configuration file, see
*note Configuration file precedence::.

If the number of HDUs is less than the number of input images,
Arithmetic will abort and notify you. However, if there are more
HDUs than FITS images, there is no problem: they will be used in
the given order (every time a FITS image comes up on the stack) and
the extra HDUs will be ignored in the end. So there is no problem
with having extra HDUs in the configuration files and by default
several HDUs with a value of ‘0’ are kept in the system-wide
configuration file when you install Gnuastro.

‘-g INT/STR’
‘--globalhdu INT/STR’
Use the value to this option as the HDU of all input FITS files.
This option is very convenient when you have many input files and
the dataset of interest is in the same HDU of all the files. When
this option is called, any values given to the ‘--hdu’ option
(explained above) are ignored and will not be used.

‘-w FITS’
‘--wcsfile FITS’
FITS Filename containing the WCS structure that must be written to
the output. The HDU/extension should be specified with ‘--wcshdu’.

When this option is used, the respective WCS will be read before
any processing is done on the command-line and directly used in the
final output. If the given file does not have any WCS, then the
default WCS (first file on the command-line with WCS) will be used
in the output.

This option will mostly be used when the default file (first of the
set of inputs) is not the one containing your desired WCS. But with
this option, you can also use Arithmetic to rewrite/change the WCS
of an existing FITS dataset from another file:

$ astarithmetic data.fits --wcsfile=other.fits -ofinal.fits

‘-W STR’
‘--wcshdu STR’
HDU/extension to read the WCS within the file given to ‘--wcsfile’.
For more, see the description of ‘--wcsfile’.

‘--envseed’
Use the environment for the random number generator settings in
operators that need them (for example, ‘mknoise-sigma’). This is
very important for obtaining reproducible results, for more see
*note Generating random numbers::.

‘-n STR’
‘--metaname=STR’
Metadata (name) of the output dataset. For a FITS image or table,
the string given to this option is written in the ‘EXTNAME’ or
‘TTYPE1’ keyword (respectively).

If this keyword is present in a FITS extension, it will be printed
in the table output of a command like ‘astfits image.fits’ (for
images) or ‘asttable table.fits -i’ (for tables). This metadata
can be very helpful for yourself in the future (when you have
forgotten the details), so it is recommended to use this option for
files that should be archived or shared with colleagues.

‘-u STR’
‘--metaunit=STR’
Metadata (units) of the output dataset. For a FITS image or table,
the string given to this option is written in the ‘BUNIT’ or
‘TTYPE1’ keyword respectively. In the case of tables, recall that
the Arithmetic program only outputs a single column, you should use
column arithmetic in Table for more than one column (see *note
Column arithmetic::). For more on the importance of metadata, see
the description of ‘--metaname’.

‘-c STR’
‘--metacomment=STR’
Metadata (comments) of the output dataset. For a FITS image or
table, the string given to this option is written in the ‘COMMENT’
or ‘TCOMM1’ keyword respectively. In the case of tables, recall
that the Arithmetic program only outputs a single column, you
should use column arithmetic in Table for more than one column (see
*note Column arithmetic::). For more on the importance of
metadata, see the description of ‘--metaname’.

‘-O’
‘--onedasimage’
Write final dataset as a FITS image/array even if it has a single
dimension. By default, if the output is 1D, it will be written as
a table, see above. If the output has more than one dimension,
this option is redundant.

‘-s’
‘--onedonstdout’
Write final dataset (only when it is 1D) to standard output, not as
a file. By default 1D datasets will be written as a table, see
above. If the output has more than one dimension, this option is
redundant.

‘-D’
‘--dontdelete’
Do not delete the output file, or files given to the ‘tofile’ or
‘tofilefree’ operators, if they already exist. Instead append the
desired datasets to the extensions that already exist in the
respective file. Note it does not matter if the final output file
name is given with the ‘--output’ option, or determined
automatically.

Arithmetic treats this option differently from its default
operation in other Gnuastro programs (see *note Input output
options::). If the output file exists, when other Gnuastro
programs are called with ‘--dontdelete’, they simply complain and
abort. But when Arithmetic is called with ‘--dontdelete’, it will
appended the dataset(s) to the existing extension(s) in the file.

‘-a’
‘--writeall’
Write all datasets on the stack as separate HDUs in the output
file. This only affects datasets with multiple dimensions (or
single-dimension datasets when the ‘--onedasimg’ is called). This
option is useful to debug Arithmetic calls: to check all the images
on the stack while you are designing your operation. The top
dataset on the stack will be on HDU number 1 of the output, the
second dataset will be on HDU number 2 and so on.

Arithmetic accepts two kinds of input: images and numbers. Images
are considered to be any of the inputs that is a file name of a
recognized type (see *note Arguments::) and has more than one
element/pixel. Numbers on the command-line will be read into the
smallest type (see *note Numeric data types::) that can store them, so
‘-2’ will be read as a ‘char’ type (which is signed on most systems and
can thus keep negative values), ‘2500’ will be read as an ‘unsigned
short’ (all positive numbers will be read as unsigned), while
‘3.1415926535897’ will be read as a ‘double’ and ‘3.14’ will be read as
a ‘float’. To force a number to be read as float, put a ‘.’ after it
(possibly followed by a zero for easier readability), or add an ‘f’
after it. Hence while ‘5’ will be read as an integer, ‘5.’, ‘5.0’ or
‘5f’ will be added to the stack as ‘float’ (see *note Reverse polish
notation::).

Unless otherwise stated (in *note Arithmetic operators::), the
operators can deal with numeric multiple data types (see *note Numeric
data types::). For example, in "‘a.fits b.fits +’", the image types can
be ‘long’ and ‘float’. In such cases, C's internal type conversion will
be used. The output type will be set to the higher-ranking type of the
two inputs. Unsigned integer types have smaller ranking than their
signed counterparts and floating point types have higher ranking than
the integer types. So the internal C type conversions done in the
example above are equivalent to this piece of C:

size_t i;
long a[100];
float b[100], out[100];
for(i=0;i<100;++i) out[i]=a[i]+b[i];

Relying on the default C type conversion significantly speeds up the
processing and also requires less RAM (when using very large images).

Some operators can only work on integer types (of any length, for
example, bitwise operators) while others only work on floating point
types, (currently only the ‘pow’ operator). In such cases, if the
operand type(s) are different, an error will be printed. Arithmetic
also comes with internal type conversion operators which you can use to
convert the data into the appropriate type, see *note Arithmetic
operators::.

The hyphen (‘-’) can be used both to specify options (see *note
Options::) and also to specify a negative number which might be
necessary in your arithmetic. In order to enable you to do this,
Arithmetic will first parse all the input strings and if the first
character after a hyphen is a digit, then that hyphen is temporarily
replaced by the vertical tab character which is not commonly used. The
arguments are then parsed and these strings will not be specified as an
option. Then the given arguments are parsed and any vertical tabs are
replaced back with a hyphen so they can be read as negative numbers.
Therefore, as long as the names of the files you want to work on, do not
start with a vertical tab followed by a digit, there is no problem. An
important consequence of this implementation is that you should not
write negative fractions like this: ‘-.3’, instead write them as ‘-0.3’.

Without any images, Arithmetic will act like a simple calculator and
print the resulting output number on the standard output like the first
example above. If you really want such calculator operations on the
command-line, AWK (GNU AWK is the most common implementation) is much
faster, easier and much more powerful. For example, the numerical
one-line example above can be done with the following command. In
general AWK is a fantastic tool and GNU AWK has a wonderful manual
(<https://www.gnu.org/software/gawk/manual/>). So if you often confront
situations like this, or have to work with large text tables/catalogs,
be sure to checkout AWK and simplify your life.

$ echo "" | awk '{print (10.32-3.84)^2.7}'
155.329

File: gnuastro.info, Node: Convolve, Next: Warp, Prev: Arithmetic, Up: Data manipulation

6.3 Convolve
============

On an image, convolution can be thought of as a process to blur or
remove the contrast in an image. If you are already familiar with the
concept and just want to run Convolve, you can jump to *note Convolution
kernel:: and *note Invoking astconvolve:: and skip the lengthy
introduction on the basic definitions and concepts of convolution.

There are generally two methods to convolve an image. The first and
more intuitive one is in the "spatial domain" or using the actual image
pixel values, see *note Spatial domain convolution::. The second method
is when we manipulate the "frequency domain", or work on the magnitudes
of the different frequencies that constitute the image, see *note
Frequency domain and Fourier operations::. Understanding convolution in
the spatial domain is more intuitive and thus recommended if you are
just starting to learn about convolution. However, getting a good grasp
of the frequency domain is a little more involved and needs some
concentration and some mathematical proofs. However, its reward is a
faster operation and more importantly a very fundamental understanding
of this very important operation.

Convolution of an image will generally result in blurring the image
because it mixes pixel values. In other words, if the image has sharp
differences in neighboring pixel values(1), those sharp differences will
become smoother. This has very good consequences in detection of signal
in noise for example. In an actual observed image, the variation in
neighboring pixel values due to noise can be very high. But after
convolution, those variations will decrease and we have a better hope in
detecting the possible underlying signal. Another case where
convolution is extensively used is in mock images and modeling in
general, convolution can be used to simulate the effect of the
atmosphere or the optical system on the mock profiles that we create,
see *note PSF::. Convolution is a very interesting and important topic
in any form of signal analysis (including astronomical observations).
So we have thoroughly(2) explained the concepts behind it in the
following sub-sections.

* Menu:

* Spatial domain convolution:: Only using the input image values.
* Frequency domain and Fourier operations:: Using frequencies in input.
* Spatial vs. Frequency domain:: When to use which?
* Convolution kernel:: How to specify the convolution kernel.
* Invoking astconvolve:: Options and argument to Convolve.

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

(1) In astronomy, the only major time we confront such sharp borders
in signal are cosmic rays. All other sources of signal in an image are
already blurred by the atmosphere or the optics of the instrument.

(2) A mathematician will certainly consider this explanation is
incomplete and inaccurate. However this text is written for an
understanding on the operations that are done on a real (not complex,
discrete and noisy) astronomical image, not any general form of abstract
function

File: gnuastro.info, Node: Spatial domain convolution, Next: Frequency domain and Fourier operations, Prev: Convolve, Up: Convolve

6.3.1 Spatial domain convolution
--------------------------------

The pixels in an input image represent different "spatial" positions,
therefore when convolution is done only using the actual input pixel
values, we name the process as being done in the "Spatial domain". In
particular this is in contrast to the "frequency domain" that we will
discuss later in *note Frequency domain and Fourier operations::. In
the spatial domain (and in realistic situations where the image and the
convolution kernel do not extend to infinity), convolution is the
process of changing the value of one pixel to the _weighted_ average of
all the pixels in its _neighborhood_.

The 'neighborhood' of each pixel (how many pixels in which direction)
and the 'weight' function (how much each neighboring pixel should
contribute depending on its position) are given through a second image
which is known as a "kernel"(1).

* Menu:

* Convolution process:: More basic explanations.
* Edges in the spatial domain:: Dealing with the edges of an image.

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

(1) Also known as filter, here we will use 'kernel'.

File: gnuastro.info, Node: Convolution process, Next: Edges in the spatial domain, Prev: Spatial domain convolution, Up: Spatial domain convolution

6.3.1.1 Convolution process
...........................

In convolution, the kernel specifies the weight and positions of the
neighbors of each pixel. To find the convolved value of a pixel, the
central pixel of the kernel is placed on that pixel. The values of each
overlapping pixel in the kernel and image are multiplied by each other
and summed for all the kernel pixels. To have one pixel in the center,
the sides of the convolution kernel have to be an odd number. This
process effectively mixes the pixel values of each pixel with its
neighbors, resulting in a blurred image compared to the sharper input
image.

Formally, convolution is one kind of linear 'spatial filtering' in
image processing texts. If we assume that the kernel has $2a+1$ and
$2b+1$ pixels on each side, the convolved value of a pixel placed at $x$
and $y$ ($C_{x,y}$) can be calculated from the neighboring pixel values
in the input image ($I$) and the kernel ($K$) from

$$C_{x,y}=\sum_{s=-a}^{a}\sum_{t=-b}^{b}K_{s,t}\times{}I_{x+s,y+t}.$$

Formally, any pixel that is outside of the image in the equation
above will be considered to be zero (although, see *note Edges in the
spatial domain::). When the kernel is symmetric about its center the
blurred image has the same orientation as the original image. However,
if the kernel is not symmetric, the image will be affected in the
opposite manner, this is a natural consequence of the definition of
spatial filtering. In order to avoid this we can rotate the kernel
about its center by 180 degrees so the convolved output can have the
same original orientation (this is done by default in the Convolve
program). Technically speaking, only if the kernel is flipped the
process is known as _Convolution_. If it is not it is known as
_Correlation_.

To be a weighted average, the sum of the weights (the pixels in the
kernel) has to be unity. This will have the consequence that the
convolved image of an object and unconvolved object will have the same
brightness (see *note Brightness flux magnitude::), which is natural,
because convolution should not eat up the object photons, it only
disperses them.

The convolution of each pixel is independent of the other pixels, and
in some cases, it may be necessary to convolve different parts of an
image separately (for example, when you have different amplifiers on the
CCD). Therefore, to speed up spatial convolution, Gnuastro first defines
a tessellation over the input; assigning each group of pixels to
"tiles". It then does the convolution in parallel on each tile. For
more on how Gnuastro's programs create the tile grid (tessellation), see
*note Tessellation::.

File: gnuastro.info, Node: Edges in the spatial domain, Prev: Convolution process, Up: Spatial domain convolution

6.3.1.2 Edges in the spatial domain
...................................

In purely 'linear' spatial filtering (convolution), there are problems
with the edges of the input image. Here we will explain the problem in
the spatial domain. For a discussion of this problem from the frequency
domain perspective, see *note Edges in the frequency domain::. The
problem originates from the fact that on the edges, in practice, the sum
of the weights we use on the actual image pixels is not unity(1). For
example, as discussed above, a profile in the center of an image will
have the same brightness before and after convolution. However, for
partially imaged profile on the edge of the image, the brightness (sum
of its pixel fluxes within the image, see *note Brightness flux
magnitude::) will not be equal, some of the flux is going to be 'eaten'
by the edges.

If you run ‘$ make check’ on the source files of Gnuastro, you can
see this effect by comparing the ‘convolve_frequency.fits’ with
‘convolve_spatial.fits’ in the ‘./tests/’ directory. In the spatial
domain, by default, no assumption will be made about pixels outside of
the image or any blank pixels in the image. The problem explained above
will also occur on the sides of blank regions (see *note Blank
pixels::). The solution to this edge effect problem is only possible in
the spatial domain. For pixels near the edge, we have to abandon the
assumption that the sum of the kernel pixels is unity during the
convolution process(2). So taking $W$ as the sum of the kernel pixels
that overlapped with non-blank and in-image pixels, the equation in
*note Convolution process:: will become:

$$C_{x,y}= { \sum_{s=-a}^{a}\sum_{t=-b}^{b}K_{s,t}\times{}I_{x+s,y+t} \over W}.$$

In this manner, objects which are near the edges of the image or blank
pixels will also have the same brightness (within the image) before and
after convolution. This correction is applied by default in Convolve
when convolving in the spatial domain. To disable it, you can use the
‘--noedgecorrection’ option. In the frequency domain, there is no way
to avoid this loss of flux near the edges of the image, see *note Edges
in the frequency domain:: for an interpretation from the frequency
domain perspective.

Note that the edge effect discussed here is different from the one in
*note If convolving afterwards::. In making mock images we want to
simulate a real observation. In a real observation, the images of the
galaxies on the sides of the CCD are first blurred by the atmosphere and
instrument, then imaged. So light from the parts of a galaxy which are
immediately outside the CCD will affect the parts of the galaxy which
are covered by the CCD. Therefore in modeling the observation, we have
to convolve an image that is larger than the input image by exactly half
of the convolution kernel. We can hence conclude that this correction
for the edges is only useful when working on actual observed images
(where we do not have any more data on the edges) and not in modeling.

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

(1) Because we assumed the overlapping pixels outside the input image
have a value of zero.

(2) Of course the sum of the kernel pixels still have to be unity in
general.

File: gnuastro.info, Node: Frequency domain and Fourier operations, Next: Spatial vs. Frequency domain, Prev: Spatial domain convolution, Up: Convolve

6.3.2 Frequency domain and Fourier operations
---------------------------------------------

Getting a good grip on the frequency domain is usually not an easy job!
So we have decided to give the issue a complete review here.
Convolution in the frequency domain (see *note Convolution theorem::)
heavily relies on the concepts of Fourier transform (*note Fourier
transform::) and Fourier series (*note Fourier series::) so we will be
investigating these important operations first. It has become something
of a cliché for people to say that the Fourier series "is a way to
represent a (wave-like) function as the sum of simple sine waves" (from
Wikipedia). However, sines themselves are abstract functions, so this
statement really adds no extra layer of physical insight.

Before jumping head-first into the equations and proofs, we will
begin with a historical background to see how the importance of
frequencies actually roots in our ancient desire to see everything in
terms of circles. A short review of how the complex plane should be
interpreted is then given. Having paved the way with these two basics,
we define the Fourier series and subsequently the Fourier transform.
The final aim is to explain discrete Fourier transform, however some
very important concepts need to be solidified first: The Dirac comb,
convolution theorem and sampling theorem. So each of these topics are
explained in their own separate sub-sub-section before going on to the
discrete Fourier transform. Finally we revisit (after *note Edges in
the spatial domain::) the problem of convolution on the edges, but this
time in the frequency domain. Understanding the sampling theorem and
the discrete Fourier transform is very important in order to be able to
pull out valuable science from the discrete image pixels. Therefore we
have included the mathematical proofs and figures so you can have a
clear understanding of these very important concepts.

* Menu:

* Fourier series historical background:: Historical background.
* Circles and the complex plane:: Interpreting complex numbers.
* Fourier series:: Fourier Series definition.
* Fourier transform:: Fourier Transform definition.
* Dirac delta and comb:: Dirac delta and Dirac comb.
* Convolution theorem:: Derivation of Convolution theorem.
* Sampling theorem:: Sampling theorem (Nyquist frequency).
* Discrete Fourier transform:: Derivation and explanation of DFT.
* Fourier operations in two dimensions:: Extend to 2D images.
* Edges in the frequency domain:: Interpretation of edge effects.

File: gnuastro.info, Node: Fourier series historical background, Next: Circles and the complex plane, Prev: Frequency domain and Fourier operations, Up: Frequency domain and Fourier operations

6.3.2.1 Fourier series historical background
............................................

Ever since the ancient times, the circle has been (and still is) the
simplest shape for abstract comprehension. All you need is a center
point and a radius and you are done. All the points on a circle are at
a fixed distance from the center. However, the moment you try to
connect this elegantly simple and beautiful abstract construct (the
circle) with the real world (for example, compute its area or its
circumference), things become really hard (ideally, impossible) because
the irrational number $\pi$ gets involved.

The key to understanding the Fourier series (thus the Fourier
transform and finally the Discrete Fourier Transform) is our ancient
desire to express everything in terms of circles or the most
exceptionally simple and elegant abstract human construct. Most people
prefer to say the same thing in a more ahistorical manner: to break a
function into sines and cosines. As the term "ancient" in the previous
sentence implies, Jean-Baptiste Joseph Fourier (1768 - 1830 A.D.) was
not the first person to do this. The main reason we know this process
by his name today is that he came up with an ingenious method to find
the necessary coefficients (radius of) and frequencies ("speed" of
rotation on) the circles for any generic (integrable) function.

[image src="gnuastro-figures/epicycles.png" alt="Middle ages epicycles along with two demonstrations of breaking a generic function using epicycles." text="../gnuastro-figures//epicycles.png"]

Figure 6.1: Epicycles and the Fourier series. Left: A demonstration of
Mercury's epicycles relative to the "center of the world" by Qutb al-Din
al-Shirazi (1236 - 1311 A.D.) retrieved from Wikipedia
(https://commons.wikimedia.org/wiki/File:Ghotb2.jpg). Middle
(https://commons.wikimedia.org/wiki/File:Fourier_series_square_wave_circles_animation.gif)
and Right: How adding more epicycles (or terms in the Fourier series)
will approximate functions. The right
(https://commons.wikimedia.org/wiki/File:Fourier_series_sawtooth_wave_circles_animation.gif)
animation is also available.

Like most aspects of mathematics, this process of interpreting
everything in terms of circles, began for astronomical purposes. When
astronomers noticed that the orbit of Mars and other outer planets, did
not appear to be a simple circle (as everything should have been in the
heavens). At some point during their orbit, the revolution of these
planets would become slower, stop, go back a little (in what is known as
the retrograde motion) and then continue going forward again.

The correction proposed by Ptolemy (90 - 168 A.D.) was the most
agreed upon. He put the planets on Epicycles or circles whose center
itself rotates on a circle whose center is the earth. Eventually, as
observations became more and more precise, it was necessary to add more
and more epicycles in order to explain the complex motions of the
planets(1). *note Figure 6.1: epicycle.(Left) shows an example
depiction of the epicycles of Mercury in the late 13th century.

Of course we now know that if they had abdicated the Earth from its
throne in the center of the heavens and allowed the Sun to take its
place, everything would become much simpler and true. But there was not
enough observational evidence for changing the "professional consensus"
of the time to this radical view suggested by a small minority(2). So
the pre-Galilean astronomers chose to keep Earth in the center and find
a correction to the models (while keeping the heavens a purely
"circular" order).

The main reason we are giving this historical background which might
appear off topic is to give historical evidence that while such
"approximations" do work and are very useful for pragmatic reasons (like
measuring the calendar from the movement of astronomical bodies). They
offer no physical insight. The astronomers who were involved with the
Ptolemaic world view had to add a huge number of epicycles during the
centuries after Ptolemy in order to explain more accurate observations.
Finally the death knell of this world-view was Galileo's observations
with his new instrument (the telescope). So the physical insight, which
is what Astronomers and Physicists are interested in (as opposed to
Mathematicians and Engineers who just like proving and optimizing or
calculating!) comes from being creative and not limiting ourselves to
such approximations. Even when they work.

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

(1) See the Wikipedia page on "Deferent and epicycle" for a more
complete historical review.

(2) Aristarchus of Samos (310 - 230 B.C.) appears to be one of the
first people to suggest the Sun being in the center of the universe.
This approach to science (that the standard model is defined by
consensus) and the fact that this consensus might be completely wrong
still applies equally well to our models of particle physics and
cosmology today.

File: gnuastro.info, Node: Circles and the complex plane, Next: Fourier series, Prev: Fourier series historical background, Up: Frequency domain and Fourier operations

6.3.2.2 Circles and the complex plane
.....................................

Before going onto the derivation, it is also useful to review how the
complex numbers and their plane relate to the circles we talked about
above. The two schematics in the middle and right of *note Figure 6.1:
epicycle. show how a 1D function of time can be made using the 2D real
and imaginary surface. Seeing the animation in Wikipedia will really
help in understanding this important concept. At each point in time, we
take the vertical coordinate of the point and use it to find the value
of the function at that point in time. *note Figure 6.2: iandtime.
shows this relation with the axes marked.

Leonhard Euler(1) (1707 - 1783 A.D.) showed that the complex
exponential ($e^{iv}$ where $v$ is real) is periodic and can be written
as: $e^{iv}=\cos{v}+isin{v}$. Therefore $e^{iv+2\pi}=e^{iv}$. Later,
Caspar Wessel (mathematician and cartographer 1745 - 1818 A.D.) showed
how complex numbers can be displayed as vectors on a plane. Euler's
identity might seem counter intuitive at first, so we will try to
explain it geometrically (for deeper physical insight). On the
real-imaginary 2D plane (like the left hand plot in each box of *note
Figure 6.2: iandtime.), multiplying a number by $i$ can be interpreted
as rotating the point by $90$ degrees (for example, the value $3$ on the
real axis becomes $3i$ on the imaginary axis). On the other hand,
$e\equiv\lim_{n\rightarrow\infty}(1+{1\over n})^n$, therefore, defining
$m\equiv nu$, we get:

$$e^{u}=\lim_{n\rightarrow\infty}\left(1+{1\over n}\right)^{nu}
=\lim_{n\rightarrow\infty}\left(1+{u\over nu}\right)^{nu}
=\lim_{m\rightarrow\infty}\left(1+{u\over m}\right)^{m}$$

Taking $u\equiv iv$ the result can be written as a generic complex
number (a function of $v$):

$$e^{iv}=\lim_{m\rightarrow\infty}\left(1+i{v\over
m}\right)^{m}=a(v)+ib(v)$$

For $v=\pi$, a nice geometric animation of going to the limit can be
seen on Wikipedia (https://commons.wikimedia.org/wiki/File:ExpIPi.gif).
We see that $\lim_{m\rightarrow\infty}a(\pi)=-1$, while
$\lim_{m\rightarrow\infty}b(\pi)=0$, which gives the famous
$e^{i\pi}=-1$ equation. The final value is the real number $-1$,
however the distance of the polygon points traversed as
$m\rightarrow\infty$ is half the circumference of a circle or $\pi$,
showing how $v$ in the equation above can be interpreted as an angle in
units of radians and therefore how $a(v)=cos(v)$ and $b(v)=sin(v)$.

Since $e^{iv}$ is periodic (let's assume with a period of $T$), it is
more clear to write it as $v\equiv{2{\pi}n\over T}t$ (where $n$ is an
integer), so $e^{iv}=e^{i{2{\pi}n\over T}t}$. The advantage of this
notation is that the period ($T$) is clearly visible and the frequency
($2{\pi}n \over T$, in units of 1/cycle) is defined through the integer
$n$. In this notation, $t$ is in units of "cycle"s.

As we see from the examples in *note Figure 6.1: epicycle. and *note
Figure 6.2: iandtime, for each constituting frequency, we need a
respective 'magnitude' or the radius of the circle in order to
accurately approximate the desired 1D function. The concepts of
"period" and "frequency" are relatively easy to grasp when using
temporal units like time because this is how we define them in every-day
life. However, in an image (astronomical data), we are dealing with
spatial units like distance. Therefore, by one "period" we mean the
_distance_ at which the signal is identical and frequency is defined as
the inverse of that spatial "period". The complex circle of *note
Figure 6.2: iandtime. can be thought of the Moon rotating about Earth
which is rotating around the Sun; so the "Real (signal)" axis shows the
Moon's position as seen by a distant observer on the Sun as time goes
by. Because of the scalar (not having any direction or vector) nature
of time, *note Figure 6.2: iandtime. is easier to understand in units of
time. When thinking about spatial units, mentally replace the "Time
(sec)" axis with "Distance (meters)". Because length has direction and
is a vector, visualizing the rotation of the imaginary circle and the
advance along the "Distance (meters)" axis is not as simple as temporal
units like time.

[image src="gnuastro-figures/iandtime.png" text="../gnuastro-figures//iandtime.eps"]

Figure 6.2: Relation between the real (signal), imaginary
($i\equiv\sqrt{-1}$) and time axes at two snapshots of time.

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

(1) Other forms of this equation were known before Euler. For
example, in 1707 A.D. (the year of Euler's birth) Abraham de Moivre
(1667 - 1754 A.D.) showed that
$(\cos{x}+i\sin{x})^n=\cos(nx)+i\sin(nx)$. In 1714 A.D., Roger Cotes
(1682 - 1716 A.D. a colleague of Newton who proofread the second edition
of Principia) showed that: $ix=\ln(\cos{x}+i\sin{x})$.

File: gnuastro.info, Node: Fourier series, Next: Fourier transform, Prev: Circles and the complex plane, Up: Frequency domain and Fourier operations

6.3.2.3 Fourier series
......................

In astronomical images, our variable (brightness, or number of
photo-electrons, or signal to be more generic) is recorded over the 2D
spatial surface of a camera pixel. However to make things easier to
understand, here we will assume that the signal is recorded in 1D
(assume one row of the 2D image pixels). Also for this section and the
next (*note Fourier transform::) we will be talking about the signal
before it is digitized or pixelated. Let's assume that we have the
continuous function $f(l)$ which is integrable in the interval $[l_0,
l_0+L]$ (always true in practical cases like images). Take $l_0$ as the
position of the first pixel in the assumed row of the image and $L$ as
the width of the image along that row. The units of $l_0$ and $L$ can
be in any spatial units (for example, meters) or an angular unit (like
radians) multiplied by a fixed distance which is more common.

To approximate $f(l)$ over this interval, we need to find a set of
frequencies and their corresponding 'magnitude's (see *note Circles and
the complex plane::). Therefore our aim is to show $f(l)$ as the
following sum of periodic functions:

$$f(l)=\displaystyle\sum_{n=-\infty}^{\infty}c_ne^{i{2{\pi}n\over L}l} $$

Note that the different frequencies ($2{\pi}n/L$, in units of cycles per
meters for example) are not arbitrary. They are all integer multiples
of the fundamental frequency of $\omega_0=2\pi/L$. Recall that $L$ was
the length of the signal we want to model. Therefore, we see that the
smallest possible frequency (or the frequency resolution) in the end,
depends on the length we observed the signal or $L$. In the case of
each dimension on an image, this is the size of the image in the
respective dimension. The frequencies have been defined in this
"harmonic" fashion to insure that the final sum is periodic outside of
the $[l_0, l_0+L]$ interval too. At this point, you might be thinking
that the sky is not periodic with the same period as my camera's view
angle. You are absolutely right! The important thing is that since
your camera's observed region is the only region we are "observing" and
will be using, the rest of the sky is irrelevant; so we can safely
assume the sky is periodic outside of it. However, this working
assumption will haunt us later in *note Edges in the frequency domain::.

The frequencies are thus determined by definition. So all we need to
do is to find the coefficients ($c_n$), or magnitudes, or radii of the
circles for each frequency which is identified with the integer $n$.
Fourier's approach was to multiply both sides with a fixed term:

$$f(l)e^{-i{2{\pi}m\over L}l}=\displaystyle\sum_{n=-\infty}^{\infty}c_ne^{i{2{\pi}(n-m)\over L}l}
$$

where $m>0$(1). We can then integrate both sides over the observation
period:

$$\int_{l_0}^{l_0+L}f(l)e^{-i{2{\pi}m\over L}l}dl
=\int_{l_0}^{l_0+L}\displaystyle\sum_{n=-\infty}^{\infty}c_ne^{i{2{\pi}(n-m)\over L}l}dl=\displaystyle\sum_{n=-\infty}^{\infty}c_n\int_{l_0}^{l_0+L}e^{i{2{\pi}(n-m)\over L}l}dl
$$

Both $n$ and $m$ are positive integers. Also, we know that a complex
exponential is periodic so after one period ($L$) it comes back to its
starting point. Therefore $\int_{l_0}^{l_0+L}e^{2{\pi}k/L}dl=0$ for any
$k>0$. However, when $k=0$, this integral becomes:
$\int_{l_0}^{l_0+T}e^0dt=\int_{l_0}^{l_0+T}dt=T$. Hence since the
integral will be zero for all $n{\neq}m$, we get:

$$\displaystyle\sum_{n=-\infty}^{\infty}c_n\int_{l_0}^{l_0+T}e^{i{2{\pi}(n-m)\over L}l}dl=Lc_m $$

The origin of the axis is fundamentally an arbitrary position. So let's
set it to the start of the image such that $l_0=0$. So we can find the
"magnitude" of the frequency $2{\pi}m/L$ within $f(l)$ through the
relation:

$$c_m={1\over L}\int_{0}^{L}f(l)e^{-i{2{\pi}m\over L}l}dl $$

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

(1) We could have assumed $m<0$ and set the exponential to positive,
but this is more clear.

File: gnuastro.info, Node: Fourier transform, Next: Dirac delta and comb, Prev: Fourier series, Up: Frequency domain and Fourier operations

6.3.2.4 Fourier transform
.........................

In *note Fourier series::, we had to assume that the function is
periodic outside of the desired interval with a period of $L$.
Therefore, assuming that $L\rightarrow\infty$ will allow us to work with
any function. However, with this approximation, the fundamental
frequency ($\omega_0$) or the frequency resolution that we discussed in
*note Fourier series:: will tend to zero: $\omega_0\rightarrow0$. In
the equation to find $c_m$, every $m$ represented a frequency (multiple
of $\omega_0$) and the integration on $l$ removes the dependence of the
right side of the equation on $l$, making it only a function of $m$ or
frequency. Let's define the following two variables:

$$\omega{\equiv}m\omega_0={2{\pi}m\over L}$$

$$F(\omega){\equiv}Lc_m$$

The equation to find the coefficients of each frequency in *note Fourier
series:: thus becomes:

$$F(\omega)=\int_{-\infty}^{\infty}f(l)e^{-i{\omega}l}dl.$$

The function $F(\omega)$ is thus the _Fourier transform_ of $f(l)$ in
the frequency domain. So through this transformation, we can find
(analyze) the magnitudes of the constituting frequencies or the value in
the frequency space(1) of our spatial input function. The great thing
is that we can also do the reverse and later synthesize the input
function from its Fourier transform. Let's do it: with the
approximations above, multiply the right side of the definition of the
Fourier Series (*note Fourier series::) with
$1=L/L=({\omega_0}L)/(2\pi)$:

$$f(l)={1\over
2\pi}\displaystyle\sum_{n=-\infty}^{\infty}Lc_ne^{{2{\pi}in\over
L}l}\omega_0={1\over
2\pi}\displaystyle\sum_{n=-\infty}^{\infty}F(\omega)e^{i{\omega}l}\Delta\omega
$$

To find the right most side of this equation, we renamed $\omega_0$ as
$\Delta\omega$ because it was our resolution, $2{\pi}n/L$ was written as
$\omega$ and finally, $Lc_n$ was written as $F(\omega)$ as we defined
above. Now, as $L\rightarrow\infty$, $\Delta\omega\rightarrow0$ so we
can write:

$$f(l)={1\over
2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i{\omega}l}d\omega $$

Together, these two equations provide us with a very powerful set of
tools that we can use to process (analyze) and recreate (synthesize) the
input signal. Through the first equation, we can break up our input
function into its constituent frequencies and analyze it, hence it is
also known as _analysis_. Using the second equation, we can synthesize
or make the input function from the known frequencies and their
magnitudes. Thus it is known as _synthesis_. Here, we symbolize the
Fourier transform (analysis) and its inverse (synthesis) of a function
$f(l)$ and its Fourier Transform $F(\omega)$ as ${\cal F}[f]$ and ${\cal
F}^{-1}[F]$.

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

(1) As we discussed before, this 'magnitude' can be interpreted as
the radius of the circle rotating at this frequency in the epicyclic
interpretation of the Fourier series, see *note Figure 6.1: epicycle.
and *note Figure 6.2: iandtime.

File: gnuastro.info, Node: Dirac delta and comb, Next: Convolution theorem, Prev: Fourier transform, Up: Frequency domain and Fourier operations

6.3.2.5 Dirac delta and comb
............................

The Dirac $\delta$ (delta) function (also known as an impulse) is the
way that we convert a continuous function into a discrete one. It is
defined to satisfy the following integral:

$$\int_{-\infty}^{\infty}\delta(l)dl=1$$

When integrated with another function, it gives that function's value at
$l=0$:

$$\int_{-\infty}^{\infty}f(l)\delta(l)dt=f(0)$$

An impulse positioned at another point (say $l_0$) is written as
$\delta(l-l_0)$:

$$\int_{-\infty}^{\infty}f(l)\delta(l-l_0)dt=f(l_0)$$

The Dirac $\delta$ function also operates similarly if we use summations
instead of integrals. The Fourier transform of the delta function is:

$${\cal F}[\delta(l)]=\int_{-\infty}^{\infty}\delta(l)e^{-i{\omega}l}dl=e^{-i{\omega}0}=1$$

$${\cal F}[\delta(l-l_0)]=\int_{-\infty}^{\infty}\delta(l-l_0)e^{-i{\omega}l}dl=e^{-i{\omega}l_0}$$

From the definition of the Dirac $\delta$ we can also define a Dirac
comb (${\rm III}_P$) or an impulse train with infinite impulses
separated by $P$:

$${\rm III}_P(l)\equiv\displaystyle\sum_{k=-\infty}^{\infty}\delta(l-kP) $$

$P$ is chosen to represent "pixel width" later in *note Sampling
theorem::. Therefore the Dirac comb is periodic with a period of $P$.
We have intentionally used a different name for the period of the Dirac
comb compared to the input signal's length of observation that we showed
with $L$ in *note Fourier series::. This difference is highlighted here
to avoid confusion later when these two periods are needed together in
*note Discrete Fourier transform::. The Fourier transform of the Dirac
comb will be necessary in *note Sampling theorem::, so let's derive it.
By its definition, it is periodic, with a period of $P$, so the Fourier
coefficients of its Fourier Series (*note Fourier series::) can be
calculated within one period:

$${\rm III}_P=\displaystyle\sum_{n=-\infty}^{\infty}c_ne^{i{2{\pi}n\over
P}l}$$

We can now find the $c_n$ from *note Fourier series:::

$$c_n={1\over P}\int_{-P/2}^{P/2}\delta(l)e^{-i{2{\pi}n\over P}l}
={1\over P}\quad\quad \rightarrow \quad\quad
{\rm III}_P={1\over P}\displaystyle\sum_{n=-\infty}^{\infty}e^{i{2{\pi}n\over P}l}
$$

So we can write the Fourier transform of the Dirac comb as:

$${\cal F}[{\rm III}_P]=\int_{-\infty}^{\infty}{\rm III}_Pe^{-i{\omega}l}dl
={1\over P}\displaystyle\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-i(\omega-{2{\pi}n\over P})l}dl={1\over P}\displaystyle\sum_{n=-\infty}^{\infty}\delta\left(\omega-{2{\pi}n\over P}\right)
$$

In the last step, we used the fact that the complex exponential is a
periodic function, that $n$ is an integer and that as we defined in
*note Fourier transform::, $\omega{\equiv}m\omega_0$, where $m$ was an
integer. The integral will be zero for any $\omega$ that is not equal
to $2{\pi}n/P$, a more complete explanation can be seen in *note Fourier
series::. Therefore, while in the spatial domain the impulses had
spacing of $P$ (meters for example), in the frequency space, the spacing
between the different impulses are $2\pi/P$ cycles per meters.

File: gnuastro.info, Node: Convolution theorem, Next: Sampling theorem, Prev: Dirac delta and comb, Up: Frequency domain and Fourier operations

6.3.2.6 Convolution theorem
...........................

The convolution (shown with the $\ast$ operator) of the two functions
$f(l)$ and $h(l)$ is defined as:

$$c(l)\equiv[f{\ast}h](l)=\int_{-\infty}^{\infty}f(\tau)h(l-\tau)d\tau
$$

See *note Convolution process:: for a more detailed physical (pixel
based) interpretation of this definition. The Fourier transform of
convolution ($C(\omega)$) can be written as:

$$ C(\omega)=\int_{-\infty}^{\infty}[f{\ast}h](l)e^{-i{\omega}l}dl=
\int_{-\infty}^{\infty}f(\tau)\left[\int_{-\infty}^{\infty}h(l-\tau)e^{-i{\omega}l}dl\right]d\tau
$$

To solve the inner integral, let's define $s{\equiv}l-\tau$, so that
$ds=dl$ and $l=s+\tau$ then the inner integral becomes:

$$\int_{-\infty}^{\infty}h(l-\tau)e^{-i{\omega}l}dl=
\int_{-\infty}^{\infty}h(s)e^{-i{\omega}(s+\tau)}ds=e^{-i{\omega}\tau}\int_{-\infty}^{\infty}h(s)e^{-i{\omega}s}ds=H(\omega)e^{-i{\omega}\tau}
$$

where $H(\omega)$ is the Fourier transform of $h(l)$. Substituting this
result for the inner integral above, we get:

$$C(\omega)=H(\omega)\int_{-\infty}^{\infty}f(\tau)e^{-i{\omega}\tau}d\tau=H(\omega)F(\omega)=F(\omega)H(\omega)
$$

where $F(\omega)$ is the Fourier transform of $f(l)$. So multiplying
the Fourier transform of two functions individually, we get the Fourier
transform of their convolution. The convolution theorem also proves a
relation between the convolutions in the frequency space. Let's define:

$$D(\omega){\equiv}F(\omega){\ast}H(\omega)$$

Applying the inverse Fourier Transform or synthesis equation (*note
Fourier transform::) to both sides and following the same steps above,
we get:

$$d(l)=f(l)h(l)$$

Where $d(l)$ is the inverse Fourier transform of $D(\omega)$. We can
therefore re-write the two equations above formally as the convolution
theorem:

$$ {\cal F}[f{\ast}h]={\cal F}[f]{\cal F}[h]
$$

$$ {\cal F}[fh]={\cal F}[f]\ast{\cal F}[h]
$$

Besides its usefulness in blurring an image by convolving it with a
given kernel, the convolution theorem also enables us to do another very
useful operation in data analysis: to match the blur (or PSF) between
two images taken with different telescopes/cameras or under different
atmospheric conditions. This process is also known as deconvolution.
Let's take $f(l)$ as the image with a narrower PSF (less blurry) and
$c(l)$ as the image with a wider PSF which appears more blurred. Also
let's take $h(l)$ to represent the kernel that should be convolved with
the sharper image to create the more blurry image. Above, we proved the
relation between these three images through the convolution theorem.
But there, we assumed that $f(l)$ and $h(l)$ are known (given) and the
convolved image is desired.

In deconvolution, we have $f(l)$ -the sharper image- and $f*h(l)$
-the more blurry image- and we want to find the kernel $h(l)$. The
solution is a direct result of the convolution theorem:

$$ {\cal F}[h]={{\cal F}[f{\ast}h]\over {\cal F}[f]}
\quad\quad
{\rm or}
\quad\quad
h(l)={\cal F}^{-1}\left[{{\cal F}[f{\ast}h]\over {\cal F}[f]}\right]
$$

While this works really nice, it has two problems:

• If ${\cal F}[f]$ has any zero values, then the inverse Fourier
transform will not be a number!

• If there is significant noise in the image, then the high
frequencies of the noise are going to significantly reduce the
quality of the final result.

A standard solution to both these problems is the Weiner
deconvolution algorithm(1).

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

(1) <https://en.wikipedia.org/wiki/Wiener_deconvolution>

File: gnuastro.info, Node: Sampling theorem, Next: Discrete Fourier transform, Prev: Convolution theorem, Up: Frequency domain and Fourier operations

6.3.2.7 Sampling theorem
........................

Our mathematical functions are continuous, however, our data collecting
and measuring tools are discrete. Here we want to give a mathematical
formulation for digitizing the continuous mathematical functions so that
later, we can retrieve the continuous function from the digitized
recorded input. Assuming that we have a continuous function $f(l)$,
then we can define $f_s(l)$ as the 'sampled' $f(l)$ through the Dirac
comb (see *note Dirac delta and comb::):

$$f_s(l)=f(l){\rm III}_P=\displaystyle\sum_{n=-\infty}^{\infty}f(l)\delta(l-nP)
$$

The discrete data-element $f_k$ (for example, a pixel in an image),
where $k$ is an integer, can thus be represented as:

$$f_k=\int_{-\infty}^{\infty}f_s(l)dl=\int_{-\infty}^{\infty}f(l)\delta(l-kP)dt=f(kP)$$

Note that in practice, our discrete data points are not found in this
fashion. Each detector pixel (in an image for example) has an area and
averages the signal it receives over that area, not a mathematical point
as the Dirac $\delta$ function defines. However, as long as the
variation in the signal over one detector pixel is not significant, this
can be a good approximation. Having put this issue to the side, we can
now try to find the relation between the Fourier transforms of the
un-sampled $f(l)$ and the sampled $f_s(l)$. For a more clear notation,
let's define:

$$F_s(\omega)\equiv{\cal F}[f_s]$$

$$D(\omega)\equiv{\cal F}[{\rm III}_P]$$

Then using the Convolution theorem (see *note Convolution theorem::),
$F_s(\omega)$ can be written as:

$$F_s(\omega)={\cal F}[f(l){\rm III}_P]=F(\omega){\ast}D(\omega)$$

Finally, from the definition of convolution and the Fourier transform of
the Dirac comb (see *note Dirac delta and comb::), we get:

$$\eqalign{
F_s(\omega) &= \int_{-\infty}^{\infty}F(\omega)D(\omega-\mu)d\mu \cr
&= {1\over P}\displaystyle\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}F(\omega)\delta\left(\omega-\mu-{2{\pi}n\over P}\right)d\mu \cr
&= {1\over P}\displaystyle\sum_{n=-\infty}^{\infty}F\left(
\omega-{2{\pi}n\over P}\right).\cr }
$$

$F(\omega)$ was only a simple function, see *note Figure 6.3:
samplingfreq.(left). However, from the sampled Fourier transform
function we see that $F_s(\omega)$ is the superposition of infinite
copies of $F(\omega)$ that have been shifted, see *note Figure 6.3:
samplingfreq.(right). From the equation, it is clear that the shift in
each copy is $2\pi/P$.

[image src="gnuastro-figures/samplingfreq.png" text="../gnuastro-figures//samplingfreq.eps"]

Figure 6.3: Sampling causes infinite repetition in the frequency domain.
FT is an abbreviation for 'Fourier transform'. $\omega_m$ represents
the maximum frequency present in the input. $F(\omega)$ is only
symmetric on both sides of 0 when the input is real (not complex). In
general $F(\omega)$ is complex and thus cannot be simply plotted like
this. Here we have assumed a real Gaussian $f(t)$ which has produced a
Gaussian $F(\omega)$.

The input $f(l)$ can have any distribution of frequencies in it. In
the example of *note Figure 6.3: samplingfreq.(left), the input
consisted of a range of frequencies equal to $\Delta\omega=2\omega_m$.
Fortunately as *note Figure 6.3: samplingfreq.(right) shows, the assumed
pixel size ($P$) we used to sample this hypothetical function was such
that $2\pi/P>\Delta\omega$. The consequence is that each copy of
$F(\omega)$ has become completely separate from the surrounding copies.
Such a digitized (sampled) data set is thus called _over-sampled_. When
$2\pi/P=\Delta\omega$, $P$ is just small enough to finely separate even
the largest frequencies in the input signal and thus it is known as
_critically-sampled_. Finally if $2\pi/P<\Delta\omega$ we are dealing
with an _under-sampled_ data set. In an under-sampled data set, the
separate copies of $F(\omega)$ are going to overlap and this will
deprive us of recovering high constituent frequencies of $f(l)$. The
effects of under-sampling in an image with high rates of change (for
example, a brick wall imaged from a distance) can clearly be visually
seen and is known as _aliasing_.

When the input $f(l)$ is composed of a finite range of frequencies,
$f(l)$ is known as a _band-limited_ function. The example in *note
Figure 6.3: samplingfreq.(left) was a nice demonstration of such a case:
for all $\omega<-\omega_m$ or $\omega>\omega_m$, we have $F(\omega)=0$.
Therefore, when the input function is band-limited and our detector's
pixels are placed such that we have critically (or over-) sampled it,
then we can exactly reproduce the continuous $f(l)$ from the discrete or
digitized samples. To do that, we just have to isolate one copy of
$F(\omega)$ from the infinite copies and take its inverse Fourier
transform.

This ability to exactly reproduce the continuous input from the
sampled or digitized data leads us to the _sampling theorem_ which
connects the inherent property of the continuous signal (its maximum
frequency) to that of the detector (the spacing between its pixels).
The sampling theorem states that the full (continuous) signal can be
recovered when the pixel size ($P$) and the maximum constituent
frequency in the signal ($\omega_m$) have the following relation(1):

$${2\pi\over P}>2\omega_m$$

This relation was first formulated by Harry Nyquist (1889 - 1976 A.D.)
in 1928 and formally proved in 1949 by Claude E. Shannon (1916 - 2001
A.D.) in what is now known as the Nyquist-Shannon sampling theorem. In
signal processing, the signal is produced (synthesized) by a transmitter
and is received and de-coded (analyzed) by a receiver. Therefore
producing a band-limited signal is necessary.

In astronomy, we do not produce the shapes of our targets, we are
only observers. Galaxies can have any shape and size, therefore
ideally, our signal is not band-limited. However, since we are always
confined to observing through an aperture, the aperture will cause a
point source (for which $\omega_m=\infty$) to be spread over several
pixels. This spread is quantitatively known as the point spread
function or PSF. This spread does blur the image which is undesirable;
however, for this analysis it produces the positive outcome that there
will be a finite $\omega_m$. Though we should caution that any detector
will have noise which will add lots of very high frequency (ideally
infinite) changes between the pixels. However, the coefficients of
those noise frequencies are usually exceedingly small.

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

(1) This equation is also shown in some places without the $2\pi$.
Whether $2\pi$ is included or not depends on how you define the
frequency

File: gnuastro.info, Node: Discrete Fourier transform, Next: Fourier operations in two dimensions, Prev: Sampling theorem, Up: Frequency domain and Fourier operations

6.3.2.8 Discrete Fourier transform
..................................

As we have stated several times so far, the input image is a digitized,
pixelated or discrete array of values ($f_s(l)$, see *note Sampling
theorem::). The input is not a continuous function. Also, all our
numerical calculations can only be done on a sampled, or discrete
Fourier transform. Note that $F_s(\omega)$ is not discrete, it is
continuous. One way would be to find the analytic $F_s(\omega)$, then
sample it at any desired "freq-pixel"(1) spacing. However, this process
would involve two steps of operations and computers in particular are
not too good at analytic operations for the first step. So here, we
will derive a method to directly find the 'freq-pixel'ated $F_s(\omega)$
from the pixelated $f_s(l)$. Let's start with the definition of the
Fourier transform (see *note Fourier transform::):

$$F_s(\omega)=\int_{-\infty}^{\infty}f_s(l)e^{-i{\omega}l}dl $$

From the definition of $f_s(\omega)$ (using $x$ instead of $n$) we get:

$$\eqalign{
F_s(\omega) &= \displaystyle\sum_{x=-\infty}^{\infty}
\int_{-\infty}^{\infty}f(l)\delta(l-xP)e^{-i{\omega}l}dl \cr
&= \displaystyle\sum_{x=-\infty}^{\infty}
f_xe^{-i{\omega}xP}
}
$$

Where $f_x$ is the value of $f(l)$ on the point $x$ or the value of the
$x$th pixel. As shown in *note Sampling theorem:: this function is
infinitely periodic with a period of $2\pi/P$. So all we need is the
values within one period: $0<\omega<2\pi/P$, see *note Figure 6.3:
samplingfreq. We want $X$ samples within this interval, so the
frequency difference between each frequency sample or freq-pixel is
$1/XP$. Hence we will evaluate the equation above on the points at:

$$\omega={u\over XP} \quad\quad u = 0, 1, 2, ..., X-1$$

Therefore the value of the freq-pixel $u$ in the frequency domain is:

$$F_u=\displaystyle\sum_{x=0}^{X-1} f_xe^{-i{ux\over X}} $$

Therefore, we see that for each freq-pixel in the frequency domain, we
are going to need all the pixels in the spatial domain(2). If the input
(spatial) pixel row is also $X$ pixels wide, then we can exactly recover
the $x$th pixel with the following summation:

$$f_x={1\over X}\displaystyle\sum_{u=0}^{X-1} F_ue^{i{ux\over X}} $$

When the input pixel row (we are still only working on 1D data) has
$X$ pixels, then it is $L=XP$ spatial units wide. $L$, or the length of
the input data was defined in *note Fourier series:: and $P$ or the
space between the pixels in the input was defined in *note Dirac delta
and comb::. As we saw in *note Sampling theorem::, the input (spatial)
pixel spacing ($P$) specifies the range of frequencies that can be
studied and in *note Fourier series:: we saw that the length of the
(spatial) input, ($L$) determines the resolution (or size of the
freq-pixels) in our discrete Fourier transformed image. Both result
from the fact that the frequency domain is the inverse of the spatial
domain.

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

(1) We are using the made-up word "freq-pixel" so they are not
confused with spatial domain "pixels".

(2) So even if one pixel is a blank pixel (see *note Blank pixels::),
all the pixels in the frequency domain will also be blank.

File: gnuastro.info, Node: Fourier operations in two dimensions, Next: Edges in the frequency domain, Prev: Discrete Fourier transform, Up: Frequency domain and Fourier operations

6.3.2.9 Fourier operations in two dimensions
............................................

Once all the relations in the previous sections have been clearly
understood in one dimension, it is very easy to generalize them to two
or even more dimensions since each dimension is by definition
independent. Previously we defined $l$ as the continuous variable in 1D
and the inverse of the period in its direction to be $\omega$. Let's
show the second spatial direction with $m$ the inverse of the period in
the second dimension with $\nu$. The Fourier transform in 2D (see *note
Fourier transform::) can be written as:

$$F(\omega, \nu)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
f(l, m)e^{-i({\omega}l+{\nu}m)}dl$$

$$f(l, m)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
F(\omega, \nu)e^{i({\omega}l+{\nu}m)}dl$$

The 2D Dirac $\delta(l,m)$ is non-zero only when $l=m=0$. The 2D
Dirac comb (or Dirac brush! See *note Dirac delta and comb::) can be
written in units of the 2D Dirac $\delta$. For most image detectors,
the sides of a pixel are equal in both dimensions. So $P$ remains
unchanged, if a specific device is used which has non-square pixels,
then for each dimension a different value should be used.

$${\rm III}_P(l, m)\equiv\displaystyle\sum_{j=-\infty}^{\infty}
\displaystyle\sum_{k=-\infty}^{\infty}
\delta(l-jP, m-kP) $$

The Two dimensional Sampling theorem (see *note Sampling theorem::)
is thus very easily derived as before since the frequencies in each
dimension are independent. Let's take $\nu_m$ as the maximum frequency
along the second dimension. Therefore the two dimensional sampling
theorem says that a 2D band-limited function can be recovered when the
following conditions hold(1):

$${2\pi\over P} > 2\omega_m \quad\quad\quad {\rm and}
\quad\quad\quad {2\pi\over P} > 2\nu_m$$

Finally, let's represent the pixel counter on the second dimension in
the spatial and frequency domains with $y$ and $v$ respectively. Also
let's assume that the input image has $Y$ pixels on the second
dimension. Then the two dimensional discrete Fourier transform and its
inverse (see *note Discrete Fourier transform::) can be written as:

$$F_{u,v}=\displaystyle\sum_{x=0}^{X-1}\displaystyle\sum_{y=0}^{Y-1}
f_{x,y}e^{-i({ux\over X}+{vy\over Y})} $$

$$f_{x,y}={1\over XY}\displaystyle\sum_{u=0}^{X-1}\displaystyle\sum_{v=0}^{Y-1}
F_{u,v}e^{i({ux\over X}+{vy\over Y})} $$

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

(1) If the pixels are not a square, then each dimension has to use
the respective pixel size, but since most detectors have square pixels,
we assume so here too

File: gnuastro.info, Node: Edges in the frequency domain, Prev: Fourier operations in two dimensions, Up: Frequency domain and Fourier operations

6.3.2.10 Edges in the frequency domain
......................................

With a good grasp of the frequency domain, we can revisit the problem of
convolution on the image edges, see *note Edges in the spatial domain::.
When we apply the convolution theorem (see *note Convolution theorem::)
to convolve an image, we first take the discrete Fourier transforms
(DFT, *note Discrete Fourier transform::) of both the input image and
the kernel, then we multiply them with each other and then take the
inverse DFT to construct the convolved image. Of course, in order to
multiply them with each other in the frequency domain, the two images
have to be the same size, so let's assume that we pad the kernel (it is
usually smaller than the input image) with zero valued pixels in both
dimensions so it becomes the same size as the input image before the
DFT.

Having multiplied the two DFTs, we now apply the inverse DFT which is
where the problem is usually created. If the DFT of the kernel only had
values of 1 (unrealistic condition!) then there would be no problem and
the inverse DFT of the multiplication would be identical with the input.
However in real situations, the kernel's DFT has a maximum of 1 (because
the sum of the kernel has to be one, see *note Convolution process::)
and decreases something like the hypothetical profile of *note Figure
6.3: samplingfreq. So when multiplied with the input image's DFT, the
coefficients or magnitudes (see *note Circles and the complex plane::)
of the smallest frequency (or the sum of the input image pixels) remains
unchanged, while the magnitudes of the higher frequencies are
significantly reduced.

As we saw in *note Sampling theorem::, the Fourier transform of a
discrete input will be infinitely repeated. In the final inverse DFT
step, the input is in the frequency domain (the multiplied DFT of the
input image and the kernel DFT). So the result (our output convolved
image) will be infinitely repeated in the spatial domain. In order to
accurately reconstruct the input image, we need all the frequencies with
the correct magnitudes. However, when the magnitudes of higher
frequencies are decreased, longer periods (shorter frequencies) will
dominate in the reconstructed pixel values. Therefore, when
constructing a pixel on the edge of the image, the newly empowered
longer periods will look beyond the input image edges and will find the
repeated input image there. So if you convolve an image in this fashion
using the convolution theorem, when a bright object exists on one edge
of the image, its blurred wings will be present on the other side of the
convolved image. This is often termed as circular convolution or cyclic
convolution.

So, as long as we are dealing with convolution in the frequency
domain, there is nothing we can do about the image edges. The least we
can do is to eliminate the ghosts of the other side of the image. So,
we add zero valued pixels to both the input image and the kernel in both
dimensions so the image that will be convolved has a size equal to the
sum of both images in each dimension. Of course, the effect of this
zero-padding is that the sides of the output convolved image will become
dark. To put it another way, the edges are going to drain the flux from
nearby objects. But at least it is consistent across all the edges of
the image and is predictable. In Convolve, you can see the padded
images when inspecting the frequency domain convolution steps with the
‘--viewfreqsteps’ option.

File: gnuastro.info, Node: Spatial vs. Frequency domain, Next: Convolution kernel, Prev: Frequency domain and Fourier operations, Up: Convolve

6.3.3 Spatial vs. Frequency domain
----------------------------------

With the discussions above it might not be clear when to choose the
spatial domain and when to choose the frequency domain. Here we will
try to list the benefits of each.

The spatial domain,
• Can correct for the edge effects of convolution, see *note Edges in
the spatial domain::.

• Can operate on blank pixels.

• Can be faster than frequency domain when the kernel is small (in
terms of the number of pixels on the sides).

The frequency domain,
• Will be much faster when the image and kernel are both large.

As a general rule of thumb, when working on an image of modeled profiles
use the frequency domain and when working on an image of real (observed)
objects use the spatial domain (corrected for the edges). The reason is
that if you apply a frequency domain convolution to a real image, you
are going to loose information on the edges and generally you do not
want large kernels. But when you have made the profiles in the image
yourself, you can just make a larger input image and crop the central
parts to completely remove the edge effect, see *note If convolving
afterwards::. Also due to oversampling, both the kernels and the images
can become very large and the speed boost of frequency domain
convolution will significantly improve the processing time, see *note
Oversampling::.

File: gnuastro.info, Node: Convolution kernel, Next: Invoking astconvolve, Prev: Spatial vs. Frequency domain, Up: Convolve

6.3.4 Convolution kernel
------------------------

All the programs that need convolution will need to be given a
convolution kernel file and extension. In most cases (other than
Convolve, see *note Convolve::) the kernel file name is optional.
However, the extension is necessary and must be specified either on the
command-line or at least one of the configuration files (see *note
Configuration files::). Within Gnuastro, there are two ways to create a
kernel image:

• MakeProfiles: You can use MakeProfiles to create a parametric
(based on a radial function) kernel, see *note MakeProfiles::. By
default MakeProfiles will make the Gaussian and Moffat profiles in
a separate file so you can feed it into any of the programs.

• ConvertType: You can write your own desired kernel into a text file
table and convert it to a FITS file with ConvertType, see *note
ConvertType::. Just be careful that the kernel has to have an odd
number of pixels along its two axes, see *note Convolution
process::. All the programs that do convolution will normalize the
kernel internally, so if you choose this option, you do not have to
worry about normalizing the kernel. Only within Convolve, there is
an option to disable normalization, see *note Invoking
astconvolve::.

The two options to specify a kernel file name and its extension are
shown below. These are common between all the programs that will do
convolution.
‘-k FITS’
‘--kernel=FITS’
The convolution kernel file name. The ‘BITPIX’ (data type) value
of this file can be any standard type and it does not necessarily
have to be normalized. Several operations will be done on the
kernel image prior to the program's processing:

• It will be converted to floating point type.

• All blank pixels (see *note Blank pixels::) will be set to
zero.

• It will be normalized so the sum of its pixels equal unity.

• It will be flipped so the convolved image has the same
orientation. This is only relevant if the kernel is not
circular. See *note Convolution process::.

‘-U STR’
‘--khdu=STR’
The convolution kernel HDU. Although the kernel file name is
optional, before running any of the programs, they need to have a
value for ‘--khdu’ even if the default kernel is to be used. So be
sure to keep its value in at least one of the configuration files
(see *note Configuration files::). By default, the system
configuration file has a value.

File: gnuastro.info, Node: Invoking astconvolve, Prev: Convolution kernel, Up: Convolve

6.3.5 Invoking Convolve
-----------------------

Convolve an input dataset (2D image or 1D spectrum for example) with a
known kernel, or make the kernel necessary to match two PSFs. The
general template for Convolve is:

$ astconvolve [OPTION...] ASTRdata

One line examples:

## Convolve mockimg.fits with psf.fits:
$ astconvolve --kernel=psf.fits mockimg.fits

## Convolve in the spatial domain:
$ astconvolve observedimg.fits --kernel=psf.fits --domain=spatial

## Convolve a 3D cube (only spatial domain is supported in 3D).
## It is also necessary to define 3D tiles and channels for
## parallelization (see the Tessellation section for more).
$ astconvolve cube.fits --kernel=kernel3d.fits --domain=spatial \
--tilesize=30,30,30 --numchannels=1,1,1

## Find the kernel to match sharper and blurry PSF images (they both
## have to have the same pixel size).
$ astconvolve --kernel=sharperimage.fits --makekernel=10 \
blurryimage.fits

## Convolve a Spectrum (column 14 in the FITS table below) with a
## custom kernel (the kernel will be normalized internally, so only
## the ratios are important). Sed is used to replace the spaces with
## new line characters so Convolve sees them as values in one column.
$ echo "1 3 10 3 1" | sed 's/ /\n/g' | astconvolve spectra.fits -c14

The only argument accepted by Convolve is an input image file. Some
of the options are the same between Convolve and some other Gnuastro
programs. Therefore, to avoid repetition, they will not be repeated
here. For the full list of options shared by all Gnuastro programs,
please see *note Common options::. In particular, in the spatial
domain, on a multi-dimensional datasets, convolve uses Gnuastro's
tessellation to speed up the run, see *note Tessellation::. Common
options related to tessellation are described in *note Processing
options::.

1-dimensional datasets (for example, spectra) are only read as
columns within a table (see *note Tables:: for more on how Gnuastro
programs read tables). Note that currently 1D convolution is only
implemented in the spatial domain and thus kernel-matching is also not
supported.

Here we will only explain the options particular to Convolve. Run
Convolve with ‘--help’ in order to see the full list of options Convolve
accepts, irrespective of where they are explained in this book.

‘--kernelcolumn’
Column containing the 1D kernel. When the input dataset is a
1-dimensional column, and the host table has more than one column,
use this option to specify which column should be used.

‘--nokernelflip’
Do not flip the kernel after reading; only for spatial domain
convolution. This can be useful if the flipping has already been
applied to the kernel. By default, the input kernel is flipped to
avoid the output getting flipped; see *note Convolution process::.

‘--nokernelnorm’
Do not normalize the kernel after reading it, such that the sum of
its pixels is unity. As described in *note Convolution process::,
the kernel is normalized by default.

‘--conv-on-blank’
Do not ignore blank pixels in the convolution. The output pixels
that were originally non-blank are not affected by this option
(they will have the same value if this option is called or not).
This option just expands/dilates the non-blank regions of your
dataset into the blank regions and only works in spatial domain
convolution. Therefore, with this option convolution can be used
as a proxy for interpolation or dilation.

By default, blank pixels are ignored during spatial domain
convolution; so the input and output have exactly the same number
of blank pixels. With this option, the blank pixels that are
sufficiently close to non-blank pixels (based on the kernel) will
be given a value based on the non-blank elements that overlap with
the kernel for that blank pixel (see *note Edges in the spatial
domain::).

‘-d STR’
‘--domain=STR’
The domain to use for the convolution. The acceptable values are
'‘spatial’' and '‘frequency’', corresponding to the respective
domain.

For large images, the frequency domain process will be more
efficient than convolving in the spatial domain. However, the
edges of the image will loose some flux (see *note Edges in the
spatial domain::) and the image must not contain any blank pixels,
see *note Spatial vs. Frequency domain::.

‘--checkfreqsteps’
With this option a file with the initial name of the output file
will be created that is suffixed with ‘_freqsteps.fits’, all the
steps done to arrive at the final convolved image are saved as
extensions in this file. The extensions in order are:

1. The padded input image. In frequency domain convolution the
two images (input and convolved) have to be the same size and
both should be padded by zeros.

2. The padded kernel, similar to the above.

3. The Fourier spectrum of the forward Fourier transform of the
input image. Note that the Fourier transform is a complex
operation (and not view able in one image!) So we either have
to show the 'Fourier spectrum' or the 'Phase angle'. For the
complex number $a+ib$, the Fourier spectrum is defined as
$\sqrt{a^2+b^2}$ while the phase angle is defined as
$\arctan(b/a)$.

4. The Fourier spectrum of the forward Fourier transform of the
kernel image.

5. The Fourier spectrum of the multiplied (through complex
arithmetic) transformed images.

6. The inverse Fourier transform of the multiplied image. If you
open it, you will see that the convolved image is now in the
center, not on one side of the image as it started with (in
the padded image of the first extension). If you are working
on a mock image which originally had pixels of precisely 0.0,
you will notice that in those parts that your convolved
profile(s) did not convert, the values are now $\sim10^{-18}$,
this is due to floating-point round off errors. Therefore in
the final step (when cropping the central parts of the image),
we also remove any pixel with a value less than $10^{-17}$.

‘--noedgecorrection’
Do not correct the edge effect in spatial domain convolution (this
correction is done in spatial domain convolution by default). For
a full discussion, please see *note Edges in the spatial domain::.

‘-m INT’
‘--makekernel=INT’
If this option is called, Convolve will do PSF-matching: the output
will be the kernel that you should convolve with the sharper image
to obtain the blurry one (see *note Convolution theorem::). The
two images must have the same size (number of pixels). This option
is not yet supported in 1-dimensional datasets. In effect, it is
only necessary to give the two PSFs of your two datasets, find the
matching kernel based on that, then apply that kernel to the
higher-resolution (sharper image).

The image given to the ‘--kernel’ option is assumed to be the
sharper (less blurry) image and the input image (with no option) is
assumed to be the more blurry image. The value given to this
option will be used as the maximum radius of the kernel. Any pixel
in the final kernel that is larger than this distance from the
center will be set to zero.

Noise has large frequencies which can make the result less reliable
for the higher frequencies of the final result. So all the
frequencies which have a spectrum smaller than the value given to
the ‘minsharpspec’ option in the sharper input image are set to
zero and not divided. This will cause the wings of the final
kernel to be flatter than they would ideally be which will make the
convolved image result unreliable if it is too high.

Some notes to on how to prepare your two input PSFs. Note that
these (and several other issues that relate to an accurate
estimation of the PSF) are practically described in the following
tutorial: *note Building the extended PSF::.

• Choose a bright (unsaturated) star and use a region box (with
Crop for example, see *note Crop::) that is sufficiently above
the noise.

• Mask all background sources that may be nearby (you can use
Segment's clumps, see *note Segment::).

• Use Warp (see *note Warp::) to warp the pixel grid so the
star's center is exactly on the center of the central pixel in
the cropped image. This will certainly slightly degrade the
result, however, it is necessary. If there are multiple good
stars, you can shift all of them, then normalize them (so the
sum of each star's pixels is one) and then take their average
to decrease this effect.

• The shifting might move the center of the star by one pixel in
any direction, so crop the central pixel of the warped image
to have a clean image for the deconvolution.

‘-c’
‘--minsharpspec’
(‘=FLT’) The minimum frequency spectrum (or coefficient, or pixel
value in the frequency domain image) to use in deconvolution, see
the explanations under the ‘--makekernel’ option for more
information.

File: gnuastro.info, Node: Warp, Prev: Convolve, Up: Data manipulation

6.4 Warp
========

Image warping is the process of mapping the pixels of one image onto a
new pixel grid. This process is sometimes known as transformation,
however following the discussion of Heckbert 1989(1) we will not be
using that term because it can be confused with only pixel value or flux
transformations. Here we specifically mean the pixel grid
transformation which is better conveyed with 'warp'.

Image warping is a very important step in astronomy, both in
observational data analysis and in simulating modeled images. In
modeling, warping an image is necessary when we want to apply grid
transformations to the initial models, for example, in simulating
gravitational lensing. Observational reasons for warping an image are
listed below:

• *Noise:* Most scientifically interesting targets are inherently
faint (have a very low Signal to noise ratio). Therefore one short
exposure is not enough to detect such objects that are drowned
deeply in the noise. We need multiple exposures so we can add them
together and increase the objects' signal to noise ratio. Keeping
the telescope fixed on one field of the sky is practically
impossible. Therefore very deep observations have to put into the
same grid before adding them.

• *Resolution:* If we have multiple images of one patch of the sky
(hopefully at multiple orientations) we can warp them to the same
grid. The multiple orientations will allow us to 'guess' the
values of pixels on an output pixel grid that has smaller pixel
sizes and thus increase the resolution of the output. This process
of merging multiple observations is known as Mosaicing.

• *Cosmic rays:* Cosmic rays can randomly fall on any part of an
image. If they collide vertically with the camera, they are going
to create a very sharp and bright spot that in most cases can be
separated easily(2). However, depending on the depth of the camera
pixels, and the angle that a cosmic rays collides with it, it can
cover a line-like larger area on the CCD which makes the detection
using their sharp edges very hard and error prone. One of the best
methods to remove cosmic rays is to compare multiple images of the
same field. To do that, we need all the images to be on the same
pixel grid.

• *Optical distortion:* In wide field images, the optical distortion
that occurs on the outer parts of the focal plane will make
accurate comparison of the objects at various locations impossible.
It is therefore necessary to warp the image and correct for those
distortions prior to the analysis.

• *Detector not on focal plane:* In some cases (like the Hubble Space
Telescope ACS and WFC3 cameras), the CCD might be tilted compared
to the focal plane, therefore the recorded CCD pixels have to be
projected onto the focal plane before further analysis.

* Menu:

* Linear warping basics:: Basics of coordinate transformation.
* Merging multiple warpings:: How to merge multiple matrices.
* Resampling:: Warping an image is re-sampling it.
* Invoking astwarp:: Arguments and options for Warp.

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

(1) Paul S. Heckbert. 1989. _Fundamentals of Texture mapping and
Image Warping_, Master's thesis at University of California, Berkeley.

(2) All astronomical targets are blurred with the PSF, see *note
PSF::, however a cosmic ray is not and so it is very sharp (it suddenly
stops at one pixel).

File: gnuastro.info, Node: Linear warping basics, Next: Merging multiple warpings, Prev: Warp, Up: Warp

6.4.1 Linear warping basics
---------------------------

Let's take $\left[\matrix{u&v}\right]$ as the coordinates of a point in
the input image and $\left[\matrix{x&y}\right]$ as the coordinates of
that same point in the output image(1). The simplest form of coordinate
transformation (or warping) is the scaling of the coordinates, let's
assume we want to scale the first axis by $M$ and the second by $N$, the
output coordinates of that point can be calculated by

$$\left[\matrix{x\cr y}\right]=
\left[\matrix{Mu\cr Nv}\right]=
\left[\matrix{M&0\cr0&N}\right]\left[\matrix{u\cr v}\right]$$

Note that these are matrix multiplications. We thus see that we can
represent any such grid warping as a matrix. Another thing we can do
with this $2\times2$ matrix is to rotate the output coordinate around
the common center of both coordinates. If the output is rotated
anticlockwise by $\theta$ degrees from the positive (to the right)
horizontal axis, then the warping matrix should become:

$$\left[\matrix{x\cr y}\right]=
\left[\matrix{ucos\theta-vsin\theta\cr usin\theta+vcos\theta}\right]=
\left[\matrix{cos\theta&-sin\theta\cr sin\theta&cos\theta}\right]
\left[\matrix{u\cr v}\right]
$$

We can also flip the coordinates around the first axis, the second axis
and the coordinate center with the following three matrices
respectively:

$$\left[\matrix{1&0\cr0&-1}\right]\quad\quad
\left[\matrix{-1&0\cr0&1}\right]\quad\quad
\left[\matrix{-1&0\cr0&-1}\right]$$

The final thing we can do with this definition of a $2\times2$ warping
matrix is shear. If we want the output to be sheared along the first
axis with $A$ and along the second with $B$, then we can use the matrix:

$$\left[\matrix{1&A\cr B&1}\right]$$

To have one matrix representing any combination of these steps, you use
matrix multiplication, see *note Merging multiple warpings::. So any
combinations of these transformations can be displayed with one
$2\times2$ matrix:

$$\left[\matrix{a&b\cr c&d}\right]$$

The transformations above can cover a lot of the needs of most
coordinate transformations. However they are limited to mapping the
point $[\matrix{0&0}]$ to $[\matrix{0&0}]$. Therefore they are useless
if you want one coordinate to be shifted compared to the other one.
They are also space invariant, meaning that all the coordinates in the
image will receive the same transformation. In other words, all the
pixels in the output image will have the same area if placed over the
input image. So transformations which require varying output pixel
sizes like projections cannot be applied through this $2\times2$ matrix
either (for example, for the tilted ACS and WFC3 camera detectors on
board the Hubble space telescope).

To add these further capabilities, namely translation and projection,
we use the homogeneous coordinates. They were defined about 200 years
ago by August Ferdinand Möbius (1790 - 1868). For simplicity, we will
only discuss points on a 2D plane and avoid the complexities of higher
dimensions. We cannot provide a deep mathematical introduction here,
interested readers can get a more detailed explanation from Wikipedia(2)
and the references therein.

By adding an extra coordinate to a point we can add the flexibility
we need. The point $[\matrix{x&y}]$ can be represented as
$[\matrix{xZ&yZ&Z}]$ in homogeneous coordinates. Therefore multiplying
all the coordinates of a point in the homogeneous coordinates with a
constant will give the same point. Put another way, the point
$[\matrix{x&y&Z}]$ corresponds to the point $[\matrix{x/Z&y/Z}]$ on the
constant $Z$ plane. Setting $Z=1$, we get the input image plane, so
$[\matrix{u&v&1}]$ corresponds to $[\matrix{u&v}]$. With this
definition, the transformations above can be generally written as:

$$\left[\matrix{x\cr y\cr 1}\right]=
\left[\matrix{a&b&0\cr c&d&0\cr 0&0&1}\right]
\left[\matrix{u\cr v\cr 1}\right]$$

We thus acquired 4 extra degrees of freedom. By giving non-zero values
to the zero valued elements of the last column we can have translation
(try the matrix multiplication!). In general, any coordinate
transformation that is represented by the matrix below is known as an
affine transformation(3):

$$\left[\matrix{a&b&c\cr d&e&f\cr 0&0&1}\right]$$

We can now consider translation, but the affine transform is still
spatially invariant. Giving non-zero values to the other two elements
in the matrix above gives us the projective transformation or
Homography(4) which is the most general type of transformation with the
$3\times3$ matrix:

$$\left[\matrix{x'\cr y'\cr w}\right]=
\left[\matrix{a&b&c\cr d&e&f\cr g&h&1}\right]
\left[\matrix{u\cr v\cr 1}\right]$$

So the output coordinates can be calculated from:

$$x={x' \over w}={au+bv+c \over gu+hv+1}\quad\quad\quad\quad
y={y' \over w}={du+ev+f \over gu+hv+1}$$

Thus with Homography we can change the sizes of the output pixels on
the input plane, giving a 'perspective'-like visual impression. This
can be quantitatively seen in the two equations above. When $g=h=0$,
the denominator is independent of $u$ or $v$ and thus we have spatial
invariance. Homography preserves lines at all orientations. A very
useful fact about Homography is that its inverse is also a Homography.
These two properties play a very important role in the implementation of
this transformation. A short but instructive and illustrated review of
affine, projective and also bi-linear mappings is provided in Heckbert
1989(5).

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

(1) These can be any real number, we are not necessarily talking
about integer pixels here.

(2) <http://en.wikipedia.org/wiki/Homogeneous_coordinates>

(3) <http://en.wikipedia.org/wiki/Affine_transformation>

(4) <http://en.wikipedia.org/wiki/Homography>

(5) Paul S. Heckbert. 1989. _Fundamentals of Texture mapping and
Image Warping_, Master's thesis at University of California, Berkeley.
Note that since points are defined as row vectors there, the matrix is
the transpose of the one discussed here.

File: gnuastro.info, Node: Merging multiple warpings, Next: Resampling, Prev: Linear warping basics, Up: Warp

6.4.2 Merging multiple warpings
-------------------------------

In *note Linear warping basics:: we saw how a basic warp/transformation
can be represented with a matrix. To make more complex warpings (for
example, to define a translation, rotation and scale as one warp) the
individual matrices have to be multiplied through matrix multiplication.
However matrix multiplication is not commutative, so the order of the
set of matrices you use for the multiplication is going to be very
important.

The first warping should be placed as the left-most matrix. The
second warping to the right of that and so on. The second
transformation is going to occur on the warped coordinates of the first.
As an example for merging a few transforms into one matrix, the
multiplication below represents the rotation of an image about a point
$[\matrix{U&V}]$ anticlockwise from the horizontal axis by an angle of
$\theta$. To do this, first we take the origin to $[\matrix{U&V}]$
through translation. Then we rotate the image, then we translate it
back to where it was initially. These three operations can be merged in
one operation by calculating the matrix multiplication below:

$$\left[\matrix{1&0&U\cr0&1&V\cr{}0&0&1}\right]
\left[\matrix{cos\theta&-sin\theta&0\cr sin\theta&cos\theta&0\cr 0&0&1}\right]
\left[\matrix{1&0&-U\cr0&1&-V\cr{}0&0&1}\right]$$

File: gnuastro.info, Node: Resampling, Next: Invoking astwarp, Prev: Merging multiple warpings, Up: Warp

6.4.3 Resampling
----------------

A digital image is composed of discrete 'picture elements' or 'pixels'.
When a real image is created from a camera or detector, each pixel's
area is used to store the number of photo-electrons that were created
when incident photons collided with that pixel's surface area. This
process is called the 'sampling' of a continuous or analog data into
digital data.

When we change the pixel grid of an image, or "warp" it, we have to
calculate the flux value of each pixel on the new grid based on the old
grid, or resample it. Because of the calculation (as opposed to
observation), any form of warping on the data is going to degrade the
image and mix the original pixel values with each other. So if an
analysis can be done on an unwarped data image, it is best to leave the
image untouched and pursue the analysis. However as discussed in *note
Warp:: this is not possible in some scenarios and re-sampling is
necessary.

When the FWHM of the PSF of the camera is much larger than the pixel
scale (see *note Sampling theorem::) we are sampling the signal in a
much higher resolution than the camera can offer. This is usually the
case in many applications of image processing (nonastronomical imaging).
In such cases, we can consider each pixel to be a point and not an area:
the PSF doesn't vary much over a single pixel.

Approximating a pixel's area to a point can significantly speed up
the resampling and also the simplicity of the code. Because resampling
becomes a problem of interpolation: points of the input grid need to be
interpolated at certain other points (over the output grid). To
increase the accuracy, you might also sample more than one point from
within a pixel giving you more points for a more accurate interpolation
in the output grid.

However, interpolation has several problems. The first one is that
it will depend on the type of function you want to assume for the
interpolation. For example, you can choose a bi-linear or bi-cubic (the
'bi's are for the 2 dimensional nature of the data) interpolation
method. For the latter there are various ways to set the constants(1).
Such parametric interpolation functions can fail seriously on the edges
of an image, or when there is a sharp change in value (for example, the
bleeding saturation of bright stars in astronomical CCDs). They will
also need normalization so that the flux of the objects before and after
the warping is comparable.

The parametric nature of these methods adds a level of subjectivity
to the data (it makes more assumptions through the functions than the
data can handle). For most applications this is fine (as discussed
above: when the PSF is over-sampled), but in scientific applications
where we push our instruments to the limit and the aim is the detection
of the faintest possible galaxies or fainter parts of bright galaxies,
we cannot afford this loss. Because of these reasons Warp will not use
parametric interpolation techniques.

Warp will do interpolation based on "pixel mixing"(2) or "area
resampling". This is also similar to what the Hubble Space Telescope
pipeline calls "Drizzling"(3). This technique requires no functions, it
is thus non-parametric. It is also the closest we can get (make least
assumptions) to what actually happens on the detector pixels.

In pixel mixing, the basic idea is that you reverse-transform each
output pixel to find which pixels of the input image it covers, and what
fraction of the area of the input pixels are covered by that output
pixel. We then multiply each input pixel's value by the fraction of its
area that overlaps with the output pixel (between 0 to 1). The output's
pixel value is derived by summing all these multiplications for the
input pixels that it covers.

Through this process, pixels are treated as an area not as a point
(which is how detectors create the image), also the brightness (see
*note Brightness flux magnitude::) of an object will be fully preserved.
Since it involves the mixing of the input's pixel values, this pixel
mixing method is a form of *note Spatial domain convolution::.
Therefore, after comparing the input and output, you will notice that
the output is slightly smoothed, thus boosting the more diffuse signal,
but creating correlated noise. In astronomical imaging the correlated
noise will be decreased later when you stack many exposures(4).

If there are very high spatial-frequency signals in the image (for
example, fringes) which vary on a scale _smaller than_ your output image
pixel size (this is rarely the case in astronomical imaging), pixel
mixing can cause ailiasing(5). Therefore, in case such fringes are
present, they have to be calculated and removed separately (which would
naturally be done in any astronomical reduction pipeline). Because of
the PSF, no astronomical target has a sharp change in their signal.
Thus this issue is less important for astronomical applications, see
*note PSF::.

To find the overlap area of the output pixel over the input pixels,
we need to define polygons and clip them (find the overlap). Usually,
it is sufficient to define a pixel with a four-vertice polygon.
However, when a non-linear distortion (for example, ‘SIP’ or ‘TPV’) is
present and the distortion is significant over an output pixel's size
(usually far from the reference point), the shadow of the output pixel
on the input grid can be curved. To account for such cases (which can
only happen when correcting for non-linear distortions), Warp has the
‘--edgesampling’ option to sample the output pixel over more vertices.
For more, see the description of this option in *note Align pixels with
WCS considering distortions::.

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

(1) see <http://entropymine.com/imageworsener/bicubic/> for a nice
introduction.

(2) For a graphic demonstration see
<http://entropymine.com/imageworsener/pixelmixing/>.

(3) <http://en.wikipedia.org/wiki/Drizzle_(image_processing)>

(4) If you are working on a single exposure image and see pronounced
Moiré patterns after Warping, check *note Moire pattern in stacking and
its correction:: for a possible way to reduce them

(5) <http://en.wikipedia.org/wiki/Aliasing>

File: gnuastro.info, Node: Invoking astwarp, Prev: Resampling, Up: Warp

6.4.4 Invoking Warp
-------------------

Warp will warp an input image into a new pixel grid by pixel mixing (see
*note Resampling::). Without any options, Warp will remove any
non-linear distortions from the image and align the output pixel
coordinates to its WCS coordinates. Any homographic warp (for example,
scaling, rotation, translation, projection, see *note Linear warping
basics::) can also be done by calling the relevant option explicitly.
The general template for invoking Warp is:

$ astwarp [OPTIONS...] InputImage

One line examples:

## Align image with celestial coordinates and remove any distortion
$ astwarp image.fits

## Align four exposures to same pixel grid and stack them with
## Arithmetic program's sigma-clipped mean operator (out of many
## stacking operators, see Arithmetic's documentation).
$ grid="--center=1.234,5.678 --width=1001,1001 --widthinpix --cdelt=0.2/3600"
$ astwarp a.fits $grid --output=A.fits
$ astwarp b.fits $grid --output=B.fits
$ astwarp c.fits $grid --output=C.fits
$ astwarp d.fits $grid --output=D.fits
$ astarithmetic A.fits B.fits C.fits D.fits 4 5 0.2 sigclip-mean \
-g1 --output=stack.fits

## Warp a previously created mock image to the same pixel grid as the
## real image (including any distortions).
$ astwarp mock.fits --gridfile=real.fits

## Rotate and then scale input image:
$ astwarp --rotate=37.92 --scale=0.8 image.fits

## Scale, then translate the input image:
$ astwarp --scale 8/3 --translate 2.1 image.fits

## Directly input a custom warping matrix (using fraction):
$ astwarp --matrix=1/5,0,4/10,0,1/5,4/10,0,0,1 image.fits

## Directly input a custom warping matrix, with final numbers:
$ astwarp --matrix="0.7071,-0.7071, 0.7071,0.7071" image.fits

If any processing is to be done, Warp needs to be given a 2D FITS
image. As in all Gnuastro programs, when an output is not explicitly
set with the ‘--output’ option, the output filename will be set
automatically based on the operation, see *note Automatic output::. For
the full list of general options to all Gnuastro programs (including
Warp), please see *note Common options::.

Warp uses pixel mixing to derive the pixel values of the output
image, see *note Resampling::. To be the most accurate, the input image
will be read as a 64-bit double precision floating point dataset and all
internal processing is done in this format. Upon writing, by default it
will be converted to 32-bit single precision floating point type (actual
observational data rarely have such precision!). In case you want a
different output type, you can use the ‘--type’ option that is common to
several Gnuastro programs. For example, if your input is a mock image
without noise, and you want to preserve the 64-bit precision, use (with
‘--type=float64’. Just note that the file size will also be double!
For more on the precision of various types, see *note Numeric data
types::.

By default (if no linear operation is requested), Warp will align the
pixel grid of the input image to the WCS coordinates it contains. This
operation and the the options that govern it are described in *note
Align pixels with WCS considering distortions::. You can Warp an input
image to the same pixel grid as a reference FITS file using the
‘--wcsfile’ option. In this case, the output image will take all the
information needed from the reference WCS file and HDU/extension
specified with ‘--wcshdu’, thus it will discard any other resampling
options given.

If you need any custom linear warping (independent of the WCS, see
*note Linear warping basics::), you need to call the respective
operation manually. These are described in *note Linear warps to be
called explicitly::. Please note that you may not use both linear and
non-linear modes simultaneously. For example, you cannot scale or
rotate the image while removing its non-linear distortions at the same
time.

The following options are shared between both modes:

‘--hstartwcs=INT’
Specify the first header keyword number (line) that should be used
to read the WCS information, see the full explanation in *note
Invoking astcrop::.

‘--hendwcs=INT’
Specify the last header keyword number (line) that should be used
to read the WCS information, see the full explanation in *note
Invoking astcrop::.

‘-C FLT’
‘--coveredfrac=FLT’
Depending on the warp, the output pixels that cover pixels on the
edge of the input image, or blank pixels in the input image, are
not going to be fully covered by input data. With this option, you
can specify the acceptable covered fraction of such pixels (any
value between 0 and 1). If you only want output pixels that are
fully covered by the input image area (and are not blank), then you
can set ‘--coveredfrac=1’ (which is the default!). Alternatively,
a value of ‘0’ will keep output pixels that are even
infinitesimally covered by the input. As a result, with
‘--coveredfrac=0’, the sum of the pixels in the input and output
images will be exactly the same.

* Menu:

* Align pixels with WCS considering distortions:: Default operation.
* Linear warps to be called explicitly:: Other warps.

File: gnuastro.info, Node: Align pixels with WCS considering distortions, Next: Linear warps to be called explicitly, Prev: Invoking astwarp, Up: Invoking astwarp

6.4.4.1 Align pixels with WCS considering distortions
.....................................................

When none of the linear warps(1) are requested, Warp will align the
input's pixel axes with it's WCS axes. In the process, any possibly
existing distortion is also removed (such as ‘TPV’ and ‘SIP’). Usually,
the WCS axes are the Right Ascension and Declination in equatorial
coordinates. The output image's pixel grid is highly customizable
through the options in this section. To learn about Warp's strategy to
build the new pixel grid, see *note Resampling::. For strong
distortions (that produce strong curvatures), you can fine-tune the
area-based resampling with ‘--edgesampling’, as described below.

On the other hand, sometimes you need to Warp an input image to the
exact same grid of an already available reference FITS image with an
existing WCS. If that image is already aligned, finding its center,
number of pixels and pixel scale can be annoying (and just increase the
complexity of your script). On the other hand, if that image is not
aligned (for example, has a certain rotation in the sky, and has a
different distortion), there are too many WCS parameters to set (some
are not yet available explicitly in the options here)! For such
scenarios, Warp has the ‘--gridfile’ option. When ‘--gridfile’ is
called, the options below that are used to define the output's WCS will
be ignored (these options: ‘--center’, ‘--widthinpix’, ‘--cdelt’,
‘--ctype’). In this case, the output's WCS and pixel grid will exactly
match the image given to ‘--gridfile’ (including any rotation, pixel
scale, or distortion or projection).

*Set ‘--cdelt’ explicitly when you plan to stack many warped images:* To
align some images and later stack them, it is necessary to be sure the
pixel sizes of all the images are the same exactly. Most of the time
the measured (during astrometry) pixel scale of the separate exposures,
will be different in the second or third digit number after the decimal
point. It is a normal/statistical error in measuring the astrometry.
On a large image, these slight differences can cause different output
sizes (of one or two pixels on a very large image).

You can fix this by explicitly setting the pixel scale of each warped
exposure with Warp's ‘--cdelt’ option that is described below. For good
strategies of setting the pixel scale, see *note Moire pattern in
stacking and its correction::.

Another problem that may arise when aligning images to new pixel
grids is the aliasing or visible Moiré patterns on the output image.
This artifact should be removed if you are stacking several exposures,
especially with a pointing pattern. If not see *note Moire pattern in
stacking and its correction:: for ways to mitigate the visible patterns.
See the description of ‘--gridfile’ below for more.

*Known issue:* Warp's WCS-based aligning works best with WCSLIB version
7.12 (released in September 2022) and above. If you have an older
version of WCSLIB, you might get a ‘wcss2p’ error otherwise.

‘-c FLT,FLT’
‘--center=FLT,FLT’
WCS coordinates of the center of the central pixel of the output
image. Since a central pixel is only defined with an odd number of
pixels along both dimensions, the output will always have an odd
number of pixels. When ‘--center’ or ‘--gridfile’ aren't given,
the output will have the same central WCS coordinate as the input.

Usually, the WCS coordinates are Right Ascension and Declination
(when the first three characters of ‘CTYPE1’ and ‘CTYPE2’ are
respectively ‘RA-’ and ‘DEC’). For more on the ‘CTYPEi’ keyword
values, see ‘--ctype’ below.

‘-w INT[,INT]’
‘--width=INT[,INT]’
Width and height of the output image in units of WCS (usually
degrees). If you want the values to be read as pixels, also call
the ‘--widthinpix’ option with ‘--width’. If a single value is
given, Warp will use the same value for the second dimension
(creating a square output). When ‘--width’ or ‘--gridfile’ aren't
given, Warp will calculate the necessary size of the output pixel
grid to fully contain the input image.

Usually the WCS coordinates are in units of degrees (defined by the
‘CUNITi’ keywords of the FITS standard). But entering a certain
number of arcseconds or arcminutes for the width can be annoying
(you will usually need to go to the calculator!). To simplify such
situations, this option also accepts division. For example
‘--width=1/60,2/60’ will make an aligned warp that is 1 arcmin
along Right Ascension and 2 arcminutes along the Declination.

With the ‘--widthinpix’ option the values will be interpreted as
numbers of pixels. In this scenario, this option should be given
_odd_ integer(s) that are greater than 1. This ensures that the
output image can have a _central_ pixel. Recall that through the
‘--center’ option, you specify the WCS coordinate of the center of
the central pixel. The central coordinate of an image with an even
number of pixels will be on the edge of two pixels, so a "central"
pixel is not well defined. If any of the given values are even,
Warp will automatically add a single pixel (to make it an odd
integer) and print a warning message.

‘--widthinpix’
When called, the values given to the ‘--width’ option will be
interpreted as the number of pixels along each dimension(s). See
the description of ‘--width’ for more.

‘-x FLT[,FLT]’
‘--cdelt=FLT[,FLT]’
Coordinate deltas or increments (‘CDELTi’ in the FITS standard), or
the pixel scale in both dimensions. If a single value is given, it
will be used for both axes. In this way, the output's pixels will
be squares on the sky at the reference point (as is usually
expected!). When ‘--cdelt’ or ‘--gridfile’ aren't given, Warp will
read the input's pixel scale and choose the larger of ‘CDELT1’ or
‘CDELT2’ so the output pixels are square.

Usually (when dealing with RA and Dec, and the ‘CUNITi’s have a
value of ‘deg’), the units of the given values are degrees/pixel.
Warp allows you to easily convert from _arcsec_ to _degrees_ by
simply appending a ‘/3600’ to the value. For example, for an
output image of pixel scale ‘0.27’ arcsec/pixel, you can use
‘--cdelt=0.27/3600’.

‘--ctype=STR,STR’
The coordinate types of the output (‘CTYPE1’ and ‘CTYPE2’ keywords
in the FITS standard), separated by a comma. By default the value
to this option is '‘RA---TAN,DEC--TAN’'. However, if ‘--gridfile’
is given, this option is ignored.

If you don't call ‘--ctype’ or ‘--gridfile’, the output WCS
coordinates will be Right Ascension and Declination, while the
output's projection will be Gnomonic
(https://en.wikipedia.org/wiki/Gnomonic_projection), also known as
Tangential (TAN). This combination is the most common in
extra-galactic imaging surveys. For other coordinates and
projections in your output use other values, as described below.

According to the FITS standard version 4.0(2): ‘CTYPEi’ is the
"type for the Intermediate-coordinate Axis $i$. Any coordinate
type that is not covered by this Standard or an officially
recognized FITS convention shall be taken to be linear. All
non-linear coordinate system names must be expressed in '4–3' form:
the first four characters specify the coordinate type, the fifth
character is a hyphen (‘-’), and the remaining three characters
specify an algorithm code for computing the world coordinate value.
Coordinate types with names of fewer than four characters are
padded on the right with hyphens, and algorithm codes with fewer
than three characters are padded on the right with SPACE. Algorithm
codes should be three characters" (see list of algorithm codes
below).

You can use any of the projection algorithms (last three characters
of each coordinate's type) provided by your host WCSLIB (a
mandatory dependency of Gnuastro; see *note WCSLIB::). For a very
elaborate and complete description of projection algorithms in the
FITS WCS standard, see Calabretta and Greisen 2002
(https://doi.org/10.1051/0004-6361:20021327). Wikipedia also has a
nice article on Map projections
(https://en.wikipedia.org/wiki/Map_projection). As an example,
WCSLIB 7.12 (released in September 2022) has the following
projection algorithms:

‘AZP’
Zenithal/azimuthal perspective
‘SZP’
Slant zenithal perspective
‘TAN’
Gnomonic (tangential)
‘STG’
Stereographic
‘SIN’
Orthographic/synthesis
‘ARC’
Zenithal/azimuthal equidistant
‘ZPN’
Zenithal/azimuthal polynomial
‘ZEA’
Zenithal/azimuthal equal area
‘AIR’
Airy
‘CYP’
Cylindrical perspective
‘CEA’
Cylindrical equal area
‘CAR’
Plate carree
‘MER’
Mercator
‘SFL’
Sanson-Flamsteed
‘PAR’
Parabolic
‘MOL’
Mollweide
‘AIT’
Hammer-Aitoff
‘COP’
Conic perspective
‘COE’
Conic equal area
‘COD’
Conic equidistant
‘COO’
Conic orthomorphic
‘BON’
Bonne
‘PCO’
Polyconic
‘TSC’
Tangential spherical cube
‘CSC’
COBE spherical cube
‘QSC’
Quadrilateralized spherical cube
‘HPX’
HEALPix
‘XPH’
HEALPix polar, aka "butterfly"

‘-G’
‘--gridfile’
FITS filename containing the final pixel grid and WCS for the
output image. The HDU/extension containing should be specified
with ‘--gridhdu’ or its short option ‘-H’. The HDU should contain
a WCS, otherwise, Warp will abort with a crash. When this option
is used, Warp will read the respective WCS and the size of the
image to resample the input. Since this WCS of this HDU contains
everything needed to construct the WCS the options above will be
ignored when ‘--gridfile’ is called: ‘--cdelt’, ‘--center’, and
‘--widthinpix’.

In the example below, let's use this option to put the image of M51
in one survey (J-PLUS) into the pixel grid of another survey (SDSS)
containing M51. The J-PLUS field of view is very large (almost
$1.5\times1.5$ deg$^2$, in $9500\times9500$ pixels), while the
field of view of SDSS in each filter is small (almost
$0.3\times0.25$ deg$^2$ in $2048\times1489$ pixels). With the
first two commands, we'll first download the two images, then we'll
extract the portion of the J-PLUS image that overlaps with the SDSS
image and align it exactly to SDSS's pixel grid. Note that these
are the two images that were used in two of Gnuastro's tutorials:
*note Building the extended PSF:: and *note Detecting large
extended targets::.

## Download the J-PLUS DR2 image of M51 in the r filter.
$ jplusbase="http://archive.cefca.es/catalogues/vo/siap"
$ wget $jplusbase/jplus-dr2/get_fits?id=67510 \
-O jplus.fits.fz

## Download the SDSS image in r filter and decompress it
## (Bzip2 is not a standard FITS compression algorithm).
$ sdssbase=https://dr12.sdss.org/sas/dr12/boss/photoObj/frames
$ wget $sdssbase/301/3716/6/frame-r-003716-6-0117.fits.bz2 \
-O sdss.fits.bz2
$ bunzip2 sdss.fits.bz2

## Warp and crop the J-PLUS image so the output exactly
## matches the SDSS pixel gid.
$ astwarp jplus.fits.fz --gridfile=sdss.fits --gridhdu=0 \
--output=jplus-on-sdss.fits

## View the two images side-by-side:
$ astscript-fits-view sdss.fits jplus-on-sdss.fits

As the example above shows, this option can therefore be very
useful when comparing images from multiple surveys. But there are
other very interesting use cases also. For example, when you are
making a mock dataset and need to add distortion to the image so it
matches the distortion of your camera. Through ‘--gridhdu’, you
can easily insert that distortion over the mock image and put the
mock image in the pixel grid of an exposure.

‘-H’
‘--gridhdu’
The HDU/extension of the reference WCS file specified with option
‘--wcsfile’ or its short version ‘-H’ (see the description of
‘--wcsfile’ for more).

‘--edgesampling=INT’
Number of extra samplings along the edge of a pixel. By default
the value is ‘0’ (the output pixel's polygon over the input will be
a quadrilateral (a polygon with four edges/vertices).

Warp uses pixel mixing to derive the output pixel values. For a
complete introduction, see *note Resampling::, and in particular
its later part on distortions. To account for this possible
curvature due to distortion, you can use this option. For example,
‘--edgesampling=1’ will add one extra vertice in the middle of each
edge of the output pixel, producing an 8-vertice polygon.
Similarly, ‘--edgesampling=5’ will put 5 extra vertices along each
edge, thus sampling the shape (and possible curvature) of the
output pixel over an input pixel with $4+5\times4=24$ vertice
polygon. Since the polygon clipping will happen for every output
pixel, a higher value to this option can significantly reduce the
running speed and increase the RAM usage of Warp; so use it with
caution: in most cases the default ‘--edgesampling=0’ is
sufficient.

To visually inspect the curvature effect on pixel area of the input
image, see option ‘--pixelareaonwcs’ in *note Pixel information
images::.

‘--checkmaxfrac’
Check each output pixel's maximum coverage on the input data and
append as the '‘MAX-FRAC’' HDU/extension to the output aligned
image. This option provides an easy visual inspection for possible
recurring patterns or fringes caused by aligning to a new pixel
grid. For more detail about the origin of these patterns and how
to mitigate them see *note Moire pattern in stacking and its
correction::.

Note that the '‘MAX-FRAC’' HDU/extension is not showing the
patterns themselves; It represents the largest area coverage on the
input data for that particular pixel. The values can be in the
range between 0 to 1, where 1 means the pixel is covering at least
one complete pixel of the input data. On the other hand, 0 means
that the pixel is not covering any pixels of the input at all.

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

(1) For linear warps, see *note Linear warps to be called
explicitly::.

(2) FITS standard version 4.0:
<https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf>

File: gnuastro.info, Node: Linear warps to be called explicitly, Prev: Align pixels with WCS considering distortions, Up: Invoking astwarp

6.4.4.2 Linear warps to be called explicitly
............................................

Linear warps include operations like rotation, scaling, sheer, etc. For
an introduction, see *note Linear warping basics::. These are warps
that don't depend on the WCS of the image and should be explicitly
requested. To align the input pixel coordinates with the WCS
coordinates, see *note Align pixels with WCS considering distortions::.

While they will correct any existing WCS based on the warp, they can
also operate on images without any WCS. For example, you have a mock
image that doesn't (yet!) have its mock WCS, and it has been created on
an over-sampled grid and convolved with an over-sampled PSF. In this
scenario, you can use the ‘--scale’ option to under-sample it to your
desired resolution. This is similar to the *note Sufi simulates a
detection:: tutorial.

Linear warps must be specified as command-line options, either as
(possibly multiple) modular warpings (for example, ‘--rotate’, or
‘--scale’), or directly as a single raw matrix (with ‘--matrix’). If
specified together, the latter (direct matrix) will take precedence and
all the modular warpings will be ignored. Any number of modular
warpings can be specified on the command-line and configuration files.
If more than one modular warping is given, all will be merged to create
one warping matrix. As described in *note Merging multiple warpings::,
matrix multiplication is not commutative, so the order of specifying the
modular warpings on the command-line, and/or configuration files makes a
difference (see *note Configuration file precedence::). The full list
of modular warpings and the other options particular to Warp are
described below.

The values to the warping options (modular warpings as well as
‘--matrix’), are a sequence of at least one number. Each number in this
sequence is separated from the next by a comma (<,>). Each number can
also be written as a single fraction (with a forward-slash </> between
the numerator and denominator). Space and Tab characters are permitted
between any two numbers, just do not forget to quote the whole value.
Otherwise, the value will not be fully passed onto the option. See the
examples above as a demonstration.

Based on the FITS standard, integer values are assigned to the center
of a pixel and the coordinate [1.0, 1.0] is the center of the first
pixel (bottom left of the image when viewed in SAO DS9). So the
coordinate center [0.0, 0.0] is half a pixel away (in each axis) from
the bottom left vertex of the first pixel. The resampling that is done
in Warp (see *note Resampling::) is done on the coordinate axes and thus
directly depends on the coordinate center. In some situations this if
fine, for example, when rotating/aligning a real image, all the edge
pixels will be similarly affected. But in other situations (for
example, when scaling an over-sampled mock image to its intended
resolution, this is not desired: you want the center of the coordinates
to be on the corner of the pixel. In such cases, you can use the
‘--centeroncorner’ option which will shift the center by $0.5$ before
the main warp, then shift it back by $-0.5$ after the main warp.

‘-r FLT’
‘--rotate=FLT’
Rotate the input image by the given angle in degrees: $\theta$ in
*note Linear warping basics::. Note that commonly, the WCS
structure of the image is set such that the RA is the inverse of
the image horizontal axis which increases towards the right in the
FITS standard and as viewed by SAO DS9. So the default center for
rotation is on the right of the image. If you want to rotate about
other points, you have to translate the warping center first (with
‘--translate’) then apply your rotation and then return the center
back to the original position (with another call to ‘--translate’,
see *note Merging multiple warpings::.

‘-s FLT[,FLT]’
‘--scale=FLT[,FLT]’
Scale the input image by the given factor(s): $M$ and $N$ in *note
Linear warping basics::. If only one value is given, then both
image axes will be scaled with the given value. When two values
are given (separated by a comma), the first will be used to scale
the first axis and the second will be used for the second axis. If
you only need to scale one axis, use ‘1’ for the axis you do not
need to scale. The value(s) can also be written (on the
command-line or in configuration files) as a fraction.

‘-f FLT[,FLT]’
‘--flip=FLT[,FLT]’
Flip the input image around the given axis(s). If only one value
is given, then both image axes are flipped. When two values are
given (separated by acomma), you can choose which axis to flip
over. ‘--flip’ only takes values ‘0’ (for no flip), or ‘1’ (for a
flip). Hence, if you want to flip by the second axis only, use
‘--flip=0,1’.

‘-e FLT[,FLT]’
‘--shear=FLT[,FLT]’
Shear the input image by the given value(s): $A$ and $B$ in *note
Linear warping basics::. If only one value is given, then both
image axes will be sheared with the given value. When two values
are given (separated by a comma), the first will be used to shear
the first axis and the second will be used for the second axis. If
you only need to shear along one axis, use ‘0’ for the axis that
must be untouched. The value(s) can also be written (on the
command-line or in configuration files) as a fraction.

‘-t FLT[,FLT]’
‘--translate=FLT[,FLT]’
Translate (move the center of coordinates) the input image by the
given value(s): $c$ and $f$ in *note Linear warping basics::. If
only one value is given, then both image axes will be translated by
the given value. When two values are given (separated by a comma),
the first will be used to translate the first axis and the second
will be used for the second axis. If you only need to translate
along one axis, use ‘0’ for the axis that must be untouched. The
value(s) can also be written (on the command-line or in
configuration files) as a fraction.

‘-p FLT[,FLT]’
‘--project=FLT[,FLT]’
Apply a projection to the input image by the given values(s): $g$
and $h$ in *note Linear warping basics::. If only one value is
given, then projection will apply to both axes with the given
value. When two values are given (separated by a comma), the first
will be used to project the first axis and the second will be used
for the second axis. If you only need to project along one axis,
use ‘0’ for the axis that must be untouched. The value(s) can also
be written (on the command-line or in configuration files) as a
fraction.

‘-m STR’
‘--matrix=STR’
The warp/transformation matrix. All the elements in this matrix
must be separated by commas(<,>) characters and as described above,
you can also use fractions (a forward-slash between two numbers).
The transformation matrix can be either a 2 by 2 (4 numbers), or a
3 by 3 (9 numbers) array. In the former case (if a 2 by 2 matrix
is given), then it is put into a 3 by 3 matrix (see *note Linear
warping basics::).

The determinant of the matrix has to be non-zero and it must not
contain any non-number values (for example, infinities or NaNs).
The elements of the matrix have to be written row by row. So for
the general Homography matrix of *note Linear warping basics::, it
should be called with ‘--matrix=a,b,c,d,e,f,g,h,1’.

The raw matrix takes precedence over all the modular warping
options listed above, so if it is called with any number of modular
warps, the latter are ignored.

‘--centeroncorner’
Put the center of coordinates on the corner of the first
(bottom-left when viewed in SAO DS9) pixel. This option is applied
after the final warping matrix has been finalized: either through
modular warpings or the raw matrix. See the explanation above for
coordinates in the FITS standard to better understand this option
and when it should be used.

‘-k’
‘--keepwcs’
Do not correct the WCS information of the input image and save it
untouched to the output image. By default the WCS (World
Coordinate System) information of the input image is going to be
corrected in the output image so the objects in the image are at
the same WCS coordinates. But in some cases it might be useful to
keep it unchanged (for example, to correct alignments).

File: gnuastro.info, Node: Data analysis, Next: Data modeling, Prev: Data manipulation, Up: Top

7 Data analysis
***************

Astronomical datasets (images or tables) contain very valuable
information, the tools in this section can help in analyzing,
extracting, and quantifying that information. For example, getting
general or specific statistics of the dataset (with *note Statistics::),
detecting signal within a noisy dataset (with *note NoiseChisel::), or
creating a catalog from an input dataset (with *note MakeCatalog::).

* Menu:

* Statistics:: Calculate dataset statistics.
* NoiseChisel:: Detect objects in an image.
* Segment:: Segment detections based on signal structure.
* MakeCatalog:: Catalog from input and labeled images.
* Match:: Match two datasets.

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

7.1 Statistics
==============

The distribution of values in a dataset can provide valuable information
about it. For example, in an image, if it is a positively skewed
distribution, we can see that there is significant data in the image.
If the distribution is roughly symmetric, we can tell that there is no
significant data in the image. In a table, when we need to select a
sample of objects, it is important to first get a general view of the
whole sample.

On the other hand, you might need to know certain statistical
parameters of the dataset. For example, if we have run a detection
algorithm on an image, and we want to see how accurate it was, one
method is to calculate the average of the undetected pixels and see how
reasonable it is (if detection is done correctly, the average of
undetected pixels should be approximately equal to the background value,
see *note Sky value::). In a table, you might have calculated the
magnitudes of a certain class of objects and want to get some general
characteristics of the distribution immediately on the command-line
(very fast!), to possibly change some parameters. The Statistics
program is designed for such situations.

* Menu:

* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
* Least squares fitting:: Fitting with various parametric functions.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.

File: gnuastro.info, Node: Histogram and Cumulative Frequency Plot, Next: 2D Histograms, Prev: Statistics, Up: Statistics

7.1.1 Histogram and Cumulative Frequency Plot
---------------------------------------------

Histograms and the cumulative frequency plots are both used to visually
study the distribution of a dataset. A histogram shows the number of
data points which lie within pre-defined intervals (bins). So on the
horizontal axis we have the bin centers and on the vertical, the number
of points that are in that bin. You can use it to get a general view of
the distribution: which values have been repeated the most? how
close/far are the most significant bins? Are there more values in the
larger part of the range of the dataset, or in the lower part?
Similarly, many very important properties about the dataset can be
deduced from a visual inspection of the histogram. In the Statistics
program, the histogram can be either output to a table to plot with your
favorite plotting program(1), or it can be shown with ASCII characters
on the command-line, which is very crude, but good enough for a fast and
on-the-go analysis, see the example in *note Invoking aststatistics::.

The width of the bins is only necessary parameter for a histogram.
In the limiting case that the bin-widths tend to zero (while assuming
the number of points in the dataset tend to infinity), then the
histogram will tend to the probability density function
(https://en.wikipedia.org/wiki/Probability_density_function) of the
distribution. When the absolute number of points in each bin is not
relevant to the study (only the shape of the histogram is important),
you can _normalize_ a histogram so like the probability density
function, the sum of all its bins will be one.

In the cumulative frequency plot of a distribution, the horizontal
axis is the sorted data values and the y axis is the index of each data
in the sorted distribution. Unlike a histogram, a cumulative frequency
plot does not involve intervals or bins. This makes it less prone to
any sort of bias or error that a given bin-width would have on the
analysis. When a larger number of the data points have roughly the same
value, then the cumulative frequency plot will become steep in that
vicinity. This occurs because on the horizontal axis, there is little
change while on the vertical axis, the indexes constantly increase.
Normalizing a cumulative frequency plot means to divide each index (y
axis) by the total number of data points (or the last value).

Unlike the histogram which has a limited number of bins, ideally the
cumulative frequency plot should have one point for every data element.
Even in small datasets (for example, a $200\times200$ image) this will
result in an unreasonably large number of points to plot (40000)! As a
result, for practical reasons, it is common to only store its value on a
certain number of points (intervals) in the input range rather than the
whole dataset, so you should determine the number of bins you want when
asking for a cumulative frequency plot. In Gnuastro (and thus the
Statistics program), the number reported for each bin is the total
number of data points until the larger interval value for that bin. You
can see an example histogram and cumulative frequency plot of a single
dataset under the ‘--asciihist’ and ‘--asciicfp’ options of *note
Invoking aststatistics::.

So as a summary, both the histogram and cumulative frequency plot in
Statistics will work with bins. Within each bin/interval, the lower
value is considered to be within then bin (it is inclusive), but its
larger value is not (it is exclusive). Formally, an interval/bin
between a and b is represented by [a, b). When the over-all range of
the dataset is specified (with the ‘--greaterequal’, ‘--lessthan’, or
‘--qrange’ options), the acceptable values of the dataset are also
defined with a similar inclusive-exclusive manner. But when the range
is determined from the actual dataset (none of these options is called),
the last element in the dataset is included in the last bin's count.

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

(1) We recommend PGFPlots (http://pgfplots.sourceforge.net/) which
generates your plots directly within TeX (the same tool that generates
your document).

File: gnuastro.info, Node: 2D Histograms, Next: Least squares fitting, Prev: Histogram and Cumulative Frequency Plot, Up: Statistics

7.1.2 2D Histograms
-------------------

In *note Histogram and Cumulative Frequency Plot:: the concept of
histograms were introduced on a single dataset. But they are only
useful for viewing the distribution of a single variable (column in a
table). In many contexts, the distribution of two variables in relation
to each other may be of interest. For example, the color-magnitude
diagrams in astronomy, where the horizontal axis is the luminosity or
magnitude of an object, and the vertical axis is the color. Scatter
plots are useful to see these relations between the objects of interest
when the number of the objects is small.

As the density of points in the scatter plot increases, the points
will fall over each other and just make a large connected region hide
potentially interesting behaviors/correlations in the densest regions.
This is where 2D histograms can become very useful. A 2D histogram is
composed of 2D bins (boxes or pixels), just as a 1D histogram consists
of 1D bins (lines). The number of points falling within each box/pixel
will then be the value of that box. Added with a color-bar, you can now
clearly see the distribution independent of the density of points (for
example, you can even visualize it in log-scale if you want).

Gnuastro's Statistics program has the ‘--histogram2d’ option for this
task. It takes a single argument (either ‘table’ or ‘image’) that
specifies the format of the output 2D histogram. The two formats will
be reviewed separately in the sub-sections below. But let's start with
the generalities that are common to both (related to the input, not the
output).

You can specify the two columns to be shown using the ‘--column’ (or
‘-c’) option. So if you want to plot the color-magnitude diagram from a
table with the ‘MAG-R’ column on the horizontal and ‘COLOR-G-R’ on the
vertical column, you can use ‘--column=MAG-r,COLOR-G-r’. The number of
bins along each dimension can be set with ‘--numbins’ (for first input
column) and ‘--numbins2’ (for second input column).

Without specifying any range, the full range of values will be used
in each dimension. If you only want to focus on a certain interval of
the values in the columns in any dimension you can use the
‘--greaterequal’ and ‘--lessthan’ options to limit the values along the
first/horizontal dimension and ‘--greaterequal2’ and ‘--lessthan2’
options for the second/vertical dimension.

* Menu:

* 2D histogram as a table for plotting:: Format and usage in table format.
* 2D histogram as an image:: Format and usage in image format

File: gnuastro.info, Node: 2D histogram as a table for plotting, Next: 2D histogram as an image, Prev: 2D Histograms, Up: 2D Histograms

7.1.2.1 2D histogram as a table for plotting
............................................

When called with the ‘--histogram=table’ option, Statistics will output
a table file with three columns that have the information of every box
as a column. If you asked for ‘--numbins=N’ and ‘--numbins2=M’, all
three columns will have $M\times N$ rows (one row for every box/pixel of
the 2D histogram). The first and second columns are the position of the
box along the first and second dimensions. The third column has the
number of input points that fall within that box/pixel.

For example, you can make high-quality plots within your paper (using
the same LaTeX engine, thus blending very nicely with your text) using
PGFPlots (https://ctan.org/pkg/pgfplots). Below you can see one such
minimal example, using your favorite text editor, save it into a file,
make the two small corrections in it, then run the commands shown at the
top. This assumes that you have LaTeX installed, if not the steps to
install a minimally sufficient LaTeX package on your system, see the
respective section in *note Bootstrapping dependencies::.

The two parts that need to be corrected are marked with '‘%% <--’':
the first one (‘XXXXXXXXX’) should be replaced by the value to the
‘--numbins’ option which is the number of bins along the first
dimension. The second one (‘FILE.txt’) should be replaced with the name
of the file generated by Statistics.

%% Replace 'XXXXXXXXX' with your selected number of bins in the first
%% dimension.
%%
%% Then run these commands to build the plot in a LaTeX command.
%% mkdir tikz
%% pdflatex --shell-escape --halt-on-error report.tex
\documentclass{article}

%% Load PGFPlots and set it to build the figure separately in a 'tikz'
%% directory (which has to exist before LaTeX is run). This
%% "externalization" is very useful to include the commands of multiple
%% plots in the middle of your paper/report, but also have the plots
%% separately to use in slides or other scenarios.
\usepackage{pgfplots}
\usetikzlibrary{external}
\tikzexternalize
\tikzsetexternalprefix{tikz/}

%% Define colormap for the PGFPlots 2D histogram
\pgfplotsset{
/pgfplots/colormap={hsvwhitestart}{
rgb255(0cm)=(255,255,255)
rgb255(0.10cm)=(128,0,128)
rgb255(0.5cm)=(0,0,230)
rgb255(1.cm)=(0,255,255)
rgb255(2.5cm)=(0,255,0)
rgb255(3.5cm)=(255,255,0)
rgb255(6cm)=(255,0,0)
}
}

%% Start the prinable document
\begin{document}

You can write a full paper here and include many figures!
Describe what the two axes are, and how you measured them.
Also, do not forget to explain what it shows and how to interpret it.
You also have separate PDFs for every figure in the `tikz' directory.
Feel free to change this text.

%% Draw the plot.
\begin{tikzpicture}
\small
\begin{axis}[
width=\linewidth,
view={0}{90},
colorbar horizontal,
xlabel=X axis,
ylabel=Y axis,
ylabel shift=-0.1cm,
colorbar style={at={(0,1.01)}, anchor=south west,
xticklabel pos=upper},
]
\addplot3[
surf,
shader=flat corner,
mesh/ordering=rowwise,
mesh/cols=XXXXXXXXX, %% <-- Number of bins in 1st column.
] file {FILE.txt}; %% <-- Name of aststatistics output.

\end{axis}
\end{tikzpicture}

%% End the printable document.
\end{document}

Let's assume you have put the LaTeX source above, into a plain-text
file called ‘report.tex’. The PGFPlots call above is configured to
build the plots as separate PDF files in a ‘tikz/’ directory(1). This
allows you to directly load those PDFs in your slides or other reports.
Therefore, before building the PDF report, you should first make a
‘tikz/’ directory:

$ mkdir tikz

To build the final PDF, you should run ‘pdflatex’ with the
‘--shell-escape’ option, so it can build the separate PDF(s) separately.
We are also adding the ‘--halt-on-error’ so it immediately aborts in the
case of an error (in the case of an error, by default LaTeX will not
abort, it will stop and ask for your input to temporarily change things
and try fixing the error, but it has a special interface which can be
hard to master).

$ pdflatex --shell-escape --halt-on-error report.tex

You can now open ‘report.pdf’ to see your very high quality 2D histogram
within your text. And if you need the plots separately (for example,
for slides), you can take the PDF inside the ‘tikz/’ directory.

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

(1) TiKZ (https://www.ctan.org/pkg/pgf) is the name of the
lower-level engine behind PGPlots.

File: gnuastro.info, Node: 2D histogram as an image, Prev: 2D histogram as a table for plotting, Up: 2D Histograms

7.1.2.2 2D histogram as an image
................................

When called with the ‘--histogram=image’ option, Statistics will output
a FITS file with an image/array extension. If you asked for
‘--numbins=N’ and ‘--numbins2=M’ the image will have a size of $N\times
M$ pixels (one pixel per 2D bin). Also, the FITS image will have a
linear WCS that is scaled to the 2D bin size along each dimension. So
when you hover your mouse over any part of the image with a FITS viewer
(for example, SAO DS9), besides the number of points in each pixel, you
can directly also see "coordinates" of the pixels along the two axes.
You can also use the optimized and fast FITS viewer features for many
aspects of visually inspecting the distributions (which we will not go
into further).

For example, let's assume you want to derive the color-magnitude
diagram (CMD) of the UVUDF survey (http://uvudf.ipac.caltech.edu). You
can run the first command below to download the table with magnitudes of
objects in many filters and run the second command to see general column
metadata after it is downloaded.

$ wget http://asd.gsfc.nasa.gov/UVUDF/uvudf_rafelski_2015.fits.gz
$ asttable uvudf_rafelski_2015.fits.gz -i

Let's assume you want to find the color to be between the ‘F606W’ and
‘F775W’ filters (roughly corresponding to the g and r filters in
ground-based imaging). However, the original table does not have color
columns (there would be too many combinations!). Therefore you can use
the *note Column arithmetic:: feature of Gnuastro's Table program for
deriving a new table with the ‘F775W’ magnitude in one column and the
difference between the ‘F606W’ and ‘F775W’ on the other column. With
the second command, you can see the actual values if you like.

$ asttable uvudf_rafelski_2015.fits.gz -cMAG_F775W \
-c'arith MAG_F606W MAG_F775W -' \
--colmetadata=ARITH_1,F606W-F775W,"AB mag" -ocmd.fits
$ asttable cmd.fits

You can now construct your 2D histogram as a $100\times100$ pixel FITS
image with this command (assuming you want ‘F775W’ magnitudes between 22
and 30, colors between -1 and 3 and 100 bins in each dimension). Note
that without the ‘--manualbinrange’ option the range of each axis will
be determined by the values within the columns (which may be larger or
smaller than your desired large).

aststatistics cmd.fits -cMAG_F775W,F606W-F775W --histogram2d=image \
--numbins=100 --greaterequal=22 --lessthan=30 \
--numbins2=100 --greaterequal2=-1 --lessthan2=3 \
--manualbinrange --output=cmd-2d-hist.fits

If you have SAO DS9, you can now open this FITS file as a normal FITS
image, for example, with the command below. Try hovering/zooming over
the pixels: not only will you see the number of objects in the UVUDF
catalog that fall in each bin, but you also see the ‘F775W’ magnitude
and color of that pixel also.

$ ds9 cmd-2d-hist.fits -cmap sls -zoom to fit

With the first command below, you can activate the grid feature of DS9
to actually see the coordinate grid, as well as values on each line.
With the second command, DS9 will even read the labels of the axes and
use them to generate an almost publication-ready plot.

$ ds9 cmd-2d-hist.fits -cmap sls -zoom to fit -grid yes
$ ds9 cmd-2d-hist.fits -cmap sls -zoom to fit -grid yes \
-grid type publication

If you are happy with the grid, coloring and the rest, you can also
use ds9 to save this as a JPEG image to directly use in your
documents/slides with these extra DS9 options (DS9 will write the image
to ‘cmd-2d.jpeg’ and quit immediately afterwards):

$ ds9 cmd-2d-hist.fits -cmap sls -zoom 4 -grid yes \
-grid type publication -saveimage cmd-2d.jpeg -quit

This is good for a fast progress update. But for your paper or more
official report, you want to show something with higher quality. For
that, you can use the PGFPlots package in LaTeX to add axes in the same
font as your text, sharp grids and many other elegant/powerful features
(like over-plotting interesting points and lines). But to load the 2D
histogram into PGFPlots first you need to convert the FITS image into a
more standard format, for example, PDF. We will use Gnuastro's *note
ConvertType:: for this, and use the ‘sls-inverse’ color map (which will
map the pixels with a value of zero to white):

$ astconvertt cmd-2d-hist.fits --colormap=sls-inverse \
--borderwidth=0 -ocmd-2d-hist.pdf

Below you can see a minimally working example of how to add axis
numbers, labels and a grid to the PDF generated above. Copy and paste
the LaTeX code below into a plain-text file called ‘cmd-report.tex’
Notice the ‘xmin’, ‘xmax’, ‘ymin’, ‘ymax’ values and how they are the
same as the range specified above.

\documentclass{article}
\usepackage{pgfplots}
\dimendef\prevdepth=0
\begin{document}

You can write all you want here...

\begin{tikzpicture}
\begin{axis}[
enlargelimits=false,
grid,
axis on top,
width=\linewidth,
height=\linewidth,
xlabel={Magnitude (F775W)},
ylabel={Color (F606W-F775W)}]

\addplot graphics[xmin=22, xmax=30, ymin=-1, ymax=3]
{cmd-2d-hist.pdf};
\end{axis}
\end{tikzpicture}
\end{document}

Run this command to build your PDF (assuming you have LaTeX and
PGFPlots).

$ pdflatex cmd-report.tex

The improved quality, blending in with the text, vector-graphics
resolution and other features make this plot pleasing to the eye, and
let your readers focus on the main point of your scientific argument.
PGFPlots can also built the PDF of the plot separately from the rest of
the paper/report, see *note 2D histogram as a table for plotting:: for
the necessary changes in the preamble.

File: gnuastro.info, Node: Least squares fitting, Next: Sky value, Prev: 2D Histograms, Up: Statistics

7.1.3 Least squares fitting
---------------------------

After completing a good observation, doing robust data reduction and
finalizing the measurements, it is commonly necessary to parameterize
the derived correlations. For example, you have derived the radial
profile of the PSF of your image (see *note Building the extended
PSF::). You now want to parameterize the radial profile to estimate the
slope. Alternatively, you may have found the star formation rate and
stellar mass of your sample of galaxies. Now, you want to derive the
star formation main sequence as a parametric relation between the two.
The fitting functions below can be used for such purposes.

Gnuastro's least squares fitting features are just wrappers over the
least squares fitting methods of the linear
(https://www.gnu.org/software/gsl/doc/html/lls.html) and nonlinear
(https://www.gnu.org/software/gsl/doc/html/nls.html) least-squares
fitting functions of the GNU Scientific Library (GSL). For the low-level
details and equations of the methods, please see the GSL documentation.
The names have been preserved here in Gnuastro to simplify the
connection with GSL and follow the details in the detailed documentation
there.

GSL is a very low-level library, designed for maximal portability to
many scenarios, and power. Therefore calling GSL's functions directly
for a fast operation requires a good knowledge of the C programming
language and many lines of code. As a low-level library, GSL is
designed to be the back-end of higher-level programs (like Gnuastro).
Through the Statistics program, in Gnuastro we provide a high-level
interface to access to GSL's very powerful least squares fitting engine
to read/write from/to standard data formats in astronomy. A fully
working example is shown below.

To activate fitting in Statistics, simply give your desired fitting
method to the ‘--fit’ option (for the full list of acceptable methods,
see *note Fitting options::). For example, with the command below,
we'll build a fake measurement table (including noise) from the
polynomial $y=1.23-4.56x+7.89x^2$. To understand how this equation
translates to the command below (part before ‘set-y’), see *note Reverse
polish notation:: and *note Column arithmetic::. We will set the X axis
to have values from 0.1 to 2, with steps of 0.01 and let's assume a
random Gaussian noise to each $y$ measurement: $\sigma_y=0.1y$. To make
the random number generation exactly reproducible, we are also setting
the seed (see *note Generating random numbers::, which also uses GSL as
a backend). To learn more about the ‘mknoise-sigma’ operator, see the
Arithmetic program's *note Random number generators::.

$ export GSL_RNG_SEED=1664015492
$ seq 0.1 0.01 2 \
| asttable --output=noisy.fits --envseed -c1 \
-c'arith 1.23 -4.56 $1 x + 7.89 $1 x $1 x + set-y \
0.1 y x set-yerr \
y yerr mknoise-sigma yerr' \
--colmetadata=1,X --colmetadata=2,Y \
--colmetadata=3,Yerr

Let's have a look at the output plot with TOPCAT using the command
below.

$ astscript-fits-view noisy.fits

To see the error-bars, after opening the scatter plot, go into the
"Form" tab for that plot. Click on the button with a green "+" sign
followed by "Forms" and select "XYError". On the side-menu, in front of
"Y Positive Error", select the ‘Yerr’ column of the input table.

As you see, the error bars do indeed increase for higher X axis
values. Since we have error bars in this example (as in any
measurement), we can use weighted fitting. Also, this isn't a linear
relation, so we'll use a polynomial to second order (a maximum power of
2 in the form of $Y=c_0+c_1X+c_2X^2$):

$ aststatistics noisy.fits -cX,Y,Yerr --fit=polynomial-weighted \
--fitmaxpower=2
Statistics (GNU Astronomy Utilities) 0.22
-------
Fitting results (remove extra info with '--quiet' or '-q)
Input file: noisy.fits (hdu: 1) with 191 non-blank rows.
X column: X
Y column: Y
Weight column: Yerr [Standard deviation of Y in each row]

Fit function: Y = c0 + (c1 * X^1) + (c2 * X^2) + ... (cN * X^N)
N: 2
c0: +1.2286211608
c1: -4.5127796636
c2: +7.8435883943

Covariance matrix:
+0.0010496001 -0.0039928488 +0.0028367390
-0.0039928488 +0.0175244127 -0.0138030778
+0.0028367390 -0.0138030778 +0.0128129806

Reduced chi^2 of fit:
+0.9740670090

As you see from the elaborate message, the weighted polynomial
fitting has found return the $c_0$, $c_1$ and $c_2$ of
$Y=c_0+c_1X+c_2X^2$ that best represents the data we inserted. Our
input values were $c_0=1.23$, $c_1=-4.56$ and $c_2=7.89$, and the fitted
values are $c_0\approx1.2286$, $c_1\approx-4.5128$ and
$c_2\approx7.8436$ (which is statistically a very good fit! given that
we knew the original values a-priori!). The covariance matrix is also
calculated, it is necessary to calculate error bars on the estimations
and contains a lot of information (e.g., possible correlations between
parameters). Finally, the reduced $\chi^2$ (or $\chi_{red}^2$) of the
fit is also printed (which was the measure to minimize). A
$\chi_{red}^2\approx1$ shows a good fit. This is good for real-world
scenarios when you don't know the original values a-priori. For more on
interpreting $\chi_{red}^2\approx1$, see Andrae et al. 2010
(https://arxiv.org/abs/1012.3754).

The comparison of fitted and input values look pretty good, but
nothing beats visual inspection! To see how this looks compared to the
data, let's open the table again:

$ astscript-fits-view noisy.fits

Repeat the steps below to show the scatter plot and error-bars.
Then, go to the "Layers" menu and select "Add Function Control". Use
the results above to fill the box in front of "Function Expression":
‘1.2286+(-4.5128*x)+(7.8436*x*x)’. You will see that the second order
polynomial falls very nicely over the points(1). But this fit is not
perfect: it also has errors (inherited from the measurement errors). We
need the covariance matrix to estimate the errors on each point, and
that can be complex to do by hand.

Fortunately GSL has the tools to easily estimate the function at any
point and also calculate its corresponding error. To access this
feature within Gnuastro's Statistics program, you should use the
‘--fitestimate’ option. You can either give an independent table file
name (with ‘--fitestimatehdu’ and ‘--fitestimatecol’ to specify the HDU
and column in that file), or just ‘self’ so it uses the same X axis
column that was used in this fit. Let's use the easier case:

$ aststatistics noisy.fits -cX,Y,Yerr --fit=polynomial-weighted \
--fitmaxpower=2 --fitestimate=self --output=est.fits

...[[truncated; same as above]]...

Requested estimation:
Written to: est.fits

The first lines of the printed text are the same as before.
Afterwards, you will see a new line printed in the output, saying that
the estimation was written in ‘est.fits’. You can now inspect the two
tables with TOPCAT again with the command below. After TOPCAT opens,
plot both scatter plots:

$ astscript-fits-view noisy.fits est.fits

It is clear that they fall nicely on top of each other. The
‘est.fits’ table also has a third column with error bars. You can
follow the same steps before and draw the error bars to see how they
compare with the scatter of the measured data. They are much smaller
than the error in each point because we had a very good sampling of the
function in our noisy data.

Another useful point with the estimated output file is that it
contains all the fitting outputs as keywords in the header:

$ astfits est.fits -h1
...[[truncated]]...

/ Fit results
FITTYPE = 'polynomial-weighted' / Functional form of the fitting.
FITMAXP = 2 / Maximum power of polynomial.
FITIN = 'noisy.fits' / Name of file with input columns.
FITINHDU= '1 ' / Name or Number of HDU with input cols.
FITXCOL = 'X ' / Name or Number of independent (X) col.
FITYCOL = 'Y ' / Name or Number of measured (Y) column.
FITWCOL = 'Yerr ' / Name or Number of weight column.
FITWNAT = 'Standard deviation' / Nature of weight column.
FRDCHISQ= 0.974067008958516 / Reduced chi^2 of fit.
FITC0 = 1.22862116084727 / C0: multiple of x^0 in polynomial
FITC1 = -4.51277966356177 / C1: multiple of x^1 in polynomial
FITC2 = 7.84358839431161 / C2: multiple of x^2 in polynomial
FCOV11 = 0.00104960011629718 / Covariance matrix element (1,1).
FCOV12 = -0.00399284880859776 / Covariance matrix element (1,2).
FCOV13 = 0.00283673901863388 / Covariance matrix element (1,3).
FCOV21 = -0.00399284880859776 / Covariance matrix element (2,1).
FCOV22 = 0.0175244126670659 / Covariance matrix element (2,2).
FCOV23 = -0.0138030778380786 / Covariance matrix element (2,3).
FCOV31 = 0.00283673901863388 / Covariance matrix element (3,1).
FCOV32 = -0.0138030778380786 / Covariance matrix element (3,2).
FCOV33 = 0.0128129806394559 / Covariance matrix element (3,3).

...[[truncated]]...

In scenarios were you don't want the estimation, but only the fitted
parameters, all that verbose, human-friendly text or FITS keywords can
be an annoying extra step. For such cases, you should use the ‘--quiet’
option like below. It will print the parameters, rows of the covariance
matrix and $\chi_{red}^2$ on separate lines with nothing extra. This
allows you to parse the values in any way that you would like.

$ aststatistics noisy.fits -cX,Y,Yerr --fit=polynomial-weighted \
--fitmaxpower=2 --quiet
+1.2286211608 -4.5127796636 +7.8435883943
+0.0010496001 -0.0039928488 +0.0028367390
-0.0039928488 +0.0175244127 -0.0138030778
+0.0028367390 -0.0138030778 +0.0128129806
+0.9740670090

As a final example, because real data usually have outliers, let's
look at the "robust" polynomial fit which has special features to remove
outliers. First, we need to add some outliers to the table. To do
this, we'll make a plain-text table with ‘echo’, and use Table's
‘--catrowfile’ to concatenate (or append) those two rows to the original
table. Finally, we'll run the same fitting step above:

$ echo "0.6 20 0.01" > outliers.txt
$ echo "0.8 20 0.01" >> outliers.txt

$ asttable noisy.fits --catrowfile=outliers.txt \
--output=with-outlier.fits

$ aststatistics with-outlier.fits -cX,Y,Yerr --fit=polynomial-weighted \
--fitmaxpower=2 --fitestimate=self \
--output=est-out.fits

Statistics (GNU Astronomy Utilities) 0.22
-------
Fitting results (remove extra info with '--quiet' or '-q)
Input file: with-outlier.fits (hdu: 1) with 193 non-blank rows.
X column: X
Y column: Y
Weight column: Yerr [Standard deviation of Y in each row]

Fit function: Y = c0 + (c1 * X^1) + (c2 * X^2) + ... (cN * X^N)
N: 2
c0: -13.6446036899
c1: +66.8463258547
c2: -30.8746303591

Covariance matrix:
+0.0007889160 -0.0027706310 +0.0022208939
-0.0027706310 +0.0113922468 -0.0100306732
+0.0022208939 -0.0100306732 +0.0094087226

Reduced chi^2 of fit:
+4501.8356719150

Requested estimation:
Written to: est-out.fit

We see that the coefficient values have changed significantly and
that $\chi_{red}^2$ has increased to $4501$! Recall that a good fit
should have $\chi_{red}^2\approx1$. These numbers clearly show that the
fit was bad, but again, nothing beats a visual inspection. To visually
see the effect of those outliers, let's plot them with the command
below. You see that those two points have clearly caused a turn in the
fitted result which is terrible.

$ astscript-fits-view with-outlier.fits est-out.fits

For such cases, GSL has Robust linear regression
(https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression).
In Gnuastro's Statistics, you can access it with
‘--fit=polynomial-robust’, like the example below. Just note that the
robust method doesn't take an error column (because it estimates the
errors internally while rejecting outliers, based on the method).

$ aststatistics with-outlier.fits -cX,Y --fit=polynomial-robust \
--fitmaxpower=2 --fitestimate=self \
--output=est-out.fits --quiet

$ astfits est-out.fits -h1 | grep ^FITC
FITC0 = 1.20422691185238 / C0: multiple of x^0 in polynomial
FITC1 = -4.4779253576348 / C1: multiple of x^1 in polynomial
FITC2 = 7.84986153686548 / C2: multiple of x^2 in polynomial

$ astscript-fits-view with-outlier.fits est-out.fits

It is clear that the coefficients are very similar to the no-outlier
scenario above and if you run the second command to view the scatter
plots on TOPCAT, you also see that the fit nicely follows the curve and
is not affected by those two points. GSL provides many methods to
reject outliers. For their full list, see the description of
‘--fitrobust’ in *note Fitting options::. For a description of the
outlier rejection methods, see the GSL manual
(https://www.gnu.org/software/gsl/doc/html/lls.html#c.gsl_multifit_robust_workspace).

You may have noticed that unlike the cases before the last Statistics
command above didn't print anything on the standard output. This is
becasue ‘--quiet’ and ‘--fitestimate’ were called together. In this
case, because all the fitting parameters are written as FITS keywords,
because of the ‘--quiet’ option, they are no longer printed on standard
output.

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

(1) After plotting, you will notice that the legend made the plot too
thin. Fortunately you have a lot of empty space within the plot. To
bring the legend in, click on the "Legend" item on the bottom-left menu,
in the "Location" tab, click on "Internal" and hold and move it to the
top-left in the box below. To make the functional fit more clear, you
can click on the "Function" item of the bottom-left menu. In the
"Style" tab, change the color and thickness.

File: gnuastro.info, Node: Sky value, Next: Invoking aststatistics, Prev: Least squares fitting, Up: Statistics

7.1.4 Sky value
---------------

One of the most important aspects of a dataset is its reference value:
the value of the dataset where there is no signal. Without knowing, and
thus removing the effect of, this value it is impossible to compare the
derived results of many high-level analyses over the dataset with other
datasets (in the attempt to associate our results with the "real"
world).

In astronomy, this reference value is known as the "Sky" value: the
value that noise fluctuates around: where there is no signal from
detectable objects or artifacts (for example, galaxies, stars, planets
or comets, star spikes or internal optical ghost). Depending on the
dataset, the Sky value maybe a fixed value over the whole dataset, or it
may vary based on location. For an example of the latter case, see
Figure 11 in Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664).

Because of the significance of the Sky value in astronomical data
analysis, we have devoted this subsection to it for a thorough review.
We start with a thorough discussion on its definition (*note Sky value
definition::). In the astronomical literature, researchers use a
variety of methods to estimate the Sky value, so in *note Sky value
misconceptions::) we review those and discuss their biases. From the
definition of the Sky value, the most accurate way to estimate the Sky
value is to run a detection algorithm (for example, *note NoiseChisel::)
over the dataset and use the undetected pixels. However, there is also
a more crude method that maybe useful when good direct detection is not
initially possible (for example, due to too many cosmic rays in a
shallow image). A more crude (but simpler method) that is usable in
such situations is discussed in *note Quantifying signal in a tile::.

* Menu:

* Sky value definition:: Definition of the Sky/reference value.
* Sky value misconceptions:: Wrong methods to estimate the Sky value.
* Quantifying signal in a tile:: Method to estimate the presence of signal.

File: gnuastro.info, Node: Sky value definition, Next: Sky value misconceptions, Prev: Sky value, Up: Sky value

7.1.4.1 Sky value definition
............................

This analysis is taken from Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664). Let's assume that all instrument
defects - bias, dark and flat - have been corrected and the magnitude
(see *note Brightness flux magnitude::) of a detected object, $O$, is
desired. The sources of flux on pixel(1) $i$ of the image can be
written as follows:

• Contribution from the target object ($O_i$).
• Contribution from other detected objects ($D_i$).
• Undetected objects or the fainter undetected regions of bright
objects ($U_i$).
• A cosmic ray ($C_i$).
• The background flux, which is defined to be the count if none of
the others exists on that pixel ($B_i$).
The total flux in this pixel ($T_i$) can thus be written as:

$$T_i=B_i+D_i+U_i+C_i+O_i.$$

By definition, $D_i$ is detected and it can be assumed that it is
correctly estimated (deblended) and subtracted, we can thus set $D_i=0$.
There are also methods to detect and remove cosmic rays, for example,
the method described in van Dokkum (2001)(2), or by comparing multiple
exposures. This allows us to set $C_i=0$. Note that in practice, $D_i$
and $U_i$ are correlated, because they both directly depend on the
detection algorithm and its input parameters. Also note that no
detection or cosmic ray removal algorithm is perfect. With these
limitations in mind, the observed Sky value for this pixel ($S_i$) can
be defined as

$$S_i\equiv{}B_i+U_i.$$

Therefore, as the detection process (algorithm and input parameters)
becomes more accurate, or $U_i\to0$, the Sky value will tend to the
background value or $S_i\to B_i$. Hence, we see that while $B_i$ is an
inherent property of the data (pixel in an image), $S_i$ depends on the
detection process. Over a group of pixels, for example, in an image or
part of an image, this equation translates to the average of undetected
pixels (Sky$=\sum{S_i}$). With this definition of Sky, the object flux
in the data can be calculated, per pixel, with

$$T_{i}=S_{i}+O_{i} \quad\rightarrow\quad
O_{i}=T_{i}-S_{i}.$$

In the fainter outskirts of an object, a very small fraction of the
photo-electrons in a pixel actually belongs to objects, the rest is
caused by random factors (noise), see Figure 1b in Akhlaghi and Ichikawa
2015 (https://arxiv.org/abs/1505.01664). Therefore even a small over
estimation of the Sky value will result in the loss of a very large
portion of most galaxies. Besides the lost area/brightness, this will
also cause an over-estimation of the Sky value and thus even more
under-estimation of the object's magnitude. It is thus very important
to detect the diffuse flux of a target, even if they are not your
primary target.

In summary, the more accurately the Sky is measured, the more
accurately the magnitude (calculated from the sum of pixel values) of
the target object can be measured (photometry). Any
under/over-estimation in the Sky will directly translate to an
over/under-estimation of the measured object's magnitude.

The *Sky value* is only correctly found when all the detected objects
($D_i$ and $C_i$) have been removed from the data.

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

(1) For this analysis the dimension of the data (image) is
irrelevant. So if the data is an image (2D) with width of $w$ pixels,
then a pixel located on column $x$ and row $y$ (where all counting
starts from zero and (0, 0) is located on the bottom left corner of the
image), would have an index: $i=x+y\times{}w$.

(2) van Dokkum, P. G. (2001). Publications of the Astronomical
Society of the Pacific. 113, 1420.

File: gnuastro.info, Node: Sky value misconceptions, Next: Quantifying signal in a tile, Prev: Sky value definition, Up: Sky value

7.1.4.2 Sky value misconceptions
................................

As defined in *note Sky value::, the sky value is only accurately
defined when the detection algorithm is not significantly reliant on the
sky value. In particular its detection threshold. However, most
signal-based detection tools(1) use the sky value as a reference to
define the detection threshold. These older techniques therefore had to
rely on approximations based on other assumptions about the data. A
review of those other techniques can be seen in Appendix A of Akhlaghi
and Ichikawa 2015 (https://arxiv.org/abs/1505.01664).

These methods were extensively used in astronomical data analysis for
several decades, therefore they have given rise to a lot of
misconceptions, ambiguities and disagreements about the sky value and
how to measure it. As a summary, the major methods used until now were
an approximation of the mode of the image pixel distribution and
$\sigma$-clipping.

• To find the mode of a distribution those methods would either have
to assume (or find) a certain probability density function (PDF) or
use the histogram. But astronomical datasets can have any
distribution, making it almost impossible to define a generic
function. Also, histogram-based results are very inaccurate (there
is a large dispersion) and it depends on the histogram bin-widths.
Generally, the mode of a distribution also shifts as signal is
added. Therefore, even if it is accurately measured, the mode is a
biased measure for the Sky value.

• Another approach was to iteratively clip the brightest pixels in
the image (which is known as $\sigma$-clipping). See *note Sigma
clipping:: for a complete explanation. $\sigma$-clipping is useful
when there are clear outliers (an object with a sharp edge in an
image for example). However, real astronomical objects have
diffuse and faint wings that penetrate deeply into the noise, see
Figure 1 in Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664).

As discussed in *note Sky value::, the sky value can only be
correctly defined as the average of undetected pixels. Therefore all
such approaches that try to approximate the sky value prior to detection
are ultimately poor approximations.

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

(1) According to Akhlaghi and Ichikawa (2015), signal-based detection
is a detection process that relies heavily on assumptions about the
to-be-detected objects. This method was the most heavily used technique
prior to the introduction of NoiseChisel in that paper.

File: gnuastro.info, Node: Quantifying signal in a tile, Prev: Sky value misconceptions, Up: Sky value

7.1.4.3 Quantifying signal in a tile
....................................

In order to define detection thresholds on the image, or calibrate it
for measurements (subtract the signal of the background sky and define
errors), we need some basic measurements. For example, the quantile
threshold in NoiseChisel (‘--qthresh’ option), or the mean of the
undetected regions (Sky) and the Sky standard deviation (Sky STD) which
are the output of NoiseChisel and Statistics. But astronomical images
will contain a lot of stars and galaxies that will bias those
measurements if not properly accounted for. Quantifying where signal is
present is thus a very important step in the usage of a dataset; for
example, if the Sky level is over-estimated, your target object's
magnitude will be under-estimated.

Let's start by clarifying some definitions: _Signal_ is defined as
the non-random source of flux in each pixel (you can think of this as
the mean in a Gaussian or Poisson distribution). In astronomical
images, signal is mostly photons coming of a star or galaxy, and counted
in each pixel. _Noise_ is defined as the random source of flux in each
pixel (or the standard deviation of a Gaussian or Poisson distribution).
Noise is mainly due to counting errors in the detector electronics upon
data collection. _Data_ is defined as the combination of signal and
noise (so a noisy image of a galaxy is one _data_set).

When a dataset does not have any signal (for example, you take an
image with a closed shutter, producing an image that only contains
noise), the mean, median and mode of the distribution are equal within
statistical errors. Signal from emitting objects, like astronomical
targets, always has a positive value and will never become negative, see
Figure 1 in Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664). Therefore, when signal is added to
the data (you take an image with an open shutter pointing to a galaxy
for example), the mean, median and mode of the dataset shift to the
positive, creating a positively skewed distribution. The shift of the
mean is the largest. The median shifts less, since it is defined after
ordering all the elements/pixels (the median is the value at a quantile
of 0.5), thus it is not affected by outliers. Finally, the mode's shift
to the positive is the least.

Inverting the argument above gives us a robust method to quantify the
significance of signal in a dataset: when the mean and median of a
distribution are approximately equal we can argue that there is no
significant signal. In other words: when the quantile of the mean
($q_{mean}$) is around 0.5. This definition of skewness through the
quantile of the mean is further introduced with a real image the
tutorials, see *note Skewness caused by signal and its measurement::.

However, in an astronomical image, some of the pixels will contain
more signal than the rest, so we cannot simply check $q_{mean}$ on the
whole dataset. For example, if we only look at the patch of pixels that
are placed under the central parts of the brightest stars in the field
of view, $q_{mean}$ will be very high. The signal in other parts of the
image will be weaker, and in some parts it will be much smaller than the
noise (for example, 1/100-th of the noise level). When the
signal-to-noise ratio is very small, we can generally assume no signal
(because its effectively impossible to measure it) and $q_{mean}$ will
be approximately 0.5.

To address this problem, we break the image into a grid of tiles(1)
(see *note Tessellation::). For example, a tile can be a square box of
size $30\times30$ pixels. By measuring $q_{mean}$ on each tile, we can
find which tiles that contain significant signal and ignore them.
Technically, if a tile's $|q_{mean}-0.5|$ is larger than the value given
to the ‘--meanmedqdiff’ option, that tile will be ignored for the next
steps. You can read this option as "mean-median-quantile-difference".

The raw dataset's pixel distribution (in each tile) is noisy, to
decrease the noise/error in estimating $q_{mean}$, we convolve the image
before tessellation (see *note Convolution process::. Convolution
decreases the range of the dataset and enhances its skewness, See
Section 3.1.1 and Figure 4 in Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664). This enhanced skewness can be
interpreted as an increase in the Signal to noise ratio of the objects
buried in the noise. Therefore, to obtain an even better measure of the
presence of signal in a tile, the mean and median discussed above are
measured on the convolved image.

There is one final hurdle: raw astronomical datasets are commonly
peppered with Cosmic rays. Images of Cosmic rays are not smoothed by
the atmosphere or telescope aperture, so they have sharp boundaries.
Also, since they do not occupy too many pixels, they do not affect the
mode and median calculation. But their very high values can greatly
bias the calculation of the mean (recall how the mean shifts the fastest
in the presence of outliers), for example, see Figure 15 in Akhlaghi and
Ichikawa 2015 (https://arxiv.org/abs/1505.01664). The effect of
outliers like cosmic rays on the mean and standard deviation can be
removed through $\sigma$-clipping, see *note Sigma clipping:: for a
complete explanation.

Therefore, after asserting that the mean and median are approximately
equal in a tile (see *note Tessellation::), the Sky and its STD are
measured on each tile after $\sigma$-clipping with the ‘--sigmaclip’
option (see *note Sigma clipping::). In the end, some of the tiles will
pass the test and will be given a value. Others (that had signal in
them) will just be assigned a NaN (not-a-number) value. But we need a
measurement over each tile (and thus pixel). We will therefore use
interpolation to assign a value to the NaN tiles.

However, prior to interpolating over the failed tiles, another point
should be considered: large and extended galaxies, or bright stars, have
wings which sink into the noise very gradually. In some cases, the
gradient over these wings can be on scales that is larger than the tiles
(for example, the pixel value changes by $0.1\sigma$ after 100 pixels,
but the tile has a width of 30 pixels).

In such cases, the $q_{mean}$ test will be successful, even though
there is signal. Recall that $q_{mean}$ is a measure of skewness. If
we do not identify (and thus set to NaN) such outlier tiles before the
interpolation, the photons of the outskirts of the objects will leak
into the detection thresholds or Sky and Sky STD measurements and bias
our result, see *note Detecting large extended targets::. Therefore,
the final step of "quantifying signal in a tile" is to look at this
distribution of successful tiles and remove the outliers.
$\sigma$-clipping is a good solution for removing a few outliers, but
the problem with outliers of this kind is that there may be many such
tiles (depending on the large/bright stars/galaxies in the image). We
therefore apply the following local outlier rejection strategy.

For each tile, we find the nearest $N_{ngb}$ tiles that had a usable
value ($N_{ngb}$ is the value given to ‘--outliernumngb’). We then sort
them and find the difference between the largest and second-to-smallest
elements (The minimum is not used because the scatter can be large).
Let's call this the tile's _slope_ (measured from its neighbors). All
the tiles that are on a region of flat noise will have similar slope
values, but if a few tiles fall on the wings of a bright star or large
galaxy, their slope will be significantly larger than the tiles with no
signal. We just have to find the smallest tile slope value that is an
outlier compared to the rest, and reject all tiles with a slope larger
than that.

To identify the smallest outlier, we will use the distribution of
distances between sorted elements. Let's assume the total number of
tiles with a good mean-median quantile difference is $N$. They are
first sorted and searching for the outlier starts on element $N/3$
(integer division). Let's take $v_i$ to be the $i$-th element of the
sorted input (with no blank values) and $m$ and $\sigma$ as the
$\sigma$-clipped median and standard deviation from the distances of the
previous $N/3-1$ elements (not including $v_i$). If the value given to
‘--outliersigma’ is displayed with $s$, the $i$-th element is considered
as an outlier when the condition below is true.

$${(v_i-v_{i-1})-m\over \sigma}>s$$

Since $i$ begins from the $N/3$-th element in the sorted array (a
quantile of $1/3=0.33$), the outlier has to be larger than the $0.33$
quantile value of the dataset (this is usually the case; otherwise, it
is hard to define it as an "outlier"!).

Once the outlying tiles have been successfully identified and set to
NaN, we use nearest-neighbor interpolation to give a value to all tiles
in the image. We do not use parametric interpolation methods (like
bicubic), because they will effectively extrapolate on the edges,
creating strong artifacts. Nearest-neighbor interpolation is very
simple: for each tile, we find the $N_{ngb}$ nearest tiles that had a
good value, the tile's value is found by estimating the median. You can
set $N_{ngb}$ through the ‘--interpnumngb’ option. Once all the tiles
are given a value, a smoothing step is implemented to remove the sharp
value contrast that can happen on the edges of tiles. The size of the
smoothing box is set with the ‘--smoothwidth’ option.

As mentioned above, the process above is used for any of the basic
measurements (for example, identifying the quantile-based thresholds in
NoiseChisel, or the Sky value in Statistics). You can use the
check-image feature of NoiseChisel or Statistics to inspect the steps
and visually see each step (all the options that start with ‘--check’).
For example, as mentioned in the *note NoiseChisel optimization::
tutorial, when given a dataset from a new instrument (with differing
noise properties), we highly recommend to use ‘--checkqthresh’ in your
first call and visually inspect how the parameters above affect the
final quantile threshold (e.g., have the wings of bright sources leaked
into the threshold?). The same goes for the ‘--checksky’ option of
Statistics or NoiseChisel.

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

(1) The options to customize the tessellation are discussed in *note
Processing options::.

File: gnuastro.info, Node: Invoking aststatistics, Prev: Sky value, Up: Statistics

7.1.5 Invoking Statistics
-------------------------

Statistics will print statistical measures of an input dataset (table
column or image). The executable name is ‘aststatistics’ with the
following general template

$ aststatistics [OPTION ...] InputImage.fits

One line examples:

## Print some general statistics of input image:
$ aststatistics image.fits

## Print some general statistics of column named MAG_F160W:
$ aststatistics catalog.fits -h1 --column=MAG_F160W

## Make the histogram of the column named MAG_F160W:
$ aststatistics table.fits -cMAG_F160W --histogram

## Find the Sky value on image with a given kernel:
$ aststatistics image.fits --sky --kernel=kernel.fits

## Print Sigma-clipped results of records with a MAG_F160W
## column value between 26 and 27:
$ aststatistics cat.fits -cMAG_F160W -g26 -l27 --sigmaclip=3,0.2

## Find the polynomial (to third order) that best fits the X and Y
## columns of 'table.fits'. Robust fitting will be used to reject
## outliers. Also, estimate the fitted polynomial on the same input
## column (with errors).
$ aststatistics table.fits --fit=polynomial-robust --fitmaxpower=3 \
-cX,Y --fitestimate=self --output=estimated.fits

## Print the median value of all records in column MAG_F160W that
## have a value larger than 3 in column PHOTO_Z:
$ aststatistics tab.txt -rPHOTO_Z -g3 -cMAG_F160W --median

## Calculate the median of the third column in the input table, but only
## for rows where the mean of the first and second columns is >5.
$ awk '($1+$2)/2 > 5 {print $3}' table.txt | aststatistics --median

Statistics can take its input dataset either from a file (image or
table) or the Standard input (see *note Standard input::). If any
output file is to be created, the value to the ‘--output’ option, is
used as the base name for the generated files. Without ‘--output’, the
input name will be used to generate an output name, see *note Automatic
output::. The options described below are particular to Statistics, but
for general operations, it shares a large collection of options with the
other Gnuastro programs, see *note Common options:: for the full list.
For more on reading from standard input, please see the description of
‘--stdintimeout’ option in *note Input output options::. Options can
also be given in configuration files, for more, please see *note
Configuration files::.

The input dataset may have blank values (see *note Blank pixels::),
in this case, all blank pixels are ignored during the calculation.
Initially, the full dataset will be read, but it is possible to select a
specific range of data elements to use in the analysis of each run. You
can either directly specify a minimum and maximum value for the range of
data elements to use (with ‘--greaterequal’ or ‘--lessthan’), or specify
the range using quantiles (with ‘--qrange’). If a range is specified,
all pixels outside of it are ignored before any processing.

When no operation is requested, Statistics will print some general
basic properties of the input dataset on the command-line like the
example below (ran on one of the output images of ‘make check’(1)).
This default behavior is designed to help give you a general feeling of
how the data are distributed and help in narrowing down your analysis.

$ aststatistics convolve_spatial_scaled_noised.fits \
--greaterequal=9500 --lessthan=11000
Statistics (GNU Astronomy Utilities) X.X
-------
Input: convolve_spatial_scaled_noised.fits (hdu: 0)
Range: from (inclusive) 9500, upto (exclusive) 11000.
Unit: counts
-------
Number of elements: 9074
Minimum: 9622.35
Maximum: 10999.7
Mode: 10055.45996
Mode quantile: 0.4001983908
Median: 10093.7
Mean: 10143.98257
Standard deviation: 221.80834
-------
Histogram:
| **
| ******
| *******
| *********
| *************
| **************
| ******************
| ********************
| *************************** *
| ***************************************** ***
|* **************************************************************
|-----------------------------------------------------------------

Gnuastro's Statistics is a very general purpose program, so to be
able to easily understand this diversity in its operations (and how to
possibly run them together), we will divided the operations into two
types: those that do not respect the position of the elements and those
that do (by tessellating the input on a tile grid, see *note
Tessellation::). The former treat the whole dataset as one and can
re-arrange all the elements (for example, sort them), but the former do
their processing on each tile independently. First, we will review the
operations that work on the whole dataset.

The group of options below can be used to get single value
measurement(s) of the whole dataset. They will print only the requested
value as one field in a line/row, like the ‘--mean’, ‘--median’ options.
These options can be called any number of times and in any order. The
outputs of all such options will be printed on one line following each
other (with a space character between them). This feature makes these
options very useful in scripts, or to redirect into programs like GNU
AWK for higher-level processing. These are some of the most basic
measures, Gnuastro is still under heavy development and this list will
grow. If you want another statistical parameter, please contact us and
we will do out best to add it to this list, see *note Suggest new
feature::.

* Menu:

* Input to Statistics:: How to specify the inputs to Statistics.
* Single value measurements:: Can be used together (like -mean, or -maximum).
* Generating histograms and cumulative frequency plots:: Histogram and CFP tables.
* Fitting options:: Least squares fitting.
* Contour options:: Table of contours.
* Statistics on tiles:: Possible to do single-valued measurements on tiles.

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

(1) You can try it by running the command in the ‘tests’ directory,
open the image with a FITS viewer and have a look at it to get a sense
of how these statistics relate to the input image/dataset.

File: gnuastro.info, Node: Input to Statistics, Next: Single value measurements, Prev: Invoking aststatistics, Up: Invoking aststatistics

7.1.5.1 Input to Statistics
...........................

The following set of options are for specifying the input/outputs of
Statistics. There are many other input/output options that are common
to all Gnuastro programs including Statistics, see *note Input output
options:: for those.

‘-c STR/INT’
‘--column=STR/INT’
The column to use when the input file is a table with more than one
column. See *note Selecting table columns:: for a full description
of how to use this option. For more on how tables are read in
Gnuastro, please see *note Tables::.

‘-g FLT’
‘--greaterequal=FLT’
Limit the range of inputs into those with values greater and equal
to what is given to this option. None of the values below this
value will be used in any of the processing steps below.

‘-l FLT’
‘--lessthan=FLT’
Limit the range of inputs into those with values less-than what is
given to this option. None of the values greater or equal to this
value will be used in any of the processing steps below.

‘-Q FLT[,FLT]’
‘--qrange=FLT[,FLT]’
Specify the range of usable inputs using the quantile. This option
can take one or two quantiles to specify the range. When only one
number is input (let's call it $Q$), the range will be those values
in the quantile range $Q$ to $1-Q$. So when only one value is
given, it must be less than 0.5. When two values are given, the
first is used as the lower quantile range and the second is used as
the larger quantile range.

The quantile of a given element in a dataset is defined by the
fraction of its index to the total number of values in the sorted
input array. So the smallest and largest values in the dataset
have a quantile of 0.0 and 1.0. The quantile is a very useful
non-parametric (making no assumptions about the input) relative
measure to specify a range. It can best be understood in terms of
the cumulative frequency plot, see *note Histogram and Cumulative
Frequency Plot::. The quantile of each horizontal axis value in
the cumulative frequency plot is the vertical axis value associate
with it.

File: gnuastro.info, Node: Single value measurements, Next: Generating histograms and cumulative frequency plots, Prev: Input to Statistics, Up: Invoking aststatistics

7.1.5.2 Single value measurements
.................................

‘-n’
‘--number’
Print the number of all used (non-blank and in range) elements.

‘--minimum’
Print the minimum value of all used elements.

‘--maximum’
Print the maximum value of all used elements.

‘--sum’
Print the sum of all used elements.

‘-m’
‘--mean’
Print the mean (average) of all used elements.

‘-t’
‘--std’
Print the standard deviation of all used elements.

‘--mad’
Print the median absolute deviation (MAD) of all used elements.

‘-E’
‘--median’
Print the median of all used elements.

‘-u FLT[,FLT[,...]]’
‘--quantile=FLT[,FLT[,...]]’
Print the values at the given quantiles of the input dataset. Any
number of quantiles may be given and one number will be printed for
each. Values can either be written as a single number or as
fractions, but must be between zero and one (inclusive). Hence, in
effect ‘--quantile=0.25 --quantile=0.75’ is equivalent to
‘--quantile=0.25,3/4’, or ‘-u1/4,3/4’.

The returned value is one of the elements from the dataset. Taking
$q$ to be your desired quantile, and $N$ to be the total number of
used (non-blank and within the given range) elements, the returned
value is at the following position in the sorted array:
$round(q\times{}N$).

‘--quantfunc=FLT[,FLT[,...]]’
Print the quantiles of the given values in the dataset. This
option is the inverse of the ‘--quantile’ and operates similarly
except that the acceptable values are within the range of the
dataset, not between 0 and 1. Formally it is known as the
"Quantile function".

Since the dataset is not continuous this function will find the
nearest element of the dataset and use its position to estimate the
quantile function.

‘--quantofmean’
Print the quantile of the mean in the dataset. This is a very good
measure of detecting skewness or outliers. The concept is used by
programs like NoiseChisel to identify the presence of signal in a
tile of the image (because signal in noise causes skewness).

For example, take this simple array: ‘1 2 20 4 5 6 3’. The mean is
‘5.85’. The nearest element to this mean is ‘6’ and the quantile
of ‘6’ in this distribution is 0.8333. Here is how we got to this:
in the sorted dataset (‘1 2 3 4 5 6 20’), ‘6’ is the 5-th element
(counting from zero, since a quantile of zero corresponds to the
minimum, by definition) and the maximum is the 6-th element (again,
counting from zero). So the quantile of the mean in this case is
$5/6=0.8333$.

In the example above, if we had ‘7’ instead of ‘20’ (which was an
outlier), then the mean would be ‘4’ and the quantile of the mean
would be 0.5 (which by definition, is the quantile of the median),
showing no outliers. As the number of elements increases, the mean
itself is less affected by a small number of outliers, but skewness
can be nicely identified by the quantile of the mean.

‘-O’
‘--mode’
Print the mode of all used elements. The mode is found through the
mirror distribution which is fully described in Appendix C of
Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664). See
that section for a full description.

This mode calculation algorithm is non-parametric, so when the
dataset is not large enough (larger than about 1000 elements
usually), or does not have a clear mode it can fail. In such
cases, this option will return a value of ‘nan’ (for the floating
point NaN value).

As described in that paper, the easiest way to assess the quality
of this mode calculation method is to use it's symmetricity (see
‘--modesym’ below). A better way would be to use the ‘--mirror’
option to generate the histogram and cumulative frequency tables
for any given mirror value (the mode in this case) as a table. If
you generate plots like those shown in Figure 21 of that paper,
then your mode is accurate.

‘--modequant’
Print the quantile of the mode. You can get the actual mode value
from the ‘--mode’ described above. In many cases, the absolute
value of the mode is irrelevant, but its position within the
distribution is important. In such cases, this option will become
handy.

‘--modesym’
Print the symmetricity of the calculated mode. See the description
of ‘--mode’ for more. This mode algorithm finds the mode based on
how symmetric it is, so if the symmetricity returned by this option
is too low, the mode is not too accurate. See Appendix C of
Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664) for a
full description. In practice, symmetricity values larger than 0.2
are mostly good.

‘--modesymvalue’
Print the value in the distribution where the mirror and input
distributions are no longer symmetric, see ‘--mode’ and Appendix C
of Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664)
for more.

‘--sigclip-std’
‘--sigclip-mad’
‘--sigclip-mean’
‘--sigclip-number’
‘--sigclip-median’
Calculate the desired statistic after applying $\sigma$-clipping
(see *note Sigma clipping::, part of the tutorial *note Clipping
outliers::). $\sigma$-clipping configuration is done with the
‘--sclipparams’ option.

Here is one scenario where this can be useful: assume you have a
table and you would like to remove the rows that are outliers (not
within the $\sigma$-clipping range). Let's assume your table is
called ‘table.fits’ and you only want to keep the rows that have a
value in ‘COLUMN’ within the $\sigma$-clipped range (to $3\sigma$,
with a tolerance of 0.1). This command will return the
$\sigma$-clipped median and standard deviation (used to define the
range later).

$ aststatistics table.fits -cCOLUMN --sclipparams=3,0.1 \
--sigclip-median --sigclip-std

You can then use the ‘--range’ option of Table (see *note Table::)
to select the proper rows. But for that, you need the actual
starting and ending values of the range ($m\pm s\sigma$; where $m$
is the median and $s$ is the multiple of sigma to define an
outlier). Therefore, the raw outputs of Statistics in the command
above are not enough.

To get the starting and ending values of the non-outlier range (and
put a '<,>' between them, ready to be used in ‘--range’), pipe the
result into AWK. But in AWK, we will also need the multiple of
$\sigma$, so we will define it as a shell variable (‘s’) before
calling Statistics (note how ‘$s’ is used two times now):

$ s=3
$ aststatistics table.fits -cCOLUMN --sclipparams=$s,0.1 \
--sigclip-median --sigclip-std \
| awk '{s='$s'; printf("%f,%f\n", $1-s*$2, $1+s*$2)}'

To pass it onto Table, we will need to keep the printed output from
the command above in another shell variable (‘r’), not print it.
In Bash, can do this by putting the whole statement within a ‘$()’:

$ s=3
$ r=$(aststatistics table.fits -cCOLUMN --sclipparams=$s,0.1 \
--sigclip-median --sigclip-std \
| awk '{s='$s'; printf("%f,%f\n", $1-s*$2, $1+s*$2)}')
$ echo $r # Just to confirm.

Now you can use Table with the ‘--range’ option to only print the
rows that have a value in ‘COLUMN’ within the desired range:

$ asttable table.fits --range=COLUMN,$r

To save the resulting table (that is clean of outliers) in another
file (for example, named ‘cleaned.fits’, it can also have a ‘.txt’
suffix), just add ‘--output=cleaned.fits’ to the command above.

‘--madclip-std’
‘--madclip-mad’
‘--madclip-mean’
‘--madclip-number’
‘--madclip-median’
Calculate the desired statistic after applying median absolute
deviation (MAD) clipping (see *note MAD clipping::, part of the
tutorial *note Clipping outliers::). MAD-clipping configuration is
done with the ‘--mclipparams’ option.

This option behaves similarly to ‘--sigclip-*’ options, read their
description for usage examples.

File: gnuastro.info, Node: Generating histograms and cumulative frequency plots, Next: Fitting options, Prev: Single value measurements, Up: Invoking aststatistics

7.1.5.3 Generating histograms and cumulative freq.
..................................................

The list of options below are for those statistical operations that
output more than one value. So while they can be called together in one
run, their outputs will be distinct (each one's output will usually be
printed in more than one line).

‘-A’
‘--asciihist’
Print an ASCII histogram of the usable values within the input
dataset along with some basic information like the example below
(from the UVUDF catalog(1)). The width and height of the histogram
(in units of character widths and heights on your command-line
terminal) can be set with the ‘--numasciibins’ (for the width) and
‘--asciiheight’ options.

For a full description of the histogram, please see *note Histogram
and Cumulative Frequency Plot::. An ASCII plot is certainly very
crude and cannot be used in any publication, but it is very useful
for getting a general feeling of the input dataset very fast and
easily on the command-line without having to take your hands off
the keyboard (which is a major distraction!). If you want to try
it out, you can write it all in one line and ignore the <\> and
extra spaces.

$ aststatistics uvudf_rafelski_2015.fits.gz --hdu=1 \
--column=MAG_F160W --lessthan=40 \
--asciihist --numasciibins=55
ASCII Histogram:
Number: 8593
Y: (linear: 0 to 660)
X: (linear: 17.7735 -- 31.4679, in 55 bins)
| ****
| *****
| ******
| ********
| *********
| ***********
| **************
| *****************
| ***********************
| ********************************
|*** ***************************************************
|-------------------------------------------------------

‘--asciicfp’
Print the cumulative frequency plot of the usable elements in the
input dataset. Please see descriptions under ‘--asciihist’ for
more, the example below is from the same input table as that
example. To better understand the cumulative frequency plot,
please see *note Histogram and Cumulative Frequency Plot::.

$ aststatistics uvudf_rafelski_2015.fits.gz --hdu=1 \
--column=MAG_F160W --lessthan=40 \
--asciicfp --numasciibins=55
ASCII Cumulative frequency plot:
Y: (linear: 0 to 8593)
X: (linear: 17.7735 -- 31.4679, in 55 bins)
| *******
| **********
| ***********
| *************
| **************
| ***************
| *****************
| *******************
| ***********************
| ******************************
|*******************************************************
|-------------------------------------------------------

‘-H’
‘--histogram’
Save the histogram of the usable values in the input dataset into a
table. The first column is the value at the center of the bin and
the second is the number of points in that bin. If the
‘--cumulative’ option is also called with this option in a run,
then the table will have three columns (the third is the cumulative
frequency plot). Through the ‘--numbins’, ‘--onebinstart’, or
‘--manualbinrange’, you can modify the first column values and with
‘--normalize’ and ‘--maxbinone’ you can modify the second columns.
See below for the description of each.

By default (when no ‘--output’ is specified) a plain text table
will be created, see *note Gnuastro text table format::. If a FITS
name is specified, you can use the common option ‘--tableformat’ to
have it as a FITS ASCII or FITS binary format, see *note Common
options::. This table can then be fed into your favorite plotting
tool and get a much more clean and nice histogram than what the raw
command-line can offer you (with the ‘--asciihist’ option).

‘--histogram2d’
Save the 2D histogram of two input columns into an output file, see
*note 2D Histograms::. The output will have three columns: the
first two are the coordinates of each box's center in the first and
second dimensions/columns. The third will be number of input
points that fall within that box.

‘-C’
‘--cumulative’
Save the cumulative frequency plot of the usable values in the
input dataset into a table, similar to ‘--histogram’.

‘--madclip’
Do median absolute deviation (MAD) clipping on the usable pixels of
the input dataset. See *note MAD clipping:: for a description on
MAD-clipping and *note Clipping outliers:: for a complete tutorial
on clipping of outliers. The MAD-clipping parameters can be set
through the ‘--mclipparams’ option (see below).

‘-s’
‘--sigmaclip’
Do $\sigma$-clipping on the usable pixels of the input dataset.
See *note Sigma clipping:: for a full description on
$\sigma$-clipping and *note Clipping outliers:: for a complete
tutorial on clipping of outliers. The $\sigma$-clipping parameters
can be set through the ‘--sclipparams’ option (see below).

‘--mirror=FLT’
Make a histogram and cumulative frequency plot of the mirror
distribution for the given dataset when the mirror is located at
the value to this option. The mirror distribution is fully
described in Appendix C of Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664) and currently it is only used to
calculate the mode (see ‘--mode’).

Just note that the mirror distribution is a discrete distribution
like the input, so while you may give any number as the value to
this option, the actual mirror value is the closest number in the
input dataset to this value. If the two numbers are different,
Statistics will warn you of the actual mirror value used.

This option will make a table as output. Depending on your
selected name for the output, it will be either a FITS table or a
plain text table (which is the default). It contains three
columns: the first is the center of the bins, the second is the
histogram (with the largest value set to 1) and the third is the
normalized cumulative frequency plot of the mirror distribution.
The bins will be positioned such that the mode is on the starting
interval of one of the bins to make it symmetric around the mirror.
With this output file and the input histogram (that you can
generate in another run of Statistics, using the ‘--onebinvalue’),
it is possible to make plots like Figure 21 of Akhlaghi and
Ichikawa 2015 (https://arxiv.org/abs/1505.01664).

The list of options below allow customization of the histogram and
cumulative frequency plots (for the ‘--histogram’, ‘--cumulative’,
‘--asciihist’, and ‘--asciicfp’ options).

‘--numbins’
The number of bins (rows) to use in the histogram and the
cumulative frequency plot tables (outputs of ‘--histogram’ and
‘--cumulative’).

‘--numasciibins’
The number of bins (characters) to use in the ASCII plots when
printing the histogram and the cumulative frequency plot (outputs
of ‘--asciihist’ and ‘--asciicfp’).

‘--asciiheight’
The number of lines to use when printing the ASCII histogram and
cumulative frequency plot on the command-line (outputs of
‘--asciihist’ and ‘--asciicfp’).

‘-n’
‘--normalize’
Normalize the histogram or cumulative frequency plot tables
(outputs of ‘--histogram’ and ‘--cumulative’). For a histogram,
the sum of all bins will become one and for a cumulative frequency
plot the last bin value will be one.

‘--maxbinone’
Divide all the histogram values by the maximum bin value so it
becomes one and the rest are similarly scaled. In some situations
(for example, if you want to plot the histogram and cumulative
frequency plot in one plot) this can be very useful.

‘--onebinstart=FLT’
Make sure that one bin starts with the value to this option. In
practice, this will shift the bins used to find the histogram and
cumulative frequency plot such that one bin's lower interval
becomes this value.

For example, when a histogram range includes negative and positive
values and zero has a special significance in your analysis, then
zero might fall somewhere in one bin. As a result that bin will
have counts of positive and negative. By setting
‘--onebinstart=0’, you can make sure that one bin will only count
negative values in the vicinity of zero and the next bin will only
count positive ones in that vicinity.

Note that by default, the first row of the histogram and cumulative
frequency plot show the central values of each bin. So in the
example above you will not see the 0.000 in the first column, you
will see two symmetric values.

If the value is not within the usable input range, this option will
be ignored. When it is, this option is the last operation before
the bins are finalized, therefore it has a higher priority than
options like ‘--manualbinrange’.

‘--manualbinrange’
Use the values given to the ‘--greaterequal’ and ‘--lessthan’ to
define the range of all bin-based calculations like the histogram.
This option itself does not take any value, but just tells the
program to use the values of those two options instead of the
minimum and maximum values of a plot. If any of the two options
are not given, then the minimum or maximum will be used
respectively. Therefore, if none of them are called calling this
option is redundant.

The ‘--onebinstart’ option has a higher priority than this option.
In other words, ‘--onebinstart’ takes effect after the range has
been finalized and the initial bins have been defined, therefore it
has the power to (possibly) shift the bins. If you want to
manually set the range of the bins _and_ have one bin on a special
value, it is thus better to avoid ‘--onebinstart’.

‘--numbins2=INT’
Similar to ‘--numbins’, but for the second column when a 2D
histogram is requested, see ‘--histogram2d’.

‘--greaterequal2=FLT’
Similar to ‘--greaterequal’, but for the second column when a 2D
histogram is requested, see ‘--histogram2d’.

‘--lessthan2=FLT’
Similar to ‘--lessthan’, but for the second column when a 2D
histogram is requested, see ‘--histogram2d’.

‘--onebinstart2=FLT’
Similar to ‘--onebinstart’, but for the second column when a 2D
histogram is requested, see ‘--histogram2d’.

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

(1) <https://asd.gsfc.nasa.gov/UVUDF/uvudf_rafelski_2015.fits.gz>

File: gnuastro.info, Node: Fitting options, Next: Contour options, Prev: Generating histograms and cumulative frequency plots, Up: Invoking aststatistics

7.1.5.4 Fitting options
.......................

With the options below, you can customize the least squares fitting
features of Statistics. For a tutorial of the usage of least squares
fitting in Statistics, please see *note Least squares fitting::. Here,
we will just review the details of each option.

To activate least squares fitting in Statistics, it is necessary to
use the ‘--fit’ option to specify the type of fit you want to do. See
the description of ‘--fit’ for the various available fitting models.
The fitting models that account for weights require three input columns,
while the non-weighted ones only take two input columns. Here is a
summary of the input columns:

1. The first input column is assumed to be the independent variable
(on the horizontal axis of a plot, or $X$ in the equations of each
fit).
2. The second input column is assumed to be the measured value (on the
vertical axis of a plot, or $Y$ in the equation above).
3. The third input column is only for fittings with a weight. It is
assumed to be the "weight" of the measurement column. The nature
of the "weight" can be set with the ‘--fitweight’ option, for
example, if you have the standard deviation of the error in $Y$,
you can use ‘--fitweight=std’ (which is the default, so unless the
default value has been changed, you will not need to set this).

If three columns are given to a model without weight, or two columns
are given to a model that requires weights, Statistics will abort and
inform you. Below you can see an example of fitting with the same
linear model, once weighted and once without weights.

$ aststatistics table.fits --column=X,Y --fit=linear
$ aststatistics table.fits --column=X,Y,Yerr --fit=linear-weighted

The output of the fitting can be in three modes listed below. For a
complete example, see the tutorial in *note Least squares fitting::).
Human friendly format
By default (for example, the commands above) the output is an
elaborate description of the model parameters. For example, $c_0$
and $c_1$ in the linear model ($Y=c_0+c_1X$). Their covariance
matrix and the reduced $\chi^2$ of the fit are also printed on the
output.
Raw numbers
If you don't need the human friendly components of the output
(which are annoying when you want to parse the outputs in some
scenarios), you can use ‘--quiet’ option. Only the raw output
numbers will be printed.
Estimate on a custom X column
Through the ‘--fitestimate’ option, you can specify an independent
table column to estimate the fit (it can also take a single value).
See the description of this option for more.

‘-f STR’
‘--fit=STR’
The name of the fitting method to use. They are based on the
linear (https://www.gnu.org/software/gsl/doc/html/lls.html) and
nonlinear (https://www.gnu.org/software/gsl/doc/html/nls.html)
least-squares fitting functions of the GNU Scientific Library
(GSL).
‘linear’
$Y=c_0+c_1X$
‘linear-weighted’
$Y=c_0+c_1X$; accounting for "weights" in $Y$.
‘linear-no-constant’
$Y=c_1X$.
‘linear-no-constant-weighted’
$Y=c_1X$; accounting for "weights" in $Y$.
‘polynomial’
$Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n$; the maximum required power
($n$) is specified by ‘--fitmaxpower’.
‘polynomial-weighted’
$Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n$; accounting for "weights" in
$Y$. The maximum required power ($n$) is specified by
‘--fitmaxpower’.
‘polynomial-robust’
$Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n$; rejects outliers. The
function to use for outlier removal can be specified with the
‘--fitrobust’ option described below. This model doesn't take
weights since they are calculated internally based on the
outlier removal function (requires two input columns). The
maximum required power ($n$) is specified by ‘--fitmaxpower’.

For a comprehensive review of "robust" fitting and the
available functions, please see the Robust linear regression
(https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression)
section of the GNU Scientific Library.

‘--fitweight=STR’
The nature of the "weight" column (when a weight is necessary for
the model). It can take one of the following values:
‘std’
Standard deviation of each $Y$ axis measurement: this is the
usual "error" associated with a measurement (for example, in
*note MakeCatalog::) and is the default value to this option.
‘var’
Variance of each $Y$ axis measurement. Assuming a Gaussian
distribution with standard deviation $\sigma$, the variance is
$\sigma^2$.
‘inv-var’
Inverse variance of each $Y$ axis measurement. Assuming a
Gaussian distribution with standard deviation $\sigma$, the
variance is $1/\sigma^2$.

‘--fitmaxpower=INT’
The maximum power (an integer) in a polynomial ($n$ in
$Y=c_0+c_1X+c_2X^2+\cdots+c_nX^n$). This is only relevant when one
of the polynomial models is given to ‘--fit’. The fit will return
$n+1$ coefficients.

‘--fitrobust=STR’
The function for rejecting outliers in the ‘polynomial-robust’
fitting model. For a comprehensive review of "robust" fitting and
the available functions, please see the Robust linear regression
(https://www.gnu.org/software/gsl/doc/html/lls.html#robust-linear-regression)
section of the GNU Scientific Library. This function can take the
following values:
‘bisquare’
Tukey’s biweight (bisquare) function, this is the default
function. According to the GSL manual, this is a good general
purpose weight function.
‘cauchy’
Cauchy’s function (also known as the Lorentzian function). It
doesn't guarantee a unique solution, so it should be used with
care.
‘fair’
The fair function. It guarantees a unique solution and has
continuous derivatives to three orders.
‘huber’
Huber's $\rho$ function. This is also a good general purpose
weight function for rejecting outliers, but can cause
difficulty in some special scenarios.
‘ols’
Ordinary Least Squares (OLS) solution with a constant weight
of unity.
‘welsch’
Welsch function which is useful when the residuals follow an
exponential distribution.

‘--fitestimate=STR/FLT’
Estimate the fitted function at a single point or a complete column
of points. The input $X$ axis positions to estimate the function
can be specified in the following ways:
• A real number: the fitted function will be estimated at that
$X$ position and the corresponding $Y$ and its error will be
printed to standard output.
• ‘self’: in this mode, the same X axis column that was used in
the fit will be used for estimating the fitted function. This
can be useful to visually/easily check the fit, see *note
Least squares fitting::.
• A file name: If the value is none of the above, Statistics
expects it to be a file name containing a table. If the file
is a FITS file, the HDU containing the table should be
specified with the ‘--fitestimatehdu’ option. The column of
the table to use for the $X$ axis points should be specified
with the ‘--fitestimatecol’ option.
The output in this mode can be customized in the following ways:
• If a single floating point value is given ‘--fitestimate’, the
fitted function will be estimated on that point and printed to
standard output.
• When nothing is given to ‘--output’, the independent column
and the estimated values and errors are printed on the
standard output.
• If a file name is given to ‘--output’, the estimated table
above is saved in that file. It can have any of the formats
in *note Recognized table formats::. As a FITS file, all the
fit outputs (coefficients, covariance matrix and reduced
$\chi^2$) are kept as FITS keywords in the same HDU of the
estimated table. For a complete example, see *note Least
squares fitting::.

When the covariance matrix (and thus the $\chi^2$) cannot be
calculated (for example if you only have two rows!), the
printed values on the terminal will be NaN. However, the FITS
standard does not allow NaN values in keyword values!
Therefore, when writing the $\chi^2$ and covariance matrix
elements into the output FITS keywords, the largest value of
the 64-bit floating point type will be written:
$1.79769313486232\times10^{308}$; see *note Numeric data
types::.

• When ‘--quiet’ is given with ‘--fitestimate’, the fitted
parameters are no longer printed on the standard output; they
are available as FITS keywords in the file given to
‘--output’.

‘--fitestimatehdu=STR/INT’
HDU name or counter (counting from zero) that contains the table to
be used for the estimating the fitted function over many points
through ‘--fitestimate’. For more on selecting a HDU, see the
description of ‘--hdu’ in *note Input output options::.

‘--fitestimatecol=STR/INT’
Column name or counter (counting from one) that contains the table
to be used for the estimating the fitted function over many points
through ‘--fitestimate’. See *note Selecting table columns::.

File: gnuastro.info, Node: Contour options, Next: Statistics on tiles, Prev: Fitting options, Up: Invoking aststatistics

7.1.5.5 Contour options
.......................

Contours are useful to highlight the 2D shape of a certain flux level
over an image. To derive contours in Statistics, you can use the option
below:

‘-R FLT[,FLT[,FLT...]]’
‘--contour=FLT[,FLT[,FLT...]]’
Write the contours for the requested levels in a file ending with
‘_contour.txt’. It will have three columns: the first two are the
coordinates of each point and the third is the level it belongs to
(one of the input values). Each disconnected contour region will
be separated by a blank line. This is the requested format for
adding contours with PGFPlots in LaTeX. If any other format can be
useful for your work please let us know so we can add it. If the
image has World Coordinate System information, the written
coordinates will be in RA and Dec, otherwise, they will be in pixel
coordinates.

Note that currently, this is a very crude/simple implementation,
please let us know if you find problematic situations so we can fix
it.

File: gnuastro.info, Node: Statistics on tiles, Prev: Contour options, Up: Invoking aststatistics

7.1.5.6 Statistics on tiles
...........................

All the options described until now were from the first class of
operations discussed above: those that treat the whole dataset as one.
However, it often happens that the relative position of the dataset
elements over the dataset is significant. For example, you do not want
one median value for the whole input image, you want to know how the
median changes over the image. For such operations, the input has to be
tessellated (see *note Tessellation::). Thus this class of options
cannot currently be called along with the options above in one run of
Statistics.

‘-t’
‘--ontile’
Do the respective single-valued calculation over one tile of the
input dataset, not the whole dataset. This option must be called
with at least one of the single valued options discussed above (for
example, ‘--mean’ or ‘--quantile’). The output will be a file in
the same format as the input. If the ‘--oneelempertile’ option is
called, then one element/pixel will be used for each tile (see
*note Processing options::). Otherwise, the output will have the
same size as the input, but each element will have the value
corresponding to that tile's value. If multiple single valued
operations are called, then for each operation there will be one
extension in the output FITS file.

‘-y’
‘--sky’
Estimate the Sky value on each tile as fully described in *note
Quantifying signal in a tile::. As described in that section,
several options are necessary to configure the Sky estimation which
are listed below. The output file will have two extensions: the
first is the Sky value and the second is the Sky standard deviation
on each tile. Similar to ‘--ontile’, if the ‘--oneelempertile’
option is called, then one element/pixel will be used for each tile
(see *note Processing options::).

The parameters for estimating the sky value can be set with the
following options, except for the ‘--sclipparams’ option (which is also
used by the ‘--sigmaclip’), the rest are only used for the Sky value
estimation.

‘-k=FITS’
‘--kernel=FITS’
File name of kernel to help in estimating the significance of
signal in a tile, see *note Quantifying signal in a tile::.

‘--khdu=STR’
Kernel HDU to help in estimating the significance of signal in a
tile, see *note Quantifying signal in a tile::.

‘--meanmedqdiff=FLT’
The maximum acceptable distance between the quantiles of the mean
and median, see *note Quantifying signal in a tile::. The initial
Sky and its standard deviation 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.

‘--sclipparams=FLT,FLT’
The $\sigma$-clipping parameters, 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 a specific number. Hence, in the latter case
the value must be an integer.

‘--mclipparams=FLT,FLT’
The MAD-clipping parameters. This is very similar to
‘--sclipparams’ above, see there for more.

‘--outliersclip=FLT,FLT’
$\sigma$-clipping parameters for the outlier rejection of the Sky
value (similar to ‘--sclipparams’).

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 positive. For more, see
the latter half of *note Quantifying signal in a tile::.

‘--outliernumngb=INT’
Number of neighboring tiles to use for outlier rejection (mostly
the wings of bright stars or galaxies). If this option is given a
value of zero, no outlier rejection will take place. For more see
the latter half of *note Quantifying signal in a tile::.

‘--outliersigma=FLT’
Multiple of sigma to define an outlier in the Sky value estimation.
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::.

‘--smoothwidth=INT’
Width of a flat kernel to convolve the interpolated tile values.
Tile interpolation is done using the median of the ‘--interpnumngb’
neighbors of each tile (see *note Processing options::). If this
option is given a value of zero or one, no smoothing will be done.
Without smoothing, strong boundaries will probably be created
between the values estimated for each tile. It is thus good to
smooth the interpolated image so strong discontinuities do not show
up in the final Sky values. The smoothing is done through
convolution (see *note Convolution process::) with a flat kernel,
so the value to this option must be an odd number.

‘--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.

‘--checksky’
Create a multi-extension FITS file showing the steps that were used
to estimate the Sky value over the input, see *note Quantifying
signal in a tile::. The file will have two extensions for each
step (one for the Sky and one for the Sky standard deviation).

‘--checkskynointerp’
Similar to ‘--checksky’, but it will stop as soon as the outlier
tiles have been identified and before it interpolates the values to
cover the whole image.

This is useful when you want the good tile values before
interpolation, and don't want to slow down your pipeline with the
extra computing that interpolation and smoothing require.

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 datacubes, 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 datacubes. 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/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/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

In the sections below, NoiseChisel's options are classified into
three general classes to help in easy navigation. *note NoiseChisel
input:: mainly discusses the options relating to input and those that
are shared in both detection and segmentation. Options to configure the
detection are described in *note Detection options:: and *note
Segmentation options:: we discuss how you can fine-tune the segmentation
of the detections. Finally, in *note NoiseChisel output:: the format of
NoiseChisel's output is discussed. The order of options here follow the
same logical order that the respective action takes place within
NoiseChisel (note that the output of ‘--help’ is sorted alphabetically).

Below, we will discuss NoiseChisel's options, classified into two
general classes, to help in easy navigation. *note NoiseChisel input::
mainly discusses the basic options relating to inputs and prior to the
detection process detection. 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 datacube, 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 datacubes.

‘--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::).