[gnuastro-commits] master dadb3f3 107/113: Imported recent work in master, minor conflicts fixed

[Mohammad Akhlaghi]

branch: master
commit dadb3f3febee052c4fc1c17a6d33d96e50a3e7b8
Merge: 5d5fd95 4c5bdff
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Imported recent work in master, minor conflicts fixed
    
    Some minor conflicts came up due to the recent checks in dimensionality
    that werer fixed with this commit.
---
 NEWS                         |  17 ++++++-
 THANKS                       |   1 +
 bin/arithmetic/arithmetic.c  |   2 +
 bin/arithmetic/operands.c    |  10 +++-
 bin/arithmetic/ui.c          |  13 ++++++
 bin/arithmetic/ui.h          |   3 ++
 bin/convertt/ui.c            |   2 +
 bin/convolve/ui.c            |  14 ++++--
 bin/fits/args.h              |  13 ++++++
 bin/fits/fits.c              |  54 ++++++++++++++++++++--
 bin/fits/main.h              |   1 +
 bin/fits/ui.h                |   1 +
 bin/mkcatalog/main.h         |   1 +
 bin/mkcatalog/ui.c           | 106 ++++++++++++++++++++++++++++++++++++++-----
 bin/mknoise/ui.c             |   2 +
 bin/mkprof/ui.c              |  18 ++++----
 bin/noisechisel/ui.c         |   9 +++-
 bin/segment/ui.c             |  17 +++++++
 bin/statistics/ui.c          |  11 ++++-
 bin/warp/ui.c                |   3 ++
 doc/announce-acknowledge.txt |   1 +
 doc/gnuastro.texi            |  92 +++++++++++++++++++++++++++----------
 lib/blank.c                  |  19 ++++++++
 lib/dimension.c              |  45 ++++++++++++++++++
 lib/fits.c                   |  21 +++++++++
 lib/gnuastro/blank.h         |   3 ++
 lib/gnuastro/dimension.h     |   8 ++++
 lib/gnuastro/fits.h          |   3 ++
 28 files changed, 431 insertions(+), 59 deletions(-)

diff --git a/NEWS b/NEWS
index d14635a..937f7c1 100644
--- a/NEWS
+++ b/NEWS
@@ -23,6 +23,10 @@ See the end of the file for license conditions.
    --usedredshift: Print the used redshift as a "Specific calculation" (in
      line with other single-valued calculations).
 
+  Fits:
+   --primaryimghdu: Copy/cut the given image HDU to the zero-th/first
+     HDU of the output file that doesn't yet exist.
+
   Statistics:
    --sigclip-number, --sigclip-median, --sigclip-mean, --sigclip-std: Do
      sigma-clipping and only print the desired value as a single-value
@@ -52,11 +56,14 @@ See the end of the file for license conditions.
      for using them in combination with their names.
    - list.h: Functions to return the last element in linked lists. For
      example `gal_list_sizet_last' or `gal_list_data_last'.
+   - gal_arithmetic_operator_string: Return operator string from code.
+   - gal_arithmetic_set_operator: Return operator code from string.
+   - gal_blank_initialize_array: Initialize an array with blank values.
+   - gal_dimension_remove_extra: Remove extra (length 1) dimensions.
    - gal_list_data_to_array_ptr: Make an array of pointers from the list.
+   - gal_fits_img_info_dim: Only return the size information of a dataset.
    - GAL_BLANK_INT: Blank value for `int' (can be 16-bit or 32-bit).
    - GAL_BLANK_UINT: Blank value for unsigned `int' (can be 16-bit or 32-bit).
-   - gal_arithmetic_operator_string: Return operator string from code.
-   - gal_arithmetic_set_operator: Return operator code from string.
 
 
 ** Removed features
@@ -77,6 +84,12 @@ See the end of the file for license conditions.
      `sigclip-mean' and `sigclip-std') return 32-bit floating point
      datasets
 
+  MakeCatalog:
+   - When a clumps catalog is requested, MakeCatalog will automatically
+     deduce the total number of clumps (at a small cost in
+     performance). Until now, it was mandatory for the clumps label dataset
+     to contain the total number of clumps in the `NUMLABS' keyword.
+
   Library:
    - gal_statistics_outlier_flat_cfp: Improved implementation with new API.
 
diff --git a/THANKS b/THANKS
index bfe6ab9..45074bc 100644
--- a/THANKS
+++ b/THANKS
@@ -23,6 +23,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Hamed Altafi                         hamed.altafi2@gmail.com
     Roland Bacon                         roland.bacon@univ-lyon1.fr
     Roberto Baena Gallé                  rbaena@iac.es
+    Zahra Bagheri                        bagheri.zahra87@gmail.com
     Karl Berry                           karl@gnu.org
     Leindert Boogaard                    boogaard@strw.leidenuniv.nl
     Nicolas Bouché                       nicolas.bouche@irap.omp.eu
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 7c85a94..f0586b4 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1211,6 +1211,8 @@ reversepolish(struct arithmeticparams *p)
           /* Read the data, note that the WCS has already been set. */
           p->operands->data=gal_array_read_one_ch(filename, hdu, NULL,
                                                   p->cp.minmapsize);
+          data=p->operands->data;
+          data->ndim=gal_dimension_remove_extra(data->ndim, data->dsize, NULL);
           if(!p->cp.quiet) printf(" - %s (hdu %s) is read.\n", filename, hdu);
         }
       else
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 577fd93..9b7cceb 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -267,6 +267,7 @@ operands_copy_named(struct arithmeticparams *p, char *name)
 void
 operands_add(struct arithmeticparams *p, char *filename, gal_data_t *data)
 {
+  size_t ndim, *dsize;
   struct operand *newnode;
 
   /* Some operators might not actually return any dataset (data=NULL), in
@@ -305,13 +306,17 @@ operands_add(struct arithmeticparams *p, char *filename, 
gal_data_t *data)
               else
                 newnode->hdu=gal_list_str_pop(&p->hdus);
 
-              /* If no WCS is set yet, use the WCS of this image. */
+              /* If no WCS is set yet, use the WCS of this image and remove
+                 possibly extra dimensions if necessary. */
               if(p->refdata.wcs==NULL)
                 {
+                  dsize=gal_fits_img_info_dim(filename, newnode->hdu, &ndim);
                   p->refdata.wcs=gal_wcs_read(filename, newnode->hdu, 0, 0,
                                               &p->refdata.nwcs);
+                  ndim=gal_dimension_remove_extra(ndim, dsize, p->refdata.wcs);
                   if(p->refdata.wcs && !p->cp.quiet)
                     printf(" - WCS: %s (hdu %s).\n", filename, newnode->hdu);
+                  free(dsize);
                 }
             }
           else newnode->hdu=NULL;
@@ -349,8 +354,9 @@ operands_pop(struct arithmeticparams *p, char *operator)
       hdu=operands->hdu;
       filename=operands->filename;
 
-      /* Read the dataset. */
+      /* Read the dataset and remove possibly extra dimensions. */
       data=gal_array_read_one_ch(filename, hdu, NULL, p->cp.minmapsize);
+      data->ndim=gal_dimension_remove_extra(data->ndim, data->dsize, NULL);
 
       /* Arithmetic changes the contents of a dataset, so the existing name
          (in the FITS `EXTNAME' keyword) should not be passed on beyond
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 0ef9a24..f1e39e5 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -357,10 +357,16 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
 static void
 ui_preparations(struct arithmeticparams *p)
 {
+  size_t ndim, *dsize;
+
   /* In case a file is specified to read the WCS from (and ignore input
      datasets), read the WCS prior to starting parsing of the arguments. */
   if(p->wcsfile)
     {
+      /* Read the number of dimensions and the size of each. */
+      dsize=gal_fits_img_info_dim(p->wcsfile, p->wcshdu, &ndim);
+
+      /* Read the WCS. */
       p->refdata.wcs=gal_wcs_read(p->wcsfile, p->wcshdu, 0, 0,
                                   &p->refdata.nwcs);
       if(p->refdata.wcs)
@@ -371,6 +377,13 @@ ui_preparations(struct arithmeticparams *p)
       else
         fprintf(stderr, "WARNING: %s (hdu %s) didn't contain a "
                 "(readable by WCSLIB) WCS.", p->wcsfile, p->wcshdu);
+
+      /* Correct the WCS dimensions if necessary. Note that we don't need
+         the `ndim' or `dsize' any more. */
+      ndim=gal_dimension_remove_extra(ndim, dsize, p->refdata.wcs);
+
+      /* Clean up. */
+      free(dsize);
     }
 }
 
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/ui.h
index ebd20f8..37bf97d 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/ui.h
@@ -51,6 +51,9 @@ enum option_keys_enum
 void
 ui_read_check_inputs_setup(int argc, char *argv[], struct arithmeticparams *p);
 
+size_t *
+ui_read_ndim_dsize(char *filename, char *hdu, size_t *ndim);
+
 void
 freeandreport(struct arithmeticparams *p, struct timeval *t1);
 
diff --git a/bin/convertt/ui.c b/bin/convertt/ui.c
index a3c3330..145ae5d 100644
--- a/bin/convertt/ui.c
+++ b/bin/convertt/ui.c
@@ -523,6 +523,8 @@ ui_make_channels_ll(struct converttparams *p)
           /* Read in the array and its WCS information. */
           data=gal_fits_img_read(name->v, hdu, p->cp.minmapsize);
           data->wcs=gal_wcs_read(name->v, hdu, 0, 0, &data->nwcs);
+          data->ndim=gal_dimension_remove_extra(data->ndim, data->dsize,
+                                                data->wcs);
           gal_list_data_add(&p->chll, data);
 
           /* A FITS file only has one channel. */
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 7d131ba..c446c0d 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -438,6 +438,9 @@ ui_read_input(struct convolveparams *p)
                                                p->cp.minmapsize);
         p->input->wcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0,
                                    &p->input->nwcs);
+        p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
+                                                  p->input->dsize,
+                                                  p->input->wcs);
       }
 
   /* The input isn't an image (wasn't read yet), so we'll read it as a
@@ -459,9 +462,14 @@ ui_read_kernel(struct convolveparams *p)
   if( p->kernelname
       && p->input->ndim>1
       && gal_array_name_recognized(p->kernelname)  )
-    p->kernel = gal_array_read_one_ch_to_type(p->kernelname, p->khdu,
-                                              NULL, INPUT_USE_TYPE,
-                                              p->cp.minmapsize);
+    {
+      p->kernel = gal_array_read_one_ch_to_type(p->kernelname, p->khdu,
+                                                NULL, INPUT_USE_TYPE,
+                                                p->cp.minmapsize);
+      p->kernel->ndim=gal_dimension_remove_extra(p->kernel->ndim,
+                                                 p->kernel->dsize,
+                                                 p->kernel->wcs);
+    }
   else
     p->kernel=ui_read_column(p, 1);
 
diff --git a/bin/fits/args.h b/bin/fits/args.h
index 5f3d899..b4c6584 100644
--- a/bin/fits/args.h
+++ b/bin/fits/args.h
@@ -89,6 +89,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "primaryimghdu",
+      UI_KEY_PRIMARYIMGHDU,
+      0,
+      0,
+      "Copy/cut image HDUs to primary/zero-th HDU.",
+      UI_GROUP_EXTENSION,
+      &p->primaryimghdu,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 432bbbd..11b3d3b 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -25,6 +25,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <errno.h>
 #include <error.h>
 #include <string.h>
+#include <unistd.h>
 
 #include <gnuastro/list.h>
 #include <gnuastro/fits.h>
@@ -54,7 +55,7 @@ fits_has_error(struct fitsparams *p, int actioncode, char 
*string, int status)
     case FITS_ACTION_UPDATE:        action="updated";      break;
     case FITS_ACTION_WRITE:         action="written";      break;
     case FITS_ACTION_COPY:          action="copied";       break;
-    case FITS_ACTION_REMOVE:        action="renoved";      break;
+    case FITS_ACTION_REMOVE:        action="removed";      break;
 
     default:
       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at `%s' so we "
@@ -299,16 +300,53 @@ fits_hdu_remove(struct fitsparams *p, int *r)
 
 
 
+/* This is similar to the library's `gal_fits_open_to_write', except that
+   it won't create an empty first extension.*/
+fitsfile *
+fits_open_to_write_no_blank(char *filename)
+{
+  int status=0;
+  fitsfile *fptr;
+
+  /* When the file exists just open it. Otherwise, create the file. But we
+     want to leave the first extension as a blank extension and put the
+     image in the next extension to be consistent between tables and
+     images. */
+  if(access(filename,F_OK) == -1 )
+    {
+      /* Create the file. */
+      if( fits_create_file(&fptr, filename, &status) )
+        gal_fits_io_error(status, NULL);
+    }
+
+  /* Open the file, ready for later steps. */
+  if( fits_open_file(&fptr, filename, READWRITE, &status) )
+    gal_fits_io_error(status, NULL);
+
+  /* Return the pointer. */
+  return fptr;
+}
+
+
+
+
+
 static void
 fits_hdu_copy(struct fitsparams *p, int cut1_copy0, int *r)
 {
   char *hdu;
-  fitsfile *in, *out;
   int status=0, hdutype;
+  fitsfile *in, *out=NULL;
   gal_list_str_t *list = cut1_copy0 ? p->cut : p->copy;
 
-  /* Open the output file. */
-  out=gal_fits_open_to_write(p->cp.output);
+  /* Open the output file.
+
+############################################# USE this after checking the
+nature of the first extension to be copied. If its a table, use the
+library's function.  #############################################
+
+*/
+
 
   /* Copy all the given extensions. */
   while(list)
@@ -320,6 +358,14 @@ fits_hdu_copy(struct fitsparams *p, int cut1_copy0, int *r)
       in=gal_fits_hdu_open(p->filename, hdu,
                            cut1_copy0 ? READWRITE : READONLY);
 
+      /* If the output isn't opened yet, open it.  */
+      if(out==NULL)
+        out = ( ( gal_fits_hdu_format(p->filename, hdu)==IMAGE_HDU
+                  && p->primaryimghdu )
+                ? fits_open_to_write_no_blank(p->cp.output)
+                : gal_fits_open_to_write(p->cp.output) );
+
+
       /* Copy to the extension. */
       if( fits_copy_hdu(in, out, 0, &status) )
         *r=fits_has_error(p, FITS_ACTION_COPY, hdu, status);
diff --git a/bin/fits/main.h b/bin/fits/main.h
index 460a067..6b50be2 100644
--- a/bin/fits/main.h
+++ b/bin/fits/main.h
@@ -61,6 +61,7 @@ struct fitsparams
   gal_list_str_t    *copy;     /* Copy extensions to output.            */
   gal_list_str_t     *cut;     /* Copy ext. to output and remove.       */
   uint8_t         numhdus;     /* Print number of HDUs in FITS file.    */
+  uint8_t   primaryimghdu;     /* Copy/cut HDU into primary HDU.        */
   uint8_t    printallkeys;     /* Print all the header keywords.        */
   uint8_t            date;     /* Set DATE to current time.             */
   gal_list_str_t    *asis;     /* Strings to write asis.                */
diff --git a/bin/fits/ui.h b/bin/fits/ui.h
index 229cc61..4249fb4 100644
--- a/bin/fits/ui.h
+++ b/bin/fits/ui.h
@@ -72,6 +72,7 @@ enum option_keys_enum
   UI_KEY_TITLE        = 1000,
   UI_KEY_OUTHDU,
   UI_KEY_COPYKEYS,
+  UI_KEY_PRIMARYIMGHDU,
 };
 
 
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 4296f72..00eba86 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -195,6 +195,7 @@ struct mkcatalogparams
   int32_t       checkuplim[2];  /* Object & clump ID to check dist.     */
 
   /* Internal. */
+  char           *relabclumps;  /* Name of new file for clump labels.   */
   time_t              rawtime;  /* Starting time of the program.        */
   gal_data_t          *values;  /* Input.                               */
   gal_data_t         *objects;  /* Object labels.                       */
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 6927a54..4aa1e97 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -532,17 +532,95 @@ ui_wcs_info(struct mkcatalogparams *p)
 
 
 
+static size_t
+ui_num_clumps(struct mkcatalogparams *p)
+{
+  char *basename;
+  size_t i, counter, numclumps=0;
+  gal_list_i32_t *tmp, **labsinobj;
+  int32_t *o=p->objects->array, *of=o+p->objects->size, *c=p->clumps->array;
+
+  /* Allocate array of lists to keep the unique labels within each object. */
+  errno=0;
+  labsinobj=calloc(p->numobjects+1, sizeof *labsinobj);
+  if(labsinobj==NULL)
+    error(EXIT_FAILURE, 0, "%s: couldn't allocate %zu bytes for "
+          "`labsinobj'", __func__, p->numobjects * sizeof *labsinobj);
+
+  /* Go over each pixel and find the unique clump labels within each
+     object. */
+  do
+    {
+      /* Do the steps if we are on a clump. */
+      if(*o>0 && *c>0)
+        {
+          /* See if the label has already been found. */
+          for(tmp=labsinobj[*o];tmp!=NULL;tmp=tmp->next) if(tmp->v==*c) break;
+
+          /* When it wasn't found, `tmp==NULL'. */
+          if(tmp==NULL)
+            {
+              ++numclumps;
+              gal_list_i32_add(&labsinobj[*o], *c);
+            }
+        }
+
+      /* Increment the clumps pointer.*/
+      ++c;
+    }
+  while(++o<of);
+
+  /* Re-write the clump values so their numbering is contiguous, since this
+     is assumed during the later steps. */
+  o=p->objects->array;
+  c=p->clumps->array;
+  do
+    {
+      /* Do the steps if we are on a clump. */
+      if(*o>0 && *c>0)
+        {
+          counter=0;
+          for(tmp=labsinobj[*o];tmp!=NULL;tmp=tmp->next)
+            { counter++; if(tmp->v==*c) {*c=counter; break;} }
+        }
+
+      /* Increment the clumps pointer.*/
+      ++c;
+    }
+  while(++o<of);
+
+  /* Write the created file into a file for the user to inspect. */
+  basename = p->cp.output ? p->cp.output : p->objectsfile;
+  p->relabclumps=gal_checkset_automatic_output(&p->cp, basename,
+                                               "-clumps-relab.fits");
+  gal_fits_img_write(p->clumps, p->relabclumps, NULL, PROGRAM_STRING);
+
+
+  /* Clean up. */
+  for(i=0;i<p->numobjects;++i)
+    if(labsinobj[i]) gal_list_i32_free(labsinobj[i]);
+  free(labsinobj);
+
+  /* Return the number of clumps. */
+  return numclumps;
+}
+
+
+
+
+
 /* The only mandatory input is the objects image, so first read that and
    make sure its type is correct. */
 static void
 ui_read_labels(struct mkcatalogparams *p)
 {
-  char *msg;
   gal_data_t *tmp, *keys=gal_data_array_calloc(2);
 
   /* Read it into memory. */
   p->objects = gal_array_read_one_ch(p->objectsfile, p->cp.hdu, NULL,
                                      p->cp.minmapsize);
+  p->objects->ndim=gal_dimension_remove_extra(p->objects->ndim,
+                                              p->objects->dsize, NULL);
 
 
   /* Make sure it has an integer type. */
@@ -611,6 +689,8 @@ ui_read_labels(struct mkcatalogparams *p)
       /* Read the clumps image. */
       p->clumps = gal_array_read_one_ch(p->usedclumpsfile, p->clumpshdu,
                                         NULL, p->cp.minmapsize);
+      p->clumps->ndim=gal_dimension_remove_extra(p->clumps->ndim,
+                                                 p->clumps->dsize, NULL);
 
       /* Check its size. */
       if( gal_dimension_is_different(p->objects, p->clumps) )
@@ -630,17 +710,7 @@ ui_read_labels(struct mkcatalogparams *p)
       keys[0].array=&p->clumpsn;            keys[1].array=&p->numclumps;
       gal_fits_key_read(p->usedclumpsfile, p->clumpshdu, keys, 0, 0);
       if(keys[0].status) p->clumpsn=NAN;
-      if(keys[1].status)
-        {
-          if( asprintf(&msg, "%s (hdu: %s): couldn't find/read `NUMLABS' in "
-                       "the header keywords, see CFITSIO error above. The "
-                       "clumps image must have the total number of clumps "
-                       "(irrespective of how many objects there are in the "
-                       "image) in this header keyword", p->usedclumpsfile,
-                       p->clumpshdu)<0 )
-            error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
-          gal_fits_io_error(keys[1].status, msg);
-        }
+      if(keys[1].status) p->numclumps=ui_num_clumps(p);
 
       /* If there were no clumps, then free the clumps array and set it to
          NULL, so for the rest of the processing, MakeCatalog things that
@@ -930,6 +1000,8 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
       p->values=gal_array_read_one_ch_to_type(p->usedvaluesfile, p->valueshdu,
                                               NULL, GAL_TYPE_FLOAT32,
                                               p->cp.minmapsize);
+      p->values->ndim=gal_dimension_remove_extra(p->values->ndim,
+                                                 p->values->dsize, NULL);
 
       /* Make sure it has the correct size. */
       if( gal_dimension_is_different(p->objects, p->values) )
@@ -977,6 +1049,8 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
           p->sky=gal_array_read_one_ch_to_type(p->usedskyfile, p->skyhdu,
                                                NULL, GAL_TYPE_FLOAT32,
                                                p->cp.minmapsize);
+          p->sky->ndim=gal_dimension_remove_extra(p->sky->ndim,
+                                                  p->sky->dsize, NULL);
 
           /* Check its size and prepare tile structure. */
           ui_preparation_check_size_read_tiles(p, p->sky, p->usedskyfile,
@@ -1006,6 +1080,8 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
       p->std=gal_array_read_one_ch_to_type(p->usedstdfile, p->stdhdu,
                                            NULL, GAL_TYPE_FLOAT32,
                                            p->cp.minmapsize);
+      p->std->ndim=gal_dimension_remove_extra(p->std->ndim,
+                                              p->std->dsize, NULL);
 
       /* Check its size and prepare tile structure. */
       ui_preparation_check_size_read_tiles(p, p->std, p->usedstdfile,
@@ -1039,6 +1115,9 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
           /* Read the mask image. */
           p->upmask = gal_array_read_one_ch(p->upmaskfile, p->upmaskhdu,
                                             NULL, p->cp.minmapsize);
+          p->upmask->ndim=gal_dimension_remove_extra(p->upmask->ndim,
+                                                     p->upmask->dsize,
+                                                     NULL);
 
           /* Check its size. */
           if( gal_dimension_is_different(p->objects, p->upmask) )
@@ -1631,6 +1710,9 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
mkcatalogparams *p)
       if(p->clumps)
         printf("  - Clumps:  %s (hdu: %s)\n", p->usedclumpsfile,
                p->clumpshdu);
+      if(p->relabclumps)
+        printf("  - RELABELED CLUMPS (no NUMLABS in original): %s\n",
+               p->relabclumps);
       if(p->values)
         printf("  - Values:  %s (hdu: %s)\n", p->usedvaluesfile,
                p->valueshdu);
diff --git a/bin/mknoise/ui.c b/bin/mknoise/ui.c
index 3347590..5c4baad 100644
--- a/bin/mknoise/ui.c
+++ b/bin/mknoise/ui.c
@@ -286,6 +286,8 @@ ui_preparations(struct mknoiseparams *p)
   p->input=gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu, NULL,
                                      GAL_TYPE_FLOAT64, p->cp.minmapsize);
   p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs);
+  p->input->ndim=gal_dimension_remove_extra(p->input->ndim, p->input->dsize,
+                                            p->input->wcs);
 
 
   /* If we are dealing with an input table, make sure the format of the
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index fccdc87..5b009da 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -1332,6 +1332,7 @@ ui_prepare_canvas(struct mkprofparams *p)
 {
   float *f, *ff;
   long width[3]={1,1,1};
+  size_t tndim, *tdsize;
   int status=0, setshift=0;
   double truncr, semiaxes[3], euler_deg[3];
   size_t i, nshift=0, *dsize=NULL, ndim_counter;
@@ -1344,7 +1345,10 @@ ui_prepare_canvas(struct mkprofparams *p)
          no merged image is desired, we just need the WCS information of
          the background image and the number of its dimensions. So
          `ndim==0' and what `dsize' points to is irrelevant. */
+      tdsize=gal_fits_img_info_dim(p->backname, p->backhdu, &tndim);
       p->wcs=gal_wcs_read(p->backname, p->backhdu, 0, 0, &p->nwcs);
+      tndim=gal_dimension_remove_extra(tndim, tdsize, p->wcs);
+      free(tdsize);
       if(p->nomerged==0)
         {
           /* If p->dsize was given as an option, free it. */
@@ -1629,8 +1633,7 @@ ui_make_log(struct mkprofparams *p)
 static void
 ui_read_ndim(struct mkprofparams *p)
 {
-  gal_data_t *keysll;
-  size_t i, ndim_counter;
+  size_t i, *dsize, ndim_counter;
 
   if(p->kernel)
     {
@@ -1658,12 +1661,9 @@ ui_read_ndim(struct mkprofparams *p)
           if(p->nomerged)
             {
               /* Get the number of the background image's dimensions. */
-              keysll=gal_data_array_calloc(1);
-              keysll->name="NAXIS";   keysll->type=GAL_TYPE_SIZE_T;
-              gal_fits_key_read(p->backname, p->backhdu, keysll, 0, 0);
-              p->ndim = *(size_t *)(keysll->array);
-              keysll->name=NULL;
-              gal_data_array_free(keysll, 1, 1);
+              dsize=gal_fits_img_info_dim(p->backname, p->backhdu, &p->ndim);
+              p->ndim=gal_dimension_remove_extra(p->ndim, dsize, NULL);
+              free(dsize);
             }
           else
             {
@@ -1671,6 +1671,8 @@ ui_read_ndim(struct mkprofparams *p)
               p->out=gal_array_read_one_ch_to_type(p->backname, p->backhdu,
                                                    NULL, GAL_TYPE_FLOAT32,
                                                    p->cp.minmapsize);
+              p->out->ndim=gal_dimension_remove_extra(p->out->ndim,
+                                                      p->out->dsize, NULL);
               p->ndim=p->out->ndim;
             }
 
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index ac259c4..1dad65e 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -577,12 +577,19 @@ ui_preparations_read_input(struct noisechiselparams *p)
   float *f;
   size_t ndim;
 
-  /* Read the input as a single precision floating point dataset. */
+  /* Read the input as a single precision floating point dataset, also load
+     the WCS and finally remove any possibly existing extra dimensions
+     (with a length of 1). */
   p->input = gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu,
                                            NULL, GAL_TYPE_FLOAT32,
                                            p->cp.minmapsize);
   p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu, 0, 0,
                                &p->input->nwcs);
+  p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
+                                            p->input->dsize,
+                                            p->input->wcs);
+
+  /* When the input doesn't have a name, use `INPUT'. */
   if(p->input->name==NULL)
     gal_checkset_allocate_copy("INPUT", &p->input->name);
 
diff --git a/bin/segment/ui.c b/bin/segment/ui.c
index e8bb046..7a01946 100644
--- a/bin/segment/ui.c
+++ b/bin/segment/ui.c
@@ -413,6 +413,11 @@ ui_prepare_inputs(struct segmentparams *p)
                                            p->cp.minmapsize);
   p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu, 0, 0,
                                &p->input->nwcs);
+  p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
+                                            p->input->dsize,
+                                            p->input->wcs);
+
+  /* Set the name. */
   if(p->input->name) free(p->input->name);
   gal_checkset_allocate_copy("INPUT", &p->input->name);
 
@@ -435,6 +440,9 @@ ui_prepare_inputs(struct segmentparams *p)
       p->conv = gal_array_read_one_ch_to_type(p->convolvedname, p->chdu,
                                               NULL, GAL_TYPE_FLOAT32,
                                               p->cp.minmapsize);
+      p->conv->ndim=gal_dimension_remove_extra(p->conv->ndim,
+                                               p->conv->dsize,
+                                               p->conv->wcs);
       p->conv->wcs=gal_wcs_copy(p->input->wcs);
 
       /* Make sure it is the same size as the input. */
@@ -453,6 +461,8 @@ ui_prepare_inputs(struct segmentparams *p)
       /* Read the dataset into memory. */
       p->olabel = gal_array_read_one_ch(p->useddetectionname, p->dhdu,
                                         NULL, p->cp.minmapsize);
+      p->olabel->ndim=gal_dimension_remove_extra(p->olabel->ndim,
+                                                 p->olabel->dsize, NULL);
       if( gal_dimension_is_different(p->input, p->olabel) )
         error(EXIT_FAILURE, 0, "`%s' (hdu: %s) and `%s' (hdu: %s) have a"
               "different dimension/size", p->useddetectionname, p->dhdu,
@@ -536,6 +546,9 @@ ui_prepare_kernel(struct segmentparams *p)
           /* Read the kernel into memory. */
           p->kernel=gal_fits_img_read_kernel(p->kernelname, p->khdu,
                                              p->cp.minmapsize);
+          p->kernel->ndim=gal_dimension_remove_extra(p->kernel->ndim,
+                                                     p->kernel->dsize,
+                                                     NULL);
 
           /* Make sure it has the same dimensions as the input. */
           if( p->kernel->ndim != p->input->ndim )
@@ -743,6 +756,8 @@ ui_read_std_and_sky(struct segmentparams *p)
       p->std=gal_array_read_one_ch_to_type(p->usedstdname, p->stdhdu,
                                            NULL, GAL_TYPE_FLOAT32,
                                            p->cp.minmapsize);
+      p->std->ndim=gal_dimension_remove_extra(p->std->ndim,
+                                              p->std->dsize, NULL);
 
       /* Make sure it has the correct size. */
       ui_check_size(p->input, p->std, tl->tottiles, p->inputname, p->cp.hdu,
@@ -801,6 +816,8 @@ ui_read_std_and_sky(struct segmentparams *p)
           sky=gal_array_read_one_ch_to_type(p->skyname, p->skyhdu,
                                             NULL, GAL_TYPE_FLOAT32,
                                             p->cp.minmapsize);
+          sky->ndim=gal_dimension_remove_extra(sky->ndim, sky->dsize,
+                                               NULL);
 
           /* Check its size. */
           ui_check_size(p->input, sky, tl->tottiles, p->inputname, p->cp.hdu,
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index a85881d..2378b31 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -848,6 +848,9 @@ ui_preparations(struct statisticsparams *p)
                                      cp->minmapsize);
       p->input->wcs=gal_wcs_read(p->inputname, cp->hdu, 0, 0,
                                  &p->input->nwcs);
+      p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
+                                                p->input->dsize,
+                                                p->input->wcs);
     }
   else
     {
@@ -857,8 +860,12 @@ ui_preparations(struct statisticsparams *p)
 
   /* Read the convolution kernel if necessary. */
   if(p->sky && p->kernelname)
-    p->kernel=gal_fits_img_read_kernel(p->kernelname, p->khdu,
-                                       cp->minmapsize);
+    {
+      p->kernel=gal_fits_img_read_kernel(p->kernelname, p->khdu,
+                                         cp->minmapsize);
+      p->kernel->ndim=gal_dimension_remove_extra(p->kernel->ndim,
+                                                 p->kernel->dsize, NULL);
+    }
 
   /* Tile and channel sanity checks and preparations. */
   if(p->ontile || p->sky)
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 0e11d19..165f926 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -348,6 +348,9 @@ ui_check_options_and_arguments(struct warpparams *p)
                                              p->cp.minmapsize);
       p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, p->hstartwcs,
                                  p->hendwcs, &p->input->nwcs);
+      p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
+                                                p->input->dsize,
+                                                p->input->wcs);
       if(p->input->wcs)
         {
           p->pixelscale=gal_wcs_pixel_scale(p->input->wcs);
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 1c554bd..431c8a0 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,6 +1,7 @@
 Alphabetically ordered list to acknowledge in the next release.
 
 Hamed Altafi
+Zahra Bagheri
 Leindert Boogaard
 Bruno Haible
 Raul Infante-Sainz
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 93b6e45..8586978 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10258,15 +10258,24 @@ output file, see explanations above.
 
 @item -R STR
 @itemx --remove=STR
-Remove the specified HDU from the input file. From CFITSIO: ``In the case
-of deleting the primary array (the first HDU in the file) then [it] will be
-replaced by a null primary array containing the minimum set of required
-keywords and no data.''. So in practice, any existing data (array) and
-meta-data in the first extension will be removed, but the number of
+Remove the specified HDU from the input file.
+
+The first (zero-th) HDU cannot be removed with this option. Consider using
+@option{--copy} or @option{--cut} in combination with
+@option{primaryimghdu} to not have an empty zero-th HDU. From CFITSIO: ``In
+the case of deleting the primary array (the first HDU in the file) then
+[it] will be replaced by a null primary array containing the minimum set of
+required keywords and no data.''. So in practice, any existing data (array)
+and meta-data in the first extension will be removed, but the number of
 extensions in the file won't change. This is because of the unique position
 the first FITS extension has in the FITS standard (for example it cannot be
 used to store tables).
 
+@item --primaryimghdu
+Copy or cut an image HDU to the zero-th HDU/extension a file that doesn't
+yet exist. This option is thus irrelevant if the output file already exists
+or the copied/cut extension is a FITS table.
+
 @end table
 
 
@@ -19784,7 +19793,7 @@ There are separate options to explicitly request a file 
name and
 HDU/extension for each of the required input datasets as fully described
 below (with the @option{--*file} format). When each dataset is in a
 separate file, these options are necessary. However, one great advantage of
-the FITS file format, that is heavily used in astronomy, is that it allows
+the FITS file format (that is heavily used in astronomy) is that it allows
 the storage of multiple datasets in one file. So in most situations (for
 example if you are using the outputs of @ref{NoiseChisel} or
 @ref{Segment}), all the necessary input datasets can be in one file.
@@ -19797,6 +19806,24 @@ necessary and the only @option{--*file} option called 
is
 default/given HDUs) in the file given to @option{--valuesfile} (before
 looking into the the main argument file).
 
+When the clumps image (necessary with the @option{--clumpscat} option) is
+used, MakeCatalog looks into the (possibly existing) @code{NUMLABS} keyword
+for the total number of clumps in the image (irrespective of how many
+objects there are). If its not present, it will count them and possibly
+re-label the clumps so the clump labels always start with 1 and finish with
+the total number of clumps in each object. The re-labeled clumps image will
+be stored with the @file{-clumps-relab.fits} suffix. This can slighly
+slow-down the run.
+
+Note that @code{NUMLABS} is automatically written by Segment in its
+outputs, so if you are feeding Segment's clump labels, you can benefit from
+the improved speed. Otherwise, if you are creating the clumps label dataset
+manually, it may be good to include the @code{NUMLABS} keyword in its
+header and also be sure that there is no gap in the clump labels. For
+example if an object has three clumps, they are labeled as 1, 2, 3. If they
+are labeled as 1, 3, 4, or any other combination of three positive integers
+that aren't an increment of the previous, you might get un-known behavior.
+
 It may happen that your labeled objects image was created with a program
 that only outputs floating point files. However, you know it only has
 integer valued pixels that are stored in a floating point container. In
@@ -20247,8 +20274,8 @@ the clump multiplied by the area of the clump) is 
removed, see
 catalog will be smaller than the total clump brightness in the
 @option{--clumpbrightness} column of the objects catalog.
 
-If no usable pixels (blank or below the threshold) are present over the
-clump or object, the stored value will be NaN (note that zero is
+If no usable pixels are present over the clump or object (for example they
+are all blank), the returned value will be NaN (note that zero is
 meaningful).
 
 @item --brightnesserr
@@ -20258,9 +20285,8 @@ clumps.
 @item --clumpbrightness
 [Objects] The total brightness of the clumps within an object. This is
 simply the sum of the pixels associated with clumps in the object. If no
-usable pixels (blank or below the threshold) are present over the clump or
-object, the stored value will be NaN, because zero (note that zero is
-meaningful).
+usable pixels are present over the clump or object (for example they are
+all blank), the stored value will be NaN (note that zero is meaningful).
 
 @item --brightnessnoriver
 [Clumps] The Sky (not river) subtracted clump brightness. By definition,
@@ -20270,20 +20296,9 @@ neighbors. In cases where you will be calculating the 
flux brightness
 difference later (one example below) the contamination will be (mostly)
 removed at that stage, which is why this column was added.
 
-One example might be this: you want to know the change in the clump flux as
-a function of threshold (see @option{--threshold}). So you will make two
-catalogs (each having this column but with different thresholds) and then
-subtract the lower threshold catalog (higher brightness) from the higher
-threshold catalog (lower brightness). The effect is most visible when the
-rivers have a high average signal-to-noise ratio. The removed contribution
-from the pixels below the threshold will be less than the river
-pixels. Therefore the river-subtracted brightness (@option{--brightness})
-for the thresholded catalog for such clumps will be larger than the
-brightness with no threshold!
-
-If no usable pixels (blank or below the possibly given threshold) are
-present over the clump or object, the stored value will be NaN (note that
-zero is meaningful).
+If no usable pixels are present over the clump or object (for example they
+are all blank), the stored value will be NaN (note that zero is
+meaningful).
 
 @item --mean
 The mean sky subtracted value of pixels within the object or clump. For
@@ -25176,6 +25191,12 @@ that corresponds to its type. If @code{input} is a 
tile over a larger
 dataset, only the region that the tile covers will be set to blank.
 @end deftypefun
 
+@deftypefun void gal_blank_initialize_array (void @code{*array}, size_t 
@code{size}, uint8_t @code{type})
+Initialize all the elements in the @code{array} to the blank value that
+corresponds to its type (identified with @code{type}), assuming the array
+has @code{size} elements.
+@end deftypefun
+
 @deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int 
@code{width})
 Write the blank value for the given data type @code{type} into a string and
 return it. The space for the string is dynamically allocated so it must be
@@ -25851,6 +25872,17 @@ have the maximum value along the given dimension and 
if it is zero, the
 minimum.
 @end deftypefun
 
+@deftypefun size_t gal_dimension_remove_extra (size_t @code{ndim}, size_t 
@code{*dsize}, struct wcsprm @code{*wcs})
+Remove extra dimensions (those that only have a length of 1) from the basic
+size information of a dataset. @code{ndim} is the number of dimensions and
+@code{dsize} is an array with @code{ndim} elements containing the size
+along each dimension in the C dimension order. When @code{wcs!=NULL}, the
+respective dimension will also be removed from the WCS.
+
+This function will return the new number of dimensions and the @code{dsize}
+elements will contain the length along each new dimension.
+@end deftypefun
+
 @deffn {Function-like macro} GAL_DIMENSION_NEIGHBOR_OP (@code{index}, 
@code{ndim}, @code{dsize}, @code{connectivity}, @code{dinc}, @code{operation})
 Parse the neighbors of the element located at @code{index} and do the
 requested operation on them. This is defined as a macro to allow easy
@@ -27544,6 +27576,12 @@ then if the image has a name and units, the respective 
string will be put
 in these pointers.
 @end deftypefun
 
+@deftypefun {size_t *} gal_fits_img_info_dim (char @code{*filename}, char 
@code{*hdu}, size_t @code{*ndim})
+Put the number of dimensions in the @code{hdu} extension of @file{filename}
+in the space that @code{ndim} points to and return the size of the dataset
+along each dimension as an allocated array with @code{*ndim} elements.
+@end deftypefun
+
 @deftypefun {gal_data_t *} gal_fits_img_read (char @code{*filename}, char 
@code{*hdu}, size_t @code{minmapsize})
 Read the contents of the @code{hdu} extension/HDU of @code{filename} into a
 Gnuastro generic data container (see @ref{Generic data container}) and
@@ -28183,6 +28221,10 @@ points to. Please see @code{gal_wcs_read_fitsptr} for 
more.
 Return a fully allocated (independent) copy of @code{wcs}.
 @end deftypefun
 
+@deftypefun void gal_wcs_remove_dimension (struct wcsprm @code{*wcs}, size_t 
@code{fitsdim})
+Remove the given FITS dimension from the given @code{wcs} structure.
+@end deftypefun
+
 @deftypefun void gal_wcs_on_tile (gal_data_t @code{*tile})
 Create a WCSLIB @code{wcsprm} structure for @code{tile} using WCS
 parameters of the tile's allocated block dataset, see @ref{Tessellation
diff --git a/lib/blank.c b/lib/blank.c
index 4445d5d..a195018 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -114,6 +114,25 @@ gal_blank_initialize(gal_data_t *input)
 
 
 
+/* Initialize an array to the given type's blank values.*/
+void
+gal_blank_initialize_array(void *array, size_t size, uint8_t type)
+{
+  size_t i, w=gal_type_sizeof(type);
+  void *b=gal_blank_alloc_write(type);
+
+  /* Set all the elements to blank. */
+  for(i=0;i<size;++i)
+    memcpy(gal_pointer_increment(array, i, type), b, w);
+
+  /* Clean up. */
+  free(b);
+}
+
+
+
+
+
 /* Print the blank value as a string. For the integer types, we'll use the
    PRIxNN keywords of `inttypes.h' (which is imported into Gnuastro from
    Gnulib, so we don't necessarily rely on the host system having it). */
diff --git a/lib/dimension.c b/lib/dimension.c
index 310b5cb..61397b8 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -795,3 +795,48 @@ gal_dimension_collapse_minmax(gal_data_t *in, size_t 
c_dim, int max1_min0)
   if(num) gal_data_free(num);
   return minmax;
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/************************************************************************/
+/********************             Other            **********************/
+/************************************************************************/
+size_t
+gal_dimension_remove_extra(size_t ndim, size_t *dsize, struct wcsprm *wcs)
+{
+  size_t i, j;
+
+  for(i=0;i<ndim;++i)
+    if(dsize[i]==1)
+      {
+        /* Correct the WCS. */
+        if(wcs) gal_wcs_remove_dimension(wcs, ndim-i);
+
+        /* Shift all subsequent dimensions to replace this one. */
+        for(j=i;j<ndim-1;++j) dsize[j]=dsize[j+1];
+
+        /* Decrement the `i' and the total number of dimension. */
+        --i;
+        --ndim;
+      }
+
+  /* Return the number of dimensions. */
+  return ndim;
+}
diff --git a/lib/fits.c b/lib/fits.c
index fbbea32..6ccdae7 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1790,6 +1790,27 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t 
*ndim, size_t **dsize,
 
 
 
+/* Get the basic array info to remove extra dimensions if necessary. */
+size_t *
+gal_fits_img_info_dim(char *filename, char *hdu, size_t *ndim)
+{
+  fitsfile *fptr;
+  size_t *dsize=NULL;
+  int status=0, type;
+
+  /* Open the given header, read the basic image information and close it
+     again. */
+  fptr=gal_fits_hdu_open(filename, hdu, READONLY);
+  gal_fits_img_info(fptr, &type, ndim, &dsize, NULL, NULL);
+  if( fits_close_file(fptr, &status) ) gal_fits_io_error(status, NULL);
+
+  return dsize;
+}
+
+
+
+
+
 /* Read a FITS image HDU into a Gnuastro data structure. */
 gal_data_t *
 gal_fits_img_read(char *filename, char *hdu, size_t minmapsize)
diff --git a/lib/gnuastro/blank.h b/lib/gnuastro/blank.h
index 0bb67da..0c52522 100644
--- a/lib/gnuastro/blank.h
+++ b/lib/gnuastro/blank.h
@@ -108,6 +108,9 @@ gal_blank_alloc_write(uint8_t type);
 void
 gal_blank_initialize(gal_data_t *input);
 
+void
+gal_blank_initialize_array(void *array, size_t size, uint8_t type);
+
 char *
 gal_blank_as_string(uint8_t type, int width);
 
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 2bdc9c4..c9bf682 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -116,6 +116,14 @@ gal_dimension_collapse_minmax(gal_data_t *in, size_t 
c_dim, int max1_min0);
 
 
 /************************************************************************/
+/********************             Other            **********************/
+/************************************************************************/
+size_t
+gal_dimension_remove_extra(size_t ndim, size_t *dsize, struct wcsprm *wcs);
+
+
+
+/************************************************************************/
 /********************          Neighbors           **********************/
 /************************************************************************/
 /* Purpose
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 9071294..f25a997 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -235,6 +235,9 @@ void
 gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, size_t **dsize,
                   char **name, char **unit);
 
+size_t *
+gal_fits_img_info_dim(char *filename, char *hdu, size_t *ndim);
+
 gal_data_t *
 gal_fits_img_read(char *filename, char *hdu, size_t minmapsize);