gnuastro-commits
[Top][All Lists]
Advanced

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

[gnuastro-commits] master 7dca196: Fits: --wcsdistortion with no data HD


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 7dca196: Fits: --wcsdistortion with no data HDUs, no crash on old WCSLIB
Date: Mon, 22 Jun 2020 21:43:56 -0400 (EDT)

branch: master
commit 7dca196b6f7f588482772f3c059866647e812689
Author: Sachin Kumar Singh <sachinkumarsingh092@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Fits: --wcsdistortion with no data HDUs, no crash on old WCSLIB
    
    Until now, when the input FITS header given to the '--wcsdistortion' option
    didn't have any data, it would crash, complaining about this! But WCS
    distortions are may be stored in FITS HDUs that don't have any data (and
    infact this is what happens with astrometry.net for example). Hence it was
    necessary to account for this scenario.
    
    To enable this, with this commit, a library function 'gal_wcs_write' was
    added to only write the WCS-related keywords in a data-free FITS header and
    the necessary checks were placed to avoid reading data when there isn't
    any.
    
    In the process (on the same branch that this commit was "rebase"d from),
    the following corrections are also made:
    
     - gcc allows 2D matrix declaration and initialiasation in the form
       `array[i][j]={0}` but clang raises warning for they same. To provide
       uniformity we used nessted loops to initialise these 2D matrices.
    
     - Until now we were assuming that the user's 'wcsprm' structure has an
       'mjdref' element, but after testing on a system with WCSLIB 5.17, we
       noticed that this element is a relatively recent addition to the
       'wcsprm' structure (it wasn't present in that version). So with this
       commit, a configure-time check has been added so if the host's WCSSLIB
       doesn't have this structure, it isn't used.
    
     - On WCSLIB 5.17, we also noticed that 'dpkeyd' also doesn't exist, thus
       causing an error during compilation. We thus don't use it any more. It
       was mainly used for extra error check and after a discussion was
       scraped.
    
     - Until now, we always assumed that WCSLIB has the 'dis.h' header for
       distortion-related operations. But this is not always the case and older
       versions of WCSLIB don't have it. This could cause problems in the
       build. With this commit, configure-time checks have been added to see if
       'dis.h' is present on the host or not. If it isn't the respective
       functions will print an informative error and abort. Also a warning will
       be printed at the end of the configuration step.
---
 NEWS                              |   3 +
 bin/fits/keywords.c               |  40 +++-
 bin/fits/ui.c                     |   1 +
 configure.ac                      |  26 +++
 doc/gnuastro.texi                 |  13 ++
 lib/Makefile.am                   |  94 +++++----
 lib/gnuastro-internal/config.h.in |   2 +
 lib/gnuastro/wcs.h                |  10 +
 lib/wcs.c                         | 109 +++++++++-
 lib/wcsdistortion.c               | 416 ++++++++++++++++++--------------------
 10 files changed, 442 insertions(+), 272 deletions(-)

diff --git a/NEWS b/NEWS
index 7773443..043c42d 100644
--- a/NEWS
+++ b/NEWS
@@ -22,7 +22,10 @@ See the end of the file for license conditions.
      file. The inverse conversion is also supported (from SIP to TPV).
 
   Library:
+   - GAL_CONFIG_HAVE_WCSLIB_DIS_H: if the host's WCSLIB supports distortions.
+   - GAL_CONFIG_HAVE_WCSLIB_MJDREF: if the host's WCSLIB reads MJDREF keyword.
    - GAL_WCS_FLTERROR: Limit to identify a floating point error for WCS.
+   - gal_wcs_write: Write the given WCS into a FITS extension with no data.
    - gal_wcs_clean_errors: clean major WCS components from errors specified
        by the FITS keyword 'CRDER', or floating point errors.
    - gal_wcs_distortion_from_string: Return distortion string/name from ID.
diff --git a/bin/fits/keywords.c b/bin/fits/keywords.c
index 84335ee..a4e47f3 100644
--- a/bin/fits/keywords.c
+++ b/bin/fits/keywords.c
@@ -441,15 +441,31 @@ static void
 keywords_distortion_wcs(struct fitsparams *p)
 {
   int nwcs;
-  struct wcsprm *inwcs;
+  size_t ndim, *insize;
+  gal_data_t *data=NULL;
   char *suffix, *output;
-  gal_data_t *data=gal_fits_img_read(p->filename, p->cp.hdu,
-                                     p->cp.minmapsize, p->cp.quietmmap);
+  struct wcsprm *inwcs, *outwcs;
+
+  /* If the extension has any data, read it, otherwise just make an empty
+     array. */
+  if(gal_fits_hdu_format(p->filename, p->cp.hdu)==IMAGE_HDU)
+    {
+      /* Read the size of the dataset (we don't need the actual size!). */
+      insize=gal_fits_img_info_dim(p->filename, p->cp.hdu, &ndim);
+      free(insize);
+
+      /* If the number of dimensions is two, then read the dataset,
+         otherwise, ignore it. */
+      if(ndim==2)
+        data=gal_fits_img_read(p->filename, p->cp.hdu, p->cp.minmapsize,
+                               p->cp.quietmmap);
+    }
 
   /* Read the input WCS structure and convert it to the desired output
      distortion. */
   inwcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &nwcs);
-  data->wcs=gal_wcs_distortion_convert(inwcs, p->distortionid, data->dsize);
+  outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid,
+                                    data?data->dsize:NULL);
 
   /* Set the output filename. */
   if(p->cp.output)
@@ -460,12 +476,26 @@ keywords_distortion_wcs(struct fitsparams *p)
         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
       output=gal_checkset_automatic_output(&p->cp, p->filename, suffix);
     }
+  gal_checkset_writable_remove(output, 0, p->cp.dontdelete);
 
   /* Write the output file. */
-  gal_fits_img_write(data, output, NULL, PROGRAM_NAME);
+  if(data)
+    {
+      /* Add the output WCS to the dataset and write it. */
+      data->wcs=outwcs;
+      gal_fits_img_write(data, output, NULL, PROGRAM_NAME);
+
+      /* Clean up, but remove the pointer first (so it doesn't free it
+         here). */
+      data->wcs=NULL;
+      gal_data_free(data);
+    }
+  else
+    gal_wcs_write(outwcs, output, NULL, PROGRAM_NAME);
 
   /* Clean up. */
   wcsfree(inwcs);
+  wcsfree(outwcs);
   if(output!=p->cp.output) free(output);
 }
 
diff --git a/bin/fits/ui.c b/bin/fits/ui.c
index 8957418..bd8ea0c 100644
--- a/bin/fits/ui.c
+++ b/bin/fits/ui.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <string.h>
 
+#include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 
 #include <gnuastro-internal/options.h>
diff --git a/configure.ac b/configure.ac
index 8aaf7ee..a04cc50 100644
--- a/configure.ac
+++ b/configure.ac
@@ -444,6 +444,23 @@ AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_WCSLIB_VERSION], 
[$has_wcslib_version],
 AC_SUBST(HAVE_WCSLIB_VERSION, [$has_wcslib_version])
 
 
+# If the WCS library supports distortion
+AC_CHECK_HEADER([wcslib/dis.h], [has_wcslib_dis_h=1],
+                [has_wcslib_dis_h=0; anywarnings=yes])
+AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_WCSLIB_DIS_H], [$has_wcslib_dis_h],
+                   [WCSLIB has distortion header in dis.h])
+AC_SUBST(HAVE_WCSLIB_DIS_H, [$has_wcslib_dis_h])
+AM_CONDITIONAL([COND_HASWCSDIS_H], [test "x$has_wcslib_dis_h" = "x1"])
+
+
+# If the WCS library has the 'mjdref' element.
+AC_CHECK_MEMBER([struct wcsprm.mjdref], [has_wcslib_mjdref=1],
+                [has_wcslib_mjdref=0], [[#include <wcslib/wcs.h>]])
+AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_WCSLIB_MJDREF], [$has_wcslib_mjdref],
+                   [WCSLIB comes with wcsprm.mjdref])
+AC_SUBST(HAVE_WCSLIB_MJDREF, [$has_wcslib_mjdref])
+
+
 # If the pthreads library has 'pthread_barrier_wait'.
 AC_CHECK_LIB([pthread], [pthread_barrier_wait], [has_pthread_barrier=1],
              [has_pthread_barrier=0])
@@ -978,6 +995,15 @@ AS_IF([test x$enable_guide_message = xyes],
                AS_ECHO(["    released in October 2015)."])
                AS_ECHO([]) ])
 
+        AS_IF([test "x$has_wcslib_dis_h" = "x0"],
+              [dependency_notice=yes
+               AS_ECHO(["  - WCSLIB 
(https://www.atnf.csiro.au/people/mcalabre/WCS) version"])
+               AS_ECHO(["    on this system doesn't support distortions (i.e., 
it doesn't"])
+               AS_ECHO(["    have 'dis.h'. This build won't crash but Gnuastro 
will not be able"])
+               AS_ECHO(["    to do distortion-related operations. If you don't 
need such"])
+               AS_ECHO(["    operations you can ignore this warning."])
+               AS_ECHO([]) ])
+
         AS_IF([test "x$has_libjpeg" = "xno"],
               [dependency_notice=yes
                AS_ECHO(["  - libjpeg (http://ijg.org), could not be linked 
with in your library"])
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8a9d1e8..4b22347 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -18475,6 +18475,14 @@ However, only more recent versions of WCSLIB also 
provide its version number.
 If the WCSLIB that is installed on the system provides its version (through 
the possibly existing @code{wcslib_version} function), this macro will have a 
value of one, otherwise it will have a value of zero.
 @end deffn
 
+@deffn Macro GAL_CONFIG_HAVE_WCSLIB_DIS_H
+This macro has a value of 1 if the host's WCSLIB has the @file{wcslib/dis.h} 
header for distortion-related operations.
+@end deffn
+
+@deffn Macro GAL_CONFIG_HAVE_WCSLIB_MJDREF
+This macro has a value of 1 if the host's WCSLIB reads and stores the 
@file{MJDREF} FITS header keyword as part of its core @code{wcsprm} structure.
+@end deffn
+
 @deffn Macro GAL_CONFIG_HAVE_PTHREAD_BARRIER
 The POSIX threads standard define barriers as an optional requirement.
 Therefore, some operating systems choose to not include it.
@@ -22184,6 +22192,11 @@ number of coordinate representations found into the 
space that @code{nwcs}
 points to. Please see @code{gal_wcs_read_fitsptr} for more.
 @end deftypefun
 
+@deftypefun gal_wcs_write (struct wcsprm @code{*wcs}, char @code{*filename}, 
gal_fits_list_key_t @code{*headers}, char @code{*program_string})
+Write the given WCS structure into the second extension of an empty FITS 
header.
+The first/primary extension will be empty like the default format of all 
Gnuastro outputs.
+@end deftypefun
+
 @deftypefun {struct wcsprm *} gal_wcs_copy (struct wcsprm @code{*wcs})
 Return a fully allocated (independent) copy of @code{wcs}.
 @end deftypefun
diff --git a/lib/Makefile.am b/lib/Makefile.am
index e1b34e8..cbd4f0b 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -50,19 +50,27 @@ libgnuastro_la_LIBADD = 
$(top_builddir)/bootstrapped/lib/libgnu.la
 
 
 
+# Conditional compilation
+if COND_HASWCSDIS_H
+  MAYBE_WCSDISTORTION = wcsdistortion.c
+endif
+
+
+
+
 
 # Specify the library .c files
-libgnuastro_la_SOURCES = arithmetic.c arithmetic-and.c arithmetic-bitand.c \
-  arithmetic-bitlsh.c arithmetic-bitor.c arithmetic-bitrsh.c               \
-  arithmetic-bitxor.c arithmetic-divide.c arithmetic-eq.c arithmetic-ge.c  \
-  arithmetic-gt.c arithmetic-le.c arithmetic-lt.c arithmetic-minus.c       \
-  arithmetic-modulo.c arithmetic-multiply.c arithmetic-ne.c                \
-  arithmetic-or.c arithmetic-plus.c array.c binary.c blank.c box.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 speclines.c        \
-  statistics.c table.c tableintern.c threads.c tiff.c tile.c               \
-  tile-internal.c timing.c txt.c type.c units.c wcs.c wcsdistortion.c
+libgnuastro_la_SOURCES = $(MAYBE_WCSDISTORTION) arithmetic.c \
+  arithmetic-and.c arithmetic-bitand.c arithmetic-bitlsh.c \
+  arithmetic-bitor.c arithmetic-bitrsh.c arithmetic-bitxor.c \
+  arithmetic-divide.c arithmetic-eq.c arithmetic-ge.c arithmetic-gt.c \
+  arithmetic-le.c arithmetic-lt.c arithmetic-minus.c arithmetic-modulo.c \
+  arithmetic-multiply.c arithmetic-ne.c arithmetic-or.c arithmetic-plus.c \
+  array.c binary.c blank.c box.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 speclines.c \
+  statistics.c table.c tableintern.c threads.c tiff.c tile.c \
+  tile-internal.c timing.c txt.c type.c units.c wcs.c
 
 
 
@@ -72,17 +80,17 @@ libgnuastro_la_SOURCES = arithmetic.c arithmetic-and.c 
arithmetic-bitand.c \
 # in the $(headersdir) directory. Some of the header files don't need to be
 # installed.
 headersdir=$(top_srcdir)/lib/gnuastro
-pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h          \
-  $(headersdir)/array.h $(headersdir)/binary.h $(headersdir)/blank.h       \
-  $(headersdir)/box.h $(headersdir)/convolve.h $(headersdir)/cosmology.h   \
-  $(headersdir)/data.h $(headersdir)/dimension.h $(headersdir)/eps.h       \
-  $(headersdir)/fits.h $(headersdir)/git.h $(headersdir)/interpolate.h     \
-  $(headersdir)/jpeg.h $(headersdir)/label.h $(headersdir)/list.h          \
-  $(headersdir)/match.h $(headersdir)/pdf.h $(headersdir)/permutation.h    \
-  $(headersdir)/pointer.h $(headersdir)/polygon.h $(headersdir)/qsort.h    \
-  $(headersdir)/speclines.h $(headersdir)/statistics.h                     \
-  $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/tiff.h       \
-  $(headersdir)/tile.h $(headersdir)/txt.h $(headersdir)/type.h            \
+pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h \
+  $(headersdir)/array.h $(headersdir)/binary.h $(headersdir)/blank.h \
+  $(headersdir)/box.h $(headersdir)/convolve.h $(headersdir)/cosmology.h \
+  $(headersdir)/data.h $(headersdir)/dimension.h $(headersdir)/eps.h \
+  $(headersdir)/fits.h $(headersdir)/git.h $(headersdir)/interpolate.h \
+  $(headersdir)/jpeg.h $(headersdir)/label.h $(headersdir)/list.h \
+  $(headersdir)/match.h $(headersdir)/pdf.h $(headersdir)/permutation.h \
+  $(headersdir)/pointer.h $(headersdir)/polygon.h $(headersdir)/qsort.h \
+  $(headersdir)/speclines.h $(headersdir)/statistics.h \
+  $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/tiff.h \
+  $(headersdir)/tile.h $(headersdir)/txt.h $(headersdir)/type.h \
   $(headersdir)/units.h $(headersdir)/wcs.h
 
 
@@ -95,21 +103,21 @@ pkginclude_HEADERS = gnuastro/config.h 
$(headersdir)/arithmetic.h          \
 # will not distributed, so we need to explicitly tell Automake to
 # distribute them here.
 internaldir=$(top_srcdir)/lib/gnuastro-internal
-EXTRA_DIST = gnuastro.pc.in $(headersdir)/README $(internaldir)/README  \
-  $(internaldir)/arithmetic-and.h $(internaldir)/arithmetic-binary.h    \
+EXTRA_DIST = gnuastro.pc.in $(headersdir)/README $(internaldir)/README \
+  $(internaldir)/arithmetic-and.h $(internaldir)/arithmetic-binary.h \
   $(internaldir)/arithmetic-bitand.h $(internaldir)/arithmetic-bitlsh.h \
-  $(internaldir)/arithmetic-bitor.h $(internaldir)/arithmetic-bitrsh.h  \
+  $(internaldir)/arithmetic-bitor.h $(internaldir)/arithmetic-bitrsh.h \
   $(internaldir)/arithmetic-bitxor.h $(internaldir)/arithmetic-divide.h \
-  $(internaldir)/arithmetic-eq.h $(internaldir)/arithmetic-ge.h         \
-  $(internaldir)/arithmetic-gt.h $(internaldir)/arithmetic-internal.h   \
-  $(internaldir)/arithmetic-le.h $(internaldir)/arithmetic-lt.h         \
-  $(internaldir)/arithmetic-minus.h $(internaldir)/arithmetic-modulo.h  \
-  $(internaldir)/arithmetic-multiply.h $(internaldir)/arithmetic-ne.h   \
-  $(internaldir)/arithmetic-or.h $(internaldir)/arithmetic-plus.h       \
-  $(internaldir)/checkset.h $(internaldir)/commonopts.h                 \
-  $(internaldir)/config.h.in $(internaldir)/fixedstringmacros.h         \
-  $(internaldir)/options.h $(internaldir)/tableintern.h                 \
-  $(internaldir)/tile-internal.h $(internaldir)/timing.h                \
+  $(internaldir)/arithmetic-eq.h $(internaldir)/arithmetic-ge.h \
+  $(internaldir)/arithmetic-gt.h $(internaldir)/arithmetic-internal.h \
+  $(internaldir)/arithmetic-le.h $(internaldir)/arithmetic-lt.h \
+  $(internaldir)/arithmetic-minus.h $(internaldir)/arithmetic-modulo.h \
+  $(internaldir)/arithmetic-multiply.h $(internaldir)/arithmetic-ne.h \
+  $(internaldir)/arithmetic-or.h $(internaldir)/arithmetic-plus.h \
+  $(internaldir)/checkset.h $(internaldir)/commonopts.h \
+  $(internaldir)/config.h.in $(internaldir)/fixedstringmacros.h \
+  $(internaldir)/options.h $(internaldir)/tableintern.h \
+  $(internaldir)/tile-internal.h $(internaldir)/timing.h \
   $(internaldir)/wcsdistortion.h
 
 
@@ -131,13 +139,15 @@ CLEANFILES = gnuastro.pc gnuastro/config.h
 gnuastro/config.h: Makefile $(internaldir)/config.h.in
        rm -f $@ $@.tmp
        $(MKDIR_P) gnuastro
-       $(SED) -e 's|@VERSION[@]|$(VERSION)|g'                            \
-              -e 's|@HAVE_LIBGIT2[@]|$(HAVE_LIBGIT2)|g'                  \
-              -e 's|@HAVE_WCSLIB_VERSION[@]|$(HAVE_WCSLIB_VERSION)|g'    \
-              -e 's|@HAVE_PTHREAD_BARRIER[@]|$(HAVE_PTHREAD_BARRIER)|g'  \
-              -e 's|@SIZEOF_LONG[@]|$(SIZEOF_LONG)|g'                    \
-              -e 's|@SIZEOF_SIZE_T[@]|$(SIZEOF_SIZE_T)|g'                \
-              -e 's|@RESTRICT_REPLACEMENT[@]|$(RESTRICT_REPLACEMENT)|g'  \
+       $(SED) -e 's|@VERSION[@]|$(VERSION)|g' \
+              -e 's|@HAVE_LIBGIT2[@]|$(HAVE_LIBGIT2)|g' \
+              -e 's|@HAVE_WCSLIB_VERSION[@]|$(HAVE_WCSLIB_VERSION)|g' \
+              -e 's|@HAVE_WCSLIB_DIS_H[@]|$(HAVE_WCSLIB_DIS_H)|g' \
+              -e 's|@HAVE_WCSLIB_MJDREF[@]|$(HAVE_WCSLIB_MJDREF)|g' \
+              -e 's|@HAVE_PTHREAD_BARRIER[@]|$(HAVE_PTHREAD_BARRIER)|g' \
+              -e 's|@SIZEOF_LONG[@]|$(SIZEOF_LONG)|g' \
+              -e 's|@SIZEOF_SIZE_T[@]|$(SIZEOF_SIZE_T)|g' \
+              -e 's|@RESTRICT_REPLACEMENT[@]|$(RESTRICT_REPLACEMENT)|g' \
                $(internaldir)/config.h.in >> $@.tmp
        chmod a-w $@.tmp
        mv $@.tmp $@
diff --git a/lib/gnuastro-internal/config.h.in 
b/lib/gnuastro-internal/config.h.in
index c3f194e..253b9a2 100644
--- a/lib/gnuastro-internal/config.h.in
+++ b/lib/gnuastro-internal/config.h.in
@@ -40,6 +40,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #define GAL_CONFIG_HAVE_FITS_IS_REENTRANT   @HAVE_FITS_IS_REENTRANT@
 #define GAL_CONFIG_HAVE_WCSLIB_VERSION      @HAVE_WCSLIB_VERSION@
+#define GAL_CONFIG_HAVE_WCSLIB_DIS_H        @HAVE_WCSLIB_DIS_H@
+#define GAL_CONFIG_HAVE_WCSLIB_MJDREF       @HAVE_WCSLIB_MJDREF@
 
 #define GAL_CONFIG_HAVE_PTHREAD_BARRIER     @HAVE_PTHREAD_BARRIER@
 
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index e66548d..5821e3c 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <wcslib/wcs.h>
 
 #include <gnuastro/data.h>
+#include <gnuastro/fits.h>
 
 /* Assumed floating point error in the WCS-related functionality. */
 #define GAL_WCS_FLTERROR 1e-12
@@ -87,6 +88,15 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
 
 
 
+/*************************************************************
+ ***********               Write WCS               ***********
+ *************************************************************/
+void
+gal_wcs_write(struct wcsprm *wcs, char *filename,
+              gal_fits_list_key_t *headers, char *program_string);
+
+
+
 
 /*************************************************************
  ***********              Distortions              ***********
diff --git a/lib/wcs.c b/lib/wcs.c
index b3f0c5a..a19d00d 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -31,8 +31,6 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <unistd.h>
 #include <assert.h>
 
-#include <wcslib/dis.h>
-
 #include <gsl/gsl_linalg.h>
 #include <wcslib/wcsmath.h>
 
@@ -43,10 +41,10 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/dimension.h>
 #include <gnuastro/permutation.h>
 
+#if GAL_CONFIG_HAVE_WCSLIB_DIS_H
+#include <wcslib/dis.h>
 #include <gnuastro-internal/wcsdistortion.h>
-
-
-
+#endif
 
 
 
@@ -274,6 +272,83 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
 
 
 /*************************************************************
+ ***********               Write WCS               ***********
+ *************************************************************/
+void
+gal_wcs_write(struct wcsprm *wcs, char *filename,
+              gal_fits_list_key_t *headers, char *program_string)
+{
+  char *wcsstr;
+  size_t ndim=0;
+  fitsfile *fptr;
+  long *naxes=NULL;
+  int status=0, nkeyrec;
+
+  /* Small sanity checks */
+  if(wcs==NULL)
+    error(EXIT_FAILURE, 0, "%s: input WCS is NULL", __func__);
+  if( gal_fits_name_is_fits(filename)==0 )
+    error(EXIT_FAILURE, 0, "%s: not a FITS suffix", filename);
+
+  /* Open the file for writing */
+  fptr=gal_fits_open_to_write(filename);
+
+  /* Create the FITS file. */
+  fits_create_img(fptr, gal_fits_type_to_bitpix(GAL_TYPE_UINT8),
+                  ndim, naxes, &status);
+  gal_fits_io_error(status, NULL);
+
+  /* Remove the two comment lines put by CFITSIO. Note that in some cases,
+     it might not exist. When this happens, the status value will be
+     non-zero. We don't care about this error, so to be safe, we will just
+     reset the status variable after these calls. */
+  fits_delete_key(fptr, "COMMENT", &status);
+  fits_delete_key(fptr, "COMMENT", &status);
+  status=0;
+
+  /* Decompose the 'PCi_j' matrix and 'CDELTi' vector. */
+  gal_wcs_decompose_pc_cdelt(wcs);
+
+  /* Clean up small errors in the PC matrix and CDELT values. */
+  gal_wcs_clean_errors(wcs);
+
+  /* Convert the WCS information to text. */
+  status=wcshdo(WCSHDO_safe, wcs, &nkeyrec, &wcsstr);
+  if(status)
+    error(0, 0, "%s: WARNING: WCSLIB error, no WCS in output.\n"
+          "wcshdu ERROR %d: %s", __func__, status,
+          wcs_errmsg[status]);
+  else
+    gal_fits_key_write_wcsstr(fptr, wcsstr, nkeyrec);
+  status=0;
+
+  /* Write all the headers and the version information. */
+  gal_fits_key_write_version_in_ptr(&headers, program_string, fptr);
+
+  /* Close the FITS file. */
+  fits_close_file(fptr, &status);
+  gal_fits_io_error(status, NULL);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************
  ***********              Distortions              ***********
  *************************************************************/
 int
@@ -320,7 +395,7 @@ gal_wcs_distortion_to_string(int distortion)
   error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
         "problem. Control should not reach the end of this function",
         __func__, PACKAGE_BUGREPORT);
-  return GAL_WCS_DISTORTION_INVALID;
+  return NULL;
 }
 
 
@@ -338,6 +413,7 @@ gal_wcs_distortion_to_string(int distortion)
 int
 gal_wcs_distortion_identify(struct wcsprm *wcs)
 {
+#if GAL_CONFIG_HAVE_WCSLIB_DIS_H
   struct disprm *dispre=NULL;
   struct disprm *disseq=NULL;
 
@@ -395,6 +471,14 @@ gal_wcs_distortion_identify(struct wcsprm *wcs)
   error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to fix "
         "the problem. Control should not reach the end of this function",
         __func__, PACKAGE_BUGREPORT);
+#else
+  /* The 'wcslib/dis.h' isn't present. */
+  error(EXIT_FAILURE, 0, "%s: the installed version of WCSLIB on this "
+        "system doesn't have the 'dis.h' header! Thus Gnuastro can't do "
+        "distortion-related operations on the world coordinate system "
+        "(WCS). To use these features, please upgrade your version "
+        "of WCSLIB and re-build Gnuastro", __func__);
+#endif
   return GAL_WCS_DISTORTION_INVALID;
 }
 
@@ -421,6 +505,7 @@ struct wcsprm *
 gal_wcs_distortion_convert(struct wcsprm *inwcs, int outdisptype,
                            size_t *fitsize)
 {
+#if GAL_CONFIG_HAVE_WCSLIB_DIS_H
   struct wcsprm *outwcs=NULL;
   int indisptype=gal_wcs_distortion_identify(inwcs);
 
@@ -458,6 +543,9 @@ gal_wcs_distortion_convert(struct wcsprm *inwcs, int 
outdisptype,
         switch(outdisptype)
           {
           case GAL_WCS_DISTORTION_SIP:
+            if(fitsize==NULL)
+              error(EXIT_FAILURE, 0, "%s: the size array is necessary "
+                    "for this conversion", __func__);
             outwcs=gal_wcsdistortion_tpv_to_sip(inwcs, fitsize);
             break;
           default:
@@ -487,6 +575,15 @@ gal_wcs_distortion_convert(struct wcsprm *inwcs, int 
outdisptype,
 
   /* Return the converted WCS. */
   return outwcs;
+#else
+  /* The 'wcslib/dis.h' isn't present. */
+  error(EXIT_FAILURE, 0, "%s: the installed version of WCSLIB on this "
+        "system doesn't have the 'dis.h' header! Thus Gnuastro can't do "
+        "distortion-related operations on the world coordinate system "
+        "(WCS). To use these features, please upgrade your version "
+        "of WCSLIB and re-build Gnuastro", __func__);
+  return NULL;
+#endif
 }
 
 
diff --git a/lib/wcsdistortion.c b/lib/wcsdistortion.c
index 2b247b0..1d2d5b9 100644
--- a/lib/wcsdistortion.c
+++ b/lib/wcsdistortion.c
@@ -81,8 +81,8 @@ wcsdistortion_get_tpvparams(struct wcsprm *wcs, double 
cd[2][2],
   if(wcs==NULL)
     error(EXIT_FAILURE, 0, "%s: input WCS structure is NULL", __func__);
 
-  for(i=0,k=0; i<2; i++)
-    for(j=0; j<2; j++)
+  for(i=0,k=0; i<2; ++i)
+    for(j=0; j<2; ++j)
       {
         /* If a value is present store it, else store 0.*/
         if((wcs->cd)[i] != 0)   cd[i][j]=(wcs->cd)[k++];
@@ -93,52 +93,52 @@ wcsdistortion_get_tpvparams(struct wcsprm *wcs, double 
cd[2][2],
         */
       }
 
-  for(j=0; j < wcs->npvmax; j++)
-  {
-    if(wcs->pv[j].i == 1)
-      {
-        /*pv_m is used to check the index in the header.*/
-        pv_m=wcs->pv[j].m;
+  for(j=0; j < wcs->npvmax; ++j)
+    {
+      if(wcs->pv[j].i == 1)
+        {
+          /*pv_m is used to check the index in the header.*/
+          pv_m=wcs->pv[j].m;
 
-        /* `index` is the index of the pv* array.*/
-        index = pv_m;
+          /* `index` is the index of the pv* array.*/
+          index = pv_m;
 
-        if( wcs->pv[pv_m].value != 0 && index == pv_m )
-          pv1[index]=wcs->pv[j].value;
-        else
-          pv1[index]=0;
+          if( wcs->pv[pv_m].value != 0 && index == pv_m )
+            pv1[index]=wcs->pv[j].value;
+          else
+            pv1[index]=0;
 
-        /* For a check
-        printf("PV1_%d\t%.10f\n", pv_m, pv1[k]);
-        */
-      }
-    else if(wcs->pv[j].i == 2)
-      {
-        /*pv_m is used to check the index in the header.*/
-        pv_m=wcs->pv[j].m;
+          /* For a check
+          printf("PV1_%d\t%.10f\n", pv_m, pv1[k]);
+          */
+        }
+      else if(wcs->pv[j].i == 2)
+        {
+          /*pv_m is used to check the index in the header.*/
+          pv_m=wcs->pv[j].m;
 
-        /* `index` is the index of the pv* array.*/
-        index = pv_m;
+          /* `index` is the index of the pv* array.*/
+          index = pv_m;
 
-        if( wcs->pv[pv_m].value != 0 && index == pv_m )
-          pv2[index]=wcs->pv[j].value;
-        else
-          pv2[index]=0;
+          if( wcs->pv[pv_m].value != 0 && index == pv_m )
+            pv2[index]=wcs->pv[j].value;
+          else
+            pv2[index]=0;
 
-        /* For a check
-        printf("PV2_%d\t%.10f\n", pv_m, pv2[k]);
-        */
-      }
-    else
-      {
-        error(EXIT_FAILURE, 0, "%s: No such axis present!", __func__);
-        break;
-      }
-  }
+          /* For a check
+          printf("PV2_%d\t%.10f\n", pv_m, pv2[k]);
+          */
+        }
+      else
+        {
+          error(EXIT_FAILURE, 0, "%s: No such axis present!", __func__);
+          break;
+        }
+    }
 
   /* To check a full array:
-  for(i = 0; i < 40; i++)
-    printf("PV2_%ld\t%.10f\n", i, pv2[i]);
+     for(i = 0; i < 40; ++i)
+     printf("PV2_%ld\t%.10f\n", i, pv2[i]);
   */
 
 }
@@ -154,7 +154,6 @@ wcsdistortion_get_sipparams(struct wcsprm *wcs, double 
cd[2][2],
                             double a_coeff[5][5], double b_coeff[5][5])
 {
   const char *cp;
-  double keyval=0;
   double *temp_cd=NULL;
   struct dpkey  *keyp=NULL;
   struct disprm *dispre=NULL;
@@ -176,8 +175,8 @@ wcsdistortion_get_sipparams(struct wcsprm *wcs, double 
cd[2][2],
      CD matix is extracted using the `gal_wcs_wrap_matrix` as a single
      allocated array (temp_cd[]), that is then used to fill cd[][]. */
   temp_cd=gal_wcs_warp_matrix(wcs);
-  for(i=0,k=0; i<2; i++)
-    for(j=0; j<2; j++)
+  for(i=0,k=0; i<2; ++i)
+    for(j=0; j<2; ++j)
       {
         /* If a value is present store it, else store 0.*/
         if(temp_cd[k] != 0)   cd[i][j]=temp_cd[k++];
@@ -190,13 +189,13 @@ wcsdistortion_get_sipparams(struct wcsprm *wcs, double 
cd[2][2],
 
   /* Extract the SIP values from the strings and numbers inside the
      wcsprm. */
-  for(i=0; i<dispre->ndp; i++, keyp++)
+  for(i=0; i<dispre->ndp; ++i, ++keyp)
     {
       /* For axis1. */
       if (keyp->j == 1)
         {
           /* Ignore any missing keyvalues. */
-          if ((keyval = dpkeyd(keyp)) == 0.0) continue;
+          if ( keyp->field == NULL ) continue;
 
           cp = strchr(keyp->field, '.') + 1;
           if (strncmp(cp, "SIP.FWD.", 8) != 0) continue;
@@ -214,10 +213,10 @@ wcsdistortion_get_sipparams(struct wcsprm *wcs, double 
cd[2][2],
       else if (keyp->j == 2)
         {
           /* Ignore any missing keyvalues. */
-          if ((keyval = dpkeyd(keyp)) == 0.0) continue;
+          if ( keyp->field == NULL ) continue;
 
           cp = strchr(keyp->field, '.') + 1;
-          if (strncmp(cp, "SIP.FWD.", 8) != 0) continue;
+          if (strncmp(cp, "SIP.REV.", 8) != 0) continue;
           cp += 8;
 
           sscanf(cp, "%ld_%ld", &m, &n);
@@ -247,15 +246,18 @@ wcsdistortion_get_sipcoeff(struct wcsprm *wcs, size_t 
*a_order,
   size_t m, n;
   double val=0;
   size_t i, j, k;
-  double cd[2][2]={0};
-  double tpvu[8][8]={0}, tpvv[8][8]={0};
+  double cd[2][2];
+  double tpvu[8][8], tpvv[8][8];
 
+  /* Initialise the 2d matrices. */
+  for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd[i][j]=0;
+  for(i=0; i<8; ++i) for(j=0; j<8; ++j) { tpvu[i][j]=0; tpvv[i][j]=0; }
 
+  /* Calculate the TPV elements. */
   wcsdistortion_calc_tpveq(wcs, cd, tpvu, tpvv);
 
-
-  for(i=0,k=0; i<2; i++)
-    for(j=0; j<2; j++)
+  for(i=0,k=0; i<2; ++i)
+    for(j=0; j<2; ++j)
       {
         /* If a value is present store it, else store 0.*/
         if((wcs->cd)[i] != 0)   cd[i][j]=(wcs->cd)[k++];
@@ -266,8 +268,8 @@ wcsdistortion_get_sipcoeff(struct wcsprm *wcs, size_t 
*a_order,
         */
       }
 
-  for(m=0; m<=4; m++)
-    for(n=0; n<=4; n++)
+  for(m=0; m<=4; ++m)
+    for(n=0; n<=4; ++n)
       {
         /*For axis = 1*/
         val=wcsdistortion_calcsip(1, m, n, tpvu, tpvv);
@@ -291,10 +293,10 @@ wcsdistortion_get_sipcoeff(struct wcsprm *wcs, size_t 
*a_order,
       }
 
   /*For a check.
-  for(m=0;m<=4;m++)
-    for(n=0;n<=4;n++)
-      printf("A_%ld_%ld = %.10E \t B_%ld_%ld = %.10E\n",
-              m, n, a_coeff[m][n], m, n, b_coeff[m][n]);
+  for(m=0;m<=4;++m)
+  for(n=0;n<=4;++n)
+  printf("A_%ld_%ld = %.10E \t B_%ld_%ld = %.10E\n",
+         m, n, a_coeff[m][n], m, n, b_coeff[m][n]);
   */
 }
 
@@ -512,8 +514,8 @@ wcsdistortion_intermidate_tpveq(double cd[2][2], double 
*pv1, double *pv2,
   /* For a check:
   {
     size i, j;
-    for(i=0; i<=4; i++)
-      for(j=0;j<=4;j++)
+    for(i=0; i<=4; ++i)
+      for(j=0;j<=4;++j)
        {
          printf("k%ld_%ld \t %.10E\n",   i, j, k[i][j]);
          printf("l%ld_%ld \t %.10E\n\n", i, j, l[i][j]);
@@ -554,7 +556,7 @@ wcsdistortion_intermidate_sipeq(double cd[2][2], double 
cd_inv[2][2],
     Grillmair, Carl; Levitan, David; Sesar, Branimir, 2012, in Software and
     Cyberinfrastructure for Astronomy II.
     Proceedings of the SPIE, Volume 8451, article id. 84511M.
-    */
+  */
 
   /* pvi_0 */
   pv1[0]=a_coeff[0][0]*cd[0][0] + b_coeff[0][0]*cd[0][1];
@@ -906,7 +908,7 @@ wcsdistortion_intermidate_sipeq(double cd[2][2], double 
cd_inv[2][2],
   /* For a check:
   {
     size_t j;
-    for(j=0; j<=16 ;j++)
+    for(j=0; j<=16 ;++j)
       {
         if (j == 3 || j == 11) continue;
         printf("pv1_%ld \t %.12E\n", j, pv1[j]);
@@ -944,18 +946,20 @@ wcsdistortion_calc_tpveq(struct wcsprm *wcs, double 
cd[2][2],
                          double tpvu[8][8], double tpvv[8][8])
 {
   /* tpvu and tpvv are u-v distortion equations in TPV convention. */
-  size_t i,j;
+  size_t i=0,j=0;
+  double cd_inv[2][2];
   double determinant=0;
-  double cd_inv[2][2]={0};
+  double k[5][5], l[5][5];
   double pv1[17]={0}, pv2[17]={0};
-  double k[5][5]={0}, l[5][5]={0};
 
+  /* Initialise the 2d matrices. */
+  for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd_inv[i][j]=0;
+  for(i=0; i<5; ++i) for(j=0; j<5; ++j) { k[i][j]=0; l[i][j]=0; }
 
+  /* Estimate the parameters. */
   wcsdistortion_get_tpvparams(wcs, cd, pv1, pv2);
-
   wcsdistortion_intermidate_tpveq(cd, pv1, pv2, k, l);
 
-
   /* We will find matrix tpvu and tpvv by finding inverse of
      CD matrix and multiplying with tpv* matrix.
      For inverse of a 2x2 matrix we use the below trasformations:
@@ -964,8 +968,6 @@ wcsdistortion_calc_tpveq(struct wcsprm *wcs, double 
cd[2][2],
       where |A| is the determinant of the matrix which is calculated by:
                           |A| = a*d-b*c.
     */
-
-
   determinant = cd[0][0]*cd[1][1] - cd[0][1]*cd[1][0];
 
   /* Inverse matrix */
@@ -989,8 +991,8 @@ wcsdistortion_calc_tpveq(struct wcsprm *wcs, double 
cd[2][2],
     to multiplycation with cd_inv.
   */
 
-  for(i=0; i<=4; i++)
-    for(j=0; j<=4; j++)
+  for(i=0; i<=4; ++i)
+    for(j=0; j<=4; ++j)
       {
         tpvu[i][j]=cd_inv[0][0]*k[i][j]+cd_inv[0][1]*l[i][j];
         tpvv[i][j]=cd_inv[1][0]*k[i][j]+cd_inv[1][1]*l[i][j];
@@ -999,7 +1001,6 @@ wcsdistortion_calc_tpveq(struct wcsprm *wcs, double 
cd[2][2],
         printf("%.10E, %.10E\n", tpvu[i][j], tpvv[i][j]);
         */
       }
-
 }
 
 
@@ -1013,10 +1014,14 @@ wcsdistortion_calc_sipeq(struct wcsprm *wcs, double 
cd[2][2],
                          double *pv1, double *pv2)
 {
   /* tpvu and tpvv are u-v distortion equations in TPV convention. */
+  size_t i=0, j=0;
+  double cd_inv[2][2];
   double determinant=0;
-  double cd_inv[2][2]={0};
-  double a_coeff[5][5]={0}, b_coeff[5][5]={0};
+  double a_coeff[5][5], b_coeff[5][5];
 
+  /* Initialise the 2d matrices. */
+  for(i=0;i<2;++i) for(j=0;j<2;++j) cd_inv[i][j]=0;
+  for(i=0;i<5;++i) for(j=0;j<5;++j) {a_coeff[i][j]=0; b_coeff[i][j]=0;}
 
   /* We will find matrix tpvu and tpvv by finding inverse of
      CD matrix and multiplying with tpv* matrix.
@@ -1026,9 +1031,7 @@ wcsdistortion_calc_sipeq(struct wcsprm *wcs, double 
cd[2][2],
       where |A| is the determinant of the matrix which is calculated by:
                           |A| = a*d-b*c.
     */
-
   wcsdistortion_get_sipparams(wcs, cd, a_coeff, b_coeff);
-
   determinant = cd[0][0]*cd[1][1] - cd[0][1]*cd[1][0];
 
   /* Inverse matrix */
@@ -1037,10 +1040,9 @@ wcsdistortion_calc_sipeq(struct wcsprm *wcs, double 
cd[2][2],
   cd_inv[1][0]=(-1*cd[1][0])/determinant; /* c */
   cd_inv[1][1]=cd[0][0]/determinant;      /* d */
 
-
+  /* Find the parameters. */
   wcsdistortion_intermidate_sipeq( cd, cd_inv, a_coeff, b_coeff, pv1, pv2);
 
-
   /*For a check.
   printf("%.10lf\t%.10lf\t%.10lf\t%.10lf\n", cd_inv[0][0], cd_inv[0][1], \
                                              cd_inv[1][0], cd_inv[1][1]);
@@ -1059,16 +1061,13 @@ wcsdistortion_calcsip(size_t axis, size_t m, size_t n, 
double tpvu[8][8],
                       double tpvv[8][8])
 {
   double sip_coeff;
-  if(axis == 1)        sip_coeff=tpvu[m][n];
-  else if(axis == 2)   sip_coeff=tpvv[m][n];
-  else
-    error(EXIT_FAILURE, 0, "%s: axis does not exists! ",
-          __func__);
+  if(     axis == 1) sip_coeff=tpvu[m][n];
+  else if(axis == 2) sip_coeff=tpvv[m][n];
+  else error(EXIT_FAILURE, 0, "%s: axis does not exists! ", __func__);
 
-  if( (axis == 1) && (m == 1) && (n == 0) )
-        sip_coeff = sip_coeff - 1.0;
-  else if( (axis == 2) && (m == 0) && (n == 1) )
-        sip_coeff = sip_coeff - 1.0;
+  if(      (axis == 1) && (m == 1) && (n == 0) ) sip_coeff = sip_coeff - 1.0;
+  else if( (axis == 2) && (m == 0) && (n == 1) ) sip_coeff = sip_coeff - 1.0;
+  else error(EXIT_FAILURE, 0, "%s: axis does not exists! ", __func__);
 
   return sip_coeff;
 
@@ -1136,7 +1135,7 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
   /* Allocate updict and vpdict and initialize them. */
   updict=malloc((ap_order+1)*sizeof(*updict));
   vpdict=malloc((bp_order+1)*sizeof(*vpdict));
-  for(i=0; i<=wcsdistortion_max(ap_order, bp_order); i++)
+  for(i=0; i<=wcsdistortion_max(ap_order, bp_order); ++i)
     {
       updict[i]=malloc(tsize*sizeof(updict));
       vpdict[i]=malloc(tsize*sizeof(vpdict));
@@ -1145,7 +1144,7 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
   /* Allocate and initialize uprime and vprime. */
   uprime=malloc(tsize*sizeof(*uprime));
   vprime=malloc(tsize*sizeof(*vprime));
-  for(i=0; i<tsize; i++)
+  for(i=0; i<tsize; ++i)
     {
       uprime[i]=u[i];
       vprime[i]=v[i];
@@ -1154,12 +1153,12 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
   /* Allocate and initialize udict and vdict.*/
   udict=malloc((a_order+1)*sizeof(*udict));
   vdict=malloc((b_order+1)*sizeof(*vdict));
-  for(i=0; i<=wcsdistortion_max(a_order, b_order); i++)
+  for(i=0; i<=wcsdistortion_max(a_order, b_order); ++i)
     {
       udict[i]=malloc(tsize*sizeof(udict));
       vdict[i]=malloc(tsize*sizeof(vdict));
     }
-  for(i=0; i<tsize; i++)
+  for(i=0; i<tsize; ++i)
     {
       udict[0][i]=1;
       vdict[0][i]=1;
@@ -1168,73 +1167,65 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
 
   /* Fill the values from the in the dicts. The rows of the
       dicts act as a key to achieve a key-value functionality. */
-  for(i=1; i<=a_order; i++)
-    for(j=0; j<tsize; j++)
+  for(i=1; i<=a_order; ++i)
+    for(j=0; j<tsize; ++j)
       udict[i][j]=udict[i-1][j]*u[j];
 
-  for(i=1; i<=b_order; i++)
-    for(j=0; j<tsize; j++)
+  for(i=1; i<=b_order; ++i)
+    for(j=0; j<tsize; ++j)
       vdict[i][j]=vdict[i-1][j]*v[j];
 
 
   /* Populating forward coefficients on a grid. */
-  for(i=0; i<=a_order; i++)
-    for(j=0; j<=a_order-i; j++)
+  for(i=0; i<=a_order; ++i)
+    for(j=0; j<=a_order-i; ++j)
       {
-        for(k=0; k<tsize; k++)
+        for(k=0; k<tsize; ++k)
           uprime[k]+=a_coeff[i][j]*udict[i][k]*vdict[j][k];
 
         /* The number of parameters for AP_* coefficients. */
-        p_ap++;
+        ++p_ap;
       }
+  for(i=0; i<=b_order; ++i)
+    for(j=0; j<=b_order-i; ++j)
+      {
+        for(k=0; k<tsize; ++k)
+          vprime[k]+=b_coeff[i][j]*udict[i][k]*vdict[j][k];
 
-
-  for(i=0; i<=b_order; i++)
-    for(j=0; j<=b_order-i; j++)
-        {
-          for(k=0; k<tsize; k++)
-            vprime[k]+=b_coeff[i][j]*udict[i][k]*vdict[j][k];
-
-          /* The number of parameters for BP_* coefficients. */
-          p_bp++;
-        }
-
+        /* The number of parameters for BP_* coefficients. */
+        ++p_bp;
+      }
 
   /*For a check.
-  for(i=0; i<=a_order; i++)
-    for(j=0; j<tsize; j++)
+  for(i=0; i<=a_order; ++i)
+    for(j=0; j<tsize; ++j)
     {
       printf("udict[%ld][%ld] = %.8lf\n", i, j, uprime[j]);
       printf("u%ld = %.10E\n", i, u[i]);
     }
   */
 
-
   /* Now we have a grid populated with forward coeffiecients.  Now we fit a
      reverse polynomial through points using multiparameter linear least
      square fittings.*/
 
-
   /* Initialize dicts. */
-  for(i=0; i<tsize; i++)
+  for(i=0; i<tsize; ++i)
     {
       updict[0][i]=1;
       vpdict[0][i]=1;
     }
 
-
   /* Fill the values from the in the dicts. The rows of the
       dicts act as a key to achieve a key-value functionality. */
-  for(i=1; i<=ap_order; i++)
-    for(j=0; j<tsize; j++)
+  for(i=1; i<=ap_order; ++i)
+    for(j=0; j<tsize; ++j)
       updict[i][j]=updict[i-1][j]*uprime[j];
 
-  for(i=1; i<=bp_order; i++)
-    for(j=0; j<tsize; j++)
+  for(i=1; i<=bp_order; ++i)
+    for(j=0; j<tsize; ++j)
       vpdict[i][j]=vpdict[i-1][j]*vprime[j];
 
-
-
   /* Allocate memory for Multi-parameter Linear Regressions. */
   X_ap = gsl_matrix_alloc (tsize, p_ap);
   X_bp = gsl_matrix_alloc (tsize, p_bp);
@@ -1248,11 +1239,10 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
   cov_ap = gsl_matrix_alloc (p_ap, p_ap);
   cov_bp = gsl_matrix_alloc (p_bp, p_bp);
 
-
   /* Allocate and initialize udiff and vdiff. */
   udiff=malloc(tsize*sizeof(*udiff));
   vdiff=malloc(tsize*sizeof(*vdiff));
-  for(i=0; i<tsize; i++)
+  for(i=0; i<tsize; ++i)
     {
       udiff[i]=u[i]-uprime[i];
       vdiff[i]=v[i]-vprime[i];
@@ -1261,12 +1251,11 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
       gsl_vector_set (y_bp, i, vdiff[i]);
     }
 
-
   /* Filling up he X_ap matrix in column-wise order. */
-  for(i=0; i<=ap_order; i++)
-    for(j=0; j<=ap_order-i; j++)
+  for(i=0; i<=ap_order; ++i)
+    for(j=0; j<=ap_order-i; ++j)
       {
-        for(k=0; k<tsize; k++)
+        for(k=0; k<tsize; ++k)
           {
             gsl_matrix_set(X_ap, k, ij, updict[i][k]*vpdict[j][k]);
 
@@ -1274,16 +1263,15 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
             printf("x_ap[%ld] = %.8lf\n", ij, updict[i][k]*vpdict[j][k]);
             */
           }
-        ij++;
+        ++ij;
       }
 
-
   ij=0;
   /* Filling up he X_bp matrix in column-wise order. */
-  for(i=0; i<=bp_order; i++)
-    for(j=0; j<=bp_order-i; j++)
+  for(i=0; i<=bp_order; ++i)
+    for(j=0; j<=bp_order-i; ++j)
       {
-        for(k=0; k<tsize; k++)
+        for(k=0; k<tsize; ++k)
           {
             gsl_matrix_set(X_bp, k, ij, updict[i][k]*vpdict[j][k]);
 
@@ -1291,7 +1279,7 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
             printf("x_bp[%ld] = %.8lf\n", ij, updict[i][k]*vpdict[j][k]);
             */
           }
-        ij++;
+        ++ij;
       }
 
   /* Initialize and do the linear least square fitting. */
@@ -1302,8 +1290,8 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
 
   /* Extract the reverse coefficients. */
   p_ap=0;
-  for(i=0; i<=ap_order; i++)
-    for(j=0; j<=ap_order-i; j++)
+  for(i=0; i<=ap_order; ++i)
+    for(j=0; j<=ap_order-i; ++j)
     {
         ap_coeff[i][j]=gsl_vector_get(c_ap, p_ap);
 
@@ -1311,11 +1299,11 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
         printf("AP_%ld_%ld = %.8E\n", i, j, ap_coeff[i][j]);
         */
 
-        p_ap++;
+        ++p_ap;
     }
   p_bp=0;
-  for(i=0; i<=bp_order; i++)
-    for(j=0; j<=bp_order-i; j++)
+  for(i=0; i<=bp_order; ++i)
+    for(j=0; j<=bp_order-i; ++j)
     {
         bp_coeff[i][j]=gsl_vector_get(c_bp, p_bp);
 
@@ -1323,12 +1311,12 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
         printf("BP_%ld_%ld = %.8E\n", i, j, bp_coeff[i][j]);
         */
 
-        p_bp++;
+        ++p_bp;
     }
 
   /*For a check.
-  for(i=0; i<p_ap; i++)
-    for(j=0; j<tsize; j++)
+  for(i=0; i<p_ap; ++i)
+    for(j=0; j<tsize; ++j)
       printf("X[%ld][%ld] = %.8lf\n", i, j, gsl_matrix_get(X_ap, j, i));
   */
 
@@ -1348,8 +1336,8 @@ wcsdistortion_fitreverse(double *u, double *v, size_t 
a_order, size_t b_order,
   free(vprime);
   free(uprime);
 
-  for(i=0; i<ap_order; i++) { free(vpdict[i]); free(updict[i]); }
-  for(i=0; i<a_order; i++) { free(vdict[i]); free(udict[i]); }
+  for(i=0; i<ap_order; ++i) { free(vpdict[i]); free(updict[i]); }
+  for(i=0; i<a_order; ++i) { free(vdict[i]); free(udict[i]); }
   free(vpdict);
   free(updict);
   free(vdict);
@@ -1366,48 +1354,35 @@ static void
 wcsdistortion_get_revkeyvalues(struct wcsprm *wcs, size_t *fitsize,
                                double ap_coeff[5][5], double bp_coeff[5][5])
 {
-  size_t tsize=0;
-  size_t i=0, j=0, k=0;
-  size_t naxis1, naxis2;
-  double crpix1, crpix2;
+  size_t tsize;
+  size_t i, j, k;
   double *u=NULL, *v=NULL;
   size_t a_order=0, b_order=0;
-  double a_coeff[5][5]={0}, b_coeff[5][5]={0};
-
-  naxis2=fitsize[0];
-  naxis1=fitsize[1];
-
-  crpix1=wcs->crpix[0];
-  crpix2=wcs->crpix[1];
+  double a_coeff[5][5], b_coeff[5][5];
+  size_t naxis1=fitsize[1], naxis2=fitsize[0];
+  double crpix1=wcs->crpix[0], crpix2=wcs->crpix[1];
 
+  /* Initialise the 2d matrices. */
   tsize=(naxis1/4)*(naxis2/4);
-
+  for(i=0;i<5;++i) for(j=0;j<5;++j) {a_coeff[i][j]=0; b_coeff[i][j]=0;}
 
   /* Allocate the size of u,v arrays. */
   u=malloc(tsize*sizeof(*u));
   v=malloc(tsize*sizeof(*v));
 
   if(u==NULL)
-        error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `u'",
-              __func__, tsize*sizeof(*u));
+    error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `u'",
+          __func__, tsize*sizeof(*u));
   if(v==NULL)
-        error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `v'",
-              __func__, tsize*sizeof(*v));
+    error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `v'",
+          __func__, tsize*sizeof(*v));
 
 
   /* Make the grid and bring it's origin to the world's origin*/
-  for(i=0; i<naxis2; i+=4)
-    for(j=0; j<naxis1; j+=4)
-      u[k++]=j-crpix1;
-
-  k=0;
-  for(i=0; i<naxis2; i+=4)
-    for(j=0; j<naxis1; j+=4)
-      v[k++]=i-crpix2;
-
-
+  k=0; for(i=0; i<naxis2; i+=4) for(j=0; j<naxis1; j+=4) u[k++]=j-crpix1;
+  k=0; for(i=0; i<naxis2; i+=4) for(j=0; j<naxis1; j+=4) v[k++]=i-crpix2;
   /*For a check.
-  for(i=0; i<tsize; i++)
+    for(i=0; i<tsize; ++i)
     printf("u%ld = %.10E\n", i, u[i]);
   */
 
@@ -1419,7 +1394,6 @@ wcsdistortion_get_revkeyvalues(struct wcsprm *wcs, size_t 
*fitsize,
   /* Free the memory allocations. */
   free(v);
   free(u);
-
 }
 
 
@@ -1459,51 +1433,49 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, 
size_t *fitsize,
   size_t a_order=0, b_order=0;
   size_t ap_order=0, bp_order=0;
   size_t m, n, num=0, numkey=100;
-  double ap_coeff[5][5]={0}, bp_coeff[5][5]={0};
+  double ap_coeff[5][5], bp_coeff[5][5];
   char *fullheader, fmt[50], sipkey[8], keyaxis[9], pcaxis[10];
 
+  /* Initialise the 2d matrices. */
   *nkeys = 0;
+  for(i=0;i<5;++i) for(j=0;j<5;++j) {ap_coeff[i][j]=0; bp_coeff[i][j]=0;}
 
   /* The format for each card. */
   sprintf(fmt, "%%-8s= %%20.12E%%50s");
 
-  // wcsdistortion_calc_tpveq(wcs, cd, tpvu, tpvv);
-
-
   /* Allcate memory for cards. */
   fullheader = malloc(numkey*80);
   if(fullheader==NULL)
-        error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `fullheader'",
-              __func__, sizeof *fullheader);
-
+    error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `fullheader'",
+          __func__, sizeof *fullheader);
 
   /* Add other necessary cards. */
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20d%50s", "WCSAXES",
           wcs->naxis, "");
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CRPIX%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.8lf%50s", keyaxis,
               wcs->crpix[i-1], "");
     }
 
-  for(i=1; i<=size; i++)
-    for(j=1; j<=size; j++)
+  for(i=1; i<=size; ++i)
+    for(j=1; j<=size; ++j)
       {
         sprintf(pcaxis, "PC%d_%d", i, j);
         sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", pcaxis,
                 wcs->pc[k++], "");
       }
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CDELT%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", keyaxis,
               wcs->cdelt[i-1], "");
     }
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CUNIT%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", keyaxis,
@@ -1515,7 +1487,7 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t 
*fitsize,
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "CTYPE2",
           "'DEC--TAN-SIP'");
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CRVAL%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.10lf%50s", keyaxis,
@@ -1527,9 +1499,11 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t 
*fitsize,
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", "LATPOLE",
           wcs->latpole, "");
 
-  for(i=1; i<=size; i++)
+#if GAL_CONFIG_HAVE_WCSLIB_MJDREF == 1
+  for(i=1; i<=size; ++i)
     sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "MJDREFI",
             wcs->mjdref[i-1], "");
+#endif
 
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "RADESYS",
           wcs->radesys);
@@ -1538,8 +1512,8 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t 
*fitsize,
           wcs->equinox, "");
 
 
-  for(m=0; m<=4; m++)
-    for(n=0; n<=4; n++)
+  for(m=0; m<=4; ++m)
+    for(n=0; n<=4; ++n)
       {
         /*For axis = 1*/
         val=wcsdistortion_calcsip(1, m, n, tpvu, tpvv);
@@ -1576,8 +1550,8 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t 
*fitsize,
 
       wcsdistortion_get_revkeyvalues(wcs, fitsize, ap_coeff, bp_coeff);
 
-      for(m=0; m<=ap_order; m++)
-        for(n=0; n<=ap_order-m; n++)
+      for(m=0; m<=ap_order; ++m)
+        for(n=0; n<=ap_order-m; ++n)
           {
             /*For axis = 1*/
             val=ap_coeff[m][n];
@@ -1600,22 +1574,18 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, 
size_t *fitsize,
               }
           }
 
-    sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20ld%50s", "AP_ORDER",
-            ap_order, "");
-    sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20ld%50s", "BP_ORDER",
-            bp_order, "");
+      sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20ld%50s", "AP_ORDER",
+              ap_order, "");
+      sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20ld%50s", "BP_ORDER",
+              bp_order, "");
 
     }
 
-
-  *nkeys = num;
-
   /*For a check.
-  printf("%s\n", fullheader);
+    printf("%s\n", fullheader);
   */
-
+  *nkeys = num;
   return fullheader;
-
 }
 
 
@@ -1638,46 +1608,45 @@ wcsdistortion_add_pvkeywords(struct wcsprm *wcs, double 
*pv1,
   size_t m, n, num=0, numkey=100;
   char *fullheader, fmt[50], pvkey[8], keyaxis[9], pcaxis[10];
 
+  /* Initialize values. */
   *nkeys = 0;
 
   /* The format for each card. */
   sprintf(fmt, "%%-8s= %%20.12E%%50s");
 
-
   /* Allcate memory for cards. */
   fullheader = malloc(numkey*80);
   if(fullheader==NULL)
         error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `fullheader'",
               __func__, sizeof *fullheader);
 
-
   /* Add other necessary cards. */
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20d%50s", "WCSAXES",
           wcs->naxis, "");
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CRPIX%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.8lf%50s", keyaxis,
               wcs->crpix[i-1], "");
     }
 
-  for(i=1; i<=size; i++)
-    for(j=1; j<=size; j++)
+  for(i=1; i<=size; ++i)
+    for(j=1; j<=size; ++j)
       {
         sprintf(pcaxis, "PC%d_%d", i, j);
         sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", pcaxis,
                 wcs->pc[k++], "");
       }
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CDELT%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", keyaxis,
               wcs->cdelt[i-1], "");
     }
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CUNIT%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", keyaxis,
@@ -1689,7 +1658,7 @@ wcsdistortion_add_pvkeywords(struct wcsprm *wcs, double 
*pv1,
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "CTYPE2",
           "'DEC--TPV'");
 
-  for(i=1; i<=size; i++)
+  for(i=1; i<=size; ++i)
     {
       sprintf(keyaxis, "CRVAL%d", i);
       sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.10lf%50s", keyaxis,
@@ -1701,9 +1670,11 @@ wcsdistortion_add_pvkeywords(struct wcsprm *wcs, double 
*pv1,
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", "LATPOLE",
           wcs->latpole, "");
 
-  for(i=1; i<=size; i++)
+#if GAL_CONFIG_HAVE_WCSLIB_MJDREF == 1
+  for(i=1; i<=size; ++i)
     sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "MJDREFI",
             wcs->mjdref[i-1], "");
+#endif
 
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "RADESYS",
           wcs->radesys);
@@ -1711,8 +1682,8 @@ wcsdistortion_add_pvkeywords(struct wcsprm *wcs, double 
*pv1,
   sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "EQUINOX",
           wcs->equinox, "");
 
-  for(m=1; m<=2; m++)
-    for(n=0; n<=16; n++)
+  for(m=1; m<=2; ++m)
+    for(n=0; n<=16; ++n)
       {
         /*For axis = 1*/
         if(m == 1)
@@ -1742,13 +1713,10 @@ wcsdistortion_add_pvkeywords(struct wcsprm *wcs, double 
*pv1,
 
       }
 
-
-  *nkeys = num;
-
   /*For a check.
   printf("%s\n", fullheader);
   */
-
+  *nkeys = num;
   return fullheader;
 
 }
@@ -1911,14 +1879,19 @@ gal_wcsdistortion_tpv_to_sip(struct wcsprm *inwcs,
 {
   int ctrl=0;                /* Don't report why a keyword wasn't used. */
   int nreject=0;             /* Number of keywords rejected for syntax. */
+  double cd[2][2];
+  size_t i=0, j=0;
   size_t fulllen=0;
   char *fullheader;
-  double cd[2][2]={0};
+  int relax=WCSHDR_all;      /* Macro: use all informal WCS extensions. */
   int nwcs, sumcheck=0;
   int nkeys=0, status=0;
-  int relax=WCSHDR_all;      /* Macro: use all informal WCS extensions. */
   struct wcsprm *outwcs=NULL;
-  double tpvu[8][8] ={0}, tpvv[8][8]={0};
+  double tpvu[8][8], tpvv[8][8];
+
+  /* Initialise the 2d matrices. */
+  for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd[i][j]=0;
+  for(i=0; i<8; ++i) for(j=0; j<8; ++j) { tpvu[i][j]=0; tpvv[i][j]=0; }
 
   /* Calculate the pv equations and extract sip coefficients from it. */
   wcsdistortion_calc_tpveq(inwcs, cd, tpvu, tpvv);
@@ -1934,6 +1907,7 @@ gal_wcsdistortion_tpv_to_sip(struct wcsprm *inwcs,
   wcsdistortion_set_internalstruct(outwcs, fullheader, fulllen, status,
                                    nwcs, sumcheck);
 
+  /* Return the output WCS. */
   return outwcs;
 
 }
@@ -1957,15 +1931,19 @@ gal_wcsdistortion_sip_to_tpv(struct wcsprm *inwcs)
 {
   int ctrl=0;                /* Don't report why a keyword wasn't used. */
   int nreject=0;             /* Number of keywords rejected for syntax. */
+  double cd[2][2];
+  size_t i=0, j=0;
   size_t fulllen=0;
   char *fullheader;
-  double cd[2][2]={0};
   int nwcs, sumcheck=0;
   int nkeys=0, status=0;
   int relax=WCSHDR_all;      /* Macro: use all informal WCS extensions. */
   struct wcsprm *outwcs=NULL;
   double pv1[17]={0}, pv2[17]={0};
 
+  /* Initialise the 2d matrices. */
+  for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd[i][j]=0;
+
   /* Calculate the sip equations and extract pv coefficients from it. */
   wcsdistortion_calc_sipeq(inwcs, cd, pv1, pv2);
 



reply via email to

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