gnuastro-commits
[Top][All Lists]
Advanced

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

[gnuastro-commits] master 3df6490 001/113: Crop works in 3D, except for


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 3df6490 001/113: Crop works in 3D, except for --polygon
Date: Fri, 16 Apr 2021 10:33:29 -0400 (EDT)

branch: master
commit 3df64909451aaab65a20e03068424fccdbedcf0a
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    Crop works in 3D, except for --polygon
    
    One of the first tools necessary to work on 3D is to the required region of
    a dataset. So as the first set of works in this branch, we are enabling
    Crop (as much as possible given the time constraints) to work on 3D. With
    this commit, all crops defined by their center (using the new `--center'
    option or catalogs) now accept 3D data. For crops defined by region, the
    `--section' option also works on 3D. However, polygons still only work on
    2D. They are not immediate for my work now, so we should do them later.
    
    To allow crop to be run in any dimensionality smoothly, the old
    dimension-specific options have been removed and new dimension-agnostic
    options are defined:
    
     - `--iwidth', `--wwidth' were removed and replaced with a generic
       `--width' option. The units to interpret the value to the option are
       specified by the `--mode' option. With the new `--width' option it is
       also possible to define a non-square crop (different widths along each
       dimension). Its units in WCS mode are no longer arcseconds but are
       exactly the same units of the WCS (degrees for angles). But now, it can
       also accept fractions. So to set a width of 5 arcseconds, you can give
       it a value of `5/3600' for the angular dimensions.
    
     - The options to define the center of the crop are now removed: `--ra',
       `--dec', `--xc', `--yc'. They have all been replaced with the single
       option `--center'. This new option can take multiple values (one value
       for each dimension). With this feature, the general behavior for 2D and
       3D data will be the same and there is no need to define new options for
       higher dimensions or have unused options, or vice-versa.
    
     - The options to specify the crop center in catalog mode are now also
       dimension-agnostic: `--racol', `--deccol', `--xcol', and `--ycol' have
       been removed and replaced with a single `--coordcol'. This new option
       can be called multiple times and the order of its calling will be used
       for the respective dimension (in FITS format).
    
    Library changes:
    
     - `gal_box_overlap' now works on data of any dimensionality and thus also
       needs the number of dimensions (elements in each input array).
    
     - `gal_box_border_from_center' now accepts an array of coordinates as one
       argument and the number of dimensions as another. This allows it to work
       on any dimensionality.
    
     - `gal_wcs_pixel_scale' now replaces the old `gal_wcs_pixel_scale_deg',
       since it doesn't only apply to degrees. The pixel scale units are
       defined by the units of the WCS.
---
 bin/crop/args.h                                    | 174 ++----
 bin/crop/astcrop.conf                              |   2 -
 bin/crop/crop.c                                    |  35 +-
 bin/crop/main.h                                    |  33 +-
 bin/crop/onecrop.c                                 | 268 +++++-----
 bin/crop/onecrop.h                                 |  47 +-
 bin/crop/ui.c                                      | 592 ++++++++++++++-------
 bin/crop/ui.h                                      |  17 +-
 bin/crop/wcsmode.c                                 | 410 +++++++++-----
 bin/crop/wcsmode.h                                 |   6 +-
 bin/mkprof/mkprof.c                                |   7 +-
 bin/warp/ui.c                                      |   2 +-
 bin/warp/warp.c                                    |   2 +-
 doc/gnuastro.texi                                  | 343 ++++++------
 lib/blank.c                                        |   2 +-
 lib/box.c                                          | 184 ++++---
 lib/gnuastro-internal/options.h                    |   3 +
 lib/gnuastro/box.h                                 |   4 +-
 lib/gnuastro/dimension.h                           |   4 +-
 lib/gnuastro/wcs.h                                 |   2 +-
 lib/options.c                                      |  21 +-
 lib/wcs.c                                          |   8 +-
 tests/Makefile.am                                  |  14 +-
 tests/crop/imgcat.sh                               |   5 +-
 tests/crop/{xcyc.sh => imgcenter.sh}               |   5 +-
 tests/crop/{xcycnoblank.sh => imgcenternoblank.sh} |   7 +-
 tests/crop/wcscat.sh                               |   5 +-
 tests/crop/{radec.sh => wcscenter.sh}              |   6 +-
 28 files changed, 1233 insertions(+), 975 deletions(-)

diff --git a/bin/crop/args.h b/bin/crop/args.h
index f973bd4..4272ae3 100644
--- a/bin/crop/args.h
+++ b/bin/crop/args.h
@@ -32,6 +32,20 @@ struct argp_option program_options[] =
   {
     /* Input. */
     {
+      "mode",
+      UI_KEY_MODE,
+      "STR",
+      0,
+      "Coordinate mode `img' or `wcs'.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->mode,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_parse_coordinate_mode
+    },
+    {
       "hstartwcs",
       UI_KEY_HSTARTWCS,
       "INT",
@@ -108,7 +122,7 @@ struct argp_option program_options[] =
 
     {
       0, 0, 0, 0,
-      "Crop by center (general settings)",
+      "Crop by center",
       ARGS_GROUP_CENTER_GENERAL
     },
     {
@@ -125,92 +139,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
-      "iwidth",
-      UI_KEY_IWIDTH,
-      "INT",
-      0,
-      "Width (pixels) when crop defined by X,Y.",
-      ARGS_GROUP_CENTER_GENERAL,
-      &p->iwidthin,
-      GAL_TYPE_SIZE_T,
-      GAL_OPTIONS_RANGE_GT_0_ODD,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-    {
-      "wwidth",
-      UI_KEY_WWIDTH,
-      "FLT",
+      "width",
+      UI_KEY_WIDTH,
+      "FLT[, ...]",
       0,
-      "Width (arcseconds) for crops defined by RA,Dec.",
+      "Width when crop is defined by its center.",
       ARGS_GROUP_CENTER_GENERAL,
-      &p->wwidth,
-      GAL_TYPE_FLOAT64,
-      GAL_OPTIONS_RANGE_GT_0,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-
-
-
-
-
-    {
-      0, 0, 0, 0,
-      "Crop by center (single crop)",
-      ARGS_GROUP_CENTER_SINGLE
-    },
-    {
-      "ra",
-      UI_KEY_RA,
-      "FLT",
-      0,
-      "Right ascension of one crop box center.",
-      ARGS_GROUP_CENTER_SINGLE,
-      &p->ra,
-      GAL_TYPE_FLOAT64,
-      GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-    {
-      "dec",
-      UI_KEY_DEC,
-      "FLT",
-      0,
-      "Declination of one crop box center.",
-      ARGS_GROUP_CENTER_SINGLE,
-      &p->dec,
-      GAL_TYPE_FLOAT64,
-      GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-    {
-      "xc",
-      UI_KEY_XC,
-      "FLT",
-      0,
-      "First axis position of one crop box center.",
-      ARGS_GROUP_CENTER_SINGLE,
-      &p->xc,
-      GAL_TYPE_FLOAT64,
+      &p->width,
+      GAL_TYPE_STRING,
       GAL_OPTIONS_RANGE_ANY,
       GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
+      GAL_OPTIONS_NOT_SET,
+      ui_parse_width_and_center
     },
     {
-      "yc",
-      UI_KEY_YC,
-      "FLT",
+      "center",
+      UI_KEY_CENTER,
+      "FLT[, ...]",
       0,
-      "Second axis position of one crop box center.",
-      ARGS_GROUP_CENTER_SINGLE,
-      &p->yc,
-      GAL_TYPE_FLOAT64,
+      "Center of a single crop.",
+      ARGS_GROUP_CENTER_GENERAL,
+      &p->center,
+      GAL_TYPE_STRING,
       GAL_OPTIONS_RANGE_ANY,
       GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
+      GAL_OPTIONS_NOT_SET,
+      ui_parse_width_and_center
     },
 
 
@@ -221,7 +175,7 @@ struct argp_option program_options[] =
 
     {
       0, 0, 0, 0,
-      "Crop by center (catalog)",
+      "Crop by center (only for catalog)",
       ARGS_GROUP_CENTER_CATALOG
     },
     {
@@ -264,53 +218,14 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
-      "racol",
-      UI_KEY_RACOL,
-      "STR/INT",
-      0,
-      "Column number/info of Right Ascension (RA).",
-      ARGS_GROUP_CENTER_CATALOG,
-      &p->racol,
-      GAL_TYPE_STRING,
-      GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-    {
-      "deccol",
-      UI_KEY_DECCOL,
+      "coordcol",
+      UI_KEY_COORDCOL,
       "STR/INT",
       0,
-      "Column number/info of Declination.",
+      "Columns no./info containing coordinates.",
       ARGS_GROUP_CENTER_CATALOG,
-      &p->deccol,
-      GAL_TYPE_STRING,
-      GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-    {
-      "xcol",
-      UI_KEY_XCOL,
-      "STR/INT",
-      0,
-      "Column number/info of X (first FITS axis).",
-      ARGS_GROUP_CENTER_CATALOG,
-      &p->xcol,
-      GAL_TYPE_STRING,
-      GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
-    {
-      "ycol",
-      UI_KEY_YCOL,
-      "STR/INT",
-      0,
-      "Column number/info of Y (second FITS axis).",
-      ARGS_GROUP_CENTER_CATALOG,
-      &p->ycol,
-      GAL_TYPE_STRING,
+      &p->coordcol,
+      GAL_TYPE_STRLL,
       GAL_OPTIONS_RANGE_ANY,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
@@ -371,19 +286,6 @@ struct argp_option program_options[] =
 
 
     /* Operating mode */
-    {
-      "mode",
-      UI_KEY_MODE,
-      "STR",
-      0,
-      "Coordinate mode `img' or `wcs'",
-      GAL_OPTIONS_GROUP_OPERATING_MODE,
-      &p->modestr,
-      GAL_TYPE_STRING,
-      GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_NOT_MANDATORY,
-      GAL_OPTIONS_NOT_SET
-    },
 
 
 
diff --git a/bin/crop/astcrop.conf b/bin/crop/astcrop.conf
index d58836a..b616160 100644
--- a/bin/crop/astcrop.conf
+++ b/bin/crop/astcrop.conf
@@ -19,8 +19,6 @@
 
 # Input image and catalog parameters:
  cathdu         1
- iwidth         201
- wwidth         3
  hstartwcs      0
  hendwcs        0
 
diff --git a/bin/crop/crop.c b/bin/crop/crop.c
index 0d12f95..9cdc9b1 100644
--- a/bin/crop/crop.c
+++ b/bin/crop/crop.c
@@ -52,7 +52,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    truncated. In the log-file there is no truncation, therefore the log
    file should be used for checking the outputs, not the outputs printed on
    the screen. */
-void
+static void
 crop_verbose_info(struct onecropparams *crp)
 {
   char *filestatus, *msg;
@@ -84,7 +84,7 @@ crop_verbose_info(struct onecropparams *crp)
 
 
 /* Print final statistics in verbose mode. */
-void
+static void
 crop_verbose_final(struct cropparams *p)
 {
   char *msg;
@@ -149,7 +149,7 @@ crop_verbose_final(struct cropparams *p)
 
 
 
-void
+static void
 crop_write_to_log(struct onecropparams *crp)
 {
   char **strarr;
@@ -185,8 +185,8 @@ crop_write_to_log(struct onecropparams *crp)
 
 
 
-void *
-imgmodecrop(void *inparam)
+static void *
+crop_mode_img(void *inparam)
 {
   struct onecropparams *crp=(struct onecropparams *)inparam;
   struct cropparams *p=crp->p;
@@ -210,7 +210,7 @@ imgmodecrop(void *inparam)
       crp->out_ind=crp->indexs[i];
       crp->outfits=NULL;
       crp->numimg=1;   /* In Image mode there is only one input image. */
-      cropname(crp);
+      onecrop_name(crp);
 
       /* Crop the image. */
       onecrop(crp);
@@ -219,7 +219,7 @@ imgmodecrop(void *inparam)
       if(crp->numimg)
         {
           /* Check if the center of the crop is filled or not. */
-          crp->centerfilled=iscenterfilled(crp);
+          crp->centerfilled=onecrop_center_filled(crp);
 
           /* Add the final headers and close output FITS image: */
           gal_fits_key_write_version(crp->outfits, NULL, PROGRAM_STRING);
@@ -261,8 +261,8 @@ imgmodecrop(void *inparam)
 
 
 
-void *
-wcsmodecrop(void *inparam)
+static void *
+crop_mode_wcs(void *inparam)
 {
   struct onecropparams *crp=(struct onecropparams *)inparam;
   struct cropparams *p=crp->p;
@@ -282,21 +282,22 @@ wcsmodecrop(void *inparam)
 
 
       /* Set the sides of the crop in RA and Dec */
-      setcsides(crp);
+      wcsmode_crop_corners(crp);
 
 
       /* Go over all the images to see if this target is within their
          range or not. */
       crp->in_ind=0;
       do
-        if(radecoverlap(crp))
+        if(wcsmode_overlap(crp))
           {
             /* Open the input FITS file. */
             crp->infits=gal_fits_hdu_open_format(p->imgs[crp->in_ind].name,
                                                  p->cp.hdu, 0);
 
             /* If a name isn't set yet, set it. */
-            if(crp->name==NULL) cropname(crp);
+            if(crp->name==NULL) onecrop_name(crp);
+
 
             /* Increment the number of images used (necessary for the
                header keywords that are written in `onecrop'). Then do the
@@ -322,7 +323,7 @@ wcsmodecrop(void *inparam)
       if(crp->numimg)
         {
           /* See if the center is filled. */
-          crp->centerfilled=iscenterfilled(crp);
+          crp->centerfilled=onecrop_center_filled(crp);
 
           /* Write all the dependency versions and close the file. */
           gal_fits_key_write_version(crp->outfits, NULL, PROGRAM_STRING);
@@ -340,7 +341,7 @@ wcsmodecrop(void *inparam)
         }
       else
         {
-          cropname(crp);
+          onecrop_name(crp);
           crp->centerfilled=0;
         }
 
@@ -399,7 +400,7 @@ crop(struct cropparams *p)
 
 
   /* Set the function to run: */
-  modefunction = p->mode==IMGCROP_MODE_IMG ? &imgmodecrop : &wcsmodecrop;
+  modefunction = p->mode==IMGCROP_MODE_IMG ? &crop_mode_img : &crop_mode_wcs;
 
 
   /* Allocate the array of structures to keep the thread and parameters for
@@ -411,8 +412,8 @@ crop(struct cropparams *p)
           __func__, nt*sizeof *crp);
 
 
-  /* Distribute the indexs into the threads (this is needed even if we
-     only have one object where p->cs0 is not defined): */
+  /* Distribute the indexs into the threads (for clarity, this is needed
+     even if we only have one object). */
   gal_threads_dist_in_threads(p->catname ? p->numout : 1, nt,
                               &indexs, &thrdcols);
 
diff --git a/bin/crop/main.h b/bin/crop/main.h
index 48fd6ed..556c8d5 100644
--- a/bin/crop/main.h
+++ b/bin/crop/main.h
@@ -40,6 +40,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* Macros */
 #define LOGFILENAME             PROGRAM_EXEC".log"
 #define FILENAME_BUFFER_IN_VERB 30
+#define MAXDIM                  3
 
 
 
@@ -56,8 +57,7 @@ enum crop_modes
 
 
 /* The sides of the image keep the celestial coordinates of the four
-   sides of this image. With respect to the pixels they are.
-*/
+   sides of this image. With respect to the pixels they are. */
 struct inputimgs
 {
   char             *name;  /* File name of input image.                   */
@@ -67,8 +67,8 @@ struct inputimgs
   struct wcsprm     *wcs;  /* WCS structure of each input image.          */
   char           *wcstxt;  /* Text output of each WCS.                    */
   int           nwcskeys;  /* Number of keywords in the header WCS.       */
-  double      corners[8];  /* RA and Dec of this image corners (within).  */
-  double        sized[2];  /* Width and height of image in degrees.       */
+  double     corners[24];  /* WCS of corners (24: for 3D, 8: for 2D).     */
+  double   sized[MAXDIM];  /* Width and height of image in degrees.       */
   double  equatorcorr[2];  /* If image crosses the equator, see wcsmode.c.*/
 };
 
@@ -81,43 +81,34 @@ struct cropparams
 {
   /* Directly from command-line */
   struct gal_options_common_params cp;  /* Common parameters.             */
-  struct gal_list_str_t       *inputs;  /* All input FITS files.          */
+  gal_list_str_t       *inputs;  /* All input FITS files.                 */
   size_t             hstartwcs;  /* Header keyword No. to start read WCS. */
   size_t               hendwcs;  /* Header keyword No. to end read WCS.   */
+  int                     mode;  /* Image or WCS mode.                    */
   uint8_t       zeroisnotblank;  /* ==1: In float or double, keep 0.0.    */
   uint8_t              noblank;  /* ==1: no blank (out of image) pixels.  */
   char                 *suffix;  /* Ending of output file name.           */
   size_t           checkcenter;  /* width of a box to check for zeros     */
-  size_t              iwidthin;  /* Image mode width (in pixels).         */
-  double                wwidth;  /* WCS mode width (in arcseconds).       */
-  double                    ra;  /* RA of one crop box center.            */
-  double                   dec;  /* Dec of one crop box center.           */
-  double                    xc;  /* Center point, one crop (FITS stnrd).  */
-  double                    yc;  /* Center point, one crop (FITS stnrd).  */
+  gal_data_t           *center;  /* Center position of crop.              */
+  gal_data_t            *width;  /* Width of crop when defined by center. */
   char                *catname;  /* Name of input catalog.                */
   char                 *cathdu;  /* HDU of catalog if its a FITS file.    */
   char                *namecol;  /* Filename (without suffix) of crop col.*/
-  char                  *racol;  /* Catalog RA column                     */
-  char                 *deccol;  /* Catalog Dec column                    */
-  char                   *xcol;  /* Catalog X column                      */
-  char                   *ycol;  /* Catalog Y column                      */
+  gal_list_str_t     *coordcol;  /* Column in cat containing coordinates. */
   char                *section;  /* Section string.                       */
   char                *polygon;  /* Input string of polygon vertices.     */
   uint8_t           outpolygon;  /* ==1: Keep the inner polygon region.   */
-  char                *modestr;  /* ==1: will use X and Y coordiates.     */
 
   /* Internal */
-  int                     mode;  /* Image or WCS mode.                    */
   size_t                 numin;  /* Number of input images.               */
   size_t                numout;  /* Number of output images.              */
-  double                   *c1;  /* First coordinate from catalog.        */
-  double                   *c2;  /* Second coordinate from catalog.       */
+  double        **centercoords;  /* A 1D array for the center position.   */
   char                  **name;  /* filename of crop in row.              */
   double             *wpolygon;  /* Array of WCS polygon vertices.        */
   double             *ipolygon;  /* Array of image polygon vertices.      */
   size_t             nvertices;  /* Number of polygon vertices.           */
-  long               iwidth[2];  /* Image mode width (in pixels).         */
-  double                   res;  /* Resolution in arcseconds              */
+  long          iwidth[MAXDIM];  /* Image mode width (in pixels).         */
+  double             *pixscale;  /* Raw resolution in each dimension.     */
   time_t               rawtime;  /* Starting time of the program.         */
   int            outnameisfile;  /* Output filename is a directory.       */
   int                     type;  /* Type of output(s).                    */
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index d9f824e..46ae330 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -71,25 +71,33 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /*******************************************************************/
 /* Read the section string and set the starting and ending pixels
    based on that. */
-void
-sectionparser(struct cropparams *p, size_t *dsize,
-              long *fpixel, long *lpixel)
+static void
+onecrop_parse_section(struct cropparams *p, size_t *dsize,
+                      long *fpixel, long *lpixel)
 {
   int add;
   long read;
-  size_t dim=0;
   char *tailptr;
+  long naxes[MAXDIM];
   char forl='f', *pt=p->section;
-  long naxes[2]={dsize[1], dsize[0]};
+  size_t i, ndim=p->imgs->ndim, dim=0;
 
-  /* If control reached here, then the cropped region is not defined by its
-     center. So it makes no sense to check if the center is blank. */
+
+  /* When the user asks for a section of the dataset, then the cropped
+     region is not defined by its center. So it makes no sense to later
+     check if the center is blank or not. Hence, we will over-write it with
+     zero. */
   p->checkcenter=0;
 
-  /* Initialize the fpixel and lpixel arrays: */
-  lpixel[0]=naxes[0];
-  lpixel[1]=naxes[1];
-  fpixel[0]=fpixel[1]=1;
+
+  /* Initialize the fpixel and lpixel arrays (note that `section' is only
+     defined in image mode, so there will only be one element in `imgs'. */
+  for(i=0;i<ndim;++i)
+    {
+      fpixel[i] = 1;
+      lpixel[i] = naxes[i] = p->imgs->dsize[ ndim - i - 1 ];
+    }
+
 
   /* Parse the string. */
   while(*pt!='\0')
@@ -99,7 +107,7 @@ sectionparser(struct cropparams *p, size_t *dsize,
         {
         case ',':
           ++dim;
-          if(dim==2)
+          if(dim>=ndim)
             error(EXIT_FAILURE, 0, "Extra `,` in `%s`", p->section);
           forl='f';
           ++pt;
@@ -138,34 +146,30 @@ sectionparser(struct cropparams *p, size_t *dsize,
           else    continue;
         }
 
-      /* Put it in the correct array. Note that for the last point, we
-         don't want to include the given pixel. Unlike CFITSIO, in Crop,
-         the given intervals are not inclusive. fpixel and lpixel will be
-         directly passed to CFITSIO. So we have to correct his here.*/
+      /* Put it in the correct array. */
       if(forl=='f')
         fpixel[dim] = add ? naxes[dim]+read : read;
       else
         lpixel[dim] = add ? naxes[dim]+read : read;
       pt=tailptr;
-
-      /* For a check:
-      printf("\n\n[%ld, %ld]: fpixel=(%ld, %ld), lpixel=(%ld, %ld)\n\n",
-             naxes[0], naxes[1],
-             fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
-      */
     }
 
+
   /* Make sure the first pixel is located before/below the last pixel. */
-  if(fpixel[0]>lpixel[0] || fpixel[1]>lpixel[1])
-    error(EXIT_FAILURE, 0, "the bottom left corner coordinates "
-          "cannot be larger or equal to the top right's! Your section "
-          "string (%s) has been read as: bottom left coordinate "
-          "(%ld, %ld) to top right coordinate (%ld, %ld)",
-          p->section, fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
-
-  /*
-  printf("\n%s\nfpixel=(%ld, %ld), lpixel=(%ld, %ld)\n\n", section,
-         fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
+  for(i=0;i<ndim;++i)
+    if(fpixel[i]>lpixel[i])
+      error(EXIT_FAILURE, 0, "the bottom left corner coordinates "
+            "cannot be larger or equal to the top right's! Your section "
+            "string (%s) has been read as: bottom left coordinate "
+            "(%ld, %ld) to top right coordinate (%ld, %ld)",
+            p->section, fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
+
+  /* For a check:
+  printf("\n%s\n", p->section);
+  printf("fpixel: ("); for(i=0;i<ndim;++i) printf("%ld, ", fpixel[i]);
+  printf("\b\b)\n");
+  printf("lpixel: ("); for(i=0;i<ndim;++i) printf("%ld, ", lpixel[i]);
+  printf("\b\b)\n\n");
   exit(0);
   */
 }
@@ -176,7 +180,7 @@ sectionparser(struct cropparams *p, size_t *dsize,
 
 
 void
-crop_polygonparser(struct cropparams *p)
+onecrop_parse_polygon(struct cropparams *p)
 {
   size_t dim=0;
   char *tailptr;
@@ -280,9 +284,9 @@ crop_polygonparser(struct cropparams *p)
 
 
 
-void
-imgpolygonflpixel(double *ipolygon, size_t nvertices, long *fpixel,
-                  long *lpixel)
+static void
+onecrop_ipolygon_fl(double *ipolygon, size_t nvertices, long *fpixel,
+                    long *lpixel)
 {
   size_t i;
   double minx=FLT_MAX, miny=FLT_MAX;
@@ -399,8 +403,8 @@ polygonmask(struct onecropparams *crp, void *array, long 
*fpixel_i,
 /*******************************************************************/
 /******************          One crop.         *********************/
 /*******************************************************************/
-void
-changezerotonan(void *array, size_t size, int type)
+static void
+onecrop_zero_to_nan(void *array, size_t size, int type)
 {
   float *fp, *ffp;
   double *dp, *fdp;
@@ -434,7 +438,7 @@ changezerotonan(void *array, size_t size, int type)
 
 /* Set the output name and image output sizes. */
 void
-cropname(struct onecropparams *crp)
+onecrop_name(struct onecropparams *crp)
 {
   char **strarr;
   struct cropparams *p=crp->p;
@@ -477,73 +481,68 @@ cropname(struct onecropparams *crp)
 
 
 
-/* Find the first and last pixel of a crop from its center point (in
-   image mode or WCS mode). */
-void
-cropflpixel(struct onecropparams *crp)
+/* Find the first and last pixel of a crop. */
+static void
+onecrop_flpixel(struct onecropparams *crp)
 {
   struct cropparams *p=crp->p;
-  int ncoord=1, nelem=2, status[2]={0,0};
-  size_t *dsize=p->imgs[crp->in_ind].dsize;
-  double pixcrd[2], imgcrd[2], phi[1], theta[1];
+  size_t ndim=p->imgs->ndim;
+
+  double center[MAXDIM];
+  int ncoord=1, status=0;
+  size_t i, *dsize=p->imgs[crp->in_ind].dsize;
   long *fpixel=crp->fpixel, *lpixel=crp->lpixel;
+  double pixcrd[MAXDIM], imgcrd[MAXDIM], phi[1], theta[1];
 
   switch(p->mode)
     {
 
     case IMGCROP_MODE_IMG:
-      if(p->catname)
-        gal_box_border_from_center(p->c1[crp->out_ind], p->c2[crp->out_ind],
-                                   p->iwidth, fpixel, lpixel);
-      else if(!isnan(p->xc))
-        gal_box_border_from_center(p->xc, p->yc, p->iwidth, fpixel, lpixel);
-      else if(p->section)
-        sectionparser(p, dsize, fpixel, lpixel);
-      else if(p->polygon)
+      if(p->section)            /* Defined by section.    */
+        onecrop_parse_section(p, dsize, fpixel, lpixel);
+      else if(p->polygon)       /* Defined by a polygon.  */
         {
           if(p->outpolygon==0)
-            imgpolygonflpixel(p->ipolygon, p->nvertices, fpixel, lpixel);
+            onecrop_ipolygon_fl(p->ipolygon, p->nvertices, fpixel, lpixel);
+        }
+      else                      /* Defined by its center. */
+        {
+          for(i=0;i<ndim;++i) center[i] = p->centercoords[i][crp->out_ind];
+          gal_box_border_from_center(center, ndim, p->iwidth, fpixel, lpixel);
         }
-      else
-        error(EXIT_FAILURE, 0, "a bug! In image mode, neither of the "
-              "following has been set: a catalog, a central pixel, "
-              "a section or a polygon in the image. Please contact us "
-              "to see how it got to this impossible place! You should "
-              "have been warned of this condition long before Crop "
-              "reaches this point");
       break;
 
+
     case IMGCROP_MODE_WCS: /* In wcsmode, crp->world is already filled.   */
       if(p->polygon)       /* Note: p->iwidth was set based on p->wwidth. */
         {
           /* Fill crp->ipolygon in wcspolygonpixel, then set flpixel*/
           fillcrpipolygon(crp);
           if(p->outpolygon==0)
-            imgpolygonflpixel(crp->ipolygon, p->nvertices, fpixel, lpixel);
+            onecrop_ipolygon_fl(crp->ipolygon, p->nvertices, fpixel, lpixel);
         }
       else
         {
-          if(wcss2p(p->imgs[crp->in_ind].wcs, ncoord, nelem, crp->world,
-                    phi, theta, imgcrd, pixcrd, status) )
-            if(status[0] || status[1])
+          /* Convert `crp->world' (in WCS) into `pixcrd' (image coord). */
+          if(wcss2p(p->imgs[crp->in_ind].wcs, ncoord, ndim, crp->world,
+                    phi, theta, imgcrd, pixcrd, &status) )
+            if(status)
               error(EXIT_FAILURE, 0, "%s: wcss2p error %d: %s", __func__,
-                    status[0] ? status[0] : status[1],
-                    wcs_errmsg[status[0] ? status[0] : status[1]]);
-          gal_box_border_from_center(pixcrd[0], pixcrd[1], p->iwidth, fpixel,
-                                     lpixel);
-          /*
-            printf("\n(%f, %f): (%ld, %ld) -- (%ld, %ld)\n\n", pixcrd[0],
-                   pixcrd[1], fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
-          */
+                    status, wcs_errmsg[status]);
+
+          /* Find the first and last pixels of this crop. */
+          gal_box_border_from_center(pixcrd, ndim, p->iwidth, fpixel, lpixel);
         }
       break;
 
+
     default:
       error(EXIT_FAILURE, 0, "%s: a bug! The domain (WCS or image) are not "
             "set. Please contact us at %s so we can see how it got to this "
             "impossible place", __func__, PACKAGE_BUGREPORT);
     }
 
+
   /* If the user only wants regions outside to the polygon, then set
      the fpixel and lpixel to cover the full input image. */
   if(p->polygon && p->outpolygon)
@@ -568,17 +567,17 @@ cropflpixel(struct onecropparams *crp)
    blank areas or not. While the fpixel_i and lpixel_i arrays keep the
    first and last pixels after the blank pixels have been removed.
 */
-void
-firstcropmakearray(struct onecropparams *crp, long *fpixel_i,
+static void
+onecrop_make_array(struct onecropparams *crp, long *fpixel_i,
                    long *lpixel_i, long *fpixel_c, long *lpixel_c)
 {
-  size_t i;
+  double crpix;
   fitsfile *ofp;
-  long naxes[2];
-  int type=crp->p->type;
-  double crpix0, crpix1;
-  int naxis=2, status=0;
+  long naxes[MAXDIM];
   char *outname=crp->name;
+  char cpname[FLEN_KEYWORD];
+  int status=0, type=crp->p->type;
+  size_t i, ndim=crp->p->imgs->ndim, totsize;
   char *cp, *cpf, blankrec[80], titlerec[80];
   char startblank[]="                      / ";
   struct inputimgs *img=&crp->p->imgs[crp->in_ind];
@@ -593,16 +592,14 @@ firstcropmakearray(struct onecropparams *crp, long 
*fpixel_i,
 
   /* Set the size of the output, in WCS mode, noblank==0. */
   if(crp->p->noblank && crp->p->mode==IMGCROP_MODE_IMG)
-    {
-      fpixel_c[0]=fpixel_c[1]=1;
-      lpixel_c[0]=naxes[0]=lpixel_i[0]-fpixel_i[0]+1;
-      lpixel_c[1]=naxes[1]=lpixel_i[1]-fpixel_i[1]+1;
-    }
+    for(i=0;i<ndim;++i)
+      {
+        fpixel_c[i] = 1;
+        lpixel_c[i]=naxes[i]=lpixel_i[i]-fpixel_i[i]+1;
+      }
   else
-    {
-      naxes[0]=crp->lpixel[0]-crp->fpixel[0]+1;
-      naxes[1]=crp->lpixel[1]-crp->fpixel[1]+1;
-    }
+    for(i=0;i<ndim;++i)
+      naxes[i]=crp->lpixel[i]-crp->fpixel[i]+1;
 
 
   /* Create the FITS file with a blank first extension, then close it, so
@@ -615,13 +612,15 @@ firstcropmakearray(struct onecropparams *crp, long 
*fpixel_i,
   fits_create_img(ofp, SHORT_IMG, 0, naxes, &status);
   fits_close_file(ofp, &status);
 
+
   /* Create the output crop image. */
   fits_open_file(&crp->outfits, outname, READWRITE, &status);
   fits_create_img(crp->outfits, gal_fits_type_to_bitpix(type),
-                  naxis, naxes, &status);
+                  ndim, naxes, &status);
   gal_fits_io_error(status, "creating image");
   ofp=crp->outfits;
 
+
   /* When CFITSIO creates a FITS extension it adds two comments linking to
      the FITS paper. Since we are mentioning the version of CFITSIO and
      only use its ruitines to read/write from/to FITS files, this is
@@ -633,12 +632,14 @@ firstcropmakearray(struct onecropparams *crp, long 
*fpixel_i,
   fits_delete_key(ofp, "COMMENT", &status);
   status=0;
 
-  /* Write the blank value if necessary. */
+
+  /* Write the blank value as a FITS keyword if necessary. */
   if( type!=GAL_TYPE_FLOAT32 && type!=GAL_TYPE_FLOAT64 )
     if(fits_write_key(ofp, gal_fits_type_to_datatype(crp->p->type), "BLANK",
                       crp->p->bitnul, "pixels with no data", &status) )
       gal_fits_io_error(status, "adding Blank");
-  if(fits_write_null_img(ofp, 1, naxes[0]*naxes[1], &status))
+  totsize = naxes[0]*naxes[1] * (ndim==3?naxes[2]:1);
+  if(fits_write_null_img(ofp, 1, totsize, &status))
     gal_fits_io_error(status, "writing null array");
 
 
@@ -646,8 +647,7 @@ firstcropmakearray(struct onecropparams *crp, long 
*fpixel_i,
      update the header keywords. */
   if(img->wcs)
     {
-      crpix0 = img->wcs->crpix[0] - (fpixel_i[0]-1) + (fpixel_c[0]-1);
-      crpix1 = img->wcs->crpix[1] - (fpixel_i[1]-1) + (fpixel_c[1]-1);
+      /* Write the WCS title and the common WCS information. */
       if(fits_write_record(ofp, blankrec, &status))
         gal_fits_io_error(status, NULL);
       sprintf(titlerec, "%sWCS information", startblank);
@@ -656,9 +656,15 @@ firstcropmakearray(struct onecropparams *crp, long 
*fpixel_i,
       fits_write_record(ofp, titlerec, &status);
       for(i=0;i<img->nwcskeys-1;++i)
         fits_write_record(ofp, &img->wcstxt[i*80], &status);
-      fits_update_key(ofp, TDOUBLE, "CRPIX1", &crpix0, NULL, &status);
-      fits_update_key(ofp, TDOUBLE, "CRPIX2", &crpix1, NULL, &status);
-      gal_fits_io_error(status, NULL);
+
+      /* Correct the CRPIX keywords. */
+      for(i=0;i<ndim;++i)
+        {
+          sprintf(cpname, "CRPIX%zu", i+1);
+          crpix = img->wcs->crpix[i] - (fpixel_i[i]-1) + (fpixel_c[i]-1);
+          fits_update_key(ofp, TDOUBLE, cpname, &crpix, NULL, &status);
+          gal_fits_io_error(status, NULL);
+        }
     }
 
 
@@ -677,12 +683,7 @@ firstcropmakearray(struct onecropparams *crp, long 
*fpixel_i,
 
 
 /* The starting and ending points are set in the onecropparams structure
-   for one crop from one image. Crop that region out.
-
-   return values are:
-   0: No crop was made (not in range of image).
-   1: The input image covered at least part of the crop image.
- */
+   for one crop from one image. Crop that region out of the input. */
 void
 onecrop(struct onecropparams *crp)
 {
@@ -690,39 +691,45 @@ onecrop(struct onecropparams *crp)
   struct inputimgs *img=&p->imgs[crp->in_ind];
 
   void *array;
-  size_t cropsize;
   int status=0, anynul=0;
   char basename[FLEN_KEYWORD];
   fitsfile *ifp=crp->infits, *ofp;
   gal_fits_list_key_t *headers=NULL;
-  long fpixel_o[2], lpixel_o[2], inc[2]={1,1};
+  size_t i, j, cropsize=1, ndim=img->ndim;
   char region[FLEN_VALUE], regionkey[FLEN_KEYWORD];
-  long naxes[]={img->dsize[1], img->dsize[0]}, fpixel_i[2] , lpixel_i[2];
+  long fpixel_o[MAXDIM], lpixel_o[MAXDIM], inc[MAXDIM];
+  long naxes[MAXDIM], fpixel_i[MAXDIM], lpixel_i[MAXDIM];
+
+
+  /* Fill the `naxes' and `inc' arrays. */
+  for(i=0;i<ndim;++i)
+    {
+      inc[ i ]   = 1;
+      naxes[ i ] = img->dsize[ ndim - i - 1 ];
+    }
 
 
   /* Find the first and last pixel of this crop box from this input
-     image. If the outer polygon region is to be kept, then set the
-     sides to the image sides.*/
-  cropflpixel(crp);
-  fpixel_i[0]=crp->fpixel[0];      fpixel_i[1]=crp->fpixel[1];
-  lpixel_i[0]=crp->lpixel[0];      lpixel_i[1]=crp->lpixel[1];
+     image. Then copy the first and last pixels into the `_i' arrays.*/
+  onecrop_flpixel(crp);
+  memcpy(fpixel_i, crp->fpixel, ndim*sizeof *fpixel_i);
+  memcpy(lpixel_i, crp->lpixel, ndim*sizeof *lpixel_i);
 
 
   /* Find the overlap and apply it if there is any overlap. */
-  if( gal_box_overlap(naxes, fpixel_i, lpixel_i, fpixel_o, lpixel_o) )
+  if( gal_box_overlap(naxes, fpixel_i, lpixel_i, fpixel_o, lpixel_o, ndim) )
     {
       /* Make the output FITS image and initialize it with an array of
-         NaN or BLANK values. Note that for FLOAT_IMG and DOUBLE_IMG,
-         it will automatically fill them with the NaN value.*/
+         NaN or BLANK values. */
       if(crp->outfits==NULL)
-        firstcropmakearray(crp, fpixel_i, lpixel_i, fpixel_o, lpixel_o);
+        onecrop_make_array(crp, fpixel_i, lpixel_i, fpixel_o, lpixel_o);
       ofp=crp->outfits;
 
 
-      /* Read the desired part of the image, then write it into this
-         array. */
+      /* Allocate an array to keep the desired crop region, then read the
+         desired pixels into it. */
       status=0;
-      cropsize=(lpixel_i[0]-fpixel_i[0]+1)*(lpixel_i[1]-fpixel_i[1]+1);
+      for(i=0;i<ndim;++i) cropsize *= ( lpixel_i[i] - fpixel_i[i] + 1 );
       array=gal_data_malloc_array(p->type, cropsize, __func__, "array");
       if(fits_read_subset(ifp, gal_fits_type_to_datatype(p->type),
                           fpixel_i, lpixel_i, inc, p->bitnul, array,
@@ -736,13 +743,18 @@ onecrop(struct onecropparams *crp)
       if(p->zeroisnotblank==0
          && (p->type==GAL_TYPE_FLOAT32
              || p->type==GAL_TYPE_FLOAT64) )
-        changezerotonan(array, cropsize, p->type);
+        onecrop_zero_to_nan(array, cropsize, p->type);
 
 
       /* If a polygon is given, remove all the pixels within or
          outside of it.*/
       if(p->polygon)
         {
+          /* A small sanity check until this part supports 3D. */
+          if(ndim!=2)
+            error(EXIT_FAILURE, 0, "%s: polygons only implemented in 2D",
+                  __func__);
+
           /* In WCS mode, crp->ipolygon was allocated and filled in
              wcspolygonflpixel (wcsmode.c). */
           if(p->mode==IMGCROP_MODE_IMG) crp->ipolygon=p->ipolygon;
@@ -759,14 +771,18 @@ onecrop(struct onecropparams *crp)
         gal_fits_io_error(status, NULL);
 
 
+      /* Write the selected region of this image as a string to include as
+         a FITS keyword. Then we want to delete the last coma `,'.*/
+      j=0;
+      for(i=0;i<ndim;++i)
+        j += sprintf(&region[j], "%ld:%ld,", fpixel_i[i], lpixel_i[i]);
+      region[j-1]='\0';
+
       /* A section has been added to the cropped image from this input
-         image, so increment crp->imgcount and save the information of
-         this image. */
+         image, so save the information of this image. */
       sprintf(basename, "ICF%zu", crp->numimg);
       gal_fits_key_write_filename(basename, img->name, &headers);
       sprintf(regionkey, "%sPIX", basename);
-      sprintf(region, "%ld:%ld,%ld:%ld", fpixel_i[0], lpixel_i[0],
-              fpixel_i[1], lpixel_i[1]);
       gal_fits_key_list_add_end(&headers, GAL_TYPE_STRING, regionkey,
                                 0, region, 0, "Range of pixels used for "
                                 "this output.", 0, NULL);
@@ -807,7 +823,7 @@ onecrop(struct onecropparams *crp)
 /******************        Check center        *********************/
 /*******************************************************************/
 int
-iscenterfilled(struct onecropparams *crp)
+onecrop_center_filled(struct onecropparams *crp)
 {
   struct cropparams *p=crp->p;
 
@@ -848,7 +864,7 @@ iscenterfilled(struct onecropparams *crp)
   array=gal_data_malloc_array(type, size, __func__, "array");
   if( fits_read_subset(ofp, gal_fits_type_to_datatype(type), fpixel, lpixel,
                        inc, p->bitnul, array, &anynul, &status) )
-    gal_fits_io_error(status, NULL);
+     gal_fits_io_error(status, NULL);
   free(array);
 
   /* CFITSIO already checks if there are any blank pixels. If there are,
diff --git a/bin/crop/onecrop.h b/bin/crop/onecrop.h
index 2384567..eba6243 100644
--- a/bin/crop/onecrop.h
+++ b/bin/crop/onecrop.h
@@ -35,46 +35,43 @@ struct onecropparams
   struct   cropparams *p;
 
   /* About input image. */
-  size_t          in_ind;  /* Index of this image in the input names.  */
-  fitsfile       *infits;  /* Pointer to the input FITS image.         */
-  long         fpixel[2];  /* Position of first pixel in input image.  */
-  long         lpixel[2];  /* Position of last pixel in input image.   */
-  double       *ipolygon;  /* Input image based polygon vertices.      */
+  size_t              in_ind;  /* Index of this image in the input names.  */
+  fitsfile           *infits;  /* Pointer to the input FITS image.         */
+  long        fpixel[MAXDIM];  /* Position of first pixel in input image.  */
+  long        lpixel[MAXDIM];  /* Position of last pixel in input image.   */
+  double           *ipolygon;  /* Input image based polygon vertices.      */
 
   /* Output (cropped) image. */
-  size_t         out_ind;  /* Index of this crop in the output list.   */
-  double        world[2];  /* World coordinates of crop center.        */
-  double        sized[2];  /* Width and height of image in degrees.    */
-  double      corners[8];  /* RA and Dec of this crop's four sides.    */
-  double  equatorcorr[2];  /* Crop crosses the equator, see wcsmode.c. */
-  fitsfile      *outfits;  /* Pointer to the output FITS image.        */
+  size_t             out_ind;  /* Index of this crop in the output list.   */
+  double       world[MAXDIM];  /* World coordinates of crop center.        */
+  double       sized[MAXDIM];  /* Width and height of image in degrees.    */
+  double         corners[24];  /* RA and Dec of this crop's corners.       */
+  double      equatorcorr[2];  /* Crop crosses the equator, see wcsmode.c. */
+  fitsfile          *outfits;  /* Pointer to the output FITS image.        */
 
   /* For log */
-  char             *name;  /* Filename of crop.                        */
-  size_t          numimg;  /* Number of images used to make this crop. */
-  unsigned char centerfilled;   /* ==1 if the center is filled.        */
+  char                 *name;  /* Filename of crop.                        */
+  size_t              numimg;  /* Number of images used to make this crop. */
+  unsigned char centerfilled;  /* ==1 if the center is filled.             */
 
   /* Thread parameters. */
-  size_t         *indexs;  /* Indexs to be used in this thread.        */
-  pthread_barrier_t   *b;  /* pthread barrier to keep threads waiting. */
+  size_t             *indexs;  /* Indexs to be used in this thread.        */
+  pthread_barrier_t       *b;  /* pthread barrier to keep threads waiting. */
 };
 
 void
-crop_polygonparser(struct cropparams *p);
+onecrop_parse_polygon(struct cropparams *p);
 
 void
-cropname(struct onecropparams *crp);
-
-void
-cropflpixel(struct onecropparams *crp);
-
-void
-onecrop(struct onecropparams *crp);
+onecrop_name(struct onecropparams *crp);
 
 int
-iscenterfilled(struct onecropparams *crp);
+onecrop_center_filled(struct onecropparams *crp);
 
 void
 crop_print_log(struct onecropparams *p);
 
+void
+onecrop(struct onecropparams *crp);
+
 #endif
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 9b9e5ac..e766bfa 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
+#include <gnuastro/dimension.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/options.h>
@@ -64,7 +65,7 @@ const char *
 argp_program_bug_address = PACKAGE_BUGREPORT;
 
 static char
-args_doc[] = "[ASCIIcatalog] ASTRdata ...";
+args_doc[] = "[Crop-Identifier] ASTRdata ...";
 
 const char
 doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will create cutouts, "
@@ -133,7 +134,6 @@ ui_initialize_options(struct cropparams *p,
 
 
   /* Initalize necessary parameters. */
-  p->xc=p->yc=p->ra=p->dec=NAN;
   p->mode         = IMGCROP_MODE_INVALID;
   cp->searchin    = GAL_TABLE_SEARCH_INVALID;
 
@@ -218,14 +218,105 @@ parse_opt(int key, char *arg, struct argp_state *state)
 
 
 
+void *
+ui_parse_coordinate_mode(struct argp_option *option, char *arg,
+                         char *filename, size_t lineno, void *junk)
+{
+  char *outstr;
 
+  /* We want to print the stored values. */
+  if(lineno==-1)
+    {
+      gal_checkset_allocate_copy( *(int *)(option->value)==IMGCROP_MODE_IMG
+                                  ? "img" : "wcs", &outstr );
+      return outstr;
+    }
+  else
+    {
+      if      (!strcmp(arg, "img")) *(int *)(option->value)=IMGCROP_MODE_IMG;
+      else if (!strcmp(arg, "wcs")) *(int *)(option->value)=IMGCROP_MODE_WCS;
+      else
+        error_at_line(EXIT_FAILURE, 0, filename, lineno, "`%s' (value to "
+                      "`--mode') not recognized as an input mode. "
+                      "Recognized values are `img' and `wcs'. This option "
+                      "is necessary to identify the nature of your input "
+                      "coordinates.\n\n"
+                      "Please run the following command for more "
+                      "information (press the `SPACE' key to go down and "
+                      "`q' to return to the command-line):\n\n"
+                      "    $ info gnuastro \"Crop modes\"\n", arg);
+      return NULL;
+    }
+}
 
 
 
 
 
+/* Parse the width and center coordinates from the comman-line or
+   configuration files. */
+void *
+ui_parse_width_and_center(struct argp_option *option, char *arg,
+                          char *filename, size_t lineno, void *junk)
+{
+  size_t i, nc;
+  double *darray;
+  gal_data_t *values;
+  char *str, sstr[GAL_OPTIONS_STATIC_MEM_FOR_VALUES];
 
+  /* We want to print the stored values. */
+  if(lineno==-1)
+    {
+      /* Set the value pointer to the correct type. */
+      values = *(gal_data_t **)(option->value);
+      darray = values->array;
 
+      /* Write the values into a string. */
+      nc=0;
+      for(i=0;i<values->size;++i)
+        {
+          if( nc > GAL_OPTIONS_STATIC_MEM_FOR_VALUES-100 )
+            error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so we "
+                  "can address the problem. The number of necessary "
+                  "characters in the statically allocated string has become "
+                  "too close to %d", __func__, PACKAGE_BUGREPORT,
+                  GAL_OPTIONS_STATIC_MEM_FOR_VALUES);
+          nc += sprintf(sstr+nc, "%g,", darray[i]);
+        }
+      sstr[nc-1]='\0';
+
+      /* Copy the string into a dynamically allocated space, because it
+         will be freed later.*/
+      gal_checkset_allocate_copy(sstr, &str);
+      return str;
+    }
+  else
+    {
+      /* If the option is already set, then ignore it. */
+      if(option->set) return NULL;
+
+      /* Read the values. */
+      values=gal_options_parse_list_of_numbers(arg, filename, lineno);
+      *(gal_data_t **)(option->value) = values;
+
+      /* If we are on the width option, then make sure none of the values
+         are negative or zero. */
+      if(!strcmp(option->name, "width"))
+        {
+          darray=values->array;
+          for(i=0;i<values->size;++i)
+            if(darray[i]<=0.0f)
+              error_at_line(EXIT_FAILURE, 0, filename, lineno, "%g is <=0. "
+                            "The values to the `--width' option must be "
+                            "larger than zero. The complete input to this "
+                            "option was `%s' (%g is input number %zu)",
+                            darray[i], arg, darray[i], i+1);
+        }
+
+      /* The return value is only for printing mode. */
+      return NULL;
+    }
+}
 
 
 
@@ -233,68 +324,41 @@ parse_opt(int key, char *arg, struct argp_state *state)
 
 
 
-/**************************************************************/
-/***************       Sanity Check         *******************/
-/**************************************************************/
-/* Read and check ONLY the options. When arguments are involved, do the
-   check in `ui_check_options_and_arguments'. */
-static void
-ui_read_check_only_options(struct cropparams *p)
-{
-  int checksum;
 
-  /* Read the mode from the string the user specified. */
-  if(p->modestr)
-    {
-      if      (!strcmp(p->modestr, "img"))   p->mode=IMGCROP_MODE_IMG;
-      else if (!strcmp(p->modestr, "wcs"))   p->mode=IMGCROP_MODE_WCS;
-      else
-        error(EXIT_FAILURE, 0, "`%s' (value to `--mode') not recognized as "
-              "an input mode. Recognized values are `img' and `wcs'. This "
-              "option is necessary when the inputs are not sufficient to "
-              "identify the nature of the coordinates.\n\n"
-              "Please run the following command for more information "
-              "(press the `SPACE' key to go down and `q' to return to the "
-              "command-line):\n\n"
-              "    $ info gnuastro \"Crop modes\"\n", p->modestr);
-    }
 
 
 
 
-  /* Check that if one of the coordinates is given, the other is also
-     given. Note that if both are not given, their sum will be 0, if both
-     are given, then the sum will be 2. If only one is given, then the sum
-     will be 1. */
-  if( (isnan(p->ra) + isnan(p->dec)) == 1 )
-    error(EXIT_FAILURE, 0, "no `--%s' given, it is mandatory with the `--%s' "
-          "option", isnan(p->ra)?"ra":"dec", isnan(p->ra)?"dec":"ra");
 
-  if( (isnan(p->xc) + isnan(p->yc)) == 1 )
-    error(EXIT_FAILURE, 0, "no `--%s' given, it is mandatory with the `--%s' "
-          "option", isnan(p->xc)?"xc":"yc", isnan(p->xc)?"yc":"xc");
 
-  if( ((p->racol!=NULL) + (p->deccol!=NULL)) == 1 )
-    error(EXIT_FAILURE, 0, "no `--%scol' given, it is mandatory with the "
-          "`--%scol' option", p->racol?"dec":"ra",
-          p->racol?"ra":"dec");
 
-  if( ((p->xcol!=NULL) + (p->ycol!=NULL)) == 1 )
-    error(EXIT_FAILURE, 0, "no `--%scol' given, it is mandatory with the "
-          "`--%scol' option", p->xcol?"y":"x", p->xcol?"x":"y");
 
 
 
-  /* Make sure that the single crop modes are not called together. */
-  checksum = ( (!isnan(p->xc)) + (!isnan(p->ra)) + (p->catname!=NULL)
-               + (p->section!=NULL) + (p->polygon!=NULL) );
+
+/**************************************************************/
+/***************       Sanity Check         *******************/
+/**************************************************************/
+/* Read and check ONLY the options. When arguments are involved, do the
+   check in `ui_check_options_and_arguments'. */
+static void
+ui_read_check_only_options(struct cropparams *p)
+{
+  int checksum;
+
+
+  /* Make sure that only one of the crop definitions is given. */
+  checksum = ( (p->center!=NULL)
+               + (p->catname!=NULL)
+               + (p->section!=NULL)
+               + (p->polygon!=NULL) );
   switch(checksum)
     {
     case 0:
       error(EXIT_FAILURE, 0, "no crop definition. You can use any of the "
-            "following options to define the crop(s): (`--xc', `--yc'), "
-            "(`--ra', `--dec'), `--catalog', `--section', `--polygon'. "
-            "Please run this command for more information:\n\n"
+            "following options to define the crop(s): `--center', "
+            "`--catalog', `--section', or `--polygon'. Please run this "
+            "command for more information:\n\n"
             "    $ info gnuastro \"Crop modes\"\n");
     case 1:
       /* Everything is ok, just ignore the switch structure. */
@@ -302,16 +366,18 @@ ui_read_check_only_options(struct cropparams *p)
     default:
       error(EXIT_FAILURE, 0, "more than one crop type specified. In each "
             "run, only one crop definition is acceptable on the "
-            "command-line, or in configuration files. You have called "
-            "%s%s%s%s%s\b\b.",
-            !isnan(p->xc) ? "(`--xc', `--yc'), " : "",
-            !isnan(p->ra) ? "(`--ra', `--dec'), " : "",
+            "command-line or in configuration files. You have called: "
+            "%s%s%s%s\b\b.",
+            p->center!=NULL  ? "`--center', " : "",
             p->catname!=NULL ? "`--catalog', " : "",
             p->section!=NULL ? "`--section', " : "",
             p->polygon!=NULL ? "`--polygon', " : "");
     }
 
 
+  /* Section is currently only defined in image mode. */
+  if(p->section) p->mode=IMGCROP_MODE_IMG;
+
 
   /* Sanity checks and mode setting based on the desired crop. */
   if(p->catname)
@@ -334,50 +400,22 @@ ui_read_check_only_options(struct cropparams *p)
       /* Atleast one of the (X,Y), and (RA,Dec) set of columns are
          necessary. Note that we have checked that they are together if
          given, so we only need to check one of the two in each couple. */
-      if(p->xcol==NULL && p->racol==NULL)
-        error(EXIT_FAILURE, 0, "no crop center position column given "
-              "to read from the input catalog (`%s'). Please use either of "
-              "(`--xcol', `--ycol') or (`--racol', `--deccol'). For more "
-              "information on how to select columns in Gnuastro, please run "
-              "the following command:\n\n"
+      if(p->coordcol==NULL)
+        error(EXIT_FAILURE, 0, "no crop center columns given to read from "
+              "the input catalog (`%s'). Please use `--coordcol' several "
+              "times (depending on dimensionality) to specify the column "
+              "keeping the center position the respective dimension.\n\n"
+              "For more information on how to select columns in Gnuastro, "
+              "please run the following command:\n\n"
               "    $ info gnuastro \"Selecting table columns\"", p->catname);
-
-      /* If only image columns are specified, then we have Image mode, if
-         only WCS columns are specified, we have WSC mode. If both are
-         specified, the mode is mandatory */
-      if(p->xcol && p->racol)
-        {
-          if(p->mode==IMGCROP_MODE_INVALID)
-            error(EXIT_FAILURE, 0, "both image and WCS coordinate columns "
-                  "are specified to read the center of the crops in the "
-                  "input catalog (%s). You can use the `--mode=img', or "
-                  "`--mode=wcs' options to specify which set of columns "
-                  "should be used", p->catname);
-        }
-      else
-        p->mode = p->xcol ? IMGCROP_MODE_IMG : IMGCROP_MODE_WCS;
     }
-  else if(p->polygon)
-    {
-      if(p->mode==IMGCROP_MODE_INVALID)
-        error(EXIT_FAILURE, 0, "the polygon option works in both image and "
-              "wcs mode, please specify how the vertices should be "
-              "interpreted with `--mode=img', or `--mode=wcs' options to "
-              "specify which set of columns should be used");
-    }
-  else if(!isnan(p->xc))
-    p->mode = IMGCROP_MODE_IMG;
-  else if(!isnan(p->ra))
-    p->mode = IMGCROP_MODE_WCS;
-  else if( p->section)
-    p->mode = IMGCROP_MODE_IMG;
 
 
   /* Parse the polygon vertices if they are given to make sure that it is
      in the proper format. */
   if(p->polygon)
     {
-      crop_polygonparser(p);
+      onecrop_parse_polygon(p);
       if(p->nvertices<3)
         error(EXIT_FAILURE, 0, "a polygon has to have 3 or more vertices, "
               "you have only given %zu (%s)", p->nvertices, p->polygon);
@@ -405,12 +443,11 @@ ui_read_check_only_options(struct cropparams *p)
 static void
 ui_check_options_and_arguments(struct cropparams *p)
 {
-  /* Make sure we do actually have inputs. */
+  /* Make sure we actually have inputs. */
   if(p->inputs==NULL)
     error(EXIT_FAILURE, 0, "no input file given");
 
-  /* Make sure an input file name was given and if it was a FITS file, that
-     a HDU is also given. */
+  /* Make sure that a HDU is also given. */
   if(p->cp.hdu==NULL )
     error(EXIT_FAILURE, 0, "no HDU specified. When the input is a FITS "
           "file, a HDU must also be specified, you can use the `--hdu' "
@@ -475,88 +512,186 @@ ui_check_options_and_arguments(struct cropparams *p)
 /**************************************************************/
 /***************       Preparations         *******************/
 /**************************************************************/
+/* When the crop is defined by its center, the final width that we need
+   must be in actual number of pixels (an integer). But the user's values
+   can be in WCS mode or even in image mode, they may be non-integers. */
+static void
+ui_set_iwidth(struct cropparams *p)
+{
+  gal_data_t *newwidth;
+  double pwidth, *warray;
+  size_t i, ndim=p->imgs->ndim;
+
+
+  /* Make sure a width value is actually given. */
+  if(p->width==NULL)
+    error(EXIT_FAILURE, 0, "no crop width specified. When crops are "
+          "defined by their center (with `--center' or `--catalog') a "
+          "width is necessary (using the `--width' option)");
+
+  /* Make sure that the width array only has one element or the same number
+     of elements as the input's dimensions. */
+  if(p->width->size!=ndim && p->width->size!=1)
+    error(EXIT_FAILURE, 0, "%zu values give to `--width', but input is %zu "
+          "dimensions. It can only take either one value (same width in all "
+          "dimensions), or the same number as the input's dimensions",
+          p->width->size, ndim);
+
+  /* If the width array has only one value, that single value should be
+     used for all dimensions. */
+  if(p->width->size==1)
+    {
+      /* Allocate the new width dataset. */
+      newwidth=gal_data_alloc(NULL, p->width->type, 1, &ndim, NULL, 0, -1,
+                              NULL, NULL, NULL);
+
+      /* Fill the new width. */
+      warray=newwidth->array;
+      for(i=0;i<ndim;++i) warray[i] = *(double *)(p->width->array);
+
+      /* Free the old (single element) width dataset and put the new one in
+         its place. */
+      gal_data_free(p->width);
+      p->width=newwidth;
+    }
+  else warray=p->width->array;
+
+  /* Fill in `p->iwidth' depending on the mode. */
+  for(i=0;i<ndim;++i)
+    {
+      /* Set iwidth. */
+      if(p->mode==IMGCROP_MODE_WCS)
+        {
+          /* Convert the width in units of the input's WCS into pixels. */
+          pwidth = warray[i]/p->pixscale[i];
+          if(pwidth<3 || pwidth>50000)
+            error(EXIT_FAILURE, 0, "%g (width along dimension %zu) "
+                  "translates to %.0f pixels. This is probably not what "
+                  "you wanted. Note that the resolution in this dimension "
+                  "is %g", warray[i], i+1, pwidth, p->pixscale[i]);
+
+          /* Write the single valued width in WCS and the image width for
+             this dimension. */
+          p->iwidth[i]=GAL_DIMENSION_FLT_TO_INT(pwidth);
+          if(p->iwidth[i]%2==0)
+            {
+              p->iwidth[i] += 1;
+              warray[i]    += p->pixscale[i];
+            }
+        }
+      else
+        {
+          p->iwidth[i]=GAL_DIMENSION_FLT_TO_INT(warray[i]);
+          if(p->iwidth[i]%2==0) p->iwidth[i] += 1;
+        }
+    }
+
+  /* For a check:
+  printf("Width: "); for(i=0;i<ndim;++i) printf("\t%ld\n\t", p->iwidth[i]);
+  exit(0);
+  */
+}
+
+
+
+
+
 static void
 ui_read_cols(struct cropparams *p)
 {
-  char *colname;
-  size_t counter=0;
-  int toomanycols=0;
+  char colname[100];
   gal_list_str_t *colstrs=NULL;
   gal_data_t *cols, *tmp, *corrtype=NULL;
-  char *ax1col = p->mode==IMGCROP_MODE_IMG ? p->xcol : p->racol;
-  char *ax2col = p->mode==IMGCROP_MODE_IMG ? p->ycol : p->deccol;
+  size_t ncoordcols, counter=0, dcounter=0, ndim=p->imgs->ndim;
 
-  /* Specify the order of columns. */
+
+  /* See if the number of columns given for coordinates corresponds to the
+     number of dimensions of the input dataset. */
+  if(p->coordcol)
+    {
+      /* Check if the number of columns given for coordinates is the same
+         as the number of dimensions in the input dataset(s). */
+      ncoordcols=gal_list_str_number(p->coordcol);
+      if( ncoordcols != ndim)
+        error(EXIT_FAILURE, 0, "`--coordcol' was called %zu times, but the "
+              "input dataset%s %zu dimensions, these values must not be "
+              "different. Recall that through `--coordcol' you are "
+              "specifying the columns containing the coordinates of the "
+              "center of the crop in a catalog", ncoordcols,
+              (p->numin==1?" has":"s have"), ndim);
+    }
+  else
+    error(EXIT_FAILURE, 0, "no coordinate columns specified. When a "
+          "is given, it is necessary to identify which columns identify "
+          "the coordinate values in which dimension.\n\n"
+          "You can do this by calling `--coordcol' multiple times, the "
+          "order must be in the same order as the input's dimensions. "
+          "For more information on how to select columns in Gnuastro, "
+          "please run the following command:\n\n"
+          "    $ info gnuastro \"Selecting table columns\"");
+
+
+  /* If a name column was also given, the put that as the first column to
+     read, otherwise just use the given set of columns (in the same
+     order). */
   if(p->namecol)
-    gal_list_str_add(&colstrs, p->namecol, 0);
-  gal_list_str_add(&colstrs, ax1col, 0);
-  gal_list_str_add(&colstrs, ax2col, 0);
+    {
+      gal_list_str_add(&colstrs, p->namecol, 0);
+      colstrs->next=p->coordcol;
+    }
+  else colstrs=p->coordcol;
+
 
   /* Read the desired columns from the file. */
   cols=gal_table_read(p->catname, p->cathdu, colstrs, p->cp.searchin,
                       p->cp.ignorecase, p->cp.minmapsize);
 
-  /* Set the number of objects. */
+
+  /* Set the number of objects (rows in each column). */
   p->numout=cols->size;
 
-  /* For a sanity check, make sure that the total number of columns read is
-     the same as those that were wanted (it might be more). */
+
+  /* Make sure more columns were not read (the name matchings might result
+     in more than one column being read from the inputs). */
+  if( gal_list_data_number(cols) != ndim + (p->namecol!=NULL) )
+    gal_tableintern_error_col_selection(p->catname, p->cathdu, "too many "
+                                        "columns were selected by the given "
+                                        "values to the options ending in "
+                                        "`col'.");
+
+
+  /* Put the information in each column in the proper place. */
   while(cols!=NULL)
     {
       /* Pop out the top node. */
       tmp=gal_list_data_pop(&cols);
 
-      /* Note that the input was a linked list, so the output order is the
-         inverse of the input order. For the position, we will store the
-         values into the `x' and `y' arrays even if they are RA/Dec. */
+      /* See which column it is. */
       switch(++counter)
         {
-        case 3:
-          colname="crop name prefix";
+        case 1:
           if(p->namecol)
             {
-              if(tmp->type==GAL_TYPE_STRING)
-                {
-                  p->name=tmp->array;
-                  tmp->array=NULL;
-                  gal_data_free(tmp);
-                }
-              else
-                {
-                  corrtype=gal_data_copy_to_new_type_free(tmp,
-                                            GAL_TYPE_STRING);
-                  p->name=corrtype->array;
-                }
+              sprintf(colname, "crop name prefix");
+              corrtype = (tmp->type==GAL_TYPE_STRING ? tmp
+                          : gal_data_copy_to_new_type_free(tmp,
+                                                            GAL_TYPE_STRING));
+              p->name=corrtype->array;
             }
           else
-            toomanycols=1;
-          break;
-
-        case 2:
-          colname="first axis position";
-          corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
-          p->c1=corrtype->array;
-          break;
-
-        case 1:
-          colname="second axis position";
-          corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
-          p->c2=corrtype->array;
+            {
+              sprintf(colname, "position in dimension %zu", dcounter+1);
+              corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+              p->centercoords[ dcounter++ ]=corrtype->array;
+            }
           break;
 
-        /* If the index isn't recognized, then it is larger, showing that
-           there was more than one match for the given criteria */
         default:
-          toomanycols=1;
+          sprintf(colname, "position in dimension %zu", dcounter+1);
+          corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+          p->centercoords[ dcounter++ ] = corrtype->array;
         }
 
-      /* Print an error if there were too many columns. */
-      if(toomanycols)
-        gal_tableintern_error_col_selection(p->catname, p->cathdu, "too many "
-                                            "columns were selected by the "
-                                            "given values to the options "
-                                            "ending in `col'.");
-
       /* Sanity check and clean up.  Note that it might happen that the
          input structure is already freed. In that case, `corrtype' will be
          NULL. */
@@ -564,8 +699,10 @@ ui_read_cols(struct cropparams *p)
         {
           /* Make sure there are no blank values in this column. */
           if( gal_blank_present(corrtype, 1) )
-            error(EXIT_FAILURE, 0, "%s column has blank values. "
-                  "Input columns cannot contain blank values", colname);
+            error(EXIT_FAILURE, 0, "%s: column with %s has blank values. "
+                  "Input columns must not contain blank values",
+                  gal_fits_name_save_as_string(p->catname, p->cathdu),
+                  colname);
 
           /* Free the unnecessary sturcture information. The correct-type
              (`corrtype') data structure's array is necessary for later
@@ -583,6 +720,46 @@ ui_read_cols(struct cropparams *p)
 
 
 
+static void
+ui_prepare_center(struct cropparams *p)
+{
+  double *carray;
+  size_t i, ndim=p->imgs->ndim;
+
+  /* Allocate space to keep the central positions. */
+  errno=0;
+  p->centercoords=malloc(ndim*sizeof *p->centercoords);
+  if( p->centercoords==NULL )
+    error(EXIT_FAILURE, 0, "%s: %zu bytes for `p->centercoords'",
+          __func__, ndim*sizeof *p->centercoords);
+
+
+  /* Set the integer widths of the crop(s) when defined by center. */
+  ui_set_iwidth(p);
+
+  /* For a catalog, we have a separate function, but for a single center
+     value, put the center values into an array. This will essentially
+     simulate a catalog with one row. So from this point on, there is no
+     difference between a catalog input and a central position input. */
+  if(p->catname)
+    ui_read_cols(p);
+  else
+    {
+      carray=p->center->array;
+      for(i=0;i<ndim;++i)
+        {
+          p->centercoords[i]=gal_data_malloc_array(GAL_TYPE_FLOAT64,
+                                                   1, __func__,
+                                                   "p->centercoords[i]");
+          p->centercoords[i][0]=carray[i];
+        }
+    }
+}
+
+
+
+
+
 /* Add all the columns of the log file. Just note that since this is a
    linked list, we have to add them in the opposite order. */
 static void
@@ -590,10 +767,10 @@ ui_make_log(struct cropparams *p)
 {
   char *comment;
 
-  /* Return if no long file is to be created. */
+  /* Return if no long file was requested. */
   if(p->cp.log==0) return;
 
-  /* If central pixels are filled. */
+  /* Column to specify if central pixels are filled. */
   asprintf(&comment, "Are the central pixels filled? (1: yes, 0: no, "
            "%u: not checked)", GAL_BLANK_UINT8);
   gal_list_data_add_alloc(&p->log, NULL, GAL_TYPE_UINT8, 1, &p->numout,
@@ -601,12 +778,12 @@ ui_make_log(struct cropparams *p)
                           "bool", comment);
   free(comment);
 
-  /* Number of images used. */
+  /* Column for number of datasets used in this crop. */
   gal_list_data_add_alloc(&p->log, NULL, GAL_TYPE_UINT16, 1, &p->numout,
                           NULL, 1, p->cp.minmapsize, "NUM_INPUTS", "count",
-                          "Number of input images used to make this crop.");
+                          "Number of input datasets used to make this crop.");
 
-  /* Row number in input catalog. */
+  /* Filename of crop. */
   gal_list_data_add_alloc(&p->log, NULL, GAL_TYPE_STRING, 1, &p->numout,
                           NULL, 1, p->cp.minmapsize, "CROP_NAME", "name",
                           "File name of crop.");
@@ -619,16 +796,10 @@ ui_make_log(struct cropparams *p)
 void
 ui_preparations(struct cropparams *p)
 {
-  char *msg;
   fitsfile *tmpfits;
-  struct timeval t1;
-  size_t input_counter;
   struct inputimgs *img;
   int status, firsttype=0;
-
-
-  /* Set the initial iwidth. */
-  p->iwidth[0]=p->iwidth[1]=p->iwidthin;
+  size_t input_counter, firstndim=0;
 
 
   /* For polygon and section, there should be no center checking. */
@@ -636,22 +807,6 @@ ui_preparations(struct cropparams *p)
     p->checkcenter=0;
 
 
-  /* Read the columns if there is a catalog, otherwise, set the number
-     of images (crops) to 1.*/
-  if(p->catname)
-    ui_read_cols(p);
-  else
-    p->numout=1;
-
-
-  /* Everything is ready, notify the user of the program starting. */
-  if(!p->cp.quiet)
-    {
-      gettimeofday(&t1, NULL);
-      printf(PROGRAM_NAME" started on %s", ctime(&p->rawtime));
-    }
-
-
   /* Allocate space for all the input images. This is done here because
      WCSLIB is unfortunately not thread-safe when reading the WCS
      information from the FITS files. In cases where the number of cropped
@@ -663,8 +818,8 @@ ui_preparations(struct cropparams *p)
   errno=0;
   p->imgs=malloc(p->numin*sizeof *p->imgs);
   if(p->imgs==NULL)
-    error(EXIT_FAILURE, errno, "ui.c: %zu bytes for p->imgs",
-          p->numin*sizeof *p->imgs);
+    error(EXIT_FAILURE, errno, "%s: %zu bytes for p->imgs",
+          __func__, p->numin*sizeof *p->imgs);
 
 
   /* Fill in the WCS information of each image. */
@@ -690,23 +845,37 @@ ui_preparations(struct cropparams *p)
       else
         if(p->mode==IMGCROP_MODE_WCS)
           error(EXIT_FAILURE, 0, "the WCS structure of %s (hdu: %s) "
-                "image is not recognized. So RA and Dec cannot be used "
-                "as input. You can try with pixel coordinates in the "
-                "Image Mode (note that the crops will lack WCS "
-                "header information)", img->name, p->cp.hdu);
+                "image is not recognized. So WCS mode cannot be used "
+                "as input coordinates. You can try with pixel coordinates "
+                "with `--mode=img'", img->name, p->cp.hdu);
       fits_close_file(tmpfits, &status);
       gal_fits_io_error(status, NULL);
 
-      /* Make sure all the images have the same type. */
+      /* Make sure all the images have the same type and dimensions. */
       if(firsttype==0)
         {
-          firsttype=p->type;
+          /* Set the basic information. */
+          firsttype = p->type;
+          firstndim = img->ndim;
           p->bitnul = gal_blank_alloc_write(p->type);
+
+          /* Make sure the number of dimensions is supported. */
+          if(firstndim>MAXDIM)
+            error(EXIT_FAILURE, 0, "%s: is as %zu dimensional dataset, Crop "
+                  "currently only supports a maximum of %d dimensions",
+                  img->name, firstndim, MAXDIM);
+
+          /* Make sure the number of coordinates given for center
+             correspond to the dimensionality of the data. */
+          if(p->center && p->center->size!=firstndim)
+            error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, but "
+                  "%zu coordinates were given to `--center'", img->name,
+                  p->cp.hdu, firstndim, p->center->size);
         }
       else
         {
           if(firsttype!=p->type)
-            error(EXIT_FAILURE, 0, "%s: type is `%s' while revious image(s) "
+            error(EXIT_FAILURE, 0, "%s: type is `%s' while previous input(s) "
                   "were `%s' type. All inputs must have the same pixel data "
                   "type.\n\nYou can use Gnuastro's Arithmetic program to "
                   "convert `%s' to `%s', please run this command for more "
@@ -716,21 +885,33 @@ ui_preparations(struct cropparams *p)
                   img->name, gal_type_name(p->type, 1),
                   gal_type_name(firsttype, 1), img->name,
                   gal_type_name(p->type, 1));
+          if(firstndim!=img->ndim)
+            error(EXIT_FAILURE, 0, "%s: type has %zu dimensions, while "
+                  "previous input(s) had %zu dimensions. All inputs must "
+                  "have the same number of dimensions", img->name, img->ndim,
+                  firstndim);
         }
 
-      /* In WCS mode, Check resolution and get the first pixel
-         positions. */
-      if(p->mode==IMGCROP_MODE_WCS) wcs_check_prepare(p, img);
+      /* In WCS mode, we need some additional preparations. */
+      if(p->mode==IMGCROP_MODE_WCS) wcsmode_check_prepare(p, img);
     }
 
+  /***************************************************/
+  /********** Until 3D is fully implemented **********/
+  if(p->imgs->ndim!=2 && p->polygon)
+    error(EXIT_FAILURE, 0, "%s: currently only 2D datasets are "
+          "usable with polygon cropping", p->imgs->name);
+  /***************************************************/
+
+
+  /* Unify central crop methods into `p->centercoords'. */
+  if(p->catname || p->center)
+    ui_prepare_center(p);
 
-  /* Report timing: */
-  if(!p->cp.quiet)
-    {
-      asprintf(&msg, "Read metadata of %zu image%s.", p->numin,
-              p->numin>1 ? "s" : "");
-      gal_timing_report(&t1, msg, 1);
-    }
+
+  /* `ui_read_cols' set the number of output crops when a catalog was
+     given, in all other cases, we only have one output. */
+  if(p->catname==NULL) p->numout=1;
 
 
   /* Prepare the log file if the user has asked for it. */
@@ -762,6 +943,8 @@ ui_preparations(struct cropparams *p)
 void
 ui_read_check_inputs_setup(int argc, char *argv[], struct cropparams *p)
 {
+  char *msg;
+  struct timeval t1;
   struct gal_options_common_params *cp=&p->cp;
 
 
@@ -807,8 +990,28 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
cropparams *p)
   ui_check_options_and_arguments(p);
 
 
+  /* To see how long it takes to read meta-data. */
+  if(!p->cp.quiet) gettimeofday(&t1, NULL);
+
+
   /* Read/allocate all the necessary starting arrays. */
   ui_preparations(p);
+
+
+  /* Report timing: */
+  if(!p->cp.quiet)
+    {
+      printf(PROGRAM_NAME" started on %s", ctime(&p->rawtime));
+      asprintf(&msg, "Read metadata of %zu dataset%s.", p->numin,
+              p->numin>1 ? "s" : "");
+      gal_timing_report(&t1, msg, 1);
+      if(p->numout>1)
+        {
+          asprintf(&msg, "Will try making %zu crops (from catalog).",
+                   p->numout);
+          gal_timing_report(NULL, msg, 1);
+        }
+    }
 }
 
 
@@ -839,8 +1042,7 @@ ui_free_report(struct cropparams *p, struct timeval *t1)
   size_t i;
 
   /* Free the simple arrays (if they were set). */
-  if(p->c1) free(p->c1);
-  if(p->c2) free(p->c2);
+  if(p->center) free(p->center);
   if(p->cp.hdu) free(p->cp.hdu);
   if(p->cathdu) free(p->cathdu);
   if(p->catname) free(p->catname);
diff --git a/bin/crop/ui.h b/bin/crop/ui.h
index fdb9a70..c59a22d 100644
--- a/bin/crop/ui.h
+++ b/bin/crop/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 /* Available letters for short options:
 
-   e k m t u v
+   a d e f g i j k m r t u v y
    A B E G H J L Q R W X Y
 */
 enum option_keys_enum
@@ -37,23 +37,15 @@ enum option_keys_enum
   /* With short-option version. */
   UI_KEY_CATALOG        = 'C',
   UI_KEY_NOBLANK        = 'b',
-  UI_KEY_CHECKCENTER    = 'c',
   UI_KEY_SUFFIX         = 'p',
   UI_KEY_NAMECOL        = 'n',
-  UI_KEY_RACOL          = 'f',
-  UI_KEY_DECCOL         = 'g',
-  UI_KEY_RA             = 'r',
-  UI_KEY_DEC            = 'd',
-  UI_KEY_XCOL           = 'i',
-  UI_KEY_YCOL           = 'j',
-  UI_KEY_XC             = 'x',
-  UI_KEY_YC             = 'y',
-  UI_KEY_IWIDTH         = 'a',
-  UI_KEY_WWIDTH         = 'w',
   UI_KEY_SECTION        = 's',
   UI_KEY_POLYGON        = 'l',
   UI_KEY_ZEROISNOTBLANK = 'z',
   UI_KEY_MODE           = 'O',
+  UI_KEY_WIDTH          = 'w',
+  UI_KEY_CENTER         = 'c',
+  UI_KEY_COORDCOL       = 'x',
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
@@ -61,6 +53,7 @@ enum option_keys_enum
   UI_KEY_HSTARTWCS,
   UI_KEY_HENDWCS,
   UI_KEY_OUTPOLYGON,
+  UI_KEY_CHECKCENTER,
 };
 
 
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index 34179a9..7943638 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <error.h>
 #include <stdio.h>
 #include <float.h>
+#include <string.h>
 #include <stdlib.h>
 
 #include <gnuastro/wcs.h>
@@ -48,12 +49,17 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* This function is called from ui.c. Its job is to check the WCS
    values of this  */
 void
-wcs_check_prepare(struct cropparams *p, struct inputimgs *img)
+wcsmode_check_prepare(struct cropparams *p, struct inputimgs *img)
 {
-  double twidth, *pixscale;
+  double *pixscale;
   struct wcsprm *wcs=img->wcs;
-  int i, status[4]={0,0,0,0}, ncoord=4, nelem=2;
-  double imgcrd[8], phi[4], theta[4], pixcrd[8];
+  size_t *dsize=img->dsize, ndim=img->ndim;
+
+  /* For two dimensions, we have four corners (2 numbers for each) and for
+     three dimensions we have 8 corners (3 numbers for each), so we'll just
+     assume the largest. */
+  int i, status[8]={0,0,0,0,0,0,0,0}, ncorners=0;
+  double imgcrd[24], phi[8], theta[8], pixcrd[24];
 
 
   /* Check if the image is aligned with the WCS coordinates. Note that
@@ -81,97 +87,143 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs 
*img)
           img->name, p->cp.hdu);
 
 
-  /* Check the image resolution (if the pixels are actually a square), then
-     compare the resolution with the other input images. Due to floating
-     point errors, some very small differences might exist in the pixel
-     scale, so break out with an error only if the pixel scales are more
-     different than 1e-6. */
-  pixscale=gal_wcs_pixel_scale_deg(wcs);
+  /* Check the nature of the coordinates, currently we can only support RA
+     and Dec, other modes haven't been checked. */
+  if( strcmp(wcs->ctype[0], "RA---TAN")
+      || strcmp(wcs->ctype[1], "DEC--TAN") )
+    error(EXIT_FAILURE, 0, "currently the only WCS types usable are "
+          "`RA---TAN' and `DEC--TAN' for the first and second axises "
+          "respectively. The WCS types of `%s' (hdu %s) are `%s' and `%s' "
+          "respectively", img->name, p->cp.hdu, wcs->ctype[0], wcs->ctype[1]);
+
+
+  /* Check if the pixels are actually a square, then compare the resolution
+     with the other input images. Due to floating point errors, some very
+     small differences might exist in the pixel scale, so break out with an
+     error only if the pixel scales are more different than 1e-6. */
+  pixscale=gal_wcs_pixel_scale(wcs);
   if( fabs(pixscale[0]-pixscale[1])/pixscale[0] > 1e-6 )
     error(EXIT_FAILURE, 0, "%s: HDU %s: The pixel scale along "
           "the two image axises is not the same. The first axis "
           "is %.15g deg/pixel, while the second is %.15g",
           img->name, p->cp.hdu, pixscale[1], pixscale[0]);
-  if(p->res==0.0f)
+  if(p->pixscale)
     {
-      /* Set the resolution of the image. */
-      p->res=pixscale[0];
+      for(i=0;i<ndim;++i)
+        if(p->pixscale[i] != pixscale[i])
+          error(EXIT_FAILURE, 0, "%s (hdu %s): has resolution of %g along "
+                "dimension %d. However, previously checked input(s) had "
+                "a resolution of %g in this dimension", img->name, p->cp.hdu,
+                pixscale[i], i+1, p->pixscale[i]);
       free(pixscale);
-
-      /* Set the widths such that iwidth and wwidth are exactly the same
-         (within their different units ofcourse). Also make sure that the
-         image size is an odd number (so the central pixel is in the
-         center). */
-      p->wwidth/=3600;                 /* Convert the width to degrees. */
-      twidth=p->wwidth/p->res;
-      if(twidth<3)
-        error(EXIT_FAILURE, 0, "--wwidth = %f (arcseconds) translates "
-              "to %.0f pixels in scale of input image(s). This is probably "
-              "not what you want", p->wwidth*3600, twidth);
-      p->iwidth[0] = (twidth-(long)twidth)>0.5 ? twidth+1 : twidth;
-      if(p->iwidth[0]%2==0)
-        {
-          p->iwidth[0]+=1;
-          p->wwidth+=p->res;
-        }
-      p->iwidth[1]=p->iwidth[0];
     }
   else
+    p->pixscale=pixscale;
+
+
+  /* Set the coordinates of the dataset's corners. Note that `dsize' is in
+     C order, while pixcrd is in FITS order.*/
+  switch(ndim)
     {
-      if(p->res!=pixscale[0])
-        error(EXIT_FAILURE, 0, "%s: HDU %s: The resolution of "
-              "this image is %f arcseconds/pixel while the "
-              "previously checked input image(s) had a resolution "
-              "of %f", img->name, p->cp.hdu, 3600*wcs->pc[3],
-              3600*p->res);
-      free(pixscale);
+    case 2:
+      ncorners=4;
+      pixcrd[0] = 1;          pixcrd[1] = 1;
+      pixcrd[2] = dsize[1];   pixcrd[3] = 1;
+      pixcrd[4] = 1;          pixcrd[5] = dsize[0];
+      pixcrd[6] = dsize[1];   pixcrd[7] = dsize[0];
+      break;
+    case 3:
+      ncorners=8;
+      pixcrd[0]  = 1;         pixcrd[1]  = 1;         pixcrd[2]  = 1;
+      pixcrd[3]  = dsize[2];  pixcrd[4]  = 1;         pixcrd[5]  = 1;
+      pixcrd[6]  = 1;         pixcrd[7]  = dsize[1];  pixcrd[8]  = 1;
+      pixcrd[9]  = dsize[2];  pixcrd[10] = dsize[1];  pixcrd[11] = 1;
+      pixcrd[12] = 1;         pixcrd[13] = 1;         pixcrd[14] = dsize[0];
+      pixcrd[15] = dsize[2];  pixcrd[16] = 1;         pixcrd[17] = dsize[0];
+      pixcrd[18] = 1;         pixcrd[19] = dsize[1];  pixcrd[20] = dsize[0];
+      pixcrd[21] = dsize[2];  pixcrd[22] = dsize[1];  pixcrd[23] = dsize[0];
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: %zu dimensional datasets not supported",
+            __func__, ndim);
     }
 
 
-  /* Get the coordinates of the first pixel in the image. Note that `dsize'
-     is in C axises, while pixcrd is in FITS axises. */
-  pixcrd[0]=1;               pixcrd[1]=1;
-  pixcrd[2]=img->dsize[1];   pixcrd[3]=1;
-  pixcrd[4]=1;               pixcrd[5]=img->dsize[0];
-  pixcrd[6]=img->dsize[1];   pixcrd[7]=img->dsize[0];
-  wcsp2s(wcs, ncoord, nelem, pixcrd, imgcrd, phi, theta,
+  /* Get the coordinates of the corners of the dataset in WCS.  */
+  wcsp2s(wcs, ncorners, ndim, pixcrd, imgcrd, phi, theta,
          img->corners, status);
 
+
   /* Check if there was no error in the conversion. */
-  for(i=0;i<4;++i)
+  for(i=0;i<ncorners;++i)
     if(status[i])
       error(EXIT_FAILURE, 0, "wcsp2s ERROR %d in row %d of pixcrd: %s",
             i, status[i], wcs_errmsg[status[i]]);
 
 
-  /* Fill in the size of the image in celestial degrees from the first
-     pixel in the image. Note that `dsize' is in C axises, while pixcrd is
-     in FITS axises. */
-  img->sized[0]=img->dsize[1]*p->res/cos(img->corners[1]*M_PI/180);
-  img->sized[1]=img->dsize[0]*p->res;
+  /* Fill in the size of the dataset in WCS from the first pixel in the
+     image. Note that `dsize' is in C axises, while the `pixscale',
+     `corners' and `sized' are in FITS axises. */
+  if(ndim==2)
+    {
+      img->sized[0] = ( img->dsize[1] * p->pixscale[0]
+                        / cos( img->corners[1] * M_PI / 180 ) );
+      img->sized[1] = img->dsize[0] * p->pixscale[1];
+    }
+  else /* 3D */
+    {
+      img->sized[0] = ( img->dsize[2] * p->pixscale[0]             /* RA  */
+                        / cos( img->corners[1] * M_PI / 180 ) );
+      img->sized[1] = img->dsize[1] * p->pixscale[1];              /* Dec */
+      img->sized[2] = img->dsize[0] * p->pixscale[2];              /* 3D  */
+    }
+
 
+  /* In case the image crosses the equator, we will calculate these values
+     here so later on, we don't have to calculate them on every check. See
+     the explanation above `point_in_dataset'.
 
-  /* In case the image crosses the equator, we will calculate these
-     values here so later on, we don't have to calculate them on every
-     check. See the explanation above radecoverlap.*/
-  if( img->corners[1]*(img->corners[1]+img->sized[1]) < 0 )
+     Note that in both 2D and 3D data, the declination is in the second
+     coordinate (index 1). */
+  if( img->corners[1] * (img->corners[1]+img->sized[1]) < 0 )
     {
-      /* re in the explanations. */
+      /* re in the comments above `point_in_dataset'. */
       img->equatorcorr[0]=img->corners[0]
         -0.5*img->sized[0]*(1-cos(img->corners[1]*M_PI/180));
 
-      /* sre in the explanations. */
+      /* sre in the comments above `point_in_dataset'. */
       img->equatorcorr[1]=img->sized[0]*cos(img->corners[1]*M_PI/180);
     }
 
 
   /* Just to check:
-  printf("\n\n%s:\n(%.10f, %.10f)\n(%.10f, %.10f)"
-         "\n(%.10f, %.10f)\n(%.10f, %.10f)\n\n", img->name,
-         img->corners[0], img->corners[1],
-         img->corners[2], img->corners[3],
-         img->corners[4], img->corners[5],
-         img->corners[6], img->corners[7]);
+  printf("\n\n%s:\n", img->name);
+  if(ndim==2)
+    printf("(%.10f, %.10f)\n"
+           "(%.10f, %.10f)\n"
+           "(%.10f, %.10f)\n"
+           "(%.10f, %.10f)\n\n",
+           img->corners[0], img->corners[1],
+           img->corners[2], img->corners[3],
+           img->corners[4], img->corners[5],
+           img->corners[6], img->corners[7]);
+  else
+    printf("(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n"
+           "(%.10f, %.10f, %.10f)\n\n",
+           img->corners[0],  img->corners[1],  img->corners[2],
+           img->corners[3],  img->corners[4],  img->corners[5],
+           img->corners[6],  img->corners[7],  img->corners[8],
+           img->corners[9],  img->corners[10], img->corners[11],
+           img->corners[12], img->corners[13], img->corners[14],
+           img->corners[15], img->corners[16], img->corners[17],
+           img->corners[18], img->corners[19], img->corners[20],
+           img->corners[21], img->corners[22], img->corners[23] );
   exit(0);
   */
 }
@@ -200,23 +252,30 @@ wcs_check_prepare(struct cropparams *p, struct inputimgs 
*img)
 /*******************************************************************/
 /* Set the four sides around the point of interest in RA and Dec.
 
-   NOTE: In this format we are working on here (where the image is
-   aligned with the celestial coordinates), the declination is
-   measured on a great circle, while the right ascension is not. So we
-   have to consider the
+   NOTE: When the image is aligned with the celestial coordinates (current
+   working paradigm), the declination is measured on a great circle, while
+   the right ascension is not. So we have to consider this in calculating
+   the difference in RA.
 */
 void
-setcsides(struct onecropparams *crp)
+wcsmode_crop_corners(struct onecropparams *crp)
 {
-  size_t i;
-  double h, hr, r, d, dr;        /* The second r is for radians. */
   struct cropparams *p=crp->p;
+
+  size_t i, ndim=p->imgs->ndim;
   double minra=FLT_MAX, mindec=FLT_MAX;
   double maxra=-FLT_MAX, maxdec=-FLT_MAX;
+  double r, d, l, dr, h[MAXDIM], hr[MAXDIM];
+  size_t rmini=-1, rmaxi=-1, dmini=-1, dmaxi=-1;
 
   /* Set the four corners of the WCS region. */
   if(p->polygon)
     {
+      /* A small sanity check. */
+      if(ndim!=2)
+        error(EXIT_FAILURE, 0, "%s: currently only supports 2D datasets, "
+              "your input dataset has %zu dimensions", __func__, ndim);
+
       /* Find their minimum and maximum values. */
       for(i=0;i<p->nvertices;++i)
         {
@@ -234,58 +293,137 @@ setcsides(struct onecropparams *crp)
     }
   else
     {
-      if(p->catname)
-        {
-          r=crp->world[0]=p->c1[crp->out_ind];
-          d=crp->world[1]=p->c2[crp->out_ind];
-        }
-      else
-        {
-          r=crp->world[0]=p->ra;
-          d=crp->world[1]=p->dec;
-        }
+      /* Set the RA and Dec to use as center. */
+      r=crp->world[0]=p->centercoords[0][crp->out_ind];
+      d=crp->world[1]=p->centercoords[1][crp->out_ind];
+      if(ndim==3) l=crp->world[2]=p->centercoords[2][crp->out_ind];
+
 
-      h=p->wwidth/2;
+      /* Calculate the declination in radians for easy readability. */
       dr=d*M_PI/180;
-      hr=h*M_PI/180;
 
-      /* Set the four corners of this crop. */
-      crp->corners[0] = r+h/cos(dr-hr);  crp->corners[1] = d-h; /* Bt Lf */
-      crp->corners[2] = r-h/cos(dr-hr);  crp->corners[3] = d-h; /* Bt Rt */
-      crp->corners[4] = r+h/cos(dr+hr);  crp->corners[5] = d+h; /* Tp Lf */
-      crp->corners[6] = r-h/cos(dr+hr);  crp->corners[7] = d+h; /* Tp Rt */
+      /* Set the half width in each dimension. For the angular dimensions,
+         also calculate it in radians. */
+      hr[0] = ( h[0] = ((double *)(p->width->array))[0] / 2 ) * M_PI / 180;
+      hr[1] = ( h[1] = ((double *)(p->width->array))[1] / 2 ) * M_PI / 180;
+      if(ndim==3)
+        h[2] = ((double *)(p->width->array))[2] / 2;
+
+      /* Set the corners of this crop. */
+      switch(ndim)
+        {
+        case 2:
+          crp->corners[0] = r+h[0]/cos(dr-hr[1]);
+          crp->corners[1] = d-h[1];                   /* Bottom left.  */
+
+          crp->corners[2] = r-h[0]/cos(dr-hr[1]);
+          crp->corners[3] = d-h[1];                   /* Bottom Right. */
+
+          crp->corners[4] = r+h[0]/cos(dr+hr[1]);
+          crp->corners[5] = d+h[1];                   /* Top Left.     */
+
+          crp->corners[6] = r-h[0]/cos(dr+hr[1]);
+          crp->corners[7] = d+h[1];                   /* Top Right.    */
+          break;
+
+        case 3:
+          /* Note that the third dimension is assumed to be independent of
+             the first two. So the first two coordinates of its corners in
+             the front and back (on the two faces in the third dimension),
+             are equal.  */
+          crp->corners[0]  = crp->corners[12] = r+h[0]/cos(dr-hr[1]);
+          crp->corners[1]  = crp->corners[13] = d-h[1];
+          crp->corners[2]  = l-h[2];                /* Bottom left front. */
+
+          crp->corners[3]  = crp->corners[15] = r-h[0]/cos(dr-hr[1]);
+          crp->corners[4]  = crp->corners[16] = d-h[1];
+          crp->corners[5]  = l-h[2];                /* Bottom right front.*/
+
+          crp->corners[6]  = crp->corners[18] = r+h[0]/cos(dr+hr[1]);
+          crp->corners[7]  = crp->corners[19] = d+h[1];
+          crp->corners[8]  = l-h[2];                /* Top Left front.    */
+
+          crp->corners[9]  = crp->corners[21] = r-h[0]/cos(dr+hr[1]);
+          crp->corners[10] = crp->corners[22] = d+h[1];
+          crp->corners[11] = l-h[2];                /* Top right front.   */
+
+          crp->corners[14] = l+h[2];                /* Bottom left back.  */
+          crp->corners[17] = l+h[2];                /* Bottom right back. */
+          crp->corners[20] = l+h[2];                /* Top left back.     */
+          crp->corners[23] = l+h[2];                /* Top right back.     */
+          break;
+        }
     }
 
-  /* Set the bottom width and height of the crop in degrees. Note that
-     the width changes as the height changes, so here we want the
-     height and the lowest declination. Note that on the bottom edge,
-     corners[0] is the maximum RA and corners[2] is the minimum RA.
-     For all the region, corners[5] is one of the maximum declinations
-     and corners[3] is one of the the minimum declinations.*/
-  crp->sized[0]=( (crp->corners[0]-crp->corners[2])
-                  / cos(crp->corners[1]*M_PI/180) );
-  crp->sized[1]=crp->corners[5]-crp->corners[3];
+
+  /* Set the bottom width and height of the crop in degrees. Note that the
+     width changes as the height changes, so here we want the height and
+     the lowest declination. Note that in 2D on the bottom edge, corners[0]
+     is the maximum RA and corners[2] is the minimum RA.  For all the 2D
+     region, corners[5] is one of the maximum declinations and corners[1]
+     is one of the the minimum declinations.
+
+     North and south hemispheres are no problem: When using the center,
+     they are set properly (in any hemisphere) and for a polygon, the
+     minimums and maximums are automatically found. */
+  rmini = ndim;                 /* First element in second corner. */
+  rmaxi = 0;                    /* First element.                  */
+  dmini = 1;                    /* Second element.                 */
+  dmaxi = ndim==2 ? 5 : 7;      /* Second element in third corner. */
+  crp->sized[0]=( (crp->corners[rmaxi]-crp->corners[rmini])
+                  / cos(crp->corners[dmini]*M_PI/180) );
+  crp->sized[1]=crp->corners[dmaxi]-crp->corners[dmini];
+  if(ndim==3)
+    crp->sized[2] = crp->corners[14] - crp->corners[2];
+
 
   /* In case the crop crosses the equator, then we need these two
-     corrections. See the complete explanations below. */
+     corrections. See the complete explanations above `point_in_dataset'. */
   if(crp->corners[1]*(crp->corners[1]+crp->sized[1]) < 0 )
     {
-      /* re in the explanations. */
+      /* re in the explanations above `point_in_dataset'. */
       crp->equatorcorr[0]=crp->corners[0]
         -0.5*crp->sized[0]*(1-cos(crp->corners[1]*M_PI/180));
 
-      /* sre in the explanations. */
+      /* sre in the explanations above `point_in_dataset'. */
       crp->equatorcorr[1]=crp->sized[0]*cos(crp->corners[1]*M_PI/180);
     }
 
+
   /* Just to check:
-  printf("\n\nCorner 1: (%.10f, %.10f)\n"
-         "Corner 2: (%.10f, %.10f)\nCorner 3: (%.10f, %.10f)\n"
-         "Corner 4: (%.10f, %.10f)\n\n",
-         crp->corners[0], crp->corners[1],
-         crp->corners[2], crp->corners[3],
-         crp->corners[4], crp->corners[5],
-         crp->corners[6], crp->corners[7]);
+  if(ndim==2)
+    {
+      printf("\n\n%g, %g:\n", r, d);
+      printf("\t(%.10f, %.10f)\n"
+             "\t(%.10f, %.10f)\n"
+             "\t(%.10f, %.10f)\n"
+             "\t(%.10f, %.10f)\n\n",
+             crp->corners[0], crp->corners[1],
+             crp->corners[2], crp->corners[3],
+             crp->corners[4], crp->corners[5],
+             crp->corners[6], crp->corners[7]);
+    }
+  else
+    {
+      printf("\n\n%g, %g, %g:\n", r, d, l);
+      printf("\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n"
+             "\t(%.10f, %.10f, %g)\n\n",
+             crp->corners[0],  crp->corners[1],  crp->corners[2],
+             crp->corners[3],  crp->corners[4],  crp->corners[5],
+             crp->corners[6],  crp->corners[7],  crp->corners[8],
+             crp->corners[9],  crp->corners[10], crp->corners[11],
+             crp->corners[12], crp->corners[13], crp->corners[14],
+             crp->corners[15], crp->corners[16], crp->corners[17],
+             crp->corners[18], crp->corners[19], crp->corners[20],
+             crp->corners[21], crp->corners[22], crp->corners[23] );
+    }
+  exit(0);
   */
 }
 
@@ -434,21 +572,36 @@ fillcrpipolygon(struct onecropparams *crp)
    (South) equations of above as before and for those points that have
    a positive declination,  we use the North formula  but replacing r1
    with re, d1 with 0 and sr with sre.
-*/
-/*
-   p[2]: RA and Dec of a point (rp and dp above).
-   i[2]: RA and Dec of first point in image. (r1 and d1 above).
-   s[2]: Size of the image in degrees (sr and sd above).
-   c[2]: Corrections if equator is passed, (se and sre above).
-*/
-int
-radecinimg(double *p, double *i, double *s, double *c)
+
+
+   INPUTS
+   ------
+   p[]: Point coordinates (rp and dp above).
+   i[]: Coordinates of first pixel in image. (r1 and d1 above).
+   s[]: Size/width of box (sr and sd above).
+   c[]: Corrections if equator is passed, (se and sre above).
+
+
+   IMPORTANT: It is assumed that the dimensions are ordered with:
+
+      0: RA
+      1: Dec
+      2: Third dimension (independent of RA and Dec).          */
+static int
+point_in_dataset(double *p, double *i, double *s, double *c, size_t ndim)
 {
   double n;
 
-  /* First check the declination. If it is not in range, you can
-     safely return 0.*/
-  if(p[1]>=i[1]  && p[1]<=i[1]+s[1])
+
+  /* If there is a third dimension, then first check that. Note that the
+     third dimension is assumed to be indendent of the first two. */
+  if(ndim==3 && ( p[2]<i[2] || p[2]>i[2]+s[2] ) )
+    return 0;
+
+
+  /* In the RA and Dec checks, first check the declination. If it is not in
+     range, you can safely return 0. */
+  if(p[1]>=i[1] && p[1]<=i[1]+s[1])
     {
       if(p[1]<=0)            /* Point is in southern hemisphere, it     */
         {                    /* doesn't matter if image passes equator! */
@@ -491,21 +644,25 @@ radecinimg(double *p, double *i, double *s, double *c)
    there is an overlap (the survey image is completely within the crop
    box). So we have to check both.  */
 int
-radecoverlap(struct onecropparams *crp)
+wcsmode_overlap(struct onecropparams *crp)
 {
   double *d, *fd;
   double *i, *s, *c;                /* for clear viewing. */
   struct cropparams *p=crp->p;
+  size_t ndim=crp->p->imgs->ndim;
 
-  /* First check if the four sides of the crop box are in the image.*/
-  fd=(d=crp->corners)+8;
+  /* First check if the corners of the crop are in the image.*/
   s=p->imgs[crp->in_ind].sized;
   i=p->imgs[crp->in_ind].corners;
   c=p->imgs[crp->in_ind].equatorcorr;
+  fd=(d=crp->corners) + (ndim==2 ? 8 : 24);
   do
     {
-      if( radecinimg(d, i, s, c) ) return 1;
-      d+=2;
+      /* As long as one of the crop corners are in the image, we know there
+         is overlap and can return a true value. We don't need to check all
+         corners. */
+      if( point_in_dataset(d, i, s, c, ndim) ) return 1;
+      d+=ndim;
     }
   while(d<fd);
 
@@ -514,13 +671,14 @@ radecoverlap(struct onecropparams *crp)
   s=crp->sized;
   i=crp->corners;
   c=crp->equatorcorr;
-  fd=(d=p->imgs[crp->in_ind].corners)+8;
+  fd=(d=p->imgs[crp->in_ind].corners) + (ndim==2 ? 8 : 24);
   do
     {
-      if( radecinimg(d, i, s, c) ) return 1;
-      d+=2;
+      if( point_in_dataset(d, i, s, c, ndim) ) return 1;
+      d+=ndim;
     }
   while(d<fd);
 
+  /* If control reaches here, there was no overlap. */
   return 0;
 }
diff --git a/bin/crop/wcsmode.h b/bin/crop/wcsmode.h
index 2866ddd..ead7fff 100644
--- a/bin/crop/wcsmode.h
+++ b/bin/crop/wcsmode.h
@@ -24,15 +24,15 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #define WCSMODE_H
 
 void
-wcs_check_prepare(struct cropparams *p, struct inputimgs *img);
+wcsmode_check_prepare(struct cropparams *p, struct inputimgs *img);
 
 void
-setcsides(struct onecropparams *crp);
+wcsmode_crop_corners(struct onecropparams *crp);
 
 void
 fillcrpipolygon(struct onecropparams *crp);
 
 int
-radecoverlap(struct onecropparams *crp);
+wcsmode_overlap(struct onecropparams *crp);
 
 #endif
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index d768b1d..ddcc89a 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -218,6 +218,7 @@ mkprof_build(void *inparam)
 
   size_t i, id;
   int lockresult;
+  double center[2];
   long lpixel_o[2];
   pthread_mutex_t *qlock=&p->qlock;
   struct builtqueue *ibq, *fbq=NULL;
@@ -252,13 +253,15 @@ mkprof_build(void *inparam)
 
       /* Get the overlapping pixels using the starting points (NOT
          oversampled). */
-      gal_box_border_from_center(p->x[id], p->y[id], mkp->width,
+      center[0]=p->x[id];
+      center[1]=p->y[id];
+      gal_box_border_from_center(center, 2, mkp->width,
                                  ibq->fpixel_i, ibq->lpixel_i);
       mkp->fpixel_i[0]=ibq->fpixel_i[0];
       mkp->fpixel_i[1]=ibq->fpixel_i[1];
       ibq->overlaps = gal_box_overlap(mkp->onaxes, ibq->fpixel_i,
                                       ibq->lpixel_i, ibq->fpixel_o,
-                                      lpixel_o);
+                                      lpixel_o, 2);
 
 
       /* Build the profile if necessary, After this, the width is
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 2f7758a..bc36fdb 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -465,7 +465,7 @@ ui_matrix_make_align(struct warpparams *p, double *tmatrix)
   /* Find the pixel scale along the two dimensions. Note that we will be
      using the scale along the image X axis for both values. */
   w=gal_wcs_warp_matrix(p->input->wcs);
-  ps=gal_wcs_pixel_scale_deg(p->input->wcs);
+  ps=gal_wcs_pixel_scale(p->input->wcs);
 
 
   /* Lets call the given WCS orientation `W', the rotation matrix we want
diff --git a/bin/warp/warp.c b/bin/warp/warp.c
index ff53451..4584ec4 100644
--- a/bin/warp/warp.c
+++ b/bin/warp/warp.c
@@ -455,7 +455,7 @@ correct_wcs_save_output(struct warpparams *p)
      signs are usually different.*/
   if( wcs->pc[1]<ABSOLUTEFLTERROR ) wcs->pc[1]=0.0f;
   if( wcs->pc[2]<ABSOLUTEFLTERROR ) wcs->pc[2]=0.0f;
-  pixelscale=gal_wcs_pixel_scale_deg(wcs);
+  pixelscale=gal_wcs_pixel_scale(wcs);
   diff=fabs(wcs->pc[0])-fabs(wcs->pc[3]);
   if( fabs(diff/pixelscale[0])<RELATIVEFLTERROR )
     wcs->pc[3] =  ( (wcs->pc[3] < 0.0f ? -1.0f : 1.0f) * fabs(wcs->pc[0]) );
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c583877..34c83b6 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -2521,7 +2521,7 @@ you have downloaded
 $ tar xf cfitsio_latest.tar.gz
 $ cd cfitsio
 $ ./configure --prefix=/usr/local --enable-sse2 --enable-reentrant
-$ make -j8                         # Replace 8 with no. CPU threads.
+$ make
 $ make utils
 $ ./testprog > testprog.lis
 $ diff testprog.lis testprog.out    # Should have no output
@@ -7697,55 +7697,67 @@ image, see @ref{Crop section syntax}.
 
 @node Crop modes, Crop section syntax, Crop, Crop
 @subsection Crop modes
-In order to be comprehensive, intuitive, and easy to use, Crop has two ways
-to define crop region: 1) From its center and (square) side length, for
-example to generate postage stamps of a given catalog. 2) The vertices of
-the crop region, this can be useful for larger crops over many targets, for
-example to crop out a uniformly deep, or contiguous, region of a large
-survey. Irrespective of how the crop region is defined, both Image/pixel,
-and World coordinate system (WCS) coordinates are acceptable and all
-coordinates are read as floating point numbers (not integers, except for
-the @option{--section} option, see below). The coordinate standards used
-are the main modes of Crop. Here, the different ways to specify the crop
-region are discussed within each standard. For the full list options,
-please see @ref{Invoking astcrop}.
+In order to be comprehensive, intuitive, and easy to use, there are two
+ways to define the crop:
+
+@enumerate
+@item
+From its center and side length. For example if you already know the
+coordinates of an object and want to inspect it in an image or to generate
+postage stamps of a catalog containing many such coordinates.
+
+@item
+The vertices of the crop region, this can be useful for larger crops over
+many targets, for example to crop out a uniformly deep, or contiguous,
+region of a large survey.
+@end enumerate
+
+Irrespective of how the crop region is defined, the coordinates to define
+the crop can be in Image (pixel) or World Coordinate System (WCS)
+standards. All coordinates are read as floating point numbers (not
+integers, except for the @option{--section} option, see below). By setting
+the @emph{mode} in Crop, you define the standard that the given coordinates
+must be interpretted. Here, the different ways to specify the crop region
+are discussed within each standard. For the full list options, please see
+@ref{Invoking astcrop}.
 
 When the crop is defined by its center, the respective (integer) central
 pixel position will be found internally according to the FITS standard. To
 have this pixel positioned in the center of the cropped region, the final
-cropped region must have (in Image mode), or will have (in WCS mode) an add
-number of pixels. Furthermore, when the crop is defined as by its center,
-Crop allows you to only keep crops what don't have any blank pixels in the
-vicinity of their center (your primary target). This can be very convenient
-when your input catalog/coordinates originated from another survey/filter
-which is not fully covered by your input image, to learn more about this
-feature, please see the description of the @option{--checkcenter} option in
-@ref{Invoking astcrop}.
+cropped region will have an add number of pixels (even if you give an even
+number to @option{--width} in image mode).
+
+Furthermore, when the crop is defined as by its center, Crop allows you to
+only keep crops what don't have any blank pixels in the vicinity of their
+center (your primary target). This can be very convenient when your input
+catalog/coordinates originated from another survey/filter which is not
+fully covered by your input image, to learn more about this feature, please
+see the description of the @option{--checkcenter} option in @ref{Invoking
+astcrop}.
 
 @table @asis
 @item Image coordinates
-In image mode, Crop uses the pixel coordinates/positions (instead of world
-coordinates). In image mode, only one image may be input. The crop(s) can
-be defined in multiple ways as listed below.
+In image mode (@option{--mode=img}), Crop interprets the pixel coordinates
+and widths in units of the input data-elements (for example pixels in an
+image, not world coordinates). In image mode, only one image may be
+input. The output crop(s) can be defined in multiple ways as listed below.
 
 @table @asis
 @item Center of multiple crops (in a catalog)
-The center of (possibly multiple) crops are read from a text file. A
-catalog can also contain WCS coordinates, so @option{--mode=img} is
-necessary to avoid ambiguity when all columns are specified. The
-@option{--xcol} and @option{--ycol} columns in the catalog are interpreted
-as the center of a square crop with a width of @option{--iwidth} pixels (an
-odd number). The columns can contain any floating point value. The value to
-@option{--output} option is seen as a directory which will host (the
-possibly multiple) separate crop files, see @ref{Crop output} for more. For
-a tutorial using this feature, please see @ref{Hubble visually checks and
-classifies his catalog}.
+The center of (possibly multiple) crops are read from a text file. In this
+mode, the columns identified with the @option{--coordcol} option are
+interpreted as the center of a crop with a width of @option{--width} pixels
+along each dimension. The columns can contain any floating point value. The
+value to @option{--output} option is seen as a directory which will host
+(the possibly multiple) separate crop files, see @ref{Crop output} for
+more. For a tutorial using this feature, please see @ref{Hubble visually
+checks and classifies his catalog}.
 
 @item Center of a single crop (on the command-line)
-The center of the crop is given on the command-line with the @option{--xc}
-and @option{--yc} options. The crop width is specified by the
-@option{--iwidth} option. The given coordinates can be any floating point
-and the width has to be an odd integer, see the explanations above.
+The center of the crop is given on the command-line with the
+@option{--center} option. The crop width is specified by the
+@option{--width} option along each dimension. The given coordinates and
+width can be any floating point number.
 
 @item Vertices of a single crop
 In Image mode there are two options to define the vertices of a region to
@@ -7757,17 +7769,17 @@ section syntax} for a full description of this method.
 The latter option (@option{--polygon}) is a higher-level method to define
 any convex polygon (with any number of vertices) with floating point
 values. Please see the description of this option in @ref{Invoking astcrop}
-for its syntax. It is available in both Image and WCS modes, hence, to read
-the vertices as pixel coordinates, @option{--mode=img} has to be called.
+for its syntax.
 @end table
 
 @item WCS coordinates
-Use the World Coordinate System (WCS) information to define the crop, not
-pixel coordinates. In this mode, the width (@option{--wwidth}) is read in
-units of arc-seconds and multiple images (tiles in a survey) can be
-input. When the cropped region (defined by center or vertices) overlaps
-with multiple of the input images/tiles, the overlapping regions will be
-taken from the respective input (they will be stitched in the crop).
+In WCS mode (@option{--mode=wcs}), the coordinates and widths are
+interpretted using the World Coordinate System (WCS, that must accompany
+the dataset), not pixel coordinates. In WCS mode, Crop accepts multiple
+datasets as input. When the cropped region (defined by its center or
+vertices) overlaps with multiple of the input images/tiles, the overlapping
+regions will be taken from the respective input (they will be stitched when
+necessary for each output crop).
 
 In this mode, the input images do not necessarily have to be the same size,
 they just need to have the same orientation and pixel resolution. Currently
@@ -7787,23 +7799,22 @@ listed below.
 @table @asis
 
 @item Center of multiple crops (in a catalog)
-Similar to catalog inputs in Image mode (above), but @option{--mode=wcs}
-has to be activated to avoid ambiguity. The central RA and Dec value for
-each crop will be read from the @option{--racol} and @option{--deccol}
-columns of the input catalog. The square cropped box will have an odd
-number of pixels.
+Similar to catalog inputs in Image mode (above), except that the values
+along each dimension are assumed to have the same units as the dataset's
+WCS information. For example, the central RA and Dec value for each crop
+will be read from the first and second calls to the @option{--coordcol}
+option. The width of the cropped box (in units of the WCS, or degrees in RA
+and Dec mode) must be specified with the @option{--width} option.
 
 @item Center of a single crop (on the command-line)
-You can specify the center of only one crop box with the @option{--ra} and
-@option{--dec} options. If it exists in the input images, it will be
-cropped similar to the catalog mode.
+You can specify the center of only one crop box with the @option{--center}
+option. If it exists in the input images, it will be cropped similar to the
+catalog mode, see above also for @code{--width}.
 
 @item Vertices of a single crop
 The @option{--polygon} option is a high-level method to define any convex
 polygon (with any number of vertices). Please see the description of this
-option in @ref{Invoking astcrop} for its syntax. It is available in both
-Image and WCS modes, hence, to read the vertices as RA and Dec coordinates,
-@option{--mode=wcs} has to be called/activated.
+option in @ref{Invoking astcrop} for its syntax.
 @end table
 
 @cartouche
@@ -7819,16 +7830,9 @@ input image with these coordinates, see @ref{Warp}.
 @end table
 
 As a summary, if you don't specify a catalog, you have to define the
-cropped region manually on the command-line. Other than the
-@option{--polygon} option, the other methods to define a single crop are
-unique to each mode, so the mode (@option{--mode}) will be ignored. When
-using a catalog to define the crop centers, if the columns to only one mode
-are given (for example only @option{--xcol} and @option{--ycol}, not the
-WCS columns) then the mode is irrelevant. However, when all image and WCS
-mode columns are given, @option{--mode} is mandatory (you can keep the
-columns in a configuration file and simply set the mode, see
-@ref{Configuration files} and @ref{Selecting table columns}). When using
-@option{--polygon} the mode is mandatory to interpret the values correctly.
+cropped region manually on the command-line.  In any case the mode is
+mandatory for Crop to be able to interpret the values given as coordinates
+or widths.
 
 
 @node Crop section syntax, Blank pixels, Crop modes, Crop
@@ -7939,32 +7943,32 @@ $ astcrop [OPTION...] [ASCIIcatalog] ASTRdata ...
 One line examples:
 
 @example
-## Crop all objects in catalog.txt from image.fits:
-$ astcrop catalog.txt image.fits
+## Crop all objects in cat.txt from image.fits:
+$ astcrop --catalog=cat.txt image.fits
 
 ## Crop all options in catalog (with RA,DEC) from all the files
 ## ending in `_drz.fits' in `/mnt/data/COSMOS/':
-$ astcrop --mode=wcs catalog.txt /mnt/data/COSMOS/*_drz.fits
+$ astcrop --mode=wcs --catalog=cat.txt /mnt/data/COSMOS/*_drz.fits
 
-## Crop-out the outer 10 border pixels of input image:
+## Crop the outer 10 border pixels of the input image:
 $ astcrop --section=10:*-10,10:*-10 --hdu=2 image.fits
 
 ## Crop region around RA and Dec of (189.16704, 62.218203):
-$ astcrop --ra=189.16704 --dec=62.218203 goodsnorth.fits
+$ astcrop --mode=wcs --center=189.16704,62.218203 goodsnorth.fits
 
 ## Crop region around pixel coordinate (568.342, 2091.719):
-$ astcrop --xc=568.342 --yc=2091.719 --iwidth=201 image.fits
+$ astcrop --mode=img --center=568.342,2091.719 --width=201 image.fits
 @end example
 
 @noindent
 Crop has one mandatory argument which is the input image name(s), shown
 above with @file{ASTRdata ...}. You can use shell expansions, for example
 @command{*} for this if you have lots of images in WCS mode. If the crop
-box centers are in a catalog, you also have to provide the catalog name as
-an argument. Alternatively, you have to provide the crop box parameters
-with command-line options. See @ref{Crop output} for how the output file
-name(s) can be specified. For the full list of general options to all
-Gnuastro programs (including Crop), please see @ref{Common options}.
+box centers are in a catalog, you can use the @option{--catalog} option. In
+other cases, you have to provide the single cropped output parameters must
+be given with command-line options. See @ref{Crop output} for how the
+output file name(s) can be specified. For the full list of general options
+to all Gnuastro programs (including Crop), please see @ref{Common options}.
 
 Floating point numbers can be used to specify the crop region (except the
 @option{--section} option, see @ref{Crop section syntax}). In such cases,
@@ -7972,11 +7976,10 @@ the floating point values will be used to find the 
desired integer pixel
 indices based on the FITS standard. Hence, Crop ultimately doesn't do any
 sub-pixel cropping (in other words, it doesn't change pixel values). If you
 need such crops, you can use @ref{Warp} to first warp the image to the a
-new pixel grid where your initial floating points can be seen as integers,
-then crop from that with Crop. For example, let's say you want a crop from
-pixels 12.982 to 80.982 along the first dimension. You should first
-translate the image by @mymath{-0.482} (note that the edge of a pixel is at
-integer multiples of @mymath{0.5}). So you should run Warp with
+new pixel grid, then crop from that. For example, let's assume you want a
+crop from pixels 12.982 to 80.982 along the first dimension. You should
+first translate the image by @mymath{-0.482} (note that the edge of a pixel
+is at integer multiples of @mymath{0.5}). So you should run Warp with
 @option{--translate=-0.482,0} and then crop the warped image with
 @option{--section=13:81}.
 
@@ -7989,12 +7992,12 @@ for more.
 
 @cindex Asynchronous thread allocation
 When in catalog mode, Crop will run in parallel unless you set
-@option{--numthreads=1}, see @ref{Multi-threaded operations}. Note that when
-multiple threads are being used, in verbose mode, the outputs will not be
-in order. This is because the threads are asynchronous and thus not started
-in order. This has no effect on the output, see @ref{Hubble visually checks
-and classifies his catalog} for a tutorial on effectively using this
-feature.
+@option{--numthreads=1}, see @ref{Multi-threaded operations}. Note that
+when multiple outputs are created with threads, the outputs will not be
+created in the same order. This is because the threads are asynchronous and
+thus not started in order. This has no effect on each output, see
+@ref{Hubble visually checks and classifies his catalog} for a tutorial on
+effectively using this feature.
 
 @menu
 * Crop options::                A list of all the options with explanation.
@@ -8018,14 +8021,6 @@ example, when using @option{--section}, the right 
ascension of the bottom
 left and top left corners will not be equal. If you only want regions
 within a given right ascension, use @option{--polygon} in WCS mode.
 
-@cartouche
-@noindent
-@strong{NOTE:} The coordinates are in the FITS format. So the first
-axis is the horizontal axis when viewed in SAO ds9 and the second axis
-is the vertical. Also in the FITS standard, counting begins from 1
-(one) not 0 (zero).
-@end cartouche
-
 @noindent
 Input image parameters:
 @table @option
@@ -8061,42 +8056,33 @@ coordinate system on the input images. See 
@option{--hstartwcs}
 Crop box parameters:
 @table @option
 
-@item -x FLT
-@itemx --xc=FLT
-X axis (first image axis, horizontal when viewed in SAO ds9) position of
-crop box center. Along with @option{--yc}, this only produces one crop in
-each run. The width of the square crop (in pixels) is the value to
-@option{--iwidth}.
-
-@item -y FLT
-@itemx --yc=FLT
-Y axis (second image axis, vertical when viewed in SAO ds9) position of
-crop box center. See @option{--xc}.
-
-@item -r FLT
-@itemx --ra=FLT
-Right Ascension of crop box center. Along with @option{--dec}, this only
-produces one crop in each run. The width of the square crop (in arcseconds)
-is the value to @option{--wwidth}.
-
-@item -d FLT
-@itemx --dec=FLT
-Declination of crop box center. see @option{--ra}.
-
-@item -a INT
-@itemx --iwidth=INT
-Width of the cropped region in units of pixels when the crop is defined by
-its center in Image mode (see @ref{Crop modes}). In order for the chosen
-central pixel to be in the center of the cropped image, this option only
-accepts an odd number (an error will be given if its even). If you want an
-even sided crop, you can run Crop afterwards with
+@item -c FLT[,FLT[,...]]
+@itemx --center=FLT[,FLT[,...]]
+The central position of the crop in the input image. The positions along
+each dimension must be separated by a comma (@key{,}) and fractions are
+also acceptable. The number of values given to this option must be the same
+as the dimensions of the input dataset. The width of the crop should be set
+with @code{--width}. The units of the coordinates are read based on the
+value to the @option{--mode} option, see below.
+
+@item -w FLT[,FLT[,...]]
+@itemx --width=FLT[,FLT[,...]]
+Width of the cropped region about its center. @option{--width} may take
+either a single value (to be used for all dimensions) or multiple values (a
+specific value for each dimension). If in WCS mode, value(s) given to this
+option will be read in the same units as the dataset's WCS information
+along this dimension. The final output will have an odd number of pixels to
+allow easy identification of the pixel which keeps your requested
+coordinate (from @option{--center} or @option{--catalog}).
+
+The @code{--width} option also accepts fractions. For example if you want
+the width of your crop to be 3 by 5 arcseconds along RA and Dec
+respectively, you can call it with: @option{--width=3/3600,5/3600}.
+
+If you want an even sided crop, you can run Crop afterwards with
 @option{--section=":*-1,:*-1"} or @option{--section=2:,2:} (depending on
 which side you don't need), see @ref{Crop section syntax}.
 
-@item -w FLT
-@itemx --wwidth=FLT
-The width of the crop box in WCS mode in units of arc-seconds.
-
 @item -l STR
 @itemx --polygon=STR
 String of crop polygon vertices. Note that currently only convex polygons
@@ -8192,31 +8178,12 @@ Section of the input image which you want to be 
cropped. See @ref{Crop
 section syntax} for a complete explanation on the syntax required for this
 input.
 
-@item -i STR/INT
-@itemx --xcol=STR/INT
-Column selection for the first FITS axis position of the box center in a
-catalog. In SAO ds9, the first FITS axis is the horizontal axis. The value
-can be either the column number (starting from 1), or a match/search in the
-table meta-data, see @ref{Selecting table columns}.
-
-@item -j STR/INT
-@itemx --ycol=STR/INT
-Column selection for the second FITS axis position of the box center. In
-SAO ds9, the second FITS axis is the vertical axis. The value can be either
-the column number (starting from 1), or a match/search in the table
-meta-data, see @ref{Selecting table columns}.
-
-@item -f STR/INT
-@itemx --racol=STR/INT
-Column selection for Right Ascension (RA) in the input catalog. The value
-can be either the column number (starting from 1), or a match/search in the
-table meta-data, see @ref{Selecting table columns}.
-
-@item -g STR/INT
-@itemx --deccol=STR/INT
-Column selection of declination in the input catalog. The value can be
-either the column number (starting from 1), or a match/search in the table
-meta-data, see @ref{Selecting table columns}.
+@item -x STR/INT
+@itemx --coordcol=STR/INT
+The column in a catalog to read as a coordinate. This option must be called
+multiple times depending on the number of dimensions in the input
+dataset. The value can be either the column number (starting from 1), or a
+match/search in the table meta-data, see @ref{Selecting table columns}.
 
 @item -n STR/INT
 @item --namecol=STR/INT
@@ -8305,18 +8272,21 @@ image or WCS. The value must either be @option{img} or 
@option{wcs}, see
 @node Crop output,  , Crop options, Invoking astcrop
 @subsubsection Crop output
 
-The value given to @option{--output} option will depend on how many crops
-were created, see @ref{Crop modes}:
+The string given to @option{--output} option will be interpretted depending
+on how many crops were requested, see @ref{Crop modes}:
 
 @itemize
 @item
 When a catalog is given, the value of the @option{--output} (see
 @ref{Common options}) will be read as the directory to store the output
 cropped images. Hence if it doesn't already exist, Crop will abort with an
-error of a ``No such file or directory'' error. The crop file names will
-consist of two parts: a variable part (the row number of each target
-starting from 1) along with a fixed string which you can set with the
-@option{--suffix} option.
+error of a ``No such file or directory'' error.
+
+The crop file names will consist of two parts: a variable part (the row
+number of each target starting from 1) along with a fixed string which you
+can set with the @option{--suffix} option. Optionally, you may also use the
+@option{--namecol} option to define a column in the input catalog to use as
+the file name instead of numbers.
 
 @item
 When only one crop is desired, the value to @option{--output} will be read
@@ -8338,9 +8308,11 @@ range from that input image in the same syntax as 
@ref{Crop section
 syntax}. So this string can be directly given to the @option{--section}
 option later.
 
-Once done, a log file will be created in the current directory named
-@file{astcrop.log}. This file will have three columns and the same
-number of rows as the number of cropped images.
+Once done, a log file can be created in the current directory with the
+@code{--log} option. This file will have three columns and the same number
+of rows as the number of cropped images. There are also comments on the top
+of the log file explaining basic information about the run and descriptions
+for the columns. A short description of the columns is also given below:
 
 @enumerate
 @item
@@ -8355,13 +8327,6 @@ given a value of 0 (see @ref{Invoking astcrop}), the 
center will not be
 checked and this column will be given a value of @code{-1}.
 @end enumerate
 
-There are also comments on the top of the log file explaining basic
-information about the run and descriptions for the columns. If the log file
-cannot be created (for example you don't have write permission in the
-directory you are running Crop in) or you have specifically asked for no
-log file (with the @option{--nolog} option), then a log file will not be
-created (unless @option{--individual} is called).
-
 
 
 
@@ -8959,7 +8924,7 @@ equivalent to this piece of C:
 @example
 size_t i;
 long a[100];
-float b[100], out;
+float b[100], out[100];
 for(i=0;i<100;++i) out[i]=a[i]+b[i];
 @end example
 
@@ -17866,7 +17831,7 @@ acceptable values to this element are defined in 
@ref{Table input
 output}. Based on C's @code{printf} standards.
 
 @item disp_precision
-Width of printing each element of the dataset to a plain text file, the
+Precision of printing each element of the dataset to a plain text file, the
 acceptable values to this element are defined in @ref{Table input
 output}. Based on C's @code{printf} standards.
 
@@ -17970,7 +17935,7 @@ memory or it is used in multiple datasets, be sure to 
set it to @code{NULL}
 @code{gal_data_free_contents}.
 @end deftypefun
 
-@deftypefun {void *} gal_data_alloc (void @code{*array}, uint8_t @code{type}, 
size_t @code{ndim}, size_t @code{*dsize}, struct wcsprm @code{*wcs}, int 
@code{clear}, size_t @code{minmapsize}, char @code{*name}, char @code{*unit}, 
char @code{*comment})
+@deftypefun {gal_data_t *} gal_data_alloc (void @code{*array}, uint8_t 
@code{type}, size_t @code{ndim}, size_t @code{*dsize}, struct wcsprm 
@code{*wcs}, int @code{clear}, size_t @code{minmapsize}, char @code{*name}, 
char @code{*unit}, char @code{*comment})
 
 Dynamically allocate a @code{gal_data_t} and initialize it will all the
 given values. See the description of @code{gal_data_initialize} and
@@ -19655,7 +19620,7 @@ which is better considering floating point errors:
 @dispmath{{\sin^2(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)}
 @end deftypefun
 
-@deftypefun {double *} gal_wcs_pixel_scale_deg (struct wcsprm @code{*wcs})
+@deftypefun {double *} gal_wcs_pixel_scale (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.
@@ -20811,19 +20776,25 @@ size in long integer type. @code{width[0]}, and 
@code{width[1]} are the
 number of pixels along the first and second FITS axis.
 @end deftypefun
 
-@deftypefun void gal_box_border_from_center (double @code{xc}, double 
@code{yc}, long @code{*width}, long @code{*fpixel}, long @code{*lpixel})
-Given the center (@code{xc} and @code{yc}) and width (two element array) of a
-box, return the coordinates of the first @code{fpixel} and last @code{lpixel}
-pixels (both are two element arrays.
+@deftypefun void gal_box_border_from_center (double @code{center}, size_t 
@code{ndim}, long @code{*width}, long @code{*fpixel}, long @code{*lpixel})
+Given the center coordinates in @code{center} and the @code{width} (along
+each dimension) of a box, return the coordinates of the first
+(@code{fpixel}) and last (@code{lpixel}) pixels. All arrays must have
+@code{ndim} elements (one for each dimension).
 @end deftypefun
 
-@deftypefun int gal_box_overlap (long @code{*naxes}, long @code{*fpixel_i}, 
long @code{*lpixel_i}, long @code{*fpixel_o}, long @code{*lpixel_o})
-We have an image of size @code{naxes} and want to get the overlap of the
-image with a box. This function will return 1 if there is an overlap and 0
-if there isn't. The input and output box are specified by their first and
-last pixels. When there is an overlap, the coordinates of the first and
-last pixels of the overlap will be put in @code{fpixel_o} and
-@code{lpixel_o}.
+@deftypefun int gal_box_overlap (long @code{*naxes}, long @code{*fpixel_i}, 
long @code{*lpixel_i}, long @code{*fpixel_o}, long @code{*lpixel_o}, size_t 
@code{ndim})
+An @code{ndim}-dimensional dataset of size @code{naxes} (along each
+dimension, in FITS order) and a box with first and last (inclusive)
+coordinate of @code{fpixel_i} and @code{lpixel_i} is given. This box
+doesn't necessarily have to lie within the dataset, it can be outside of
+it, or only patially overlap. This function will change the values of
+@code{fpixel_i} and @code{lpixel_i} to exactly cover the overlap in the
+input dataset's coordinates.
+
+This function will return 1 if there is an overlap and 0 if there
+isn't. When there is an overlap, the coordinates of the first and last
+pixels of the overlap will be put in @code{fpixel_o} and @code{lpixel_o}.
 @end deftypefun
 
 
diff --git a/lib/blank.c b/lib/blank.c
index 487f77f..44dee32 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -195,7 +195,7 @@ gal_blank_present(gal_data_t *input, int updateflag)
         error(EXIT_FAILURE, 0, "%s: tile mode is currently not supported for "
               "strings", __func__);
       strf = (str=input->array) + input->size;
-      do if(!strcmp(*str++,GAL_BLANK_STRING)) return 1; while(str<strf);
+      do if(!strcmp(*str,GAL_BLANK_STRING)) return 1; while(++str<strf);
       break;
 
     /* Complex types */
diff --git a/lib/box.c b/lib/box.c
index 365d027..1b085fd 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -24,6 +24,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <math.h>
 #include <stdio.h>
+#include <errno.h>
+#include <error.h>
 #include <stdlib.h>
 
 #include <gnuastro/box.h>
@@ -86,29 +88,40 @@ gal_box_ellipse_in_box(double a, double b, double 
theta_rad, long *width)
 
 
 /* We have the central pixel and box size of the crop box, find the
-   starting and ending pixels: */
+   starting (first) and ending (last) pixels: */
 void
-gal_box_border_from_center(double xc, double yc, long *width,
+gal_box_border_from_center(double *center, size_t ndim, long *width,
                            long *fpixel, long *lpixel)
 {
-  long lxc, lyc;
+  size_t i;
+  long tmp;
   double intpart;
 
-  /* Round the double values in a the long values: */
-  lxc=xc;
-  if (fabs(modf(xc, &intpart))>=0.5)
-    ++lxc;
-  lyc=yc;
-  if (fabs(modf(yc, &intpart))>=0.5)
-    ++lyc;
-
-  /* Set the initial values for the actual image: */
-  fpixel[0]=lxc-width[0]/2;      fpixel[1]=lyc-width[1]/2;
-  lpixel[0]=lxc+width[0]/2;      lpixel[1]=lyc+width[1]/2;
-  /*
-  printf("\n\nCenter is on: %ld, %ld\n", lxc, lyc);
-  printf("Starting and ending pixels: (%ld, %ld) -- (%ld, %ld)\n\n\n",
-         fpixel[0], fpixel[1], lpixel[0], lpixel[1]);
+  /* Go over the dimensions. */
+  for(i=0;i<ndim;++i)
+    {
+      /* When the decimal fraction of the center (in floating point) is
+         larger than 0.5, then it is must be incremented. */
+      tmp = center[i] + ( fabs(modf(center[i], &intpart))>=0.5 ? 1 : 0 );
+
+      /* Set the first and last pixels in this dimension. */
+      fpixel[i]=tmp-width[i]/2;
+      lpixel[i]=tmp+width[i]/2;
+    }
+
+  /* For a check:
+  if(ndim==2)
+    {
+      printf("center: %g, %g\n", center[0], center[1]);
+      printf("fpixel: %ld, %ld\n", fpixel[0], fpixel[1]);
+      printf("lpixel: %ld, %ld\n", lpixel[0], lpixel[1]);
+    }
+  else if(ndim==3)
+    {
+      printf("center: %g, %g, %g\n", center[0], center[1], center[2]);
+      printf("fpixel: %ld, %ld, %ld\n", fpixel[0], fpixel[1], fpixel[2]);
+      printf("lpixel: %ld, %ld, %ld\n", lpixel[0], lpixel[1], lpixel[2]);
+    }
   */
 }
 
@@ -167,7 +180,7 @@ gal_box_border_from_center(double xc, double yc, long 
*width,
            -----------------------------
 
 
-   So in short the arrays we are dealing with here are:
+   So, in short the arrays we are dealing with here are:
 
    fpixel_i:    Coordinates of the first pixel in input image.
    lpixel_i:    Coordinates of the last pixel in input image.
@@ -188,80 +201,85 @@ gal_box_border_from_center(double xc, double yc, long 
*width,
 */
 int
 gal_box_overlap(long *naxes, long *fpixel_i, long *lpixel_i,
-                long *fpixel_o, long *lpixel_o)
+                long *fpixel_o, long *lpixel_o, size_t ndim)
 {
-  long width[2];
+  size_t i;
+  long width;
 
   /* In case you want to see how things are going:
-  printf("\n\nImage size: [%ld,  %ld]\n", naxes[0], naxes[1]);
-  printf("fpixel_i -- lpixel_i: (%ld, %ld) -- (%ld, %ld)\n\n", fpixel_i[0],
-         fpixel_i[1], lpixel_i[0], lpixel_i[1]);
+  printf("\n\nImage size: [");
+  for(i=0;i<ndim;++i) {printf("%ld, ", naxes[i]);} printf("\b\b]\n");
+  printf("fpixel_i -- lpixel_i: (");
+  for(i=0;i<ndim;++i) {printf("%ld, ", fpixel_i[i]);} printf("\b\b) -- (");
+  for(i=0;i<ndim;++i) {printf("%ld, ", lpixel_i[i]);} printf("\b\b)\n");
   */
 
-  width[0]=lpixel_i[0]-fpixel_i[0]+1;
-  width[1]=lpixel_i[1]-fpixel_i[1]+1;
-
-  /* Set the initial values for the cropped array: */
-  fpixel_o[0]=1;          fpixel_o[1]=1;
-  lpixel_o[0]=width[0];   lpixel_o[1]=width[1];
-
-
-  /* Check the four corners to see if they should be adjusted: To
-     understand the first, look at this, suppose | separates pixels
-     and the * shows where the image actually begins.
-
-     |-2|-1| 0* 1| 2| 3| 4|        (survey image)
-     *1 | 2| 3| 4| 5| 6| 7|        (crop image)
-
-     So when fpixel_i is negative, e.g., fpixel_i=-2, then the index
-     of the pixel in the cropped image we want to begin with, that
-     corresponds to 1 in the survey image is: fpixel_o= 4 = 2 +
-     -1*fpixel_i.*/
-  if(fpixel_i[0]<1)
-    {
-      fpixel_o[0]=-1*fpixel_i[0]+2;
-      fpixel_i[0]=1;
-    }
-  if(fpixel_i[1]<1)
+  /* Do the calculations for each dimension: */
+  for(i=0;i<ndim;++i)
     {
-      fpixel_o[1]=-1*fpixel_i[1]+2;
-      fpixel_i[1]=1;
+      /* Set the width and initialize the overlap values. */
+      fpixel_o[i] = 1;
+      lpixel_o[i] = width = lpixel_i[i] - fpixel_i[i] + 1;
+
+
+      /* Check the four corners to see if they should be adjusted: To
+         understand the first, look at this, suppose | separates pixels and
+         the * shows where the image actually begins.
+
+             |-2|-1| 0* 1| 2| 3| 4|        ( input image )
+             *1 | 2| 3| 4| 5| 6| 7|        ( crop image  )
+
+         So when fpixel_i is negative, e.g., fpixel_i=-2, then the index of
+         the pixel in the cropped image we want to begin with, that
+         corresponds to 1 in the survey image is:
+
+             fpixel_o = 4 = 2 + -1*fpixel_i
+      */
+      if(fpixel_i[i]<1)
+        {
+          fpixel_o[i] = -1*fpixel_i[i]+2;
+          fpixel_i[i] = 1;
+        }
+
+
+      /*The same principle applies to the end of an image. Take "s" is the
+        maximum size along a specific axis in the survey image and "c" is
+        the size along the same axis on the cropped image.  Assume the the
+        cropped region's last pixel in that axis will be 2 pixels larger
+        than s:
+
+             |s-1|   s* s+1| s+2|           (survey image)
+             |c-3| c-2| c-1|   c*           (crop image)
+
+        So you see that if the outer pixel is n pixels away then in the
+        cropped image we should only fill upto c-n.*/
+      if(lpixel_i[i]>naxes[i])
+        {
+          lpixel_o[i] = width - (lpixel_i[i]-naxes[i]);
+          lpixel_i[i] = naxes[i];
+        }
     }
 
+  /* In case you want to see the results.
+  printf("\nAfter correction:\n");
+  printf("Input image: (");
+  for(i=0;i<ndim;++i) {printf("%ld, ", fpixel_i[i]);} printf("\b\b) -- (");
+  for(i=0;i<ndim;++i) {printf("%ld, ", lpixel_i[i]);} printf("\b\b)\n");
+  printf("output image: (");
+  for(i=0;i<ndim;++i) {printf("%ld, ", fpixel_o[i]);} printf("\b\b) -- (");
+  for(i=0;i<ndim;++i) {printf("%ld, ", lpixel_o[i]);} printf("\b\b)\n");
+  */
 
-  /*The same principle applies to the end of an image. Take "s"
-    is the maximum size along a specific axis in the survey image
-    and "c" is the size along the same axis on the cropped image.
-    Assume the the cropped region's last pixel in that axis will
-    be 2 pixels larger than s:
-
-    |s-1|   s* s+1| s+2|           (survey image)
-    |c-3| c-2| c-1|   c*           (crop image)
 
-    So you see that if the outer pixel is n pixels away then in
-    the cropped image we should only fill upto c-n.*/
-  if (lpixel_i[0]>naxes[0])
-    {
-      lpixel_o[0]=width[0]-(lpixel_i[0]-naxes[0]);
-      lpixel_i[0]=naxes[0];
-    }
-  if (lpixel_i[1]>naxes[1])
-    {
-      lpixel_o[1]=width[1]-(lpixel_i[1]-naxes[1]);
-      lpixel_i[1]=naxes[1];
-    }
+  /* If the first image pixel in every dimension is larger than the input
+     image's width or the last image pixel is before the input image, then
+     there is no overlap and we must return 0. */
+  for(i=0;i<ndim;++i)
+    if( fpixel_i[i]>naxes[i] || lpixel_i[i]<1 )
+      return 0;
 
-  /* In case you wish to see the results.
-  printf("\nAfter correction:\n");
-  printf("Input image: (%ld, %ld) -- (%ld, %ld)\n", fpixel_i[0],
-         fpixel_i[1], lpixel_i[0], lpixel_i[1]);
-  printf("output image:(%ld, %ld) -- (%ld, %ld)\n\n\n", fpixel_o[0],
-         fpixel_o[1], lpixel_o[0], lpixel_o[1]);
-  */
 
-  if(fpixel_i[0]>naxes[0] || fpixel_i[1]>naxes[1]
-     || lpixel_i[0]<1 || lpixel_i[1]<1)
-    return 0;
-  else
-    return 1;
+  /* The first and last image pixels are within the image, so there is
+     overlap. */
+  return 1;
 }
diff --git a/lib/gnuastro-internal/options.h b/lib/gnuastro-internal/options.h
index caddc38..84b3978 100644
--- a/lib/gnuastro-internal/options.h
+++ b/lib/gnuastro-internal/options.h
@@ -49,6 +49,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #define GAL_OPTIONS_MAX_VALUE_LEN 10
 
 
+/* Statically allocated space for printing option values. */
+#define GAL_OPTIONS_STATIC_MEM_FOR_VALUES 2000
+
 
 
 /* Standard groups of options. From the Argp manual (in GNU C Library): "In
diff --git a/lib/gnuastro/box.h b/lib/gnuastro/box.h
index 5ee522a..a1aa764 100644
--- a/lib/gnuastro/box.h
+++ b/lib/gnuastro/box.h
@@ -55,12 +55,12 @@ void
 gal_box_ellipse_in_box(double a, double b, double theta_rad, long *width);
 
 void
-gal_box_border_from_center(double xc, double yc, long *width,
+gal_box_border_from_center(double *center, size_t ndim, long *width,
                            long *fpixel, long *lpixel);
 
 int
 gal_box_overlap(long *naxes, long *fpixel_i, long *lpixel_i,
-                long *fpixel_o, long *lpixel_o);
+                long *fpixel_o, long *lpixel_o, size_t ndim);
 
 
 
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 226d534..314891e 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -65,8 +65,8 @@ gal_dimension_num_neighbors(size_t ndim);
 /************************************************************************/
 /********************          Coordinates         **********************/
 /************************************************************************/
-#define GAL_DIMENSION_FLT_TO_INT(FLT) ( (FLT)-(int)(FLT) > 0.5f  \
-                                        ? (int)(FLT)+1 : (int)(FLT) )
+#define GAL_DIMENSION_FLT_TO_INT(FLT) ( (FLT)-(long)(FLT) > 0.5f  \
+                                        ? (long)(FLT)+1 : (long)(FLT) )
 
 void
 gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 8a2a70e..78c6dec 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -86,7 +86,7 @@ double
 gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
 
 double *
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs);
+gal_wcs_pixel_scale(struct wcsprm *wcs);
 
 double
 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
diff --git a/lib/options.c b/lib/options.c
index bc57997..d3978a4 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -514,7 +514,6 @@ gal_options_read_tableformat(struct argp_option *option, 
char *arg,
    element of the option. The last element of the array is `-1' to allow
    finding the number of elements within it later (similar to a string
    which terminates with a '\0' element). */
-#define PARSE_SIZES_STATICSTR_LEN 2000
 void *
 gal_options_parse_sizes_reverse(struct argp_option *option, char *arg,
                                 char *filename, size_t lineno, void *junk)
@@ -523,7 +522,7 @@ gal_options_parse_sizes_reverse(struct argp_option *option, 
char *arg,
   double *v;
   gal_data_t *values;
   size_t nc, num, *array;
-  char *str, sstr[PARSE_SIZES_STATICSTR_LEN];
+  char *str, sstr[GAL_OPTIONS_STATIC_MEM_FOR_VALUES];
 
   /* We want to print the stored values. */
   if(lineno==-1)
@@ -537,12 +536,12 @@ gal_options_parse_sizes_reverse(struct argp_option 
*option, char *arg,
       nc=0;
       for(i=num-1;i>=0;--i)
         {
-          if( nc > PARSE_SIZES_STATICSTR_LEN-100 )
+          if( nc > GAL_OPTIONS_STATIC_MEM_FOR_VALUES-100 )
             error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so we "
                   "can address the problem. The number of necessary "
                   "characters in the statically allocated string has become "
                   "too close to %d", __func__, PACKAGE_BUGREPORT,
-                  PARSE_SIZES_STATICSTR_LEN);
+                  GAL_OPTIONS_STATIC_MEM_FOR_VALUES);
           nc += sprintf(sstr+nc, "%zu,", array[i]);
         }
       sstr[nc-1]='\0';
@@ -567,14 +566,16 @@ gal_options_parse_sizes_reverse(struct argp_option 
*option, char *arg,
       for(i=0;i<values->size;++i)
         {
           if(v[i]<0)
-            error(EXIT_FAILURE, 0, "the given value in `%s' (%g) is not 0 "
-                  "or positive. The values to the `--%s' option must be "
-                  "positive", arg, v[i], option->name);
+            error_at_line(EXIT_FAILURE, 0, filename, lineno, "the given "
+                          "value in `%s' (%g) is not 0 or positive. The "
+                          "values to the `--%s' option must be positive",
+                          arg, v[i], option->name);
 
           if(ceil(v[i]) != v[i])
-            error(EXIT_FAILURE, 0, "the given value in `%s' (%g) is not an "
-                  "integer. The values to the `--%s' option must be "
-                  "integers", arg, v[i], option->name);
+            error_at_line(EXIT_FAILURE, 0, filename, lineno, "the given "
+                          "value in `%s' (%g) is not an integer. The "
+                          "values to the `--%s' option must be integers",
+                          arg, v[i], option->name);
         }
 
       /* Write the values into an allocated size_t array and finish it with
diff --git a/lib/wcs.c b/lib/wcs.c
index eb566e9..bb321bf 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -339,7 +339,7 @@ gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
   if(wcs->pc)
     {
       /* Get the pixel scale. */
-      ps=gal_wcs_pixel_scale_deg(wcs);
+      ps=gal_wcs_pixel_scale(wcs);
 
       /* The PC matrix and the CDELT elements might both contain scale
          elements (during processing the scalings might be added only to PC
@@ -407,9 +407,9 @@ gal_wcs_angular_distance_deg(double r1, double d1, double 
r2, double d2)
 
 
 
-/* Return the pixel scale of the image in both dimentions in degrees. */
+/* Return the pixel scale of the image (in the units of the input dataset). */
 double *
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs)
+gal_wcs_pixel_scale(struct wcsprm *wcs)
 {
   gsl_vector S;
   gsl_matrix A, V;
@@ -466,7 +466,7 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
           "The input WCS has %d dimensions", __func__, wcs->naxis);
 
   /* Get the pixel scales along each axis in degrees, then multiply. */
-  pixscale=gal_wcs_pixel_scale_deg(wcs);
+  pixscale=gal_wcs_pixel_scale(wcs);
 
   /* Clean up and return the result. */
   out = pixscale[0] * pixscale[1] * 3600.0f * 3600.0f;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 73190ca..11ab1ee 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -87,18 +87,18 @@ if COND_COSMICCAL
   cosmiccal/simpletest.sh: prepconf.sh.log
 endif
 if COND_CROP
-  MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/xcyc.sh                
\
-  crop/xcycnoblank.sh crop/section.sh crop/radec.sh crop/imgpolygon.sh \
-  crop/imgoutpolygon.sh crop/wcspolygon.sh
+  MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/imgcenter.sh    \
+  crop/imgcenternoblank.sh crop/section.sh crop/wcscenter.sh            \
+  crop/imgpolygon.sh crop/imgoutpolygon.sh crop/wcspolygon.sh
 
   crop/imgcat.sh: mkprof/mosaic1.sh.log
   crop/wcscat.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log     \
                   mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
-  crop/xcyc.sh: mkprof/mosaic1.sh.log
-  crop/xcycnoblank.sh: mkprof/mosaic1.sh.log
+  crop/imgcenter.sh: mkprof/mosaic1.sh.log
+  crop/imgcenternoblank.sh: mkprof/mosaic1.sh.log
   crop/section.sh: mkprof/mosaic1.sh.log
-  crop/radec.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log      \
-                 mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
+  crop/wcscenter.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log      \
+                     mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
   crop/imgpolygon.sh: mkprof/mosaic1.sh.log
   crop/imgoutpolygon.sh: mkprof/mosaic1.sh.log
   crop/wcspolygon.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
diff --git a/tests/crop/imgcat.sh b/tests/crop/imgcat.sh
index b782fcd..824db0d 100755
--- a/tests/crop/imgcat.sh
+++ b/tests/crop/imgcat.sh
@@ -54,5 +54,6 @@ if [ ! -f $img      ]; then echo "$img does not exist.";   
exit 77; fi
 # enable multithreaded access to files, the tests pass. It is the
 # users choice to enable this feature.
 cat=$topsrc/tests/$prog/cat.txt
-$execname $img --catalog=$cat --mode=img --suffix=_imgcat.fits --numthreads=1 \
-          --zeroisnotblank --xcol=X_CENTER --ycol=Y_CENTER --namecol=NAME
+$execname $img --catalog=$cat --suffix=_imgcat.fits --numthreads=1      \
+          --zeroisnotblank --coordcol=X_CENTER --coordcol=Y_CENTER      \
+          --namecol=NAME --mode=img --width=201
diff --git a/tests/crop/xcyc.sh b/tests/crop/imgcenter.sh
similarity index 89%
rename from tests/crop/xcyc.sh
rename to tests/crop/imgcenter.sh
index f472ee0..f56d149 100755
--- a/tests/crop/xcyc.sh
+++ b/tests/crop/imgcenter.sh
@@ -1,4 +1,4 @@
-# Crop a box based on --xc and --yc.
+# Make a single crop with center defined in Image mode.
 #
 # See the Tests subsection of the manual for a complete explanation
 # (in the Installing gnuastro section).
@@ -53,4 +53,5 @@ if [ ! -f $img      ]; then echo "$img does not exist.";   
exit 77; fi
 # The number of threads is one so if CFITSIO does is not configured to
 # enable multithreaded access to files, the tests pass. It is the
 # users choice to enable this feature.
-$execname $img --xc=251 --yc=251 --output=crop_xcyc.fits --numthreads=1
+$execname $img --center=251,251 --output=crop_imgcenter.fits    \
+          --numthreads=1 --mode=img --width=201
diff --git a/tests/crop/xcycnoblank.sh b/tests/crop/imgcenternoblank.sh
similarity index 86%
rename from tests/crop/xcycnoblank.sh
rename to tests/crop/imgcenternoblank.sh
index 394ab4a..e23826b 100755
--- a/tests/crop/xcycnoblank.sh
+++ b/tests/crop/imgcenternoblank.sh
@@ -1,4 +1,5 @@
-# Crop a box based on --xc and --yc, with no blank pixels.
+# Crop a region defined by its center in image coordinates, with no blank
+# pixels.
 #
 # See the Tests subsection of the manual for a complete explanation
 # (in the Installing gnuastro section).
@@ -53,5 +54,5 @@ if [ ! -f $img      ]; then echo "$img does not exist.";   
exit 77; fi
 # The number of threads is one so if CFITSIO does is not configured to
 # enable multithreaded access to files, the tests pass. It is the
 # users choice to enable this feature.
-$execname $img --xc=500 --yc=500 --noblank --output=crop_xcycnb.fits \
-          --numthreads=1
+$execname $img --center=500,500 --noblank --numthreads=1              \
+          --output=crop_imgcenternoblank.fits --mode=img --width=201
diff --git a/tests/crop/wcscat.sh b/tests/crop/wcscat.sh
index 011194b..48a1923 100755
--- a/tests/crop/wcscat.sh
+++ b/tests/crop/wcscat.sh
@@ -55,5 +55,6 @@ done
 # enable multithreaded access to files, the tests pass. It is the
 # users choice to enable this feature.
 cat=$topsrc/tests/$prog/cat.txt
-$execname $img --catalog=$cat --mode=wcs --suffix=_wcscat.fits        \
-          --zeroisnotblank --racol=4 --deccol=DEC_CENTER --numthreads=1
+$execname $img --catalog=$cat --suffix=_wcscat.fits            \
+          --zeroisnotblank --coordcol=4 --coordcol=DEC_CENTER  \
+          --numthreads=1 --mode=wcs --width=3/3600
diff --git a/tests/crop/radec.sh b/tests/crop/wcscenter.sh
similarity index 88%
rename from tests/crop/radec.sh
rename to tests/crop/wcscenter.sh
index 001f961..5d62411 100755
--- a/tests/crop/radec.sh
+++ b/tests/crop/wcscenter.sh
@@ -1,4 +1,4 @@
-# Crop a box based on the --ra and --dec options.
+# Crop a box based on center in WCS coordinates.
 #
 # See the Tests subsection of the manual for a complete explanation
 # (in the Installing gnuastro section).
@@ -49,5 +49,5 @@ done
 
 # Actual test script
 # ==================
-$execname $img --ra=0.99917157 --dec=1.0008283 --wwidth=0.3 \
-          --output=crop_radec.fits
+$execname $img --center=0.99917157,1.0008283 --output=crop_wcscenter.fits \
+          --mode=wcs --width=0.3/3600



reply via email to

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