[gnuastro-commits] master d57e0e6 087/113: New spectrum measurement from 3D inputs in MakeCatalog

[Mohammad Akhlaghi]

branch: master
commit d57e0e637598b721c51b55ec48c377a2e748e533
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    New spectrum measurement from 3D inputs in MakeCatalog
    
    MakeCatalog can now estimate the spectra of 3D detections using the
    `--spectrum' option. This was implemented by defining an extra step after
    parsing of objects (and finding the label's 2D projection map). The output
    spectra are written into separate tables (one for each label).
    
    The spectra for each label are writte into a list of datasets which are
    held together using an array of datasets during the processing. Once
    processing is finished, they are written into the output tables. Maybe
    later we can add an option to write the tables using a separate thread
    using CFITSIO's low-level functionality) to speed it up (similar to how we
    write the separate profiles of MakeProfiles into the final image: with each
    thread giving their output to a separate thread that is devoted to the
    writing). But for now, I don't have time to implement that.
---
 bin/mkcatalog/args.h      |  13 ++
 bin/mkcatalog/main.h      |   3 +
 bin/mkcatalog/mkcatalog.c | 134 ++++++++++++-----
 bin/mkcatalog/mkcatalog.h |   1 +
 bin/mkcatalog/parse.c     | 375 ++++++++++++++++++++++++++++++++++++++++++----
 bin/mkcatalog/ui.c        | 130 +++++++++++++++-
 bin/mkcatalog/ui.h        |   1 +
 doc/gnuastro.texi         | 236 +++++++++++++++++++----------
 lib/data.c                |   8 +
 9 files changed, 756 insertions(+), 145 deletions(-)

diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 18fdaea..bd57e58 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -231,6 +231,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "spectrum",
+      UI_KEY_SPECTRUM,
+      0,
+      0,
+      "Object spectrum for cube (3D) datasets.",
+      GAL_OPTIONS_GROUP_OUTPUT,
+      &p->spectrum,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 4fa7a2c..b2f8e28 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -183,6 +183,7 @@ struct mkcatalogparams
   uint8_t         subtractsky;  /* ==1: subtract the Sky from values.   */
   float           sfmagnsigma;  /* Surface brightness multiple of sigma.*/
   float             sfmagarea;  /* Surface brightness area (arcsec^2).  */
+  uint8_t            spectrum;  /* Object spectrum for 3D datasets.     */
 
   char            *upmaskfile;  /* Name of upper limit mask file.       */
   char             *upmaskhdu;  /* HDU of upper limit mask file.        */
@@ -224,6 +225,8 @@ struct mkcatalogparams
   uint8_t      uprangewarning;  /* A warning must be printed.           */
   size_t         *hostobjid_c;  /* To sort the clumps table by Obj.ID.  */
   size_t         *numclumps_c;  /* To sort the clumps table by Obj.ID.  */
+  gal_data_t   *specsliceinfo;  /* Slice information for spectra.       */
+  gal_data_t         *spectra;  /* Array of datasets containing spectra.*/
 
   char        *usedvaluesfile;  /* Ptr to final name used for values.   */
   char        *usedclumpsfile;  /* Ptr to final name used for clumps.   */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index cfa385d..50d6764 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -42,6 +42,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/permutation.h>
 
 #include <gnuastro-internal/timing.h>
+#include <gnuastro-internal/checkset.h>
 
 #include "main.h"
 #include "mkcatalog.h"
@@ -147,9 +148,10 @@ mkcatalog_single_object(void *in_prm)
     {
       /* For easy reading. Note that the object IDs start from one while
          the array positions start from 0. */
-      pp.ci     = NULL;
-      pp.object = tprm->indexs[i] + 1;
-      pp.tile   = &p->tiles[ tprm->indexs[i] ];
+      pp.ci       = NULL;
+      pp.object   = tprm->indexs[i] + 1;
+      pp.tile     = &p->tiles[   tprm->indexs[i] ];
+      pp.spectrum = &p->spectra[ tprm->indexs[i] ];
 
       /* Initialize the parameters for this object/tile. */
       parse_initialize(&pp);
@@ -612,41 +614,84 @@ sort_clumps_by_objid(struct mkcatalogparams *p)
 static void
 mkcatalog_write_outputs(struct mkcatalogparams *p)
 {
-  /*char *str;*/
+  size_t i, scounter;
+  char str[200], *fname;
   gal_list_str_t *comments;
+  int outisfits=gal_fits_name_is_fits(p->objectsout);
 
-
-  /* OBJECT catalog */
-  comments=mkcatalog_outputs_same_start(p, 0, "Detection");
-
-  /* Reverse the comments list (so it is printed in the same order here),
-     write the objects catalog and free the comments. */
-  gal_list_str_reverse(&comments);
-  gal_table_write(p->objectcols, comments, p->cp.tableformat, p->objectsout,
-                  "OBJECTS", 0);
-  gal_list_str_free(comments, 1);
-
-
-  /* CLUMPS catalog */
-  if(p->clumps)
+  /* If a catalog is to be generated. */
+  if(p->objectcols)
     {
-      /* Make the comments. */
-      comments=mkcatalog_outputs_same_start(p, 1, "Clumps");
-
-      /* Write objects catalog
-         ---------------------
+      /* OBJECT catalog */
+      comments=mkcatalog_outputs_same_start(p, 0, "Detection");
 
-         Reverse the comments list (so it is printed in the same order here),
-         write the objects catalog and free the comments. */
+      /* Reverse the comments list (so it is printed in the same order
+         here), write the objects catalog and free the comments. */
       gal_list_str_reverse(&comments);
-      gal_table_write(p->clumpcols, comments, p->cp.tableformat, p->clumpsout,
-                      "CLUMPS", 0);
+      gal_table_write(p->objectcols, comments, p->cp.tableformat,
+                      p->objectsout, "OBJECTS", 0);
       gal_list_str_free(comments, 1);
+
+
+      /* CLUMPS catalog */
+      if(p->clumps)
+        {
+          /* Make the comments. */
+          comments=mkcatalog_outputs_same_start(p, 1, "Clumps");
+
+          /* Write objects catalog
+             ---------------------
+
+             Reverse the comments list (so it is printed in the same order
+             here), write the objects catalog and free the comments. */
+          gal_list_str_reverse(&comments);
+          gal_table_write(p->clumpcols, comments, p->cp.tableformat,
+                          p->clumpsout, "CLUMPS", 0);
+          gal_list_str_free(comments, 1);
+        }
     }
 
+  /* Spectra. */
+  if(p->spectra)
+    {
+      /* Inform the user (Writing many FITS extensions can take
+         long). */
+      if(p->objectcols && outisfits)
+        printf("  - Catalog(s) complete, writing spectra.\n");
+
+      /* Start counting and writing the files. Note that due to some
+         conditions (for example in debugging), a `p->spectra[i]' may not
+         actually contain any data. So we'll also count the number of
+         spectra that are written. */
+      scounter=0;
+      for(i=0;i<p->numobjects;++i)
+        if(p->spectra[i].ndim>0)
+          {
+            /* Increment the written spectra-counter. */
+            ++scounter;
+
+            /* Write the spectra based on the requested format. */
+            if(outisfits)
+              {
+                /* Write the table. */
+                sprintf(str, "SPECTRUM_%zu", i+1);
+                gal_table_write(&p->spectra[i], NULL, GAL_TABLE_FORMAT_BFITS,
+                                p->objectsout, str, 0);
+              }
+            else
+              {
+                sprintf(str, "-spec-%zu.txt", i+1);
+                fname=gal_checkset_automatic_output(&p->cp, p->objectsout,
+                                                    str);
+                gal_table_write(&p->spectra[i], NULL, GAL_TABLE_FORMAT_TXT,
+                                fname, NULL, 0);
+                free(fname);
+              }
+          }
+    }
 
   /* Configuration information. */
-  if(gal_fits_name_is_fits(p->objectsout))
+  if(outisfits)
     {
       gal_fits_key_write_filename("input", p->objectsfile, &p->cp.okeys, 1);
       gal_fits_key_write_config(&p->cp.okeys, "MakeCatalog configuration",
@@ -657,14 +702,35 @@ mkcatalog_write_outputs(struct mkcatalogparams *p)
   /* Inform the user */
   if(!p->cp.quiet)
     {
-      if(p->clumpsout==p->objectsout)
-        printf("  - Output catalog: %s\n", p->objectsout);
-      else
+      if(p->objectcols)
+        {
+          if(p->clumpsout && strcmp(p->clumpsout,p->objectsout))
+            {
+              printf("  - Output objects catalog: %s\n", p->objectsout);
+              if(p->clumps)
+                printf("  - Output clumps catalog: %s\n", p->clumpsout);
+            }
+          else
+            printf("  - Catalog written to %s\n", p->objectsout);
+
+        }
+
+      if(p->spectra)
         {
-          printf("  - Output objects catalog: %s\n", p->objectsout);
-          if(p->clumps)
-            printf("  - Output clumps catalog: %s\n", p->clumpsout);
+          if(outisfits)
+            {
+              if(p->objectcols)
+                printf("  - Spectra in %zu extensions named `SPECTRUM_NN'.\n",
+                       p->numobjects);
+              else
+                printf("  - Output: %s (Spectra in %zu extensions named "
+                       "`SPECTRUM_NN').\n)", p->objectsout, p->numobjects);
+            }
+          else
+            printf("  - Spectra in %zu files with `-spec-NN.txt' suffix.\n",
+                   p->numobjects);
         }
+
     }
 }
 
diff --git a/bin/mkcatalog/mkcatalog.h b/bin/mkcatalog/mkcatalog.h
index be94551..dd0c7bc 100644
--- a/bin/mkcatalog/mkcatalog.h
+++ b/bin/mkcatalog/mkcatalog.h
@@ -41,6 +41,7 @@ struct mkcatalog_passparams
   gsl_rng              *rng;    /* Random number generator.             */
   size_t    clumpstartindex;    /* Clump starting row in final catalog. */
   gal_data_t       *up_vals;    /* Container for upper-limit values.    */
+  gal_data_t      *spectrum;    /* Spectrum of each object.             */
 };
 
 void
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index 52225f5..4c4cc55 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -106,6 +106,322 @@ parse_initialize(struct mkcatalog_passparams *pp)
 
 
 
+static size_t *
+parse_spectrum_pepare(struct mkcatalog_passparams *pp, size_t *start_end_inc,
+                      int32_t **st_o, float **st_v, float **st_std)
+{
+  size_t *tsize;
+  gal_data_t *spectile;
+  struct mkcatalogparams *p=pp->p;
+  size_t coord[3], minmax[6], numslices=p->objects->dsize[0];
+  gal_data_t *area, *sum, *esum, *proj, *eproj, *oarea, *osum, *oesum;
+
+  /* Get the coordinates of the spectral tile's starting element, then make
+     the tile. */
+  gal_dimension_index_to_coord(gal_pointer_num_between(p->objects->array,
+                                                       pp->tile->array,
+                                                       p->objects->type),
+                               p->objects->ndim, p->objects->dsize, coord);
+  minmax[0]=0;                                   /* Changed to first slice.*/
+  minmax[1]=coord[1];
+  minmax[2]=coord[2];
+  minmax[3]=p->objects->dsize[0]-1;              /* Changed to last slice. */
+  minmax[4]=coord[1]+pp->tile->dsize[1]-1;
+  minmax[5]=coord[2]+pp->tile->dsize[2]-1;
+  spectile=gal_tile_series_from_minmax(p->objects, minmax, 1);
+
+  /* Find the starting (and ending) pointers on each of the datasets. */
+  *st_o   = gal_tile_start_end_ind_inclusive(spectile, p->objects,
+                                             start_end_inc);
+  *st_v   = (float *)(p->values->array) + start_end_inc[0];
+  *st_std = ( p->std
+                 ? ( p->std->size==p->objects->size
+                     ? (float *)(p->std->array) + start_end_inc[0]
+                     : NULL )
+                 : NULL );
+
+  /* Allocate the columns. */
+  area  = gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "AREA", "counter",
+                         "Area of object in a slice");
+  sum   = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "SUM", p->values->unit,
+                         "Sum of values with this label.");
+  esum  = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "SUM_ERR", p->values->unit,
+                         "Error in SUM column.");
+  proj  = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "SUM_PROJECTED", p->values->unit,
+                         "Sum of full projected 2D area on a slice.");
+  eproj = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "SUM_PROJECTED_ERR",
+                         p->values->unit, "Error in SUM_PROJECTED column.");
+  oarea = gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "AREA_OTHER", "counter",
+                         "Area covered by other labels in a slice.");
+  osum  = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "SUM_OTHER", p->values->unit,
+                         "Sum of values in other labels on a slice.");
+  oesum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+                         p->cp.minmapsize, "SUM_OTHER_ERR", p->values->unit,
+                         "Error in SUM_OTHER column.");
+
+  /* Fill up the contents of the first element (note that the first
+     `gal_data_t' is actually in an array, so the skeleton is already
+     allocated, we just have to allocate its contents. */
+  gal_data_initialize(pp->spectrum, NULL, p->specsliceinfo->type, 1,
+                      &numslices, NULL, 0, p->cp.minmapsize, NULL, NULL,
+                      NULL);
+  gal_data_copy_to_allocated(p->specsliceinfo, pp->spectrum);
+  pp->spectrum->next=gal_data_copy(p->specsliceinfo->next);
+
+
+  /* Add all the other columns in the final spectrum table. */
+  pp->spectrum->next->next                       = area;
+  area->next                                     = sum;
+  area->next->next                               = esum;
+  area->next->next->next                         = proj;
+  area->next->next->next->next                   = eproj;
+  area->next->next->next->next->next             = oarea;
+  area->next->next->next->next->next->next       = osum;
+  area->next->next->next->next->next->next->next = oesum;
+
+  /* Clean up and return. */
+  tsize=spectile->dsize;
+  spectile->dsize=NULL;
+  gal_data_free(spectile);
+  return tsize;
+}
+
+
+
+
+
+/* Since spectra will be a large table for many objects, it is very
+   important to not consume too much space in for columns that don't need
+   it. This function will check the integer columns and if they are smaller
+   than the maximum values of smaller types, */
+static void
+parse_spectrum_uint32_to_best_type(gal_data_t **input)
+{
+  gal_data_t *tmp=gal_statistics_maximum(*input);
+
+  /* If maximum is smaller than UINT8_MAX, convert it to uint8_t */
+  if( *(uint32_t *)(tmp->array) < UINT8_MAX )
+    *input = gal_data_copy_to_new_type_free(*input, GAL_TYPE_UINT8);
+
+  /* otherwise, if it is smaller than UINT16_MAX, convert it to uint16_t */
+  else if( *(uint32_t *)(tmp->array)      < UINT16_MAX )
+    *input = gal_data_copy_to_new_type_free(*input, GAL_TYPE_UINT16);
+
+  /* Clean up. */
+  gal_data_free(tmp);
+}
+
+
+
+
+
+static void
+parse_spectrum_end(struct mkcatalog_passparams *pp, gal_data_t *xybin)
+{
+  size_t i;
+  double *searr, *pearr, *osearr;
+  struct mkcatalogparams *p=pp->p;
+
+  /* The datasets and their pointers. */
+  gal_data_t *area  = pp->spectrum->next->next;
+  gal_data_t *sum   = area->next;
+  gal_data_t *esum  = area->next->next;
+  gal_data_t *proj  = area->next->next->next;
+  gal_data_t *eproj = area->next->next->next->next;
+  gal_data_t *oarea = area->next->next->next->next->next;
+  gal_data_t *osum  = area->next->next->next->next->next->next;
+  gal_data_t *oesum = area->next->next->next->next->next->next->next;
+
+  /* Apply corrections to the columns that need it. */
+  searr  = esum->array;
+  pearr  = eproj->array;
+  osearr = oesum->array;
+  for(i=0; i<p->objects->dsize[0]; ++i)
+    {
+      searr[i]  = sqrt( searr[i]  );
+      pearr[i]  = sqrt( pearr[i]  );
+      osearr[i] = sqrt( osearr[i] );
+    }
+
+  /* Convert the `double' type columns to `float'. The extra precision of
+     `double' was necessary when we were summing values in each slice. But
+     afterwards, it is not necessary at all (the measurement error is much
+     larger than a double-precision floating point number (15
+     decimals). But the extra space gained (double) is very useful in not
+     wasting too much memory and hard-disk space or online transfer time.*/
+  sum   = gal_data_copy_to_new_type_free(sum,   GAL_TYPE_FLOAT32);
+  esum  = gal_data_copy_to_new_type_free(esum,  GAL_TYPE_FLOAT32);
+  proj  = gal_data_copy_to_new_type_free(proj,  GAL_TYPE_FLOAT32);
+  eproj = gal_data_copy_to_new_type_free(eproj, GAL_TYPE_FLOAT32);
+  osum  = gal_data_copy_to_new_type_free(osum,  GAL_TYPE_FLOAT32);
+  oesum = gal_data_copy_to_new_type_free(oesum, GAL_TYPE_FLOAT32);
+
+  /* For the two area columns, find their maximum value and convert the
+     dataset to the smallest type that can hold them. */
+  parse_spectrum_uint32_to_best_type(&area);
+  parse_spectrum_uint32_to_best_type(&oarea);
+
+  /* List the datasets and write them into the pointer for this object
+     (exact copy of the statement in `parse_spectrum_pepare'). */
+  pp->spectrum->next->next                       = area;
+  area->next                                     = sum;
+  area->next->next                               = esum;
+  area->next->next->next                         = proj;
+  area->next->next->next->next                   = eproj;
+  area->next->next->next->next->next             = oarea;
+  area->next->next->next->next->next->next       = osum;
+  area->next->next->next->next->next->next->next = oesum;
+}
+
+
+
+
+/* Each spectrum is a multi-column table (note that the slice counter and
+   wavelength are written in the end):
+
+     Column 3:  Number of object pixels.
+     Column 4:  Sum of object pixel values.
+     Column 5:  Error in Column 2.
+     Column 6:  Sum over all 2D projection over whole specturm.
+     Column 7:  Error in Column 4.
+     Column 8:  Area of other labels in this slice.
+     Column 9:  Flux by other objects in projected area.
+     Column 10: Error in Column 9.
+ */
+static void
+parse_spectrum(struct mkcatalog_passparams *pp, gal_data_t *xybin)
+{
+  struct mkcatalogparams *p=pp->p;
+
+  gal_data_t *area;
+  float *st_v, *st_std;
+  uint32_t *narr, *oarr;
+  size_t *tsize, start_end_inc[2];
+  uint8_t *xybinarr = xybin ? xybin->array : NULL;
+  int32_t *O, *OO, *st_o, *objarr=p->objects->array;
+  size_t tid, *dsize=p->objects->dsize, num_increment=1;
+  double var, *sarr, *searr, *parr, *pearr, *osarr, *osearr;
+  size_t increment=0, pind=0, sind=0, ndim=p->objects->ndim, c[3];
+  float st, sval, *V=NULL, *ST=NULL, *std=p->std?p->std->array:NULL;
+
+  /* Prepare the columns to write in. */
+  tsize  = parse_spectrum_pepare(pp, start_end_inc, &st_o, &st_v, &st_std);
+  area   = pp->spectrum->next->next;
+  narr   = area->array;
+  sarr   = area->next->array;
+  searr  = area->next->next->array;
+  parr   = area->next->next->next->array;
+  pearr  = area->next->next->next->next->array;
+  oarr   = area->next->next->next->next->next->array;
+  osarr  = area->next->next->next->next->next->next->array;
+  osearr = area->next->next->next->next->next->next->next->array;
+
+  /* If tile-id isn't necessary, set `tid' to a blank value. */
+  tid = (p->std && p->std->size>1 && st_std == NULL) ? 0 : GAL_BLANK_SIZE_T;
+
+  /* Parse each contiguous patch of memory covered by this object. */
+  while( start_end_inc[0] + increment <= start_end_inc[1] )
+    {
+      /* Set the contiguous range to parse. The pixel-to-pixel counting
+         along the fastest dimension will be done over the `O' pointer. */
+      if( p->values        ) V  = st_v   + increment;
+      if( p->std && st_std ) ST = st_std + increment;
+      OO = ( O = st_o + increment ) + pp->tile->dsize[ndim-1];
+
+      /* Parse the tile. */
+      do
+        {
+          /* Only continue if this voxel is useful: it isn't NaN, or its
+             covered by the projected area or object's label. */
+          if( !isnan(*V) && (xybin && xybinarr[pind]==2) )
+            {
+              /* Get the error in measuing this pixel's flux. */
+              if(p->std)
+                {
+                  /* If the standard deviation is given on a tile
+                     structure, estimate the tile ID. */
+                  if(tid != GAL_BLANK_SIZE_T)
+                    {
+                      gal_dimension_index_to_coord(O-objarr, ndim, dsize, c);
+                      tid=gal_tile_full_id_from_coord(&p->cp.tl, c);
+                    }
+
+                  /* Get the error associated with this voxel. Note that if
+                     we are given a variance dataset already, there is no
+                     need to use `st*st', we can directly use `sval'. */
+                  sval = st_std ? *ST : (p->std->size>1?std[tid]:std[0]);
+                  st = p->variance ? sqrt(sval) : sval;
+                  var = (p->variance ? sval : st*st) + fabs(*V);
+                }
+              else var = NAN;
+
+
+              /* Projected spectra: see if we have a value of `2' in the
+                 `xybin' array (showing that there is atleast one non-blank
+                 element there over the whole spectrum.  */
+              parr [ sind ] += *V;
+              pearr[ sind ] += var;
+
+              /* Calculate the number of labeled/detected pixels that
+                 don't belong to this object. */
+              if(*O>0)
+                {
+                  if(*O==pp->object)
+                    {
+                      ++narr[ sind ];
+                      sarr  [ sind ] += *V;
+                      searr [ sind ] += var;
+                    }
+                  else
+                    {
+                      ++oarr [ sind ];
+                      osarr  [ sind ] += *V;
+                      osearr [ sind ] += var;
+                    }
+                }
+            }
+
+          /* Increment the pointers. */
+          if( xybin            ) ++pind;
+          if( p->values        ) ++V;
+          if( p->std && st_std ) ++ST;
+        }
+      while(++O<OO);
+
+      /* Increment to the next contiguous region of this tile. */
+      increment += ( gal_tile_block_increment(p->objects, tsize,
+                                              num_increment++, NULL) );
+
+      /* Increment the slice number, `sind', and reset the projection (2D)
+         index `pind' if we have just finished parsing a slice. */
+      if( (num_increment-1)%pp->tile->dsize[1]==0 )
+        {
+          pind=0;
+          ++sind;
+        }
+    }
+
+  /* Finalize the spectrum generation and clean up. */
+  parse_spectrum_end(pp, xybin);
+  free(tsize);
+
+  /* For a check.
+  gal_table_write(pp->spectrum, NULL, GAL_TABLE_FORMAT_BFITS,
+                  "spectrum.fits", "SPECTRUM", 0);
+  */
+}
+
+
+
+
+
 void
 parse_objects(struct mkcatalog_passparams *pp)
 {
@@ -115,10 +431,11 @@ parse_objects(struct mkcatalog_passparams *pp)
 
   double *oi=pp->oi;
   gal_data_t *xybin=NULL;
+  size_t *tsize=pp->tile->dsize;
   uint8_t *xybinarr=NULL, *u, *uf;
   float st, sval, *V=NULL, *SK=NULL, *ST=NULL;
   size_t d, pind=0, increment=0, num_increment=1;
-  int32_t *O, *OO, *C=NULL, *objects=p->objects->array;
+  int32_t *O, *OO, *C=NULL, *objarr=p->objects->array;
   float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
 
   /* If tile processing isn't necessary, set `tid' to a blank value. */
@@ -153,12 +470,13 @@ parse_objects(struct mkcatalog_passparams *pp)
                : NULL );
 
 
-  /* If an XY projection area is requested, we'll need to allocate an array
-     to keep the projected space.*/
-  if( oif[    OCOL_NUMALLXY ]
+  /* If an XY projection area is necessary, we'll need to allocate an array
+     to keep the projected space. */
+  if( p->spectrum
+      || oif[ OCOL_NUMALLXY ]
       || oif[ OCOL_NUMXY    ] )
     {
-      xybin=gal_data_alloc(NULL, GAL_TYPE_UINT8, 2, &pp->tile->dsize[1], NULL,
+      xybin=gal_data_alloc(NULL, GAL_TYPE_UINT8, 2, &tsize[1], NULL,
                            1, p->cp.minmapsize, NULL, NULL, NULL);
       xybinarr=xybin->array;
     }
@@ -173,7 +491,7 @@ parse_objects(struct mkcatalog_passparams *pp)
       if( p->values            ) V  = pp->st_v   + increment;
       if( p->sky && pp->st_sky ) SK = pp->st_sky + increment;
       if( p->std && pp->st_std ) ST = pp->st_std + increment;
-      OO = ( O = pp->st_o + increment ) + pp->tile->dsize[ndim-1];
+      OO = ( O = pp->st_o + increment ) + tsize[ndim-1];
 
       /* Parse the tile. */
       do
@@ -189,15 +507,15 @@ parse_objects(struct mkcatalog_passparams *pp)
 
 
               /* Add to the area of this object. */
+              if(xybin) xybinarr[ pind ]=1;
               if(oif[ OCOL_NUMALL   ]) oi[ OCOL_NUMALL ]++;
-              if(oif[ OCOL_NUMALLXY ]) xybinarr[ pind ]=1;
 
 
               /* Geometric coordinate measurements. */
               if(c)
                 {
                   /* Convert the index to coordinate. */
-                  gal_dimension_index_to_coord(O-objects, ndim, dsize, c);
+                  gal_dimension_index_to_coord(O-objarr, ndim, dsize, c);
 
                   /* If we need tile-ID, get the tile ID now. */
                   if(tid!=GAL_BLANK_SIZE_T)
@@ -238,9 +556,9 @@ parse_objects(struct mkcatalog_passparams *pp)
               if( p->values && !( p->hasblank && isnan(*V) ) )
                 {
                   /* General flux summations. */
+                  if(xybin) xybinarr[ pind ]=2;
                   if(oif[ OCOL_NUM    ]) oi[ OCOL_NUM     ]++;
                   if(oif[ OCOL_SUM    ]) oi[ OCOL_SUM     ] += *V;
-                  if(oif[ OCOL_NUMXY  ]) xybinarr[ pind ]=2;
 
                   /* Get the necessary clump information. */
                   if(p->clumps && *C>0)
@@ -321,12 +639,12 @@ parse_objects(struct mkcatalog_passparams *pp)
       while(++O<OO);
 
       /* Increment to the next contiguous region of this tile. */
-      increment += ( gal_tile_block_increment(p->objects, dsize,
+      increment += ( gal_tile_block_increment(p->objects, tsize,
                                               num_increment++, NULL) );
 
       /* If a 2D projection is requested, see if we should initialize (set
          to zero) the projection-index (`pind') not. */
-      if(xybin && (num_increment-1)%p->objects->dsize[1]==0 )
+      if(xybin && (num_increment-1)%tsize[1]==0 )
         pind=0;
     }
 
@@ -352,6 +670,10 @@ parse_objects(struct mkcatalog_passparams *pp)
       */
     }
 
+  /* Generate the Spectrum. */
+  if(p->spectrum)
+    parse_spectrum(pp, xybin);
+
   /* Clean up. */
   if(c)     free(c);
   if(sc)    free(sc);
@@ -384,6 +706,7 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
   double *ci, *cir;
   gal_data_t *xybin=NULL;
+  size_t *tsize=pp->tile->dsize;
   int32_t *O, *OO, *C=NULL, nlab;
   uint8_t *u, *uf, *cif=p->ciflag;
   float st, sval, *V=NULL, *SK=NULL, *ST=NULL;
@@ -438,9 +761,8 @@ parse_clumps(struct mkcatalog_passparams *pp)
     {
       xybin=gal_data_array_calloc(pp->clumpsinobj);
       for(i=0;i<pp->clumpsinobj;++i)
-        gal_data_initialize(&xybin[i], NULL, GAL_TYPE_UINT8, 2,
-                            &pp->tile->dsize[1], NULL, 1, p->cp.minmapsize,
-                            NULL, NULL, NULL);
+        gal_data_initialize(&xybin[i], NULL, GAL_TYPE_UINT8, 2, &tsize[1],
+                            NULL, 1, p->cp.minmapsize, NULL, NULL, NULL);
     }
 
 
@@ -453,7 +775,7 @@ parse_clumps(struct mkcatalog_passparams *pp)
       if( p->values            ) V  = pp->st_v   + increment;
       if( p->sky && pp->st_sky ) SK = pp->st_sky + increment;
       if( p->std && pp->st_std ) ST = pp->st_std + increment;
-      OO = ( O = pp->st_o + increment ) + pp->tile->dsize[ndim-1];
+      OO = ( O = pp->st_o + increment ) + tsize[ndim-1];
 
       /* Parse the tile */
       do
@@ -649,12 +971,12 @@ parse_clumps(struct mkcatalog_passparams *pp)
       while(++O<OO);
 
       /* Increment to the next contiguous region of this tile. */
-      increment += ( gal_tile_block_increment(p->objects, dsize,
+      increment += ( gal_tile_block_increment(p->objects, tsize,
                                               num_increment++, NULL) );
 
       /* If a 2D projection is requested, see if we should initialize (set
          to zero) the projection-index (`pind') not. */
-      if(xybin && (num_increment-1)%p->objects->dsize[1]==0 )
+      if(xybin && (num_increment-1) % tsize[1]==0 )
         pind=0;
     }
 
@@ -699,7 +1021,6 @@ void
 parse_median(struct mkcatalog_passparams *pp)
 {
   struct mkcatalogparams *p=pp->p;
-  size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
 
   float *V;
   double *ci;
@@ -707,9 +1028,10 @@ parse_median(struct mkcatalog_passparams *pp)
   int32_t *O, *OO, *C=NULL;
   gal_data_t **clumpsmed=NULL;
   size_t i, increment=0, num_increment=1;
-  size_t counter=0, *ccounter=NULL, tsize=pp->oi[OCOL_NUM];
-  gal_data_t *objmed=gal_data_alloc(NULL, p->values->type, 1, &tsize, NULL, 0,
-                                    p->cp.minmapsize, NULL, NULL, NULL);
+  size_t *tsize=pp->tile->dsize, ndim=p->objects->ndim;
+  size_t counter=0, *ccounter=NULL, tmpsize=pp->oi[OCOL_NUM];
+  gal_data_t *objmed=gal_data_alloc(NULL, p->values->type, 1, &tmpsize, NULL,
+                                    0, p->cp.minmapsize, NULL, NULL, NULL);
 
   /* Allocate space for the clump medians. */
   if(p->clumps)
@@ -726,9 +1048,10 @@ parse_median(struct mkcatalog_passparams *pp)
                                     __func__, "ccounter");
       for(i=0;i<pp->clumpsinobj;++i)
         {
-          tsize=pp->ci[ i * CCOL_NUMCOLS + CCOL_NUM ];
-          clumpsmed[i]=gal_data_alloc(NULL, p->values->type, 1, &tsize, NULL,
-                                      0, p->cp.minmapsize, NULL, NULL, NULL);
+          tmpsize=pp->ci[ i * CCOL_NUMCOLS + CCOL_NUM ];
+          clumpsmed[i]=gal_data_alloc(NULL, p->values->type, 1, &tmpsize,
+                                      NULL, 0, p->cp.minmapsize, NULL, NULL,
+                                      NULL);
         }
     }
 
@@ -740,7 +1063,7 @@ parse_median(struct mkcatalog_passparams *pp)
          along the fastest dimension will be done over the `O' pointer. */
       V = pp->st_v + increment;
       if(p->clumps) C = pp->st_c + increment;
-      OO = ( O = pp->st_o + increment ) + pp->tile->dsize[ndim-1];
+      OO = ( O = pp->st_o + increment ) + tsize[ndim-1];
 
       /* Parse the next contiguous region of this tile. */
       do
@@ -770,7 +1093,7 @@ parse_median(struct mkcatalog_passparams *pp)
       while(++O<OO);
 
       /* Increment to the next contiguous region of this tile. */
-      increment += ( gal_tile_block_increment(p->objects, dsize,
+      increment += ( gal_tile_block_increment(p->objects, tsize,
                                               num_increment++, NULL) );
     }
 
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 836d32f..090cfe2 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -552,12 +552,18 @@ ui_read_labels(struct mkcatalogparams *p)
   p->objects=gal_data_copy_to_new_type_free(p->objects, GAL_TYPE_INT32);
 
 
-  /* Currently MakeCatalog is only implemented for 2D images. */
+  /* Currently MakeCatalog is only implemented for 2D images or 3D cubes. */
   if(p->objects->ndim!=2 && p->objects->ndim!=3)
     error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, MakeCatalog "
           "currently only supports 2D or 3D datasets", p->objectsfile,
           p->cp.hdu, p->objects->ndim);
 
+  /* Make sure the `--spectrum' option is not given on a 2D image.  */
+  if(p->spectrum && p->objects->ndim!=3)
+    error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, but `--spectrum' "
+          "is currently only defined on 3D datasets", p->objectsfile,
+          p->cp.hdu, p->objects->ndim);
+
   /* See if the total number of objects is given in the header keywords. */
   keys[0].name="NUMLABS";
   keys[0].type=GAL_TYPE_SIZE_T;
@@ -689,6 +695,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
   size_t i;
 
   if(p->upperlimit) *values=1;
+  if(p->spectrum)   *values=*std=1;
 
   /* Go over all the object columns. Note that the objects and clumps (if
      the `--clumpcat' option is given) inputs are mandatory and it is not
@@ -883,6 +890,10 @@ ui_subtract_sky(struct mkcatalogparams *p)
     error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
           "the problem. For some reason, the size doesn't match", __func__,
           PACKAGE_BUGREPORT);
+
+  /* Inform the user that this operation is done (if necessary). */
+  if(!p->cp.quiet)
+    printf("  - Sky subtracted from input values.\n");
 }
 
 
@@ -1316,6 +1327,90 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
 
 
 
+/* When a spectrum is requested, the slice information (slice number and
+   slice WCS) is common to all different spectra. So instead of calculating
+   it every time, we'll just make it once here, then copy it for every
+   object.
+
+   The Slice information is going to be written in every spectrum. So we
+   don't want it to take too much space. Therefore, only when the number of
+   slices is less than 65535 (2^16-1), will we actually use a 32-bit
+   integer type for the slice number column.
+*/
+static void
+ui_preparations_spectrum_wcs(struct mkcatalogparams *p)
+{
+  double *xarr, *yarr, *zarr;
+  gal_data_t *x, *y, *z, *coords;
+  size_t i, numslices=p->objects->dsize[0];
+  size_t slicenumtype=numslices>=65535 ? GAL_TYPE_UINT32 : GAL_TYPE_UINT16;
+
+  /* A small sanity check. */
+  if(p->objects->ndim!=3)
+    error(EXIT_FAILURE, 0, "%s (hdu %s) is a %zuD dataset, but `--spectrum' "
+          "is currently only defined on 3D datasets", p->objectsfile,
+          p->cp.hdu, p->objects->ndim);
+
+  /* Allocate space for the slice number as well as the X and Y positions
+     for WCS conversion. Note that the `z' axis is going to be converted to
+     WCS later, so we'll just give it the basic information now.*/
+  x=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 0,
+                   p->cp.minmapsize, NULL, NULL, NULL);
+  y=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 0,
+                   p->cp.minmapsize, NULL, NULL, NULL);
+  z=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 0,
+                   p->cp.minmapsize, p->ctype[2], p->objects->wcs->cunit[2],
+                   "Slice WCS coordinates.");
+
+  /* Write values into the 3 coordinates. */
+  xarr=x->array; yarr=y->array; zarr=z->array;
+  for(i=0;i<numslices;++i) { zarr[i]=i+1; xarr[i]=yarr[i]=1; }
+
+
+  /* Convert the coordinates to WCS. We are doing this inplace to avoid too
+     much memory/speed consumption. */
+  coords=x;
+  coords->next=y;
+  coords->next->next=z;
+  gal_wcs_img_to_world(coords, p->objects->wcs, 1);
+
+  /* For a check.
+  for(i=0;i<numslices;++i)
+    printf("%g, %g, %g\n", xarr[i], yarr[i], zarr[i]);
+  exit(0);
+  */
+
+  /* Allocate the slice counter array (we are doing it again because we
+     want it to be in integer type now). */
+  p->specsliceinfo=gal_data_alloc(NULL, slicenumtype, 1, &numslices, NULL, 0,
+                                  p->cp.minmapsize, "SLICE", "counter",
+                                  "Slice number in cube (counting from 1).");
+  if(p->specsliceinfo->type==GAL_TYPE_UINT16)
+    for(i=0;i<numslices;++i) ((uint16_t *)(p->specsliceinfo->array))[i]=i+1;
+  else
+    for(i=0;i<numslices;++i) ((uint32_t *)(p->specsliceinfo->array))[i]=i+1;
+
+  /* Set the slice WCS column information. Note that `z' is now the WCS
+     coordinate value of the third dimension, and to avoid wasting extra
+     space (this column is repeated one very object's spectrum), we'll
+     convert it to a 32-bit floating point dataset. */
+  p->specsliceinfo->next=gal_data_copy_to_new_type(z, GAL_TYPE_FLOAT32);
+
+  /* For a final check.
+  gal_table_write(p->specsliceinfo, NULL, GAL_TABLE_FORMAT_BFITS,
+                  "specsliceinfo.fits", "test-debug",0);
+  */
+
+  /* Clean up. */
+  gal_data_free(x);
+  gal_data_free(y);
+  gal_data_free(z);
+}
+
+
+
+
+
 /* Sanity checks and preparations for the upper-limit magnitude. */
 static void
 ui_preparations_upperlimit(struct mkcatalogparams *p)
@@ -1372,9 +1467,9 @@ void
 ui_preparations(struct mkcatalogparams *p)
 {
   /* If no columns are requested, then inform the user. */
-  if(p->columnids==NULL)
-    error(EXIT_FAILURE, 0, "no columns requested, please run again with "
-          "`--help' for the full list of columns you can ask for");
+  if(p->columnids==NULL && p->spectrum==0)
+    error(EXIT_FAILURE, 0, "no measurements requested! Please run again "
+          "with `--help' for the possible list of measurements");
 
 
   /* Set the actual filenames to use. */
@@ -1406,6 +1501,14 @@ ui_preparations(struct mkcatalogparams *p)
   ui_one_tile_per_object(p);
 
 
+  /* If a spectrum is requested, generate the two WCS columns. */
+  if(p->spectrum)
+    {
+      ui_preparations_spectrum_wcs(p);
+      p->spectra=gal_data_array_calloc(p->numobjects);
+    }
+
+
   /* Allocate the reference random number generator and seed values. It
      will be cloned once for every thread. If the user hasn't called
      `envseed', then we want it to be different for every run, so we need
@@ -1585,7 +1688,7 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
mkcatalogparams *p)
 void
 ui_free_report(struct mkcatalogparams *p, struct timeval *t1)
 {
-  size_t d;
+  size_t d, i;
 
   /* The temporary arrays for WCS coordinates. */
   if(p->wcs_vo ) gal_list_data_free(p->wcs_vo);
@@ -1631,9 +1734,26 @@ ui_free_report(struct mkcatalogparams *p, struct timeval 
*t1)
   gal_data_free(p->upmask);
   gal_data_free(p->clumps);
   gal_data_free(p->objects);
+  gal_list_data_free(p->specsliceinfo);
   if(p->upcheckout) free(p->upcheckout);
   gal_data_array_free(p->tiles, p->numobjects, 0);
 
+  /* Clean up the spectra. */
+  if(p->spectra)
+    {
+      /* Note that each element of the array is the first node in a list of
+         datasets. So we can't free the first one with
+         `gal_list_data_free', we'll delete all the nodes after it in the
+         loop. */
+      for(i=0;i<p->numobjects;++i)
+        {
+          gal_list_data_free( p->spectra[i].next );
+          p->spectra[i].next=NULL;
+          gal_data_free_contents(&p->spectra[i]);
+        }
+      gal_data_array_free(p->spectra, p->numobjects, 0);
+    }
+
   /* If the Sky or its STD image were given in tiles, then we defined a
      tile structure to deal with them. The initialization of the tile
      structure is checked with its `ndim' element. */
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 1fae21e..596df75 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -91,6 +91,7 @@ enum option_keys_enum
   UI_KEY_SUBTRACTSKY,
   UI_KEY_SFMAGNSIGMA,
   UI_KEY_SFMAGAREA,
+  UI_KEY_SPECTRUM,
   UI_KEY_UPMASKFILE,
   UI_KEY_UPMASKHDU,
   UI_KEY_UPNUM,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8014ace..ee791da 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -483,8 +483,8 @@ Invoking MakeCatalog
 
 * MakeCatalog inputs and basic settings::  Input files and basic settings.
 * Upper-limit settings::        Settings for upper-limit measurements.
-* MakeCatalog output columns::  Available columns in MakeCatalog's output 
table.
-* MakeCatalog output file::     File names of MakeCatalog's output table.
+* MakeCatalog measurements::    Available measurements in MakeCatalog.
+* MakeCatalog output::          File names of MakeCatalog's output table.
 
 Match
 
@@ -11891,11 +11891,11 @@ dimension in its WCS matrix. Therefore, when the WCS 
is important for later
 processing, be sure that the input is aligned with the respective axises:
 all non-diagonal elements in the WCS matrix are zero.
 
-@cindex IFU
 @cindex Data cubes
 @cindex 3D data-cubes
 @cindex Cubes (3D data)
 @cindex Narrow-band image
+@cindex IFU: Integral Field Unit
 @cindex Integral field unit (IFU)
 One common application of this operator is the creation of pseudo
 broad-band or narrow-band 2D images from 3D data cubes. For example
@@ -17589,21 +17589,27 @@ it. For example the magnitudes, positions and 
elliptical properties of the
 galaxies that are in the image.
 
 MakeCatalog is Gnuastro's program for localized measurements over a
-dataset. The role of MakeCatalog in a scientific analysis and the benefits
-of this model of data analysis (were detection/segmentation is separated
-from measurement) is discussed in @url{https://arxiv.org/abs/1611.06387v1,
-Akhlaghi [2016]}. We strongly recommend reading this short paper for a
-better understanding of this methodology thus effective usage of
-MakeCatalog, in combination with the other Gnuastro's programs. However,
-that paper cannot undergo any more change, so this manual is the definitive
-guide.
-
-It is important to define your regions of interest @emph{before} running
-MakeCatalog. MakeCatalog is specialized in doing measurements accurately
-and efficiently. Therefore MakeCatalog will not do detection, segmentation,
-or defining apertures on requested positions in your dataset. Following
-Gnuastro's modularity principle, There are separate and highly specialized
-and customizable programs in Gnuastro for these other jobs:
+dataset. In other words, MakeCatalog is Gnuastro's program to convert
+low-level datasets (like images), to high level catalogs. The role of
+MakeCatalog in a scientific analysis and the benefits of its model (where
+detection/segmentation is separated from measurement) is discussed in
+@url{https://arxiv.org/abs/1611.06387v1, Akhlaghi [2016]}@footnote{A
+published paper cannot undergo any more change, so this manual is the
+definitive guide.} and summarized in @ref{Detection and catalog
+production}. We strongly recommend reading this short paper for a better
+understanding of this methodology. Understanding the effective usage of
+MakeCatalog, will thus also help effective use of other (lower-level)
+Gnuastro's programs like @ref{NoiseChisel} or @ref{Segment}.
+
+It is important to define your regions of interest for measurements
+@emph{before} running MakeCatalog. MakeCatalog is specialized in doing
+measurements accurately and efficiently. Therefore MakeCatalog will not do
+detection, segmentation, or defining apertures on requested positions in
+your dataset. Following Gnuastro's modularity principle, there are separate
+and highly specialized and customizable programs in Gnuastro for these
+other jobs as shown below (for a usage example in a real-world analysis,
+see @ref{General program usage tutorial} and @ref{Detecting large extended
+targets}).
 
 @itemize
 @item
@@ -17620,19 +17626,23 @@ and customizable programs in Gnuastro for these other 
jobs:
 @end itemize
 
 These programs will/can return labeled dataset(s) to be fed into
-MakeCatalog. The labeled dataset must have the same size/dimensions as the
-input, but only with integer valued pixels that have the label/counter for
-the feature the pixel belongs to.
-
-These labels are then directly used to make the necessary measurements. For
-example the flux weighted average position of all the pixels with a label
-of 42 will be considered as the central position@footnote{See
+MakeCatalog. A labeled dataset for measurement has the same size/dimensions
+as the input, but with integer valued pixels that have the label/counter
+for each sub-set of pixels that must be measured together. For example all
+the pixels covering one galaxy in an image, get the same label.
+
+The requested measurements are then done on similarly labeled pixels. The
+final result is a catalog where each row corresponds to the measurements on
+pixels with a specfic label. For example the flux weighted average position
+of all the pixels with a label of 42 will be written into the 42nd row of
+the output catalog/table's central position column@footnote{See
 @ref{Measuring elliptical parameters} for a discussion on this and the
-derivation of positional parameters.} of the 42nd row of the output
-catalog. Similarly, the sum of all these pixels will be the 42nd row in the
-brightness column and etc. Pixels with labels equal to or smaller than zero
-will be ignored by MakeCatalog. In other words, the number of rows in
-MakeCatalog's output is already known before running it.
+derivation of positional parameters, which includes the
+center.}. Similarly, the sum of all these pixels will be the 42nd row in
+the brightness column and etc. Pixels with labels equal to, or smaller than
+zero will be ignored by MakeCatalog. In other words, the number of rows in
+MakeCatalog's output is already known before running it (the maximum value
+of the labeled dataset).
 
 Before getting into the details of running MakeCatalog (in @ref{Invoking
 astmkcatalog}, we'll start with a discussion on the basics of its approach
@@ -17658,23 +17668,26 @@ steps is provided for readily adding your own new 
measurements/columns.
 @node Detection and catalog production, Quantifying measurement limits, 
MakeCatalog, MakeCatalog
 @subsection Detection and catalog production
 
-Most other common tools in low-level astronomical data-analysis (for
+Most existing common tools in low-level astronomical data-analysis (for
 example
 SExtractor@footnote{@url{https://www.astromatic.net/software/sextractor}})
-merge the two processes of detection and measurement into one. Gnuastro's
-modularized methodology to separating detection from measurements is
-therefore new to many experienced astronomers and deserves a short review
-here. Further discussion on the benefits of this methodology can be seen in
+merge the two processes of detection and measurement (catalog production)
+in one program. However, in light of Gnuastro's modularized approach
+(modeled on the Unix system) detection is separated from measurements and
+catalog production. This modularity is therefore new to many experienced
+astronomers and deserves a short review here. Further discussion on the
+benefits of this methodology can be seen in
 @url{https://arxiv.org/abs/1611.06387v1, Akhlaghi [2016]}.
 
-The raw output of any of the programs mentioned in @ref{MakeCatalog} can be
-directly fed into MakeCatalog to get to a fast catalog. This is good when
-no further customization is necessary and you want a fast/simple
-catalog. But the modular approach taken by Gnuastro has many benefits that
-will become clear as you get more experienced in astronomical data analysis
-and want to be more creative in using your valuable data for the exciting
-scientific project you are working on. In short the reasons for this
-modularity can be classified as below:
+As discussed in the introduction of @ref{MakeCatalog}, detection
+(identifying which pixels to do measurements on) can be done with different
+programs. Their outputs (a labeled dataset) can be directly fed into
+MakeCatalog to do the measurements and write the result as a
+catalog/table. Beyond that, Gnuastro's modular approach has many benefits
+that will become clear as you get more experienced in astronomical data
+analysis and want to be more creative in using your valuable data for the
+exciting scientific project you are working on. In short the reasons for
+this modularity can be classified as below:
 
 @itemize
 
@@ -17693,7 +17706,7 @@ same scale. However, it is imperative that the same 
pixels be used in
 measuring the colors of galaxies.
 
 To solve the problem, NoiseChisel can be run on the reference image to
-generate the labeled detection image. After wards, the labeled image can be
+generate the labeled detection image. Afterwards, the labeled image can be
 warped into the grid of the other color (using @ref{Warp}). MakeCatalog
 will then generate the same catalog for both colors (with the different
 labeled images). It is currently customary to warp the images to the same
@@ -18295,18 +18308,18 @@ including MakeCatalog.
 The various measurements/columns of MakeCatalog are requested as options,
 either on the command-line or in configuration files, see
 @ref{Configuration files}. The full list of available columns is available
-in @ref{MakeCatalog output columns}. Depending on the requested columns,
+in @ref{MakeCatalog measurements}. Depending on the requested columns,
 MakeCatalog needs more than one input dataset, for more details, please see
 @ref{MakeCatalog inputs and basic settings}. The upper-limit measurements
 in particular need several configuration options which are thoroughly
 discussed in @ref{Upper-limit settings}. Finally, in @ref{MakeCatalog
-output file} the output file(s) created by MakeCatalog are discussed.
+output} the output file(s) created by MakeCatalog are discussed.
 
 @menu
 * MakeCatalog inputs and basic settings::  Input files and basic settings.
 * Upper-limit settings::        Settings for upper-limit measurements.
-* MakeCatalog output columns::  Available columns in MakeCatalog's output 
table.
-* MakeCatalog output file::     File names of MakeCatalog's output table.
+* MakeCatalog measurements::    Available measurements in MakeCatalog.
+* MakeCatalog output::          File names of MakeCatalog's output table.
 @end menu
 
 @node MakeCatalog inputs and basic settings, Upper-limit settings, Invoking 
astmkcatalog, Invoking astmkcatalog
@@ -18333,7 +18346,7 @@ error). All such auxiliary input files have to have the 
same size (number
 of pixels in each dimension) as the input labeled image. Their numeric data
 type is irrelevant (they will be converted to 32-bit floating point
 internally). For the full list of available measurements, see
-@ref{MakeCatalog output columns}.
+@ref{MakeCatalog measurements}.
 
 The ``values'' dataset is used for measurements like brightness/magnitude,
 or flux-weighted positions. If it is a real image, by default it is assumed
@@ -18395,7 +18408,7 @@ options and general setting options are described below.
 @item -l STR
 @itemx --clumpsfile=STR
 The file containing the labeled clumps dataset when @option{--clumpscat} is
-called (see @ref{MakeCatalog output file}). When @option{--clumpscat} is
+called (see @ref{MakeCatalog output}). When @option{--clumpscat} is
 called, but this option isn't, MakeCatalog will look into the main input
 file (given as an argument) for the required extension/HDU (value to
 @option{--clumpshdu}).
@@ -18478,7 +18491,7 @@ magnitude}.
 @end table
 
 
-@node Upper-limit settings, MakeCatalog output columns, MakeCatalog inputs and 
basic settings, Invoking astmkcatalog
+@node Upper-limit settings, MakeCatalog measurements, MakeCatalog inputs and 
basic settings, Invoking astmkcatalog
 @subsubsection Upper-limit settings
 
 The upper-limit magnitude was discussed in @ref{Quantifying measurement
@@ -18618,17 +18631,34 @@ option.
 @end table
 
 
-@node MakeCatalog output columns, MakeCatalog output file, Upper-limit 
settings, Invoking astmkcatalog
-@subsubsection MakeCatalog output columns
+@node MakeCatalog measurements, MakeCatalog output, Upper-limit settings, 
Invoking astmkcatalog
+@subsubsection MakeCatalog measurements
 
 The final group of options particular to MakeCatalog are those that specify
-which columns should be written into the final output table. For each
-column there is an option, if it has been called on the command line or in
-any of the configuration files, it will included as a column in the output
-catalog in the same order (see @ref{Configuration file precedence}). Some
-of the columns apply to both objects and clumps and some are particular to
-only one of them. The latter cases are explicitly marked with [Objects] or
-[Clumps] to specify the catalog they will be placed in.
+which measurements/columns should be written into the final output
+table. The most basic measurement is one which only produces one final
+value for each label (for example its total brightness: a single
+number).
+
+In this case, all the different label's measurements can be written as one
+column in a final table/catalog that contains other columns for other
+similar single-number measurements. The majority of this section is devoted
+to MakeCatalog's single-valued measurements. However, MakeCatalog can also
+do measurements that produce more than one value for each label. Currently
+the only such measurement is generation of spectra from 3D cubes with the
+@option{--spectrum} option and it is discussed in the end of this section.
+
+Command-line options are used to identify which single-valued measurements
+you want in the final catalog(s) and in what order. If any of the options
+below is called on the command line or in any of the configuration files,
+it will be included as a column in the output catalog. The order of the
+columns is in the same order as the options were seen by MakeCatalog (see
+@ref{Configuration file precedence}). Some of the columns apply to both
+``objects'' and ``clumps'' and some are particular to only one of them (for
+the definition of ``objects'' and ``clumps'', see
+@ref{Segment}). Columns/options that are unique to one catalog (only
+objects, or only clumps), are explicitly marked with [Objects] or [Clumps]
+to specify the catalog they will be placed in.
 
 @table @option
 
@@ -19024,13 +19054,51 @@ with the first FITS axis in degrees.
 
 @end table
 
+@cindex 3D data-cubes
+@cindex Cubes (3D data)
+@cindex IFU: Integral Field Unit
+@cindex Integral field unit (IFU)
+@cindex Spectrum (of astronomical source)
+Above, all of MakeCatalog's single-valued measurements were listed. As
+mentioned in the start of this section, MakeCatalog can also do
+multi-valued measurements per label. Currently the only such measurement is
+the creation of spectra from 3D data cubes as discussed below:
+
+@table @option
+@item --spectrum
+Generate a spectrum (measurement along the first two FITS dimensions) for
+each label when the input dataset is a 3D data cube. With this option, a
+seprate table/spectrum will be generated for every label. If the output is
+a FITS file, each label's spectrum will be written into an extension of
+that file with a standard name of @code{SPECTRUM_NN} (the label will be
+replaced with @code{NN}). If the output is a plain text file, each label's
+spectrum will be written into a separate file with the suffix
+@file{spec-NN.txt}. See @ref{MakeCatalog output} for more on specifying
+MakeCatalog's output file.
+
+The spectra will contain one row for every slice (third FITS dimension) of
+the cube. Since the physical nature of the third dimension is different,
+two types of spectra (along with their errors) are measured: 1) Sum of
+values in each slice that only have the requested label. 2) Sum of values
+on the 2D projection of the whole label (the area of this projection can be
+requested with the @option{--areaxy} column above).
+
+Labels can overlap when they are projected onto the first two FITS
+dimensions (the spatial domain). To help separate them, MakeCatalog does a
+third measurement on each slice: the area, sum of values and error of all
+pixels that belong to other labels but overlap with the 2D projection. This
+can be used to see how reliable the emission line measurement is (on the
+projected spectra) and also if multiple lines (labeled regions) belong to
+the same physical object.
+@end table
 
 
 
-@node MakeCatalog output file,  , MakeCatalog output columns, Invoking 
astmkcatalog
-@subsubsection MakeCatalog output file
+@node MakeCatalog output,  , MakeCatalog measurements, Invoking astmkcatalog
+@subsubsection MakeCatalog output
 
-When complete, MakeCatalog will store its measurements as a table. If an
+After it has completed all the requested measurements (see @ref{MakeCatalog
+measurements}), MakeCatalog will store its measurements in table(s). If an
 output filename is given (see @option{--output} in @ref{Input output
 options}), the format of the table will be deduced from the name. When it
 isn't given, the input name will be appended with a @file{_cat} suffix (see
@@ -19040,12 +19108,19 @@ options}. @option{--tableformat} is also necessary 
when the requested
 output name is a FITS table (recall that FITS can accept ASCII and binary
 tables, see @ref{Table}).
 
-By default only a single catalog/table will be created for ``objects'',
-however, if @option{--clumpscat} is called, a secondary catalog/table will
-also be created. For more on ``objects'' and ``clumps'', see
-@ref{Segment}. In short, if you only have one set of labeled images, you
-don't have to worry about clumps (they are deactivated by default). Here is
-a full description of MakeCatalog's output options:
+By default (when @option{--spectrum} isn't called) only a single
+catalog/table will be created for ``objects'', however, if
+@option{--clumpscat} is called, a secondary catalog/table will also be
+created. For more on ``objects'' and ``clumps'', see @ref{Segment}. In
+short, if you only have one set of labeled images, you don't have to worry
+about clumps (they are deactivated by default).
+
+When @option{--spectrum} is called, it is not mandatory to specify any
+single-valued measurement columns. In this case, the output will only be
+the spectra of each labeled region. See the description of
+@option{--spectrum} in @ref{MakeCatalog measurements}.
+
+The full list of MakeCatalog's output options are elaborated below.
 
 @table @option
 @item -C
@@ -23864,10 +23939,11 @@ not have any zero values (a dimension of length zero 
is not defined).
 @end deftypefun
 
 @deftypefun void gal_data_free_contents (gal_data_t @code{*data})
-Free all the non-@code{NULL} pointers in @code{gal_data_t}. If @code{data}
-is actually a tile (@code{data->block!=NULL}, see @ref{Tessellation
-library}), then @code{data->array} is not freed. For a complete description
-of @code{gal_data_t} and its contents, see @ref{Generic data container}.
+Free all the non-@code{NULL} pointers in @code{gal_data_t} except for
+@code{next} and @code{block}. If @code{data} is actually a tile
+(@code{data->block!=NULL}, see @ref{Tessellation library}), then
+@code{data->array} is not freed. For a complete description of
+@code{gal_data_t} and its contents, see @ref{Generic data container}.
 @end deftypefun
 
 @deftypefun void gal_data_free (gal_data_t @code{*data})
@@ -26717,13 +26793,13 @@ function-like macro that is described below.
 
 @deftypefun {gal_data_t *} gal_tile_series_from_minmax (gal_data_t 
@code{*block}, size_t @code{*minmax}, size_t @code{number})
 Construct a list of tile(s) given coordinates of the minimum and maximum of
-each tile. The minimum and maximums are assumed to be inclusive. The
-returned pointer is an allocated @code{gal_data_t} array that can later be
-freed with @code{gal_data_array_free} (see @ref{Arrays of
-datasets}). Internally, each element of the output array points to the next
-element, so the output may also be treated as a list of datasets (see
-@ref{List of gal_data_t}) and passed onto the other functions described in
-this section.
+each tile. The minimum and maximums are assumed to be inclusive and in C
+order (slowest dimension first). The returned pointer is an allocated
+@code{gal_data_t} array that can later be freed with
+@code{gal_data_array_free} (see @ref{Arrays of datasets}). Internally, each
+element of the output array points to the next element, so the output may
+also be treated as a list of datasets (see @ref{List of gal_data_t}) and
+passed onto the other functions described in this section.
 
 The array keeping the minimum and maximum coordinates for each tile must
 have the following format. So in total @code{minmax} must have
diff --git a/lib/data.c b/lib/data.c
index 45b4762..96f7f37 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -777,12 +777,20 @@ gal_data_copy_to_allocated(gal_data_t *in, gal_data_t 
*out)
           "of dimensions, the dimensions are %zu and %zu respectively",
           __func__, out->ndim, in->ndim);
 
+  /* Free possibly allocated meta-data strings. */
+  if(out->name)    free(out->name);
+  if(out->unit)    free(out->unit);
+  if(out->comment) free(out->comment);
+
   /* Write the basic meta-data. */
   out->flag           = in->flag;
   out->next           = in->next;
   out->status         = in->status;
   out->disp_width     = in->disp_width;
   out->disp_precision = in->disp_precision;
+  gal_checkset_allocate_copy(in->name,    &out->name);
+  gal_checkset_allocate_copy(in->unit,    &out->unit);
+  gal_checkset_allocate_copy(in->comment, &out->comment);
 
   /* Do the copying. */
   switch(out->type)