gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 9e260a0: All libraries documented


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 9e260a0: All libraries documented
Date: Thu, 4 May 2017 14:14:07 -0400 (EDT)

branch: master
commit 9e260a02486d5a384627a3d892def191098f73c0
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    All libraries documented
    
    All the libraries are now documented. In the process some other small
    corrections were made in some remaining function APIs to make them more
    consistent. Since their corrected version is now documented for the first
    time with this commit, there is no need to go over them. But from this
    commit onwards, any change in the library APIs must be documented in the
    commit message (as well as the manual), so they can be included in the NEWS
    file on every release.
---
 bin/arithmetic/operands.c   |   3 +-
 bin/convertt/ui.c           |   2 +-
 bin/convolve/ui.c           |   2 +-
 bin/crop/ui.c               |   4 +-
 bin/mkcatalog/ui.c          |   3 +-
 bin/mknoise/ui.c            |   3 +-
 bin/mkprof/ui.c             |   3 +-
 bin/noisechisel/detection.c |   8 +-
 bin/noisechisel/ui.c        |   3 +-
 bin/statistics/ui.c         |   4 +-
 bin/warp/ui.c               |   4 +-
 doc/gnuastro.texi           | 556 ++++++++++++++++++++++++++++++++++++++------
 lib/binary.c                |  33 +--
 lib/data.c                  |   2 +-
 lib/gnuastro/wcs.h          |  16 +-
 lib/interpolate.c           |   8 +-
 lib/permutation.c           |   2 +-
 lib/wcs.c                   |  43 ++--
 18 files changed, 544 insertions(+), 155 deletions(-)

diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index e8c7dd0..b4e025e 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -120,8 +120,7 @@ pop_operand(struct imgarithparams *p, char *operator)
       /* In case this is the first image that is read, then read the WCS
          information.*/
       if(p->popcounter==0)
-        gal_wcs_read(filename, hdu, 0, 0, &p->refdata.nwcs,
-                     &p->refdata.wcs);
+        p->refdata.wcs=gal_wcs_read(filename, hdu, 0, 0, &p->refdata.nwcs);
 
       /* Read the input image as a double type array. */
       data=gal_fits_img_read(filename, hdu, p->cp.minmapsize);
diff --git a/bin/convertt/ui.c b/bin/convertt/ui.c
index 053947f..96f9f7e 100644
--- a/bin/convertt/ui.c
+++ b/bin/convertt/ui.c
@@ -433,7 +433,7 @@ ui_make_channels_ll(struct converttparams *p)
 
           /* Read in the array and its WCS information. */
           data=gal_fits_img_read(name->v, hdu, p->cp.minmapsize);
-          gal_wcs_read(name->v, hdu, 0, 0, &data->nwcs, &data->wcs);
+          data->wcs=gal_wcs_read(name->v, hdu, 0, 0, &data->nwcs);
           gal_list_data_add(&p->chll, data);
 
           /* A FITS file only has one channel. */
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 803f385..eaa0c65 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -346,7 +346,7 @@ ui_preparations(struct convolveparams *p)
   /* Read the input image as a float64 array and its WCS info. */
   p->input=gal_fits_img_read_to_type(p->filename, cp->hdu,
                                      GAL_TYPE_FLOAT32, cp->minmapsize);
-  gal_wcs_read(p->filename, cp->hdu, 0, 0, &p->input->nwcs, &p->input->wcs);
+  p->input->wcs=gal_wcs_read(p->filename, cp->hdu, 0, 0, &p->input->nwcs);
 
 
   /* See if there are any blank values. */
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 04ced63..9323e3e 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -678,8 +678,8 @@ ui_preparations(struct cropparams *p)
       img->name=gal_list_str_pop(&p->inputs);
       tmpfits=gal_fits_hdu_open_format(img->name, p->cp.hdu, 0);
       gal_fits_img_info(tmpfits, &p->type, &img->ndim, &img->dsize);
-      gal_wcs_read_from_fitsptr(tmpfits, &img->nwcs, &img->wcs,
-                                p->hstartwcs, p->hendwcs);
+      img->wcs=gal_wcs_read_fitsptr(tmpfits, p->hstartwcs, p->hendwcs,
+                                    &img->nwcs);
       if(img->wcs)
         {
           gal_wcs_decompose_pc_cdelt(img->wcs);
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index e6db9ff..a4a7a1f 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -501,8 +501,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
           break;
         }
   if(readwcs)
-    gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs,
-                 &p->input->wcs);
+    p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs);
 
 
   /* Clean up. */
diff --git a/bin/mknoise/ui.c b/bin/mknoise/ui.c
index 5aa402a..bc88d11 100644
--- a/bin/mknoise/ui.c
+++ b/bin/mknoise/ui.c
@@ -269,8 +269,7 @@ ui_preparations(struct mknoiseparams *p)
 
 
   /* Read the WSC structure. */
-  gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs,
-               &p->input->wcs);
+  p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs);
 
 
   /* If we are dealing with an input table, make sure the format of the
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 453e78c..a1d2f37 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -627,8 +627,7 @@ ui_prepare_canvas(struct mkprofparams *p)
       p->shift[0]=p->shift[1]=0;
 
       /* Read the WCS structure of the background image. */
-      gal_wcs_read(p->backname, p->backhdu, 0, 0, &p->out->nwcs,
-                   &p->out->wcs);
+      p->out->wcs=gal_wcs_read(p->backname, p->backhdu, 0, 0, &p->out->nwcs);
 
     }
   else
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 417b158..ccedd64 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -77,7 +77,7 @@ detection_initial(struct noisechiselparams *p)
 
   /* Erode the image. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  gal_binary_erode(p->binary, p->erode, p->erodengb, 1);
+  gal_binary_erode(p->binary, p->erode, p->erodengb==4 ? 1 : 2, 1);
   if(!p->cp.quiet)
     {
       asprintf(&msg, "Eroded %zu time%s (%zu-connectivity).", p->erode,
@@ -100,7 +100,7 @@ detection_initial(struct noisechiselparams *p)
 
   /* Do the opening. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  gal_binary_open(p->binary, p->opening, p->openingngb, 1);
+  gal_binary_open(p->binary, p->opening, p->openingngb==4 ? 1 : 2, 1);
   if(!p->cp.quiet)
     {
       asprintf(&msg, "Opened (depth: %zu, %s connectivity).",
@@ -251,7 +251,7 @@ detection_fill_holes_open(void *in_prm)
         }
 
       /* Open all the regions. */
-      gal_binary_open(copy, 1, 4, 1);
+      gal_binary_open(copy, 1, 1, 1);
 
       /* Write the copied region back into the large input and AFTERWARDS,
          correct the tile's pointers, the pointers must not be corrected
@@ -956,7 +956,7 @@ detection(struct noisechiselparams *p)
   /* If the user asked for dilation, then apply it. */
   if(p->dilate)
     {
-      gal_binary_dilate(workbin, p->dilate, 8, 1);
+      gal_binary_dilate(workbin, p->dilate, p->input->ndim, 1);
       num_true_initial = gal_binary_connected_components(workbin, &p->olabel,
                                                          8);
     }
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 7848a8e..1559089 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -559,8 +559,7 @@ ui_preparations(struct noisechiselparams *p)
   p->input = gal_fits_img_read_to_type(p->inputname, p->cp.hdu,
                                        GAL_TYPE_FLOAT32,
                                        p->cp.minmapsize);
-  gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs,
-               &p->input->wcs);
+  p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs);
   if(p->input->name==NULL)
     gal_checkset_allocate_copy("INPUT", &p->input->name);
 
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index fd97f56..5e0f0ff 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -792,8 +792,8 @@ ui_preparations(struct statisticsparams *p)
       p->inputformat=INPUT_FORMAT_IMAGE;
       p->input=gal_fits_img_read(p->inputname, cp->hdu, cp->minmapsize);
       if(p->ontile)
-        gal_wcs_read(p->inputname, cp->hdu, 0, 0, &p->input->nwcs,
-                     &p->input->wcs);
+        p->input->wcs=gal_wcs_read(p->inputname, cp->hdu, 0, 0,
+                                   &p->input->nwcs);
     }
   else
     {
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 755fd66..9bec0bb 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -350,8 +350,8 @@ ui_check_options_and_arguments(struct warpparams *p)
       p->input=gal_fits_img_read_to_type(p->inputname, p->cp.hdu,
                                          GAL_TYPE_FLOAT64,
                                          p->cp.minmapsize);
-      gal_wcs_read(p->inputname, p->cp.hdu, p->hstartwcs,
-                   p->hendwcs, &p->input->nwcs, &p->input->wcs);
+      p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, p->hstartwcs,
+                                 p->hendwcs, &p->input->nwcs);
     }
   else
     error(EXIT_FAILURE, 0, "no input file is specified");
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index d1ff72f..92e1337 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -542,6 +542,7 @@ Gnuastro library
 * Dimensions::                  Dealing with coordinates and dimensions.
 * Linked lists::                Various types of linked lists.
 * FITS files::                  Working with FITS data.
+* World Coordinate System::     Dealing with the world coordinate system.
 * Text files::                  Functions to work on Text files.
 * Table input output::          Reading and writing table columns.
 * Arithmetic on datasets::      Arithmetic operations on a dataset.
@@ -549,8 +550,11 @@ Gnuastro library
 * Bounding box::                Finding the bounding box.
 * Polygons::                    Working with the vertices of a polygon.
 * Qsort functions::             Helper functions for Qsort.
+* Permutations::                Re-order (or permute) the values in a dataset.
 * Statistical operations::      Functions for basic statistics.
-* World Coordinate System::     Dealing with the world coordinate system.
+* Binary datasets::             Datasets that can only have values of 0 or 1.
+* Convolution functions::       Library functions to do convolution.
+* Interpolation::               Interpolate (over blank values possibly).
 * Git wrappers::                Wrappers for functions in libgit2.
 
 Multithreaded programming (@file{threads.h})
@@ -16635,6 +16639,7 @@ documentation will correspond to your installed version.
 * Dimensions::                  Dealing with coordinates and dimensions.
 * Linked lists::                Various types of linked lists.
 * FITS files::                  Working with FITS data.
+* World Coordinate System::     Dealing with the world coordinate system.
 * Text files::                  Functions to work on Text files.
 * Table input output::          Reading and writing table columns.
 * Arithmetic on datasets::      Arithmetic operations on a dataset.
@@ -16642,8 +16647,11 @@ documentation will correspond to your installed 
version.
 * Bounding box::                Finding the bounding box.
 * Polygons::                    Working with the vertices of a polygon.
 * Qsort functions::             Helper functions for Qsort.
+* Permutations::                Re-order (or permute) the values in a dataset.
 * Statistical operations::      Functions for basic statistics.
-* World Coordinate System::     Dealing with the world coordinate system.
+* Binary datasets::             Datasets that can only have values of 0 or 1.
+* Convolution functions::       Library functions to do convolution.
+* Interpolation::               Interpolate (over blank values possibly).
 * Git wrappers::                Wrappers for functions in libgit2.
 @end menu
 
@@ -18784,7 +18792,7 @@ Free all the datasets in @code{list} along with all the 
allocated spaces in
 each.
 @end deftypefun
 
address@hidden FITS files, Text files, Linked lists, Gnuastro library
address@hidden FITS files, World Coordinate System, Linked lists, Gnuastro 
library
 @subsection FITS files (@file{fits.h})
 
 @cindex FITS
@@ -19319,7 +19327,148 @@ FITS file (see @ref{List of strings}.
 @end deftypefun
 
 
address@hidden Text files, Table input output, FITS files, Gnuastro library
+
address@hidden World Coordinate System, Text files, FITS files, Gnuastro library
address@hidden World Coordinate System (@file{wcs.h})
+
+The FITS standard defines the world coordinate system (WCS) as a mechanism
+to associate physical values to positions within a dataset. For example, it
+can be used to convert pixel coordinates in an image to celestial
+coordinates like the right ascension and declination. The functions in this
+section are mainly just wrappers over CFITSIO, WCSLIB and GSL library
+functions to help in common applications.
+
+
address@hidden {struct wcsprm *} gal_wcs_read_fitsptr (fitsfile @code{*fptr}, 
size_t @code{hstartwcs}, size_t @code{hendwcs}, int @code{*nwcs})
address@hidden thread-safe}] Return the WCSLIB @code{wcsprm} structure that
+is read from the CFITSIO @code{fptr} pointer to an opened FITS file. Also
+put the number of coordinate representations found into the space that
address@hidden points to. To read the WCS structure directly from a filename,
+see @code{gal_wcs_read} below. After processing has finished, you can free
+the returned structure with WCSLIB's @code{wcsvfree} keyword:
+
address@hidden
+status = wcsvfree(&nwcs,&wcs);
address@hidden example
+
+If you don't want to search the full FITS header for WCS-related FITS
+keywords (for example due to conflicting keywords), but only a specific
+range of the header keywords you can use the @code{hstartwcs} and
address@hidden arguments to specify the keyword number range (counting from
+zero). If @code{hendwcs} is larger than @code{hstartwcs}, then only
+keywords in the given range will be checked. Hence, to ignore this feature
+(and search the full FITS header), give both these arguments the same
+value.
+
+If the WCS information couldn't be read from the FITS file, this function
+will return a @code{NULL} pointer and put a zero in @code{nwcs}. A WCSLIB
+error message will also be printed in @code{stderr} if there was an error.
+
+This function is just a wrapper over WCSLIB's @code{wcspih} function which
+is not thread-safe. Therefore, be sure to not call this function
+simultaneously (over multiple threads).
address@hidden deftypefun
+
address@hidden {struct wcsprm *} gal_wcs_read (char @code{*filename}, char 
@code{*hdu}, size_t @code{hstartwcs}, size_t @code{hendwcs}, int @code{*nwcs})
address@hidden thread-safe}] Return the WCSLIB structure that is read from
+the HDU/extension @code{hdu} of the file @code{filename}. Also put the
+number of coordinate representations found into the space that @code{nwcs}
+points to. Please see @code{gal_wcs_read_fitsptr} for more.
address@hidden deftypefun
+
address@hidden {struct wcsprm *} gal_wcs_copy (struct wcsprm @code{*wcs})
+Return a fully allocated (independent) copy of @code{wcs}.
address@hidden deftypefun
+
address@hidden void gal_wcs_on_tile (gal_data_t @code{*tile})
+Create a WCSLIB @code{wcsprm} structure for @code{tile} using WCS
+parameters of the tile's allocated block dataset, see @ref{Tessellation
+library} for the definition of tiles. If @code{tile} already has a WCS
+structure, this function won't do anything.
+
+In many cases, tiles are created for internal/low-level processing. Hence
+for preformance reasons, when creating the tiles they don't have any WCS
+structure. When needed, this function can be used to add a WCS structure to
+each tile tile by copying the WCS structure of its block and correcting the
+reference point's coordinates within the tile.
address@hidden deftypefun
+
address@hidden {double *} gal_wcs_warp_matrix (struct wcsprm @code{*wcs})
+Return the Warping matrix of the given WCS structure as an array of double
+precision floating points. This will be the final matrix, irrespective of
+the type of storage in the WCS structure. Recall that the FITS standard has
+several methods to store the matrix. The output is an allocated square
+matrix with each side equal to the number of dimensions.
address@hidden deftypefun
+
address@hidden void gal_wcs_decompose_pc_cdelt (struct wcsprm @code{*wcs})
+Decompose the @code{PCi_j} and @code{CDELTi} elements of
address@hidden According to the FITS standard, in the @code{PCi_j} WCS
+formalism, the rotation matrix elements @mymath{m_{ij}} are encoded in the
address@hidden keywords and the scale factors are encoded in the
address@hidden keywords. There is also another formalism (the @code{CDi_j}
+formalism) which merges the two into one matrix.
+
+However, WCSLIB's internal operations are apparently done in the
address@hidden formalism. So its outputs are also all in that format by
+default. When the input is a @code{CDi_j}, WCSLIB will still read the
+matrix directly into the @code{PCi_j} matrix and the @code{CDELTi} values
+are set to @code{1} (one). This function is designed to correct such
+issues: after it is finished, the @code{CDELTi} values in @code{wcs} will
+correspond to the pixel scale, and the @code{PCi_j} will correction show
+the rotation.
address@hidden deftypefun
+
address@hidden double gal_wcs_angular_distance_deg (double @code{r1}, double 
@code{d1}, double @code{r2}, double @code{d2})
+Return the angular distance (in degrees) between a point located at
+(@code{r1}, @code{d1}) to (@code{r2}, @code{d2}). All input coordinates are
+in degrees. The distance (along a great circle) on a sphere between two
+points is calculated with the equation below.
+
address@hidden {\cos(d)=\sin(d_1)\sin(d_2)+\cos(d_1)\cos(d_2)\cos(r_1-r_2)}
+
+However, since the the pixel scales are usually very small numbers, this
+function won't use that direct formula. It will be use the
address@hidden://en.wikipedia.org/wiki/Haversine_formula, Haversine formula}
+which is better considering floating point errors:
+
address@hidden(d)\over 2}=\sin^2\left( {d_1-d_2\over 2} 
\right)+\cos(d_1)\cos(d_2)\sin^2\left( {r_1-r_2\over 2} \right)}
address@hidden deftypefun
+
address@hidden {double *} gal_wcs_pixel_scale_deg (struct wcsprm @code{*wcs})
+Return the pixel scale for each dimension of @code{wcs} in degrees. The
+output is an array of double precision floating point type with one element
+for each dimension.
address@hidden deftypefun
+
address@hidden double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
+Return the pixel area of @code{wcs} in arcsecond squared.
address@hidden deftypefun
+
address@hidden void gal_wcs_world_to_img (struct wcsprm @code{*wcs}, double 
@code{*ra}, double @code{*dec}, double @code{**x}, double @code{**y}, size_t 
@code{size})
+Convert the arrays of input world coordinates (@code{ra} and @code{dec})
+into arrays of image coordinates (@code{x} and @code{y}). Each is assumed
+to be a separate one-dimensional array of @code{size} elements. If
address@hidden or @code{*y==NULL}, then space will be allocated for them
+by this function, otherwise, it is assumed that they already contain the
+space necessary to write the values.
+
+If you don't need the input values after this conversion any more, you can
+pass the pointers of the @code{ra} and @code{dec} arrays and the outputs
+will be written into them. This can help to avoid extra allocations and
+freeing.
address@hidden deftypefun
+
address@hidden void gal_wcs_img_to_world (struct wcsprm @code{*wcs}, double 
@code{*x}, double @code{*y}, double @code{**ra}, double @code{**dec}, size_t 
@code{size})
+Convert the arrays of input image coordinates (@code{x} and @code{y}) into
+arrays of world coordinates (@code{ra} and @code{dec}). Each is assumed to
+be a separate one-dimensional array of @code{size} elements. See
address@hidden for more.
address@hidden deftypefun
+
+
+
address@hidden Text files, Table input output, World Coordinate System, 
Gnuastro library
 @subsection Text files (@file{txt.h})
 
 FITS files are the primary data container in astronomy. FITS indeed as many
@@ -20433,10 +20582,10 @@ declared in @file{gnuastro/box.h}. All coordinates in 
this header are in
 the FITS format (first axis is the horizontal and the second axis is
 vertical).
 
address@hidden void gal_box_ellipse_in_box (double @code{A}, double @code{B}, 
double @code{theta_rad}, long @code{*width})
address@hidden void gal_box_ellipse_in_box (double @code{a}, double @code{b}, 
double @code{theta_rad}, long @code{*width})
 Any ellipse can be enclosed into a rectangular box. The purpose of this
-function is to give the height and width of that box. @code{A} is the
-ellipse major axis, @code{B} is the minor axis, @code{theta_rad} is the
+function is to give the height and width of that box. @code{a} is the
+ellipse major axis, @code{b} is the minor axis, @code{theta_rad} is the
 position angle in radians. The @code{width} array will contain the output
 size in long integer type. @code{width[0]}, and @code{width[1]} are the
 number of pixels along the first and second FITS axis.
@@ -20573,7 +20722,7 @@ both polygons have to be sorted in an anti-clock-wise 
manner.
 
 
 
address@hidden Qsort functions, Statistical operations, Polygons, Gnuastro 
library
address@hidden Qsort functions, Permutations, Polygons, Gnuastro library
 @subsection Qsort functions (@file{qsort.h})
 
 The C programming language comes with the @code{qsort} (Quick sort)
@@ -20639,12 +20788,63 @@ types}, for example 
@code{gal_qsort_int32_decreasing}, or
 
 
 
address@hidden Permutations, Statistical operations, Qsort functions, Gnuastro 
library
address@hidden Permutations (@file{permutation.h})
address@hidden permutation
+Permutation is the technical name for re-ordering of values. The need for
+permutations occurs a lot during (mainly low-level) processing. To do
+permutation, you must provide two inputs: an array of values (that you want
+to re-order inplace) and a permutation array which contains the new index
+of each element (let's call it @code{perm}). The diagram below shows the
+input array before and after the re-ordering.
 
address@hidden
+permute:    AFTER[ i       ] = BEFORE[ perm[i] ]     i = 0 .. N-1
+inverse:    AFTER[ perm[i] ] = BEFORE[ i       ]     i = 0 .. N-1
address@hidden example
 
address@hidden GNU Scientific Library
+The functions here are a re-implementation of the GNU Scientific Library's
address@hidden function. The reason we didn't use that function was
+that it uses system-specific types (like @code{long} and @code{int}) which
+can have different widths on different systems, hence are not easily
+convertable to Gnuastro's fixed width types (see @ref{Numeric data
+types}). There is also a separate function for each type, heavily using
+macros to allow a @code{base} function to work on all the types. Thus it is
+hard to read/understand. Hence, Gnuastro contains a re-write of their steps
+in a new type-agnostic method which is a single function that can work on
+any type.
+
+As described in GSL's source code and manual, this implementation comes
+from Donald Knuth's @emph{Art of computer programming} book, in the
+"Sorting and Searching" chapter of Volume 3 (3rd ed). Exercise 10 of
+Section 5.2 defines the problem and in the answers, Knuth describes the
+solution. So if you are interested, please have a look there for more.
+
+We are in contact with the GSL developers and in the
address@hidden's @url{http://savannah.gnu.org/task/?14497, Task
+14497}. If this task is still ``postponed'' when you are reading this and
+you are interested to help, your help would be very welcome. Both Gnuastro
+and GSL developers are very busy, hence both would appreciate your help.}
+we will submit these implementations to GSL. If they are finally
+incorporated there, we will delete this section in future versions.
+
address@hidden void gal_permutation_check (size_t @code{*permutation}, size_t 
@code{size})
+Print how @code{permutation} will re-order an array that has @code{size}
+elements for each element in one one line.
address@hidden deftypefun
 
address@hidden void gal_permutation_apply (gal_data_t @code{*input}, size_t 
@code{*permutation})
+Apply @code{permutation} on the @code{input} dataset (can have any type),
+see above for the definition of permutation.
address@hidden deftypefun
 
address@hidden void gal_permutation_apply_inverse (gal_data_t @code{*input}, 
size_t @code{*permutation})
+Apply the inverse of @code{permutation} on the @code{input} dataset (can
+have any type), see above for the definition of permutation.
address@hidden deftypefun
 
address@hidden Statistical operations, World Coordinate System, Qsort 
functions, Gnuastro library
address@hidden Statistical operations, Binary datasets, Permutations, Gnuastro 
library
 @subsection Statistical operations (@file{statistics.h})
 
 After reading a dataset into memory from a file or fully simulating it with
@@ -20925,91 +21125,309 @@ array[3]: Standard deviation.
 
 
 
address@hidden Binary datasets, Convolution functions, Statistical operations, 
Gnuastro library
address@hidden Binary datasets (@file{binary.h})
+
address@hidden Thresholding
address@hidden Binary datasets
address@hidden Dataset: binary
+Binary datasets only have two (usable) values: 0 (also known as background)
+or 1 (also known as foreground). They are created after some binary
+classificataion is applied to the dataset. The most common is thresholding:
+in an image, pixels with a value above the threshold is commonly assigned a
+value of 1 and those that have a value less than the threshold are assigned
+a value of 0.
+
address@hidden Connectivity
address@hidden Immediate neighbors
address@hidden Neighbors, immediate
+Since there is only two values, in the processing of binary images, you are
+usually concered with the positioning of an element and its vicinity
+(neighbors). When a dataset has more than one dimension, multiple classes
+of immediate neighbors (that are touching the element) can be defined for
+each data-element. To separate these different classes of immediate
+neighbors, we define @emph{connectivity}.
+
+The classification is done by the distance from element center to the
+neighbor's center. The nearest immediate neighbors have a connectivity of
+1, the second nearest class of neighbors have a connectivity of 2 and so
+on. In total, the largest possible connectivity for data with @code{ndim}
+dimensions is @code{ndim}. For example in a 2D dataset, 4-connected
+neighbors (that share an edge and have a distance of 1 pixel) have a
+connectivity of 1. The other 4 neighbors that only share a vertice (with a
+distance of @mymath{\sqrt{2}} pixels) have a connectivity of
+2. Conventionally, the class of connectivity-2 neighbors also includes the
+connectivity 1 neighbors, so for example we call them 8-connected neighbors
+in 2D datasets.
+
+Ideally, one bit is sufficient for each element of a binary
+dataset. However, CPUs are not designed to work on individual bits, the
+smallest unit of memory addresses is a byte (containing 8 bits on modern
+CPUs). Therefore, in Gnuastro, the type used for binary dataset is
address@hidden (see @ref{Numeric data types}). Although it does take
+8-times more memory, this choice offers much better performance and the
+some extra (useful) features.
+
+The advantage of using a full byte for each element of a binary dataset is
+that you can also have other values (that will be ignored in the
+processing). One such common ``other'' value in real datasets is a blank
+value (to mark regions that should not be processed because there is no
+data). The constant @code{GAL_BLANK_UINT8} value must be used in these
+cases (see @ref{Library blank values}). Another is some temporary value(s)
+that can be given to a processed pixel to avoid having another copy of the
+dataset as in @code{GAL_BINARY_TMP_VALUE} that is described below.
+
address@hidden Macro GAL_BINARY_TMP_VALUE
+The functions described below work on a @code{uint8_t} type dataset with
+values of 1 or 0 (no other pixel will be touched). However, in some cases,
+it is necessary to put temporary values in each element during the
+processing of the functions. This temporary value has a special meaning for
+the operation and will be operated on. So if your input datasets have
+values other than 0 and 1 that you don't want these functions to work on,
+be sure they are not equal to this macro's value. Note that this value is
+also different from @code{GAL_BLANK_UINT8}, so your input datasets may also
+contain blank elements.
address@hidden deffn
 
 
address@hidden World Coordinate System, Git wrappers, Statistical operations, 
Gnuastro library
address@hidden World Coordinate System (@file{wcs.h})
address@hidden {gal_data_t *} gal_binary_erode (gal_data_t @code{*input}, 
size_t @code{num}, int @code{connectivity}, int @code{inplace})
+Do @code{num} erosions on the @code{connectivity}-connected neighbors of
address@hidden (see above for the defintion of connectivity).
+
+If @code{inplace} is non-zero @emph{and} the input's type is
address@hidden, then the erosion will be done within the input
+dataset and the returned pointer will be @code{input}. Otherwise,
address@hidden is copied (and converted if necessary) to
address@hidden and erosion will be done on this new dataset which
+will also be returned. This function will only work on the elements with a
+value of 1 or 0. It will leave all the rest unchanged.
 
-The FITS world coordinate system is a mechanism to give physical values to
-the data points. In the case of images, it can be used to convert pixel
-images to celestial coordinates like the right ascension and
-declination. Gnuastro doesn't have any functions to actually do the
-transformations, the standard tool for doing this in astronomy is WCSLIB
-(see @ref{WCSLIB}). The functions introduced here are just wrappers around
-the lower-level functions provided by WCSLIB.
-
address@hidden void gal_wcs_xy_array_to_radec (struct wcsprm @code{*wcs}, 
double @code{*xy}, double @code{*radec}, size_t @code{number}, size_t 
@code{stride})
-Convert the image coordinates @code{xy} (a two column array with
address@hidden rows) to world coordinates @code{wcs} (also a two column
-array, which has to have been allocated prior to this function) using the
-conversion parameters in @code{wcs}. @code{stride} is the number of columns
-in a larger array, see below.
-
-It might happen that you have one array which contains both the image and
-world coordinates along with many more columns. @code{stride} is the width
-of that larger array. If such an array doesn't exist and @code{xy} and
address@hidden are separately allocated arrays with
address@hidden@code{number} elements, then just set the value of
address@hidden to @code{2}.
address@hidden Erosion
address@hidden Mathematical morphology
+Erosion (inverse of dilation) is an operation in mathematical morphology
+where each foreground pixel that is touching a background pixel is flipped
+(changed to background). The @code{connectivity} value determines the
+definition of ``touching''. Erosion will thus decrease the area of the
+foregound regions by one layer of pixels.
 @end deftypefun
 
address@hidden void gal_wcs_radec_array_to_xy (struct wcsprm @code{*wcs}, 
double @code{*radec}, double @code{*xy}, size_t @code{number}, size_t 
@code{stride})
-Similar to @code{gal_wcs_xy_array_to_radec} but convert the world
-coordinates in @code{radec} to image coordinates in @code{xy}.
address@hidden {gal_data_t *} gal_binary_dilate (gal_data_t @code{*input}, 
size_t @code{num}, int @code{connectivity}, int @code{inplace})
+Do @code{num} dilations on the @code{connectivity}-connected neighbors of
address@hidden (see above for the defintion of connectivity). For more on
address@hidden and the output, see @code{gal_binary_erode}.
+
address@hidden Dilation
+Dilation (inverse of erosion) is an operation in mathematical morphology
+where each background pixel that is touching a foreground pixel is flipped
+(changed to foreground). The @code{connectivity} value determines the
+definition of ``touching''. Dilation will thus increase the area of the
+foregound regions by one layer of pixels.
 @end deftypefun
 
address@hidden double gal_wcs_angular_distance_deg (double @code{r1}, double 
@code{d1}, double @code{r2}, double @code{d2})
-Return the angular distance (in degrees) between a point located at
-(@code{r1}, @code{d1}) to (@code{r2}, @code{d2}). All input coordinates are
-in degrees. The distance (along a great circle) on a sphere between two
-points is calculated with the equation below. Since the pixel sides are
-usually very small, we won't be using the direct formula:
address@hidden {gal_data_t *} gal_binary_open (gal_data_t @code{*input}, size_t 
@code{num}, int @code{connectivity}, int @code{inplace})
+Do @code{num} openings on the @code{connectivity}-connected neighbors of
address@hidden (see above for the defintion of connectivity). For more on
address@hidden and the output, see @code{gal_binary_erode}.
+
address@hidden Opening (Mathematical morphology)
+Opening is an operation in mathematical morphology which is defined as
+erosion followed by dilation (see above for the definitions of erosion and
+dilation). Opening will thus remove the outer structure of the
+foreground. In this implementation, @code{num} erosions are going to be
+applied on the dataset, then @code{num} dilations.
address@hidden deftypefun
 
address@hidden {\cos(d)=\sin(d_1)\sin(d_2)+\cos(d_1)\cos(d_2)\cos(r_1-r_2)}
 
-Here, we will be using the
address@hidden://en.wikipedia.org/wiki/Haversine_formula, Haversine formula}
-which is better considering floating point errors:
address@hidden size_t gal_binary_connected_components (gal_data_t 
@code{*binary}, gal_data_t @code{**out}, int @code{connectivity})
address@hidden Connected components
address@hidden Breadth first search
+Return the number of connected components in @code{binary} through the
+breadth first search algorithm (finding all pixels belonging to one
+component before going on to the next). Connection between two pixels is
+defined based on the value to @code{connectivity}. @code{out} is a dataset
+with the same size as @code{binary} with @code{GAL_TYPE_INT32} type. Every
+pixel in @code{out} will have the label of the connected component it
+belongs to. The labeling of connected components starts from 1, so a label
+of zero is given to the input's background pixels.
+
+When @code{*out!=NULL} (its space is already allocated), it will be cleared
+(to zero) at the start of this function. Otherwise, when @code{*out==NULL},
+the necessary dataset to keep the output will be allocated by this
+function.
 
address@hidden(d)\over 2}=\sin^2\left( {d_1-d_2\over 2} 
\right)+\cos(d_1)\cos(d_2)\sin^2\left( {r_1-r_2\over 2} \right)}
address@hidden must have a type of @code{GAL_TYPE_UINT8}, otherwise this
+function will abort with an error. Other than blank pixels (with a value of
address@hidden defined in @ref{Library blank values}), all other
+non-zero pixels in @code{binary} will be considered as foreground (and will
+be labeled). Blank pixels in the input will also be blank in the output.
address@hidden deftypefun
+
address@hidden {gal_data_t *} gal_binary_connected_adjacency_matrix (gal_data_t 
@code{*adjacency}, size_t @code{*numconnected})
address@hidden Adjacency matrix
address@hidden Matrix, adjacency
+Find the number of connected labels and new labels based on an an adjacency
+matrix (which must be binary array of type @code{GAL_TYPE_UINT8}). The
+returned dataset is a list of new labels for each old label. In other
+words, this function will find the objects that are connected (possibly
+through a third object) and in the output array, the respective elements
+for all input labels is going to have the same value. The total number of
+connected labels is put into the space that @code{numconnected} points to.
+
+An adjacency matrix defines connection between two labels. For example,
+let's assume we have 5 labels and we know that labels 1 and 5 are connected
+to label 3, but are not connected with each other. Also, labels 2 and 4 are
+not touching any other label. So in total we have 3 final labels: one
+combined object (merged from labels 1, 3, and 5) and the initial labels 2
+and 4. The input adjacency matrix would look like this (note the extra row
+and column for a label 0 which is ignored):
+
address@hidden
+            INPUT                             OUTPUT
+            =====                             ======
+          in_lab  1  2  3  4  5   |
+                                  |       numconnected = 3
+               0  0  0  0  0  0   |
+in_lab 1 -->   0  0  0  1  0  0   |
+in_lab 2 -->   0  0  0  0  0  0   |  Returned: new labels for the
+in_lab 3 -->   0  1  0  0  0  1   |       5 initial objects
+in_lab 4 -->   0  0  0  0  0  0   |   | 0 | 1 | 2 | 1 | 3 | 1 |
+in_lab 5 -->   0  0  0  1  0  0   |
address@hidden example
+
+Although the adjacency matrix as used here is symmetric, currently this
+function assumes that it is filled on both sides of the diagonal.
 @end deftypefun
 
address@hidden void gal_wcs_pixel_scale_deg (struct wcsprm @code{*wcs}, double 
@code{*dx}, double @code{*dy})
-Return the pixel scale of the WCS structure @code{wcs} for the two image
-dimensions @code{dx} and @code{dy} in degrees. Here, the X-axis is the
-first FITS axis, or the horizontal axis when viewed in SAO ds9.
address@hidden void gal_binary_fill_holes (gal_data_t @code{*input})
+Fill all the holes that are bounded within a 4-connected region of the
+binary @code{input} dataset. This function currently only works on a 2D
+dataset.
 @end deftypefun
 
address@hidden double gal_wcs_pixel_area_arcsec2 (struct wcsprm @code{*wcs})
-Given the WCS structure @code{wcs} calculate the pixel area in
-arcsec-squared.
address@hidden Convolution functions, Interpolation, Binary datasets, Gnuastro 
library
address@hidden Convolution functions (@file{convolve.h})
+
+Convolution is a very common operation during data analysis and is
+thoroughly described as part of Gnuastro's @ref{Convolve} program which is
+fully devoted to this job. Because of the complete introduction that was
+presented there, we will directly skip onto the currently available
+convolution functions in Gnuastro's library.
+
+As of this version, only spatial domain convolution is available in
+Gnuastro's libraries. We haven't had the time to liberate the frequency
+domain function convolution and de-convolution functions that are available
+in the Convolve address@hidden any help would be greatly
+appreciated.}.
+
address@hidden {gal_data_t *} gal_convolve_spatial (gal_data_t @code{*tiles}, 
gal_data_t @code{*kernel}, size_t @code{numthreads}, int @code{edgecorrection}, 
int @code{convoverch})
+Convolve each node of the list of @code{tiles} (see @ref{List of
+gal_data_t} and @ref{Tessellation library}) with @code{kernel} using
address@hidden When @code{edgecorrection} is non-zero, it will correct
+for the edge dimming effects as discussed in @ref{Edges in the spatial
+domain}.
+
+To create a tessellation that fully covers an input image, you may use
address@hidden, or @code{gal_tile_full_two_layers} to also define
+channels over your input dataset. These functions are discussed in
address@hidden grid}. You may then pass the list of tiles to this function. This
+is the recommended way to call this function because spatial domain
+convolution is slow and breaking the job into many small tiles and working
+on simultaneously on several threads can greatly speed up the processing.
+
+If the tiles are defined within a channel (a larger tile), by default
+convolution will be done within the channel, so pixels on the edge of a
+channel will not be affected by their neighbors that are in another
+channel. See @ref{Tessellation} for the necessity of channels in
+astronomical data analysis. This behavior may be disabled when
address@hidden is non-zero. In this case, it will ignore channel borders
+(if they exist) and mix all pixels that cover the kernel within the
+dataset.
 @end deftypefun
 
address@hidden Git wrappers,  , World Coordinate System, Gnuastro library
address@hidden void gal_convolve_spatial_correct_ch_edge (gal_data_t 
@code{*tiles}, gal_data_t @code{*kernel}, size_t @code{numthreads}, int 
@code{edgecorrection}, gal_data_t @code{*tocorrect})
+Correct the edges of channels in an already convolved image when it was
+initially convolved with @code{gal_convolve_spatial} and
address@hidden In that case, strong boundaries might exist on the
+channel edges. So if you later need to remove those boundaries at later
+steps of your processing, you can call this function. It will only do
+convolution on the tiles that are near the edge and were effected by the
+channel borders. Other pixels in the image will not be touched. Hence, it
+is much faster.
address@hidden deftypefun
+
address@hidden Interpolation, Git wrappers, Convolution functions, Gnuastro 
library
address@hidden Interpolation (@file{interpolate.h})
+
+During data anlysis, it often happens that parts of the data cannot be
+given a value. For example your image was saturated due to a very bright
+start and you have to mask that star's footprint. One other common
+situation in Gnuastro is when we do processing on tiles (for example to
+estimate the Sky value and its Standard deviation, see @ref{Sky
+value}). Some tiles must not be used for the estimation of the Sky value,
+for example because they cover a large galaxy. So we need to fill them in
+with blank values. But ultimately, we need a Sky value for every
+pixel. This job (assigning a value to blank element(s) based on their
+nearby neighbors with a value) is known as interpolation.
+
+There are many ways to do interpolation, but (mainly due to lack of time),
+currently Gnuastro only contains the (median of) nearest-neighbor
+method. The power of this method of interpolation is its non-parametric
+nature. The produced values are also always within the range of the known
+values and strong outliers do not get created. We will hopefully implement
+other methods too (wrappers around the GNU Scientific Library's very
+complete set of functions), but currently the developers are too busy. So
+if you do have the chance to help your contribution would be very welcome
+and appreciated.
+
address@hidden {gal_data_t *} gal_interpolate_close_neighbors (gal_data_t 
@code{*input}, struct gal_tile_two_layer_params @code{*tl}, size_t 
@code{numneighbors}, size_t @code{numthreads}, int @code{onlyblank}, int 
@code{aslinkedlist})
+Interpolate the values in the image using the median value of their
address@hidden closest neighbors. If @code{onlyblank} is non-zero,
+then only blank elements will be interpolated and pixels that already have
+a value will be left untouched. This function is multi-threaded and will
+run on @code{numthreads} threads (see @code{gal_threads_number} in
address@hidden programming}).
+
address@hidden is Gnuastro's two later tessellation structure used to define
+tiles over an image and is fully described in @ref{Tile grid}. When
address@hidden, then it is assumed that the @code{input->array} contains
+one value per tile and interpolation will respect certain tessellation
+properties, for example to not interpolate over channel borders.
+
+If several datasets have the same set of blank values, you don't need to
+call this function multiple times. When @code{aslinkedlist} is non-zero,
+then @code{input} will be seen as a @ref{List of gal_data_t} and for all
+the same neighbor position checking will be done for all the datasets in
+the list. Ofcourse, the values for each dataset will be different, so a
+different value will be written in the each dataset, but the neighbor
+checking that is the most CPU intensive part will only be done once.
address@hidden deftypefun
+
+
address@hidden Git wrappers,  , Interpolation, Gnuastro library
 @subsection Git wrappers (@file{git.h})
 
 @cindex Git
 @cindex libgit2
 Git is one of the most common tools for version control and it can often be
 useful during development, for example see @code{COMMIT} keyword in
address@hidden headers}. The functions introduced here are described in the
address@hidden/git.h} header. At installation time, Gnuastro will also
-check for the existence of libgit2 and store the value in the
address@hidden, see @ref{Configuration
-information}. @file{gnuastro/git.h} includes @file{gnuastro/config.h}
-internally, so won't have to include both for this macro.
address@hidden headers}. At installation time, Gnuastro will also check for
+the existence of libgit2, and store the value in the
address@hidden, see @ref{Configuration information} and
address@hidden dependencies}. @file{gnuastro/git.h} includes
address@hidden/config.h} internally, so you won't have to include both for
+this macro.
 
 @deftypefun {char *} gal_git_describe ( )
 When libgit2 is present and the program is called within a directory that
-is version controlled, then this function will return a string containing
-the commit description (similar to Gnuastro's unofficial version number,
-see @ref{Version numbering}). If there are uncommitted changes in the
-running directory, it will add a address@hidden' prefix to the
-description. When there is no tagged point in the previous commit, this
-function will return a uniquely abbreviated commit object as fallback. This
-function is used for generating the value of the @code{COMMIT} keyword in
address@hidden headers}. The output string is similar to the output of the
-following command:
+is version controlled, this function will return a string containing the
+commit description (similar to Gnuastro's unofficial version number, see
address@hidden numbering}). If there are uncommitted changes in the running
+directory, it will add a address@hidden' prefix to the description. When
+there is no tagged point in the previous commit, this function will return
+a uniquely abbreviated commit object as fallback. This function is used for
+generating the value of the @code{COMMIT} keyword in @ref{Output
+headers}. The output string is similar to the output of the following
+command:
 
 @example
 $ git describe --dirty --always
diff --git a/lib/binary.c b/lib/binary.c
index ee3ced4..da884d0 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -256,8 +256,8 @@ binary_erode_dilate(gal_data_t *input, size_t num, int 
connectivity,
       for(counter=0;counter<num;++counter)
         switch(connectivity)
           {
-          case 4: binary_erode_dilate_2d_4con(binary, d0e1); break;
-          case 8: binary_erode_dilate_2d_8con(binary, d0e1); break;
+          case 1: binary_erode_dilate_2d_4con(binary, d0e1); break;
+          case 2: binary_erode_dilate_2d_8con(binary, d0e1); break;
           default:
             error(EXIT_FAILURE, 0, "%s: %d not acceptable for connectivity "
                   "in a 2D dataset", __func__, connectivity);
@@ -340,14 +340,7 @@ gal_binary_open(gal_data_t *input, size_t num, int 
connectivity,
 /*********************************************************************/
 /*****************      Connected components      ********************/
 /*********************************************************************/
-/* Find the connected components in the input binary dataset `binary'
-   through the breadth first search algorithm. `binary' has to have an
-   `uint8' datatype and only zero and non-zero values in it will be
-   distinguished. The output dataset (which will contain a label on each
-   pixel) maybe already allocated (with type `int32'). If `*out!=NULL', the
-   labels will be reset to zero before the start and the labels will be
-   written into it. If `*out==NULL', the necessary dataset will be
-   allocated here and put into it. */
+/* Find connected components in an intput dataset. */
 size_t
 gal_binary_connected_components(gal_data_t *binary, gal_data_t **out,
                                 int connectivity)
@@ -456,23 +449,7 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
 
 /* Given an adjacency matrix (which should be binary), find the number of
    connected objects and return an array of new labels for each old
-   label. In other words, this function will find the objects that are
-   connected (possibly through a third object) and in the output array, the
-   respective elements for all these objects is going to have the same
-   value. The total number of connected labels is put into the place
-   pointed to by `numconnected'.
-
-   Labels begin from 1 (0 is kept for non-labeled regions usually). So if
-   you have 3 initial objects/labels, the input matrix to this function
-   should have a size of 4x4. The first (label 0) row and column are not
-   going to be parsed/checked.
-
-   The adjacency matrix needs to be completely filled (on both sides of the
-   diagonal) for this function.
-
-   If the input adjacency matrix has a size of `amsize * amsize', the
-   output will have a size of `amsize' with each index having a new label
-   in its place. */
+   label.  */
 gal_data_t *
 gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
                                       size_t *numconnected)
@@ -527,7 +504,7 @@ gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
                 /* Give it the new label. */
                 newlabs[p]=curlab;
 
-                /* Go over the adjacecny matrix row for this touching
+                /* Go over the adjacency matrix row for this touching
                    object and see if there are any not-yet-labeled objects
                    that are touching it. */
                 for(j=1;j<num;++j)
diff --git a/lib/data.c b/lib/data.c
index e06c243..eaf59b4 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -266,7 +266,7 @@ gal_data_initialize(gal_data_t *data, void *array, uint8_t 
type,
 
 
   /* Copy the WCS structure. */
-  data->wcs=gal_wcs_copy(wcs, ndim);
+  data->wcs=gal_wcs_copy(wcs);
 
 
   /* Allocate space for the dsize array, only if the data are to have any
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 278666d..8a2a70e 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -55,13 +55,13 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 /*************************************************************
  ***********               Read WCS                ***********
  *************************************************************/
-void
-gal_wcs_read_from_fitsptr(fitsfile *fptr, int *nwcs, struct wcsprm **wcs,
-                          size_t hstartwcs, size_t hendwcs);
+struct wcsprm *
+gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, size_t hendwcs,
+                     int *nwcs);
 
-void
+struct wcsprm *
 gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
-             size_t hendwcs, int *nwcs, struct wcsprm **wcs);
+             size_t hendwcs, int *nwcs);
 
 
 
@@ -71,7 +71,7 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
 /**********              Utilities                 ************/
 /**************************************************************/
 struct wcsprm *
-gal_wcs_copy(struct wcsprm *in, size_t ndim);
+gal_wcs_copy(struct wcsprm *wcs);
 
 void
 gal_wcs_on_tile(gal_data_t *tile);
@@ -99,10 +99,6 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
 /**********              Conversion                ************/
 /**************************************************************/
 void
-gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy, double *radec,
-                          size_t number, size_t width);
-
-void
 gal_wcs_world_to_img(struct wcsprm *wcs, double *ra, double *dec,
                      double **x, double **y, size_t size);
 
diff --git a/lib/interpolate.c b/lib/interpolate.c
index ba3b30f..a3a98d4 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -45,16 +45,16 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /*********************************************************************/
 /********************      Nearest neighbor       ********************/
 /*********************************************************************/
-/* We want the flags to be powers of two so we can use bit-wise checks. */
+/* These are bit-flags, so we're using hexadecimal notation. */
 #define INTERPOLATE_FLAGS_NO       0
-#define INTERPOLATE_FLAGS_CHECKED  1<<0
-#define INTERPOLATE_FLAGS_BLANK    1<<1
+#define INTERPOLATE_FLAGS_CHECKED  0x1
+#define INTERPOLATE_FLAGS_BLANK    0x2
 
 
 
 
 
-/* Paramters for interpolation on threads. */
+/* Parameters for interpolation on threads. */
 struct interpolate_params
 {
   gal_data_t                    *input;
diff --git a/lib/permutation.c b/lib/permutation.c
index d3e532d..892f1b2 100644
--- a/lib/permutation.c
+++ b/lib/permutation.c
@@ -38,7 +38,7 @@ gal_permutation_check(size_t *permutation, size_t size)
 {
   size_t i;
   for(i=0; i<size; ++i)
-    printf("out[ %-4zu ]    =     in [ %-4zu ]\n", i, permutation[i]);
+    printf("after[ %-5zu ]    =   before [ %-5zu ]\n", i, permutation[i]);
 }
 
 
diff --git a/lib/wcs.c b/lib/wcs.c
index 6c02d84..c0d9b08 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -65,11 +65,12 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    ===================================
    Don't call this function within a thread or use a mutex.
 */
-void
-gal_wcs_read_from_fitsptr(fitsfile *fptr, int *nwcs, struct wcsprm **wcs,
-                          size_t hstartwcs, size_t hendwcs)
+struct wcsprm *
+gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, size_t hendwcs,
+                     int *nwcs)
 {
   /* Declaratins: */
+  struct wcsprm *wcs;
   int nkeys=0, status=0;
   char *fullheader, *to, *from;
   int relax    = WCSHDR_all; /* Macro: use all informal WCS extensions. */
@@ -111,28 +112,28 @@ gal_wcs_read_from_fitsptr(fitsfile *fptr, int *nwcs, 
struct wcsprm **wcs,
     }
 
   /* WCSlib function */
-  status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, nwcs, wcs);
+  status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, nwcs, &wcs);
   if(status)
     {
       fprintf(stderr, "\n##################\n"
               "WCSLIB Warning: wcspih ERROR %d: %s.\n"
               "##################\n",
               status, wcs_errmsg[status]);
-      *wcs=NULL; *nwcs=0;
+      wcs=NULL; *nwcs=0;
     }
   if (fits_free_memory(fullheader, &status) )
     gal_fits_io_error(status, "problem in fitsarrayvv.c for freeing "
                            "the memory used to keep all the headers");
 
   /* Set the internal structure: */
-  status=wcsset(*wcs);
+  status=wcsset(wcs);
   if(status)
     {
       fprintf(stderr, "\n##################\n"
               "WCSLIB Warning: wcsset ERROR %d: %s.\n"
               "##################\n",
             status, wcs_errmsg[status]);
-      *wcs=NULL; *nwcs=0;
+      wcs=NULL; *nwcs=0;
     }
 
   /* Initialize the wcsprm struct
@@ -140,28 +141,33 @@ gal_wcs_read_from_fitsptr(fitsfile *fptr, int *nwcs, 
struct wcsprm **wcs,
     error(EXIT_FAILURE, 0, "%s: wcsset ERROR %d: %s.\n", __func__,
           status, wcs_errmsg[status]);
   */
+
+  /* Return the WCS structure. */
+  return wcs;
 }
 
 
 
 
 
-void
+struct wcsprm *
 gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
-             size_t hendwcs, int *nwcs, struct wcsprm **wcs)
+             size_t hendwcs, int *nwcs)
 {
   int status=0;
   fitsfile *fptr;
+  struct wcsprm *wcs;
 
   /* Check HDU for realistic conditions: */
   fptr=gal_fits_hdu_open_format(filename, hdu, 0);
 
   /* Read the WCS information: */
-  gal_wcs_read_from_fitsptr(fptr, nwcs, wcs, hstartwcs, hendwcs);
+  wcs=gal_wcs_read_fitsptr(fptr, hstartwcs, hendwcs, nwcs);
 
-  /* Close the FITS file: */
+  /* Close the FITS file and return. */
   fits_close_file(fptr, &status);
   gal_fits_io_error(status, NULL);
+  return wcs;
 }
 
 
@@ -188,12 +194,12 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
 /**************************************************************/
 /* Copy a given WSC structure into another one. */
 struct wcsprm *
-gal_wcs_copy(struct wcsprm *in, size_t ndim)
+gal_wcs_copy(struct wcsprm *wcs)
 {
   struct wcsprm *out;
 
   /* If the input WCS is NULL, return a NULL WCS. */
-  if(in)
+  if(wcs)
     {
       /* Allocate the output WCS structure. */
       errno=0;
@@ -206,10 +212,10 @@ gal_wcs_copy(struct wcsprm *in, size_t ndim)
          the first invokation, and only the first invokation, wcsprm::flag
          must be set to -1 to initialize memory management"*/
       out->flag=-1;
-      wcsini(1, ndim, out);
+      wcsini(1, wcs->naxis, out);
 
       /* Copy the input WCS to the output WSC structure. */
-      wcscopy(1, in, out);
+      wcscopy(1, wcs, out);
     }
   else
     out=NULL;
@@ -241,7 +247,7 @@ gal_wcs_on_tile(gal_data_t *tile)
   else
     {
       /* Copy the block's WCS into the tile. */
-      tile->wcs=gal_wcs_copy(block->wcs, block->ndim);
+      tile->wcs=gal_wcs_copy(block->wcs);
 
       /* Find the coordinates of the tile's starting index. */
       start_ind=gal_data_ptr_dist(block->array, tile->array, block->type);
@@ -444,10 +450,7 @@ gal_wcs_pixel_scale_deg(struct wcsprm *wcs)
 
 
 /* Report the arcsec^2 area of the pixels in the image based on the
-   WCS information in that image. We first use the angular distance of
-   two edges of one pixel in radians. Then the radians are multiplied
-   to give stradians and finally, the stradians are converted to
-   arcsec^2. */
+   WCS information in that image. */
 double
 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
 {



reply via email to

[Prev in Thread] Current Thread [Next in Thread]