The current stable release
 is Gnu=
astro
 0.6 (June 4th, 2018).
 Use =
a
+ is Gnu=
astro
+ 0.7 (August 8th, 2018).
+ Use =
a
mirror if possible.
=20
0.0511864
@end example
=20
@noindent
The outer wings where therefore nonparametrically detected until
@mymath{\rm{S/N}\approx0.05}.
+@mymath{\rm{S/N}\approx0.05}!
+
+This is very good, but how much is this in units of standard magnitudes =
per
+arcseconds squared? To find that, we'll first need to calcaulate how man=
y
+pixels of this image are in one arcsecondsquared. Fortunately the WCS
+headers of Gnuastro's output FITS files (and in particular @code{CDELT})
+give us this information.
+
+@example
+$ n=3D$(astfits r_detected.fits h1 \
+  awk '/CDELT1/ @{p=3D1/($3*3600); print p*p@}')
+$ echo $n
+@end example
+
+@noindent
+Now, let's calculate the average skysubtracted flux in the border regio=
n
+per pixel.
+
+@example
+$ f=3D$(astarithmetic r_detected.fits boundary.fits not nan where seti =
\
+ i sumvalue i numvalue / q hINPUTNOSKY)
+$ echo $f
+@end example
+
+@noindent
+We can just multiply the two to get the average flux on this border in o=
ne
+arcsecond squared. We also have the rband SDSS zeropoint
+magnitude@footnote{From
+@url{http://classic.sdss.org/dr7/algorithms/fluxcal.html}} to be
+24.80. Therefore we can get the surface brightness limit (in magnitudes =
per
+arcsecond squared) using the following command. Just note that @code{log=
}
+in AWK is in base2 (not 10) and it doesn't have a @code{log10}
+operator. So we'll do an extra division by @code{log(10)} to correct for
+this.
+
+@example
+$ z=3D24.80
+$ echo "$n $f $z"  awk '@{print 2.5*log($1*$2)/log(10)+$3@}'
+> 30.0646
+@end example
+
+This shows that on a singleexposure SDSS image, we have reached a surfa=
ce
+brightness limit of roughly 30 magnitudes per arcseconds squared!
=20
In interpreting this value, you should just have in mind that NoiseChise=
l
works based on the contiguity of signal in the pixels. Therefore the lar=
ger
@@ 12046,18 +12288,26 @@ $ astarithmetic img1.fits img2.fits img3.fits m=
edian \
h0 h1 h2 out=3Dmedian.fits
@end example
=20
If the output is an image, and the @option{output} option is not given=
,
automatic output will use the name of the first FITS image encountered t=
o
generate an output file name, see @ref{Automatic output}. Also, output W=
CS
information will be taken from the first input image encountered. When t=
he
output is a single number, that number will be printed in the standard
output and no output file will be created. Arithmetic's notation for giv=
ing
operands to operators is described in @ref{Reverse polish notation}. To
ignore certain pixels, set them as blank, see @ref{Blank pixels}, for
example with the @command{where} operator (see @ref{Arithmetic
operators}). See @ref{Common options} for a review of the options in all
Gnuastro programs. Arithmetic just redefines the @option{hdu} option a=
s
explained below:
+If the output is not a single number, and the @option{output} option i=
s
+not given, automatic output will use the name of the first FITS image
+encountered to generate an output file name, see @ref{Automatic
+output}. Also, output WCS information will be taken from the first input
+image encountered. When the output is a single number, that number will =
be
+printed in the standard output and no output file will be
+created. Arithmetic's notation for giving operands to operators is
+described in @ref{Reverse polish notation}. To ignore certain pixels, se=
t
+them as blank, see @ref{Blank pixels}, for example with the @command{whe=
re}
+operator (see @ref{Arithmetic operators}). See @ref{Common options} for =
a
+review of the options in all Gnuastro programs. Arithmetic just redefine=
s
+the @option{hdu} option as explained below.
+
+Through operators like those starting with @code{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 @option{tableformat} option, see @ref{Inpu=
t
+output options}). You can disable this feature (write 1D arrays as FITS
+images/arrays) with the @option{onedasimage} option.
=20
@table @option
=20
@@ 12094,6 +12344,12 @@ option is very convenient when you have many inp=
ut 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 @option{hdu} option (explained above) are igno=
red
and will not be used.
+
+@item O
+@itemx onedasimage
+When final dataset to write as output only has one dimension, write it a=
s a
+FITS image/array. By default, if the output is 1D, it will be written as=
a
+table, see above.
@end table
=20
Arithmetic accepts two kinds of input: images and numbers. Images are
@@ 14422,17 +14678,19 @@ 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 deriv=
ed
results of many highlevel 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
where there is no signal from objects (for example galaxies, stars, plan=
ets
or comets). 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 @url{https://arxiv.org/abs/1505.01664,
Akhlaghi and Ichikawa (2015)}.
+(in the attempt to associate our results with the ``real'' world).
+
+In astronomy, this reference value is known as the ``Sky'' value: the va=
lue
+that noise fluctuates around: where there is no signal from detectable
+objects or artifacts (for example galaxies, stars, planets or comets, st=
ar
+spikes or internal optical ghost). Depending on the dataset, the Sky val=
ue
+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
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
=20
Because of the significance of the Sky value in astronomical data analys=
is,
we have devoted this subsection to it for a thorough review. We start wi=
th
a thorough discussion on its definition (@ref{Sky value definition}). I=
n
+a thorough discussion on its definition (@ref{Sky value definition}). In
the astronomical literature, researchers use a variety of methods to
estimate the Sky value, so in @ref{Sky value misconceptions}) we review
those and discuss their biases. From the definition of the Sky value, th=
e
@@ 14458,28 +14716,29 @@ This analysis is taken from @url{https://arxiv.=
org/abs/1505.01664, Akhlaghi
and Ichikawa (2015)}. Let's assume that all instrument defects  bias,
dark and flat  have been corrected and the brightness (see @ref{Flux
Brightness and magnitude}) of a detected object, @mymath{O}, is
desired. The sources of flux on pixel @mymath{i}@footnote{For this analy=
sis
the dimension of the data (image) is irrelevant. So if the data is an im=
age
+desired. The sources of flux on pixel@footnote{For this analysis the
+dimension of the data (image) is irrelevant. So if the data is an image
(2D) with width of @mymath{w} pixels, then a pixel located on column
@mymath{x} and row @mymath{y} (where all counting starts from zero and (=
0,
0) is located on the bottom left corner of the image), would have an ind=
ex:
@mymath{i=3Dx+y\times{}w}.} of the image can be written as follows:
+@mymath{i=3Dx+y\times{}w}.} @mymath{i} of the image can be written as
+follows:
=20
@itemize
@item
Contribution from the target object, (@mymath{O_i}).
+Contribution from the target object (@mymath{O_i}).
@item
Contribution from other detected objects, (@mymath{D_i}).
+Contribution from other detected objects (@mymath{D_i}).
@item
Undetected objects or the fainter undetected regions of bright
objects, (@mymath{U_i}).
+objects (@mymath{U_i}).
@item
@cindex Cosmic rays
A cosmic ray, (@mymath{C_i}).
+A cosmic ray (@mymath{C_i}).
@item
@cindex Background flux
The background flux, which is defined to be the count if none of the
others exists on that pixel, (@mymath{B_i}).
+others exists on that pixel (@mymath{B_i}).
@end itemize
@noindent
The total flux in this pixel (@mymath{T_i}) can thus be written as:
@@ 14489,9 +14748,9 @@ The total flux in this pixel (@mymath{T_i}) can t=
hus be written as:
@cindex Cosmic ray removal
@noindent
By definition, @mymath{D_i} is detected and it can be assumed that it is
correctly estimated (deblended) and subtracted, thus @mymath{D_i=3D0}. T=
here
are also methods to detect and remove cosmic rays, for example the metho=
d
described in van Dokkum (2001)@footnote{van Dokkum,
+correctly estimated (deblended) and subtracted, we can thus set
+@mymath{D_i=3D0}. There are also methods to detect and remove cosmic ray=
s,
+for example the method described in van Dokkum (2001)@footnote{van Dokku=
m,
P. G. (2001). Publications of the Astronomical Society of the Pacific.
113, 1420.}, or by comparing multiple exposures. This allows us to set
@mymath{C_i=3D0}. Note that in practice, @mymath{D_i} and @mymath{U_i} a=
re
@@ 14501,37 +14760,39 @@ algorithm is perfect. With these limitations in=
mind, the observed Sky
value for this pixel (@mymath{S_i}) can be defined as
=20
@cindex Sky value
@dispmath{S_i=3DB_i+U_i.}
+@dispmath{S_i\equiv{}B_i+U_i.}
=20
@noindent
Therefore, as the detection process (algorithm and input parameters)
becomes more accurate, or @mymath{U_i\to0}, the sky value will tend to
the background value or @mymath{S_i\to B_i}. Therefore, while
+becomes more accurate, or @mymath{U_i\to0}, the Sky value will tend to t=
he
+background value or @mymath{S_i\to B_i}. Hence, we see that while
@mymath{B_i} is an inherent property of the data (pixel in an image),
@mymath{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. With this definition of sky, the
object flux in the data can be calculated with
+@mymath{S_i} depends on the detection process. Over a group of pixels, f=
or
+example in an image or part of an image, this equation translates to the
+average of undetected pixels (Sky@mymath{=3D\sum{S_i}}). With this
+definition of Sky, the object flux in the data can be calculated, per
+pixel, with
=20
@dispmath{ T_{i}=3DS_{i}+O_{i} \quad\rightarrow\quad
O_{i}=3DT_{i}S_{i}.}
=20
@noindent
@cindex photoelectrons
Hence, the more accurately @mymath{S_i} is measured, the more accurately
the brightness (sum of pixel values) of the target object can be measure=
d
(photometry). Any under(over)estimation in the sky will directly
translate to an over(under)estimation of the measured object's
brightness. In the fainter outskirts of an object a very small fraction =
of
the photoelectrons in the pixels actually belong to objects (see Figure=
1b
in @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
(2015)}). Therefore even a small over estimation of the sky value will
+In the fainter outskirts of an object, a very small fraction of the
+photoelectrons in a pixel actually belongs to objects, the rest is caus=
ed
+by random factors (noise), see Figure 1b in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+(2015)}. 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 overestimation of the Sky
value and thus even more underestimation of the object's brightness. It=
is
thus very important to detect the diffuse flux of a target, even if they
are not your primary target.
=20
+In summary, the more accurately the Sky is measured, the more accurately
+the brightness (sum of pixel values) of the target object can be measure=
d
+(photometry). Any under/overestimation in the Sky will directly transla=
te
+to an over/underestimation of the measured object's brightness.
+
@cartouche
@noindent
The @strong{Sky value} is only correctly found when all the detected
@@ 14552,15 +14813,16 @@ signalbased detection is a detection process t=
hat relies heavily on
assumptions about the tobedetected objects. This method was the most
heavily used technique prior to the introduction of NoiseChisel in that
paper.} use the sky value as a reference to define the detection
threshold. So these old techniques had to rely on approximations based o=
n
other assumptions about the data. A review of those other techniques can=
be
seen in Appendix A of Akhlaghi and Ichikawa (2015)@footnote{Akhlaghi M.,
Ichikawa. T. (2015). Astrophysical Journal Supplement Series.}. Since th=
ey
were extensively used in astronomical data analysis for several decades,
such approximations have given rise to a lot of misconceptions, ambiguit=
ies
and disagreements about the sky value and how to measure it. As a summar=
y,
the major methods used until now were an approximation of the mode of th=
e
image pixel distribution and @mymath{\sigma}clipping.
+threshold. These older techniques therefore had to rely on approximation=
s
+based on other assumptions about the data. A review of those other
+techniques can be seen in Appendix A of
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
+
+These methods were extensively used in astronomical data analysis for
+several decades, therefore they have given rise to a lot of misconceptio=
ns,
+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 @mymath{\sigma}clipping.
=20
@itemize
@cindex Histogram
@@ 14573,19 +14835,19 @@ assume (or find) a certain probability density =
function (PDF) or use the
histogram. But astronomical datasets can have any distribution, making i=
t
almost impossible to define a generic function. Also, histogrambased
results are very inaccurate (there is a large dispersion) and it depends=
on
the histogram binwidths.
+the histogram binwidths. Generally, the mode of a distribution also shi=
fts
+as signal is added. Therefore, even if it is accurately measured, the mo=
de
+is a biased measure for the Sky value.
=20
@cindex sigmaclipping
@item
Another approach was to iteratively clip the brightest pixels in the
image (which is known as @mymath{\sigma}clipping, since the reference
was found from the image mean and its standard deviation or
@mymath{\sigma}). See @ref{Sigma clipping} for a complete
explanation. The problem with @mymath{\sigma}clipping was that real
astronomical objects have diffuse and faint wings that penetrate
deeply into the noise. So only removing their brightest parts is
completely useless in removing the systematic bias an object's fainter
parts cause in the sky value.
+Another approach was to iteratively clip the brightest pixels in the ima=
ge
+(which is known as @mymath{\sigma}clipping). See @ref{Sigma clipping} f=
or
+a complete explanation. @mymath{\sigma}clipping is useful when there ar=
e
+clear outliers (an object with a sharp edge in an image for
+example). However, real astronomical objects have diffuse and faint wing=
s
+that penetrate deeply into the noise, see Figure 1 in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
@end itemize
=20
As discussed in @ref{Sky value}, the sky value can only be correctly
@@ 14602,76 +14864,78 @@ are ultimately poor approximations.
@cindex Noise
@cindex Signal
@cindex Gaussian distribution
Put simply, noise can characterized with a certain spread about a
characteristic value. In the Gaussian distribution (most commonly used t=
o
model noise) the spread is defined by the standard deviation about the
characteristic mean. Before continuing let's clarify some definitions
first: @emph{Data} is defined as the combination of signal and noise (so=
a
noisy image is one @emph{data}set). @emph{Signal} is defined as the mea=
n
of the noise on each element (after sky subtraction, see @ref{Sky value
definition}).

Let's assume that the @emph{background} (see @ref{Sky value definition})=
is
subtracted and is zero. When a data set doesn't have any signal (only
noise), the mean, median and mode of the distribution are equal within
statistical errors and approximately equal to the background value. Sig=
nal
always has a positive value and will never become negative, see Figure 1=
in
+Put simply, noise can be characterized with a certain spread about the
+measured value. In the Gaussian distribution (most commonly used to mode=
l
+noise) the spread is defined by the standard deviation about the
+characteristic mean.
+
+Let's start by clarifying some definitions first: @emph{Data} is defined=
as
+the combination of signal and noise (so a noisy image is one
+@emph{data}set). @emph{Signal} is defined as the mean of the noise on ea=
ch
+element. We'll also assume that the @emph{background} (see @ref{Sky valu=
e
+definition}) is subtracted and is zero.
+
+When a data set doesn't have any signal (only noise), the mean, median a=
nd
+mode of the distribution are equal within statistical errors and
+approximately equal to the background value. Signal always has a positi=
ve
+value and will never become negative, see Figure 1 in
@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
(2015)}. Therefore, as more signal is added to the raw noise, the mean,
median and mode of the dataset (which has both signal and noise) shift t=
o
the positive. The mean's shift is the largest. The median shifts less,
since it is defined based on an ordered distribution and so is not affec=
ted
by a small number of outliers. The distribution's mode shifts the least =
to
the positive.
+(2015)}. Therefore, as more signal is added, the mean, median and mode o=
f
+the dataset shift to the positive. The mean's shift is the largest. The
+median shifts less, since it is defined based on an ordered distribution
+and so is not affected by a small number of outliers. The distribution's
+mode shifts the least to the positive.
=20
@cindex Mode
@cindex Median
Inverting the argument above gives us a robust method to quantify the
significance of signal in a dataset. Namely, when the mode and median of=
a
distribution are approximately equal, we can argue that there is no
significant signal. To allow for gradients (which are commonly present i=
n
groundbased images), we can consider the image to be made of a grid of
tiles (see @ref{Tessellation}@footnote{The options to customize the
tessellation are discussed in @ref{Processing options}.}). Hence, from t=
he
difference of the mode and median on each tile, we can `detect' the
significance of signal in it. The median of a distribution is defined to=
be
the value of the distribution's middle point after sorting (or 0.5
quantile). Thus, to estimate the presence of signal, we'll compare with =
the
quantile of the mode with 0.5, if the difference is larger than the valu=
e
given to the @option{modmedqdiff} option, this tile will be ignored. Y=
ou
can read this option as ``modemedianquantilediff''.

This method to use the input's skewness is possible because of a new
algorithm to find the mode of a distribution that was defined in Appendi=
x C
of Akhlaghi and Ichikawa (2015). However, the raw dataset's distribution=
is
noisy (noise also affects the sorting), so using the argument above on t=
he
+significance of signal in a dataset. Namely, when the mean and median of=
a
+distribution are approximately equal, or the mean's quantile is around 0=
.5,
+we can argue that there is no significant signal.
+
+To allow for gradients (which are commonly present in groundbased image=
s),
+we can consider the image to be made of a grid of tiles (see
+@ref{Tessellation}@footnote{The options to customize the tessellation ar=
e
+discussed in @ref{Processing options}.}). Hence, from the difference of =
the
+mean and median on each tile, we can estimate the significance of signal=
in
+it. The median of a distribution is defined to be the value of the
+distribution's middle point after sorting (or 0.5 quantile). Thus, to
+estimate the presence of signal, we'll compare with the quantile of the
+mean with 0.5. If the absolute difference in a tile is larger than the
+value given to the @option{meanmedqdiff} option, that tile will be
+ignored. You can read this option as ``meanmedianquantiledifference''=
.
+
+The raw dataset's distribution is noisy, so using the argument above on =
the
raw input will give a noisy result. To decrease the noise/error in
estimating the mode, we will use convolution (see @ref{Convolution
process}). Convolution decreases the range of the dataset and enhances i=
ts
skewness, See Section 3.1.1 and Figure 4 in Akhlaghi and Ichikawa
(2015). 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 mesh, the
image can be convolved with a given kernel first.

Note that through the difference of the mode and median we have actually
`detected' data in the distribution. However this ``detection'' was only
based on the total distribution of the data in each tile (a much lower
resolution). This is the main limitation of this technique. The best
approach is thus to do detection over the dataset, mask all the detected
pixels and use the undetected regions to estimate the sky and its standa=
rd
deviation.
+skewness, See Section 3.1.1 and Figure 4 in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}. Th=
is
+enhanced skewness can be interpreted as an increase in the Signal to noi=
se
+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.
+
+Through the difference of the mean and median we have actually `detected=
'
+data in the distribution. However this ``detection'' was only based on t=
he
+total distribution of the data in each tile (a much lower resolution). T=
his
+is the main limitation of this technique. The best approach is thus to d=
o
+detection over the dataset, mask all the detected pixels and use the
+undetected regions to estimate the sky and its standard deviation (possi=
bly
+over a tessellation). This is how NoiseChisel works: it uses the argumen=
t
+above to find tiles that are used to find its thresholds. Several
+higherlevel steps are done on the thresholded pixels to define the
+higherlevel detections (see @ref{NoiseChisel}).
=20
@cindex Cosmic rays
The mean value of the tiles that have an approximately equal mode and
median will be the Sky value. However there is one final hurdle:
astronomical datasets are commonly plagued with Cosmic rays. Images of
Cosmic rays aren't smoothed by the atmosphere or telescope aperture, so
they have sharp boundaries. Also, since they don't occupy too many pixel=
s,
they don't 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), see Figure 15 in Akhlag=
hi
and Ichikawa (2015) for one example.
+There is one final hurdle: raw astronomical datasets are commonly pepper=
ed
+with Cosmic rays. Images of Cosmic rays aren't smoothed by the atmospher=
e
+or telescope aperture, so they have sharp boundaries. Also, since they
+don't occupy too many pixels, they don't 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
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
=20
The effect of outliers like cosmic rays on the mean and standard deviati=
on
can be removed through @mymath{\sigma}clipping, see @ref{Sigma clipping=
}
@@ 14680,8 +14944,56 @@ median are approximately equal in a tile (see @r=
ef{Tessellation}), the
final Sky value and its standard deviation are determined after
@mymath{\sigma}clipping with the @option{sigmaclip} option.
=20
+In the end, some of the tiles will pass the mean and median quantile
+difference test. However, prior to interpolating over the failed tiles,
+another point should be considered: large and extended galaxies, or brig=
ht
+stars, have wings which sink into the noise very gradually. In some case=
s,
+the gradient over these wings can be on scales that is larger than the
+tiles. The meanmedian distance test will pass on such tiles and will ca=
use
+a strong peak in the interpolated tile grid, see @ref{Detecting large
+extended targets}.
+
+The tiles that exist over the wings of large galaxies or bright stars ar=
e
+outliers in the distribution of tiles that passed the meanmedian quanti=
le
+distance test. Therefore, the final step of ``quantifying signal in a
+tile'' is to look at this distribution and remove the
+outliers. @mymath{\sigma}clipping is a good solution for removing a few
+outliers, but the problem with outliers of this kind is that there may b=
e
+many such tiles (depending on the large/bright stars/galaxies in the
+image). Therefore a novel outlier rejection algorithm will be used.
+
+@cindex Outliers
+@cindex Identifying outliers
+To identify the first outlier, we'll use the distribution of distances
+between sorted elements. If there are @mymath{N} successful tiles, for
+every tile, the distance between the adjacent @mymath{N/2} previous
+elements is found, giving a distribution containing @mymath{N/21}
+points. The @mymath{\sigma}clipped median and standard deviation of thi=
s
+distribution is then found (@mymath{\sigma}clipping is configured with
+@option{outliersclip}). Finally, if the distance between the element a=
nd
+its previous element is more than @option{outliersigma} multiples of t=
he
+@mymath{\sigma}clipped standard deviation added with the
+@mymath{\sigma}clipped median, that element is considered an outlier an=
d
+all tiles larger than that value are ignored.
+
+Formally, if we assume there are @mymath{N} elements. They are first
+sorted. Searching for the outlier starts on element @mymath{N/2} (intege=
r
+division). Let's take @mymath{v_i} to be the @mymath{i}th element of th=
e
+sorted input (with no blank values) and @mymath{m} and @mymath{\sigma} a=
s
+the @mymath{\sigma}clipped median and standard deviation from the
+distances of the previous @mymath{N/21} elements (not including
+@mymath{v_i}). If the value given to @option{outliersigma} is displaye=
d
+with @mymath{s}, the @mymath{i}th element is considered as an outlier w=
hen
+the condition below is true.
=20
+@dispmath{{(v_iv_{i1})m\over \sigma}>s}
=20
+Since @mymath{i} begins from the median, the outlier has to be larger th=
an
+the median. You can use the check images (for example @option{checksky=
}
+in the Statistics program or @option{checkqthresh},
+@option{checkdetsky} and @option{checksky} options in NoiseChisel fo=
r
+any of its steps that uses this outlier rejection) to inspect the steps =
and
+see which tiles have been discarded as outliers prior to interpolation.
=20
=20
=20
@@ 15232,15 +15544,14 @@ tile, see @ref{Quantifying signal in a tile}.
Kernel HDU to help in estimating the significance of signal in a tile, s=
ee
@ref{Quantifying signal in a tile}.
=20
@item mirrordist=3DFLT
Maximum distance (as a multiple of error) to estimate the difference
between the input and mirror distributions in finding the mode, see
Appendix C of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichika=
wa
2015}, also see @ref{Quantifying signal in a tile}.

@item modmedqdiff=3DFLT
The maximum acceptable distance between the mode and median, see
@ref{Quantifying signal in a tile}.
+@item meanmedqdiff=3DFLT
+The maximum acceptable distance between the quantiles of the mean and
+median, see @ref{Quantifying signal in a tile}. The initial Sky and its
+standard deviation estimates are measured on tiles where the quantiles o=
f
+their mean and median are less distant than the value given to this
+option. For example @option{meanmedqdiff=3D0.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.
=20
@item sclipparams=3DFLT,FLT
The @mymath{\sigma}clipping parameters, see @ref{Sigma clipping}. This
@@ 15252,6 +15563,23 @@ 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.
=20
+@item outliersclip=3DFLT,FLT
+Sigmaclipping parameters for the outlier rejection of the Sky value
+(similar to @option{sclipparams}).
+
+Outlier rejection is useful when the dataset contains a large and diffus=
e
+(almost flat within each tile) signal. The flatness of the profile will
+cause it to successfully pass the meanmedian quantile difference test, =
so
+we'll need to use the distribution of successful tiles for removing thes=
e
+false positive. For more, see the latter half of @ref{Quantifying signal=
in
+a tile}.
+
+@item outliersigma=3DFLT
+Multiple of sigma to define an outlier in the Sky value estimation. If t=
his
+option is given a value of zero, no outlier rejection will take place. F=
or
+more see @option{outliersclip} and the latter half of @ref{Quantifying
+signal in a tile}.
+
@item smoothwidth=3DINT
Width of a flat kernel to convolve the interpolated tile values. Tile
interpolation is done using the median of the @option{interpnumngb}
@@ 15492,11 +15820,11 @@ calculated on the dataset that is convolved wit=
h the wider kernel, then the
quantiles are estimated on the image convolved with the sharper kernel.
=20
@item
@option{noerodequant}: to specify a quantile threshold where erosion
will not apply. This is useful to detect sharper pointlike sources that
will be missed due to too much erosion. To completely ignore this featur=
es
give this option a value of 1 (only the largest valued pixel in the inpu=
t
will not be eroded).
+The quantile difference to identify tiles with no significant signal is
+measured between the @emph{mean} and median. In the published paper, it =
was
+between the @emph{mode} and median. The quantile of the mean is more
+sensitive to skewness (the presence of signal), so it is preferable to t=
he
+quantile of the mode. For more see @ref{Quantifying signal in a tile}.
=20
@item
Outlier rejection in quantile thresholds: When there are large galaxies =
or
@@ 15504,9 +15832,9 @@ bright stars in the image, their gradient may be =
on a smaller scale than
the selected tile size. In such cases, those tiles will be identified as
tiles with no signal and thus preserved. An outlier identification
algorithm has been added to NoiseChisel and can be configured with the
following three options: @option{qthreshoutnum},
@option{qthreshoutsigma} and @option{qthreshoutsclip}. See their
description in @ref{Detection options} for more.
+following options: @option{outliersigma} and @option{outliersclip}. =
For
+a more complete description, see the latter half of @ref{Quantifying sig=
nal
+in a tile}.
=20
@item
@option{detgrowquant}: is used to grow the final true detections until=
a
@@ 15797,26 +16125,52 @@ that is discussed in that section.
@subsubsection Detection options
=20
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
+groups: 1) Signal, and 2) Noise. Through the parameters below, you can
customize the detection process in NoiseChisel. Recall that you can alwa=
ys
see the full list of Gnuastro's options with the @option{help} (see
+see the full list of NoiseChisel's options with the @option{help} (see
@ref{Getting help}), or @option{printparams} (or @option{P}) to see
their values (see @ref{Operating mode options}).
=20
@table @option
=20
@item r FLT
@itemx mirrordist=3DFLT
Maximum distance (as a multiple of error) to estimate the difference
between the input and mirror distributions in finding the mode, see
Appendix C of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichika=
wa
2015}, also see @ref{Quantifying signal in a tile}.

@item Q FLT
@itemx modmedqdiff=3DFLT
The maximum acceptable distance between the mode and median, see
@ref{Quantifying signal in a tile}. The quantile threshold will be found=
on
tiles that satisfy this mode and median difference.
+@itemx meanmedqdiff=3DFLT
+The maximum acceptable distance between the quantiles of the mean and
+median in each tile, see @ref{Quantifying signal in a tile}. The quantil=
e
+threshold estimates are measured on tiles where the quantiles of their m=
ean
+and median are less distant than the value given to this option. For
+example @option{meanmedqdiff=3D0.01} means that only tiles where the m=
ean's
+quantile is between 0.49 and 0.51 (recall that the median's quantile is
+0.5) will be used.
+
+@item outliersclip=3DFLT,FLT
+Sigmaclipping parameters for the outlier rejection of the quantile
+threshold. The format of the given values is similar to
+@option{sigmaclip} below. In NoiseChisel, outlier rejection is used in
+three stages:
+@itemize
+@item
+Identifying the quantile thresholds (@option{qthresh},
+@option{noerodequant}, and @option{detgrowquant}).
+@item
+Identifying the first estimate of the Sky and its standard deviation val=
ues
+for pseudodetections (@option{dthresh}).
+@item
+Identifying the final estimate of the Sky and its standard deviation.
+@end itemize
+
+Outlier rejection is useful when the dataset contains a large and diffus=
e
+(almost flat within each tile) signal. The flatness of the profile will
+cause it to successfully pass the meanmedian quantile difference test, =
so
+we'll need to use the distribution of successful tiles for removing thes=
e
+false positive. For more, see the latter half of @ref{Quantifying signal=
in
+a tile}.
+
+@item outliersigma=3DFLT
+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
+@option{outliersclip} and the latter half of @ref{Quantifying signal i=
n a
+tile}.
=20
@item t FLT
@itemx qthresh=3DFLT
@@ 15840,57 +16194,6 @@ 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=
).
=20
@item qthreshoutnum=3DINT
Number of tiles (in sorted array) to use as a moving window when estimat=
ing
the outlier tiles. If it is given a value of zero, no outlier rejection
will take place.

This option is useful when the dataset contains a large and diffuse (alm=
ost
flat within each tile) signal. The flatness of the profile will cause it=
to
successfully pass the tests of @ref{Quantifying signal in a tile}. As a
result, when outlier rejection is disabled, the flat and diffuse signal
will be interpreted as sky. As a result, the quantile threshold in those
tiles (and the ones around them) can be strongly overestimated, causing=
a
failure in their detection.

Outlier identification proceeds as follows: all nonblank tile values ar=
e
sorted by flux. A sliding window, with width equal to the value given to
this option, is parsed over the dataset and is used as a reference to
identify the first outlier element. The first element where the distance=
to
the previous (sorted) element is sigma units (@option{qthreshoutsigma}=
)
away from the distribution of distances of the element in the window is
considered an outlier and all tiles larger than that value are ignored f=
or
the later step (interpolation).

Formally, Let's take @mymath{v_i} to be the @mymath{i}th element of the
sorted input (with no blank values) and @mymath{m} and @mymath{\sigma} a=
s
the @mymath{\sigma}clipped median and standard deviation from the
distances of the previous @mymath{N} elements (@mymath{N} is the value t=
o
this option). The @option{qthreshoutsclip} option can be used to
configure the @mymath{\sigma}clipping. If the value given to
@option{qthreshoutsigma} is displayed with @mymath{s}, the returned
outlier is the first (sorted) elemented where the following condition
becomes true:

@dispmath{(v_i=E2=88=92v_{i=E2=88=921})=E2=88=92m>s\times\sigma}

Note that by this definition, the outlier cannot be any of the first
@mymath{N} elements, since they constitute the window that is used for t=
he
first checked element. It is therefore important that while @mymath{N}
isn't too small, it also isn't larger than (approximately) half of the
initially measured tiles. A warning will be printed if this condition
doesn't hold. You can use @option{checkqthresh} to inspect the steps.

@item qthreshoutsigma=3DFLT
Multiple of sigma to define an outlier, see description of
@option{qthreshoutnum} for more.

@item qthreshoutsclip=3DFLT,FLT
Sigmaclipping parameters for the outlier rejection of the quantile
threshold. The format of the given values is similar to
@option{sigmaclip} below. see description of @option{qthreshoutnum} =
for
a description of the outlier rejection algorithm.

@item smoothwidth=3DINT
Width of flat kernel used to smooth the interpolated quantile thresholds=
,
see @option{qthresh} for more.
@@ 15968,18 +16271,6 @@ them. Once opening is complete, we have @emph{in=
itial} detections.
The structuring element used for opening, see @option{erodengb} for mo=
re
information about a structuring element.
=20
@item s FLT,FLT
@itemx sigmaclip=3DFLT,FLT
The @mymath{\sigma}clipping parameters, see @ref{Sigma clipping}. This
option takes two values which are separated by a comma (@key{,}). Each
value can either be written as a single number or as a fraction of two
numbers (for example @code{3,1/10}). The first value to this option is t=
he
multiple of @mymath{\sigma} that will be clipped (@mymath{\alpha} in tha=
t
section). The second value is the exit criteria. If it is less than 1, t=
hen
it is interpreted as tolerance and if it is larger than one it is assume=
d
to be the fixed number of iterations. Hence, in the latter case the valu=
e
must be an integer.

@item B FLT
@itemx minskyfrac=3DFLT
Minimum fraction (value between 0 and 1) of Sky (undetected) areas in a
@@ 16008,10 +16299,24 @@ to the final quantile threshold, this behavior =
can be disabled with
pixel size as the input, but with the @option{oneelempertile} option,
only one pixel will be used for each tile (see @ref{Processing options})=
.
=20
+@item s FLT,FLT
+@itemx sigmaclip=3DFLT,FLT
+The @mymath{\sigma}clipping parameters for measuing the initial and fin=
al
+Sky values from the undetected pixels, see @ref{Sigma clipping}.
+
+This option takes two values which are separated by a comma (@key{,}). E=
ach
+value can either be written as a single number or as a fraction of two
+numbers (for example @code{3,1/10}). The first value to this option is t=
he
+multiple of @mymath{\sigma} that will be clipped (@mymath{\alpha} in tha=
t
+section). The second value is the exit criteria. If it is less than 1, t=
hen
+it is interpreted as tolerance and if it is larger than one it is assume=
d
+to be the fixed number of iterations. Hence, in the latter case the valu=
e
+must be an integer.
+
@item R FLT
@itemx dthresh=3DFLT
The detection threshold: a multiple of the initial sky standard deviatio=
n
added with the initial sky approximation (which you can inspect with
+The detection threshold: a multiple of the initial Sky standard deviatio=
n
+added with the initial Sky approximation (which you can inspect with
@option{checkdetsky}). This flux threshold is applied to the initially
undetected regions on the unconvolved image. The background pixels that =
are
completely engulfed in a 4connected foreground region are converted to
@@ 22397,8 +22702,7 @@ function uses C's @code{sizeof} operator to measu=
re the size of each type.
@end deftypefun
=20
@deftypefun {char *} gal_type_name (uint8_t @code{type}, int @code{long_=
name})
Return a string that contains the name of @code{type}. This can be used =
in
messages to the users when your function/program accepts many types. It =
can
+Return a string literal that contains the name of @code{type}. It can
return both short and long formats of the type names (for example
@code{f32} and @code{float32}). If @code{long_name} is nonzero, the lon=
g
format will be returned, otherwise the short name will be returned. The
@@ 22728,6 +23032,14 @@ that corresponds to its type. If @code{input} is=
a tile over a larger
dataset, only the region that the tile covers will be set to blank.
@end deftypefun
=20
+@deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int @code=
{width})
+Write the blank value for the given data type @code{type} into a string =
and
+return it. The space for the string is dynamically allocated so it must =
be
+freed after you are done with it. If @code{width!=3D0}, then the final s=
tring
+will be padded with white space characters to have the requested width i=
f
+it is smaller.
+@end deftypefun
+
@deftypefun int gal_blank_is (void @code{*pointer}, uint8_t @code{type})
Return 1 if the contents of @code{pointer} (assuming a type of @code{typ=
e})
is blank. Otherwise, return 0. Note that this function only works on one
@@ 22771,6 +23083,12 @@ Create a dataset of the the same size as the inp=
ut, but with an
those that aren't.
@end deftypefun
=20
+@deftypefun void gal_blank_flag_apply (gal_data_t @code{*input}, gal_dat=
a_t @code{*flag})
+Set all nonzero and nonblank elements of @code{flag} to blank in
+@code{input}. @code{flag} has to have an unsigned 8bit type and be the
+same size as @code{input}.
+@end deftypefun
+
=20
@deftypefun void gal_blank_remove (gal_data_t @code{*input})
Remove blank elements from a dataset, convert it to a 1D dataset, adjust
@@ 22787,11 +23105,10 @@ appropriate measure. This check is highly recom=
mended because it will avoid
strange bugs in later steps.
@end deftypefun
=20
@deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int @code=
{width})
Write the blank value for the given data type @code{type} into a string =
and
return it. The space for the string is dynamically allocated so it must =
be
freed after you are done with it.
@end deftypefun
+
+
+
+
=20
@node Library data container, Dimensions, Library blank values, Gnuastro=
library
@subsection Data container (@file{data.h})
@@ 27209,36 +27526,40 @@ above.
@end deftypefun
=20
=20
@deftypefun {gal_data_t *} gal_statistics_outlier_positive (gal_data_t @=
code{*input}, size_t @code{window_size}, float @code{sigma}, float @code{=
sigclip_multip}, float @code{sigclip_param}, int @code{inplace}, int @cod=
e{quiet})
+@deftypefun {gal_data_t *} gal_statistics_outlier_positive (gal_data_t @=
code{*input}, float @code{sigma}, float @code{sigclip_multip}, float @cod=
e{sigclip_param}, int @code{inplace}, int @code{quiet})
Find the first positive outlier in the @code{input} distribution. The
returned dataset contains a single element: the first positive outlier. =
It
is one of the dataset's elements, in the same type as the input. If the
process fails for any reason (for example @code{window_size} is larger t=
han
the number of @code{input}'s nonblank elements, or no outlier was found=
),
a @code{NULL} pointer will be returned.
+process fails for any reason (for example no outlier was found), a
+@code{NULL} pointer will be returned.
=20
All (possibly existing) blank elements are first removed from the input
dataset, then it is sorted. A sliding window of @code{window_size} eleme=
nts
is parsed over the dataset and is used as a reference to identify the
outlier. The first element where the distance to the previous (sorted)
element is @code{sigma} units away from the distribution of distances of
the @code{window_size} previous elements is considered an outlier and
returned by this function.

Formally, Let's take @mymath{v_i} to be the @mymath{i}th element of the
sorted input (with no blank values) and @mymath{m} and @mymath{\sigma} a=
s
the @mymath{\sigma}clipped median and standard deviation from the
distances of the previous @code{window_size} elements. If the value give=
n
to @code{sigma} is displayed with @mymath{s}, the returned outlier is th=
e
first (sorted) elemented where the following condition becomes true:
+dataset, then it is sorted. A sliding window equal to half the width of =
the
+dataset elements is parsed over the dataset. Starting from the middle of
+the dataset (median) in the direction of increasing values. This window =
is
+used as a reference. The first element where the distance to the previou=
s
+(sorted) element is @code{sigma} units away from the distribution of
+distances in its window is considered an outlier and returned by this
+function.
+
+Formally, if we assume there are @mymath{N} nonblank elements. They are
+first sorted. Searching for the outlier starts on element @mymath{N/2}
+(integer division). Let's take @mymath{v_i} to be the @mymath{i}th elem=
ent
+of the sorted input (with no blank values) and @mymath{m} and
+@mymath{\sigma} as the @mymath{\sigma}clipped median and standard
+deviation from the distances of the previous @mymath{N/2} elements (not
+including @mymath{v_i}). If the value given to @code{sigma} is displayed
+with @mymath{s}, the @mymath{i}th element is considered as an outlier w=
hen
+the condition below is true.
=20
@dispmath{{(v_iv_{i1})m\over \sigma}>s}
=20
The @code{sigclip_multip} and @code{sigclip_param} arguments specify the
properties of the @mymath{\sigma}clipping. You see that by this
definition, the outlier cannot be any of the first @code{window_size}
elements, since they constitute the window that is used for the first
checked element.
+properties of the @mymath{\sigma}clipping (see @ref{Sigma clipping} for
+more). You see that by this definition, the outlier cannot be any of the
+lower half elements. The advantage of this algorithm compared to
+@mymath{\sigma}clippign is that it only looks backwards (in the sorted
+array) and parses it in one direction.
=20
If @code{inplace!=3D0}, the removing of blank elements and sorting will =
be
done within the input dataset's allocated space. Otherwise, this functio=
n
diff git a/lib/Makefile.am b/lib/Makefile.am
index a6ee75a..8a36723 100644
 a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ 60,8 +60,8 @@ libgnuastro_la_SOURCES =3D arithmetic.c arithmeticand.=
c arithmeticbitand.c \
checkset.c convolve.c cosmology.c data.c eps.c fits.c git.c =
\
interpolate.c jpeg.c label.c list.c match.c options.c pdf.c =
\
permutation.c pointer.c polygon.c qsort.c dimension.c statistics.c =
\
 table.c tableintern.c threads.c tiff.c tile.c timing.c txt.c type.c =
\
 wcs.c
+ table.c tableintern.c threads.c tiff.c tile.c tileinternal.c timing.c=
\
+ txt.c type.c wcs.c
=20
=20
=20
@@ 107,7 +107,7 @@ EXTRA_DIST =3D gnuastro.pc.in $(headersdir)/README $(=
internaldir)/README \
$(internaldir)/checkset.h $(internaldir)/commonopts.h =
\
$(internaldir)/config.h.in $(internaldir)/fixedstringmacros.h =
\
$(internaldir)/options.h $(internaldir)/tableintern.h =
\
 $(internaldir)/timing.h
+ $(internaldir)/tileinternal.h $(internaldir)/timing.h
=20
=20
=20
diff git a/lib/blank.c b/lib/blank.c
index b84a5ce..b3c2bb7 100644
 a/lib/blank.c
+++ b/lib/blank.c
@@ 113,6 +113,184 @@ gal_blank_initialize(gal_data_t *input)
=20
=20
=20
+/* Print the blank value as a string. For the integer types, we'll use t=
he
+ PRIxNN keywords of `inttypes.h' (which is imported into Gnuastro from
+ Gnulib, so we don't necessarily rely on the host system having it). *=
/
+char *
+gal_blank_as_string(uint8_t type, int width)
+{
+ char *blank=3DNULL, *fmt;
+
+ /* Print the given value. */
+ switch(type)
+ {
+ case GAL_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "%s: bit types are not implemented yet",
+ __func__);
+ break;
+
+ case GAL_TYPE_STRING:
+ if(width)
+ {
+ if( asprintf(&blank, "%*s", width, GAL_BLANK_STRING)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, "%s", GAL_BLANK_STRING)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT8:
+ fmt =3D width ? "%*"PRIu8 : "%"PRIu8;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint8_t)GAL_BLANK_UINT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint8_t)GAL_BLANK_UINT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT8:
+ fmt =3D width ? "%*"PRId8 : "%"PRId8;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int8_t)GAL_BLANK_INT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int8_t)GAL_BLANK_INT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT16:
+ fmt =3D width ? "%*"PRIu16 : "%"PRIu16;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint16_t)GAL_BLANK_UINT16)<0=
)
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint16_t)GAL_BLANK_UINT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT16:
+ fmt =3D width ? "%*"PRId16 : "%"PRId16;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int16_t)GAL_BLANK_INT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int16_t)GAL_BLANK_INT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT32:
+ fmt =3D width ? "%*"PRIu32 : "%"PRIu32;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint32_t)GAL_BLANK_UINT32)<0=
)
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint32_t)GAL_BLANK_UINT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT32:
+ fmt =3D width ? "%*"PRId32 : "%"PRId32;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int32_t)GAL_BLANK_INT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int32_t)GAL_BLANK_INT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT64:
+ fmt =3D width ? "%*"PRIu64 : "%"PRIu64;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint64_t)GAL_BLANK_UINT64)<0=
)
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint64_t)GAL_BLANK_UINT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT64:
+ fmt =3D width ? "%*"PRId64 : "%"PRId64;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int64_t)GAL_BLANK_INT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int64_t)GAL_BLANK_INT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_FLOAT32:
+ if(width)
+ {
+ if( asprintf(&blank, "%*f", width, (float)GAL_BLANK_FLOAT32)<=
0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, "%f", (float)GAL_BLANK_FLOAT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_FLOAT64:
+ if(width)
+ {
+ if( asprintf(&blank, "%*f", width, (double)GAL_BLANK_FLOAT64)=
<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, "%f", (double)GAL_BLANK_FLOAT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized", __func__=
,
+ type);
+ }
+ return blank;
+}
+
+
+
+
+
=20
/* Return 1 if the contents of the pointer (with the given type) is
blank. */
@@ 318,7 +496,7 @@ gal_blank_flag(gal_data_t *input)
/* Go over the pixels and set the output values. */
switch(input>type)
{
 /* Numeric types */
+ /* Numeric types */
case GAL_TYPE_UINT8: FLAG_BLANK( uint8_t ); break;
case GAL_TYPE_INT8: FLAG_BLANK( int8_t ); break;
case GAL_TYPE_UINT16: FLAG_BLANK( uint16_t ); break;
@@ 330,19 +508,19 @@ gal_blank_flag(gal_data_t *input)
case GAL_TYPE_FLOAT32: FLAG_BLANK( float ); break;
case GAL_TYPE_FLOAT64: FLAG_BLANK( double ); break;
=20
 /* String. */
+ /* String. */
case GAL_TYPE_STRING:
do *o++ =3D !strcmp(*str,GAL_BLANK_STRING); while(++str