Fixes taks#13767.  Relying on CFITSIO, Gnuastro was already able to read
files compressed with fpack, it was just that it did not accept the
standard extension.

Fixes bug#45380: if optimization flag -O2 or higher is used, a bus error
is signaled when trying to access "p->cp.nolog" after calling the
function pointed to by modefunction.  Not using modefunction and calling
the appropriate function "by hand" avoids this crash.

Fixes bug#46212.

Until now the maximum standard deviation in the meshes was used as the
depth. However, the maximum has a strong scatter such that one outling
mesh can completely shift to unreasonable values. So, although the
maximum standrard deviation is still reported by NoiseChisel, but
Makecatalog uses the median standard deviation.

I also added a median function to `statistics.[ch]'.

A tiny spelling mistake was also fixed in MakeCatalog's `--help'
output.

Since ImageWarp is also a starting point in reading images with
non-standard WCS information, I also added this option to
ImageWarp. If it is needed in all applications, these options can be
added to the common arguments easily now. The main readfitswcs
function in lib/fitsarrayvv.[ch] was changed to allow these two values
for easy generalization if needed.

Previously the surface brightness limit and magnitude limit were
defined in terms of a fixed projected area (arcsec^2) and the S/N of
the objects or clumps. But in the first case, the projected area
should not be used since noise is only dependent on the pixel and the
pixel is the unit of our data-collection. Its physical or projected
size are irrelevant. In the latter case, the area used to find the S/N
was not clear and thus this value was very unrealistic. Therefore I
tried to explain how to use number counts (from the outputs of
MakeCatalog) in order to estimate the limiting magnitude.

When the scripts are run in parallel (for example within a Makefile
called with the -j option) then the log files can cause problems when
you want to clean them up in the end. One option was to set the log
file name to be based on the output file name, but I have not seen
that commonly and it would require a lot of extra code. But now, when
the log file doesn't matter, the user can simply igore it being made.

Convolve was not reading the WCS information of the input. A line is
now added in the ui.c to do the job.

The individual outputs of MakeProfiles were simply named as, for
example, 0.fits, 1.fits and so on. But this would not work when
MakeProfiles is to be run in parallel and is also ambiguous. So now
the row numbers are suffixed by the output file name and put into the
output directory.

The output directory and file name fixing in MakeProfiles was also
bugy with problems I hadn't noticed before, so I fixed it.

In order to do this job, two functions were added in checkset.c:
`dirpart' and `notdirpart'.

A new subsection was also added to the manual explaining the two basic
methods to run processes simultaneously.

It was necessary that the branch for the NoiseChisel paper
reproduction (version 0.0) would be separate from all the work that
was done on the programs after the last proof submission of the
paper.

After the paper was published, I completed the reproduction system and
in the process, some issues needed to be added or fixed in the
programs. They were all applied to that branch. But they are all also
necessary for the future versions of Gnuastro, so they were done
independently here.

The --makekernel option in Convlve, now gets an argument that is the
maximum radius to keep non-zero kernel values. The output kernel will
also be cropped to only be large enough to have no only-zero rows or
columns.

This was done because the high frequencies were set to zero during the
division and so the larger radii would become too flat. So the user
can specify how far out the final kernel should be.

Also a minor bug fix in NoiseChisel: The check to see if the number of
clumps or detections is more than the minimum specified by the user is
now done before anything else. In the case of zero good S/N values,
this would cause a problem in the next steps.

Convolve now has a `--makekernel' option which will find the kernel
which can be used to match the PSFs of two images.

Two new sections were added to explain the Discrete Fourier transform
and how they change in 2 dimensions. I also went through all the other
explanations and made the whole frequency space explanations more
cogent.

There were two bugs that are now fixed:

  - MakeCatalog: In a grown clump label image, the detections that
    have one or no clumps are fully filled. So there is no river
    pixels for them. I had not corrected the clump S/N calculations in
    MakeCatalog also consider this. Now this issue is corrected.

  - MakeCatalog: Although brightness was correctly calculated, in the
    final printing, it would be divided by the area! It is now fixed.

  - NoiseChisel: When finding the S/N quantile for clumps, it was
    mistakenly using the detection S/N quantile for both clumps and
    detections! Now it uses the correct value. Until now, these two
    were always the same, so I hadn't noticed this!

With this new option, the user can ask to save the grown clumps image
instead of the original clumps image in the output.

Additions
=========

- ImageCrop: Added the two options `--hstartwcs' and `--hendwcs' to
  read specific portions of the input file's header. As explained in
  the manual this is necessary now that WCSLIB includes distortion in
  the coordinates but old FITS files were written with the assumption
  that WCSLIB does not support them.

Bug fixes
=========

- NoiseChisel: The flux weighted center of the detections (which are
  used to find the sky standard deviation on the mesh grid) used the
  flux of all the pixels. But a weight can only have a positive value!
  So the negative values would make the result go wrong. Now, the xys
  array has a third column to keep the sum of positive flux pixels
  which were used as weights in finding the flux weighted center.

To get the river information of each clump, I had mistakenly just
copied from the getclumpinfo function in NoiseChisel's
clumps.c. However, over there all the clumps were within one
detection. So the IDs of neighboring clumps (to a river pixel) were
unique. But in MakeCatalog, two clumps can be separated by one river
pixel and belong to separate objects with different object IDs.

Another consequence was to use the object ID from the neighbor of the
river pixel and not the object on which the river lies.

The total brightness of the clump is also stored internally. Only when
we want to write it to the catalog, do we subtract the average river
flux multiplied by the number of clump pixels.

In the fitsarrayvv.c functions, I had mistakenly used a variable
called `numblank'. Some time in the past, this variable probably
actually did mean the number of blanks. But now it was only a boolean
value specifying if there are any blanks in the input or not. So I
renamed all wrong occurances of `numblank' (as relates to the input
image and FITS blank pixels) to `anyblank' which correctly represents
what this variable holds.

Note that NoiseChisel's segmentation does make use of a similarly
termed variable for the number of pixels left to grow. That is a
different usage and those variable names were left unchanged.

Until now, in labeling the dilated image, I used 4-connected
components. However, in the segmentation step, the 8-connected
neighbors of each river pixel are checked. When the rivers of the
edges of the detections are checked, these two different neighbor
definitions would cause a problem: two separate 4-connected regions
can be 8-connected.

The fix of this bug was a simple change of `4' to `8' in the main
NoiseChisel function. However it took me several (very frustrating)
hours to find it! While Valgrind's messages are usually very useful,
in this case (because they resulted from my own convention), they
caused even more confusion. The valgrind error said something like
this: "You are trying to reference an internal array of another
function that was allocated and freed in another funciton". In short
the array it was complaining about didn't exist when the function
giving the error was called!!! This really confused me! I finally
found it through individual checks of all the detections without using
Valgrind.

In the meantime, I also removed all the checks for NaN (masked) pixels
in oversegment (clumps.c). Since NaN pixels are no longer in the
indexs that are fed to this function.

Throughout Gnuastro, until now, I had used flux and brightness
interchangeably while this is only true for one pixel. So for example
using terms like "total flux" for the sum of pixels is wrong. I have
explained everything completely in the "Flux Brightness and magnitude"
subsection of the manual. All the programs (mainly variable names and
option explanations) and the manual were corrected to make the proper
distinction between these two terms.

Also the other MathJax scripts that might be necessary for the HTML
webpages to be run with LibreJS were also added. However, LibreJS
still doesn't recognize the scripts. I have asked
help-librejs@gnu.org.

The detection limit of the image for clumps and objects is now
reported in the output of MakeCatalog. In order to do this, I had to
save the S/N thresholds in NoiseChisel to header keywords in the
NoiseChisel outputs. MakeCatalog then reads those headers and makes
the information. The manual also contains a complete explanation.

I also had to update the old "readkeyword" function in fitsarrayvv.c
to "readkeywords". Now it can read as many keywords as the user wants.

Until now, MakeCatalog was not prepared to handle blank pixels in the
labeled images. But now, the labeled pixels are easily ignored.

I also added the names of the input files in the catalogs and also a
column specifying the 5 sigma magnitude in the input image.

NoiseChisel underwent to large changes. Mainly as a result of the
previous addition to ImageCrop where an arbitary polygon could be
cropped. So I had to find a way to treat large blank areas in a the
image.

  1. The S/N value was previously found on each large mesh separately
     and would be interpolated and smoothed! On each mesh there
     usually isn't enough area to provide the right number of
     points. Unlike the sky and its standard deviation (noise) which
     can vary over the image, the Signal to noise ratio will only
     change when the data have different noise properties, for example
     more or less correlated noise. This can happen in an image,
     especially in surveys but it doesn't have to be on the same grid,
     so there will be problems.

     Now that ImageCrop can crop out an arbitary polygon, we can crop
     out the parts that have similar noise properties and work on
     them. Over the image the number statistics will be much better
     (for example in estimating the quantile). So NoiseChisel will now
     calculate the detection and segmentation Signal to noise ratios
     over the full image, not within one mesh.

  2. Since I wanted to work on a random polygon a large area of the
     image would be blank. Until now, I had not really put too much
     emphasis on blank pixels, but since their number significantly
     increased I had to find a solution. Now NoiseChisel's functions
     also account for unsigned char and long type blank values too.

I had forgot to confirm if the old `--inpolygon' option actually
worked! So I tried it out and saw that it is much better if it is
`--outpolygon'! Since by default only the region inner to the polygon
is shown, so the old method was mostly redundant. In the few cases
that the user would want to keep the outer parts, they would have to
go and change all their configuration files! Now this option is not in
configuration files any more.

A test was also written for this option.

The polygon mode was tested in the WCS mode and some corrections
needed to be made for the steps I wrote last night. It now works very
nicely both in image mode and WCS mode.

Two tests were also added to accompany this new feature.

MakeCatalog's ui.c is complete for now (all the inputs are read). In
order to make things clear and elegant, I made a clear function to
check for data in other extensions of an image or in another file
depending on the user input.

Until now, there was a macro in checkset.h and also a function in
fitsarrayvv.c which done the job almost similarly. The mask checking
functions of all the programs were updated.

Until now, for each string option (for example `--hdu'), I would check
and allocate the space, individually. Also when printing, I had to
check individually to see if it has a space character or not! So I
wrote a simple function in checkset.c and a macro in checkset.h to do
the job and now for all the options in all the programs, we can simply
call those functions. This significantly increases the readability and
decreases the potential for bugs in the code.

Until now there was only a small number of cases in each each program,
so it wasn't bothering me. But while making MakeProfiles (which has
lots of such functions), I noticed how stupidly repetative the checks
are! So I wrote this function and the macro. When I saw how convenient
they are, I applied them to all the programs, to enhance them too.

The fast writing of the NoiseChisel options was very buggy, so I went
through it again and made corrections.

While doing the corrections, I noticed an improvement: I could create
a function to write the mesh values to a FITS file instead of
separately writing checks everytime I wanted to do so. This idea
actually came to my mind when I wanted to add the --meshbasedcheck
option. With this functionality added, several checks would have to be
done multiple times in the course of one function which would make
things very hard to read for an outside reader. They were all also
repetative and could have caused bugs! These parts of the functions
are now much more cleaner and elegant than before.

All the options that are particular to NoiseChisel are now explained
in the manual. Some of the options also underwent small corrections:

  - --minbfrac and --minnumfalse were moved under the Input category
    from the detection category in the output of --help. This was
    because these two options are used by the detection as well as
    segmentation steps.

  - --checkthresh was renamed to --checkthreshold. The user can even
    use --checkt to get the functionality. It is just more clear like
    this.

  - --numerosion was changed to --erode. This makes it much more
    easier and comparable to the --opening and --dilate options.

One test was added for NoiseChisel. I also modified the default
background and zeropoint magnitudes in MakeNoise to fit with the
default zeropoint magnitude in MakeProfiles.

The false clump removal process is now implemented in NoiseChisel. In
order to do it on a real image (with NaN values), I had to make some
modifications to the oversegmentation function. The list of top indexs
for each clump would not be accurate when there were NaN pixels in the
clump. So by setting some standards for that list and some checks in
the oversegmentation, it now stores the top non-NaN pixel of each
clump.

Work is progressing in the segmentation of detections in NoiseChisel,
the detections are now over-segmented and their signal to noise is
calculated. I am now busy bringing in the removal of false clumps
function to find the true clumps.

While working on segmentation, I noticed that there is a problem in
the imgindextomeshid function. This function would get an image index
and output the ID of the mesh that pixel belongs to. The problem was
this: when the mesh size is not an exact multiple of the image size
(in any dimension), the last mesh on all dimensions will have a
different size. If it is larger, then before this correction, this
function would give one larger mesh index for the respective
dimension, which would result in reading values outside of the
garrays. Now there is a nice check.

Until now, for each step of detection and segmentation, it was
necessary to make a temporary sky subtracted image. When the input
image is large, this will result in a large chunk of memory being
allocated and deallocated. It would also make the code a little
confusing.

Now a new member has been put in the noisechiselparams structure to
keep an array the size of the input image. This array will be used by
both the detection and segmentation processes and is allocated in ui.c
before any actual processing. Also thanks to this array (with a
different name) it is much more clear which one of the input image or
sky subtracted image we are dealing with.

Note that in particular for segmentation, using the sky subtracted
image to measure the total flux would make us add more noise (due to
the subtraction of the sky). So we use the actual input image.

There were serveral corrections made in this commit:

 1. In MakeProfiles and MakeNoise that used the random number
    generator, everytime a random number was needed, a generator would
    be allocated and a number pulled out of it. Now I have put the
    main allocation of the random number generator inside ui.c and
    then the functions that need it either directly use it (MakeNoise)
    or get a clone of it for each thread and for each profile it is
    memcpy'd into the clone. This makes things much more cleaner for
    the human reader and the computer!

 2. In MakeNoise I also removed the backgroundinmean option. Since it
    just made things more complicated and I couldn't see how it would
    be useful. Adding noise without adding the background (the default
    behavior prior to this commit) was just useless.

 3. A bug in Convolve was fixed. The floating point error results
    after convolution were not correctly set to zero. Now it is fixed.

There were two main updates in this commit:

1. The mesh IDs (that are sometimes on the whole image and sometimes
   on channels) were very confusing for me and so the functions were
   ambiguous and thus caused a lot of mis-understanding. After playing
   a lot with them and thinking more about them, I got a much better
   understanding of how to clearly separate the use cases. I have
   explained everything extensively in mesh.c so someone new (and me
   in the future!) can easily get going.

2. Thanks to these new functions, there was no more need to store an
   array the size of the image for the standard deviation of every
   pixel. Which ever function needs the standard deviation for pixel
   can easily find it from the garray in the mesh structure. This is
   very good for memory usages, because for large images, a full copy
   was very heavy.

3. The NoiseChisel clump S/N calculations are now patially done. They
   still need some work to do, but things are on the right track.

The over-segmentation algorithm is now in place and the segmentation
process is starting. So far only the clumps in the noise regions can
be found. Now, I have to add the functions to find the signal to noise
of the clumps, then find the S/N limit of the clumps and finally do
the same job on the detections to detect true clumps within them.

The garrays in the mesh structure play a very important role: They
keep the values to be used on each mesh. However, the method I was
dealing with them before was very troublesome. Because in some cases,
it was necessary to use the full image and in other cases to use only
meshs on a channel.

All the functions in mesh.c were corrected to accurately account for
such issues. Now the caller doesn't have to know what the user asked
for (to use the meshes over a grid or in channels). All this is done
internally in mesh.c. The caller only has to use garray1 or garray2.

NoiseChisel's detection algorithm is nearly fully in place now. Both
the background regions and the detection regions are now processed and
the Signal to noise ratio is found for all their connected
components. The Signal to noise ratio limit is also found and smoothed
for all meshs. We just have to remove all the false detections in the
detected regions.

The Histogram and Cumulative frequency functions in statistics were
also corrected to work without any bugs. Until now, there was a
problem when a data element would fall on the last interval.

Detection is nearly complete. Currently it goes up to the point to
calculate the Signal to noise ratio of the noise and detected objects.

The false detection removal threshold is now applied to the image. We
just have to do the filling-holes, opening and finding the signal to
noise steps to find and remove false detections.