[gnuastro-commits] master 7e7ca228: Library (arithmetic.h): stitch and noblank are two new operators

[Mohammad Akhlaghi]

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

    Library (arithmetic.h): stitch and noblank are two new operators
    
    Until now, there was no easy way to stitch multiple datasets together in
    Gnuastro (for example when adding the separate CCD amplifiers into one
    image), or a way remove blank elements of a dataset.
    
    With this commit the two new 'stitch' and 'noblank' operators have been
    added to Gnuastro's arithmetic library (and thus the Arithmetic program)
    for the two operations above.
    
    In the process, I noticed that the 'unique' operator wasn't defined in the
    library, but only in the Arithmetic program! With this commit, it has been
    moved into the library to be generally available.
---
 NEWS                        |  13 +++
 bin/arithmetic/arithmetic.c |  27 ++----
 bin/arithmetic/arithmetic.h |   1 -
 doc/gnuastro.texi           | 117 +++++++++++++++++++------
 lib/arithmetic.c            | 203 ++++++++++++++++++++++++++++++++++++++++++++
 lib/gnuastro/arithmetic.h   |   4 +
 6 files changed, 315 insertions(+), 50 deletions(-)

diff --git a/NEWS b/NEWS
index 7f498358..b78c18fa 100644
--- a/NEWS
+++ b/NEWS
@@ -16,6 +16,14 @@ See the end of the file for license conditions.
      written with the help of S. Zahra Hosseini Shahisavandi and Sepideh
      Eskandarlou.
 
+  Arithmetic:
+   - New operators:
+     - stitch: connect any number of input images along the given dimension
+       (for example to stitch the images of separate amplifiers of a CCD
+       into one image).
+     - noblank: Remove blank elements from the input and return a 1D output
+       (irrespective of how many dimensions the input had).
+
   Crop:
    --oneelemstdout: when a crop has a single pixel and this option is
      called, the single pixel's value will be printed on the standard
@@ -61,6 +69,11 @@ See the end of the file for license conditions.
      the same name). This was suggested by Ignacio Ruiz Cejudo and
      implemented by Sepideh Eskandarlou.
 
+  Library:
+   - GAL_ARITHMETIC_OP_UNIQUE: to find unique elements in dataset.
+   - GAL_ARITHMETIC_OP_NOBLANK: to remove blank elements from dataset.
+   - GAL_ARITHMETIC_OP_STITCH: to stitch multiple datasets.
+
 ** Removed features
 
 ** Changed features
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index d89c1d38..b4003836 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -79,6 +79,7 @@ pop_number_of_operands(struct arithmeticparams *p, int op, 
char *token_string,
   /* See if this operator needs any parameters. If so, pop them. */
   switch(op)
     {
+    case GAL_ARITHMETIC_OP_STITCH:
     case GAL_ARITHMETIC_OP_QUANTILE:
       numparams=1;
       break;
@@ -1043,18 +1044,6 @@ arithmetic_tofile(struct arithmeticparams *p, char 
*token, int freeflag)
 
 
 
-void
-arithmetic_unique(struct arithmeticparams *p, char *token, int operator)
-{
-  /* Pass the popped operand to the statistics library. */
-  gal_data_t *input = gal_statistics_unique(operands_pop(p, token), 1);
-  operands_add(p, NULL, input);
-}
-
-
-
-
-
 void
 arithmetic_add_dimension(struct arithmeticparams *p, char *token,
                          int operator)
@@ -1279,8 +1268,6 @@ arithmetic_set_operator(char *string, size_t 
*num_operands)
         { op=ARITHMETIC_OP_COLLAPSE_MEAN;         *num_operands=0; }
       else if (!strcmp(string, "collapse-number"))
         { op=ARITHMETIC_OP_COLLAPSE_NUMBER;       *num_operands=0; }
-      else if (!strcmp(string, "unique"))
-        { op=ARITHMETIC_OP_UNIQUE;                *num_operands=0; }
       else if (!strcmp(string, "add-dimension-slow"))
         { op=ARITHMETIC_OP_ADD_DIMENSION_SLOW;    *num_operands=0; }
       else if (!strcmp(string, "add-dimension-fast"))
@@ -1341,10 +1328,10 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
 
         case -1:
           /* This case is when the number of operands is itself an
-             operand. So except for sigma-clipping (that has other
-             parameters), the first popped operand must be an
-             integer number, we will use that to construct a linked
-             list of any number of operands within the single 'd1'
+             operand. So except for operators that have high-level
+             parameters (like sigma-clipping), the first popped operand
+             must be an integer number, we will use that to construct a
+             linked list of any number of operands within the single 'd1'
              pointer. */
           numop=pop_number_of_operands(p, operator, operator_string, &d2);
           for(i=0;i<numop;++i)
@@ -1416,10 +1403,6 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
           arithmetic_collapse(p, operator_string, operator);
           break;
 
-        case ARITHMETIC_OP_UNIQUE:
-          arithmetic_unique(p, operator_string, operator);
-          break;
-
         case ARITHMETIC_OP_ADD_DIMENSION_SLOW:
         case ARITHMETIC_OP_ADD_DIMENSION_FAST:
           arithmetic_add_dimension(p, operator_string, operator);
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 3df9d097..f1f68614 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -50,7 +50,6 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_COLLAPSE_MAX,
   ARITHMETIC_OP_COLLAPSE_MEAN,
   ARITHMETIC_OP_COLLAPSE_NUMBER,
-  ARITHMETIC_OP_UNIQUE,
   ARITHMETIC_OP_ADD_DIMENSION_SLOW,
   ARITHMETIC_OP_ADD_DIMENSION_FAST,
   ARITHMETIC_OP_REPEAT,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 1851ac7c..257e6347 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15357,6 +15357,19 @@ With this operator, this job is trivial:
 $ astarithmetic seg-crop.fits unique
 @end example
 
+@item noblank
+Remove all blank elements from the first popped operand.
+Since the blank pixels are being removed, the output dataset will always be 
single-dimensional, independent of the dimensionality of the input.
+
+Recall that by default, single-dimensional datasets are stored as a table 
column in the output.
+But you can use @option{--onedasimage} or @option{--onedonstdout} to 
respectively store them as a single-dimensional FITS array/image, or to print 
them on the standard output.
+
+For example with the command below, the non-blank pixel values of 
@file{cropped.fits} are printed on the command-line (the @option{--quiet} 
option is used to remove the extra information that Arithmetic prints as it 
reads the inputs, its version and its running time).
+
+@example
+$ astarithmetic cropped.fits nonblank --onedonstdout --quiet
+@end example
+
 @item size
 Size of the dataset along a given FITS (or FORTRAN) dimension (counting from 
1).
 The desired dimension should be the first popped operand and the dataset must 
be the second popped operand.
@@ -15659,6 +15672,54 @@ Every pixel in this 2D image will have the flux of the 
sum of the 100 slices.
 
 @table @command
 
+@item stitch
+Stitch (connect) any number of given images together along the given dimension.
+The output has the same number of dimensions as the input, but the number of 
pixels along the requested dimension will be different from the inputs.
+The @code{stitch} operator takes atleast three operands:
+@itemize
+@item
+The first popped operand (placed just before @code{stitch}) is the direction 
(dimension) that the images should be stitched along.
+The first FITS dimension is along the horizontal, therefore a value of 
@code{1} will stitch them horizontally.
+Similarly, giving a value of @code{2} will result in a vertical stitch.
+
+@item
+The second popped operand is the number of images that should be stitched.
+
+@item
+Depending on the value given to the second popped operand, @code{stitch} will 
pop the given number of datasets from the stack and stitch them along the given 
dimension.
+The popped images have to have the same number of pixels along the other 
dimension.
+The order of the stitching is defined by how they are placed in the 
command-line, not how they are popped (after being popped, they are placed in a 
list in the same order).
+@end itemize
+
+For example, in the commands below, we'll first crop out fixed sized regions 
of @mymath{100\times300} pixels of a larger image (@file{large.fits} first.
+In the first call of Arithmetic below, we will stitch the bottom set of crops 
together along the first (horizontal) axis.
+In the second Arithmetic call, we will stitch all 6 along both dimensions.
+
+@example
+## Crop the fixed-size regions of a larger image ('-O' is the
+## short form of the '--mode' option).
+$ astcrop large.fits -Oimg --section=1:100,1:300     -oa.fits
+$ astcrop large.fits -Oimg --section=101:200,1:300   -ob.fits
+$ astcrop large.fits -Oimg --section=201:300,1:300   -oc.fits
+$ astcrop large.fits -Oimg --section=1:100,301:600   -od.fits
+$ astcrop large.fits -Oimg --section=101:200,301:600 -oe.fits
+$ astcrop large.fits -Oimg --section=201:300,301:600 -of.fits
+
+## Stitch the bottom three crops into one image.
+$ astarithmetic a.fits b.fits c.fits 3 1 stitch -obottom.fits
+
+# Stitch all the 6 crops along both dimensions
+$ astarithmetic a.fits b.fits c.fits 3 1 stitch \
+                d.fits e.fits f.fits 3 1 stitch \
+                2 2 stitch -g1 -oall.fits
+@end example
+
+The start of the last command is like the one before it (stitching the bottom 
three crops along the first FITS dimension, producing a @mymath{300\times300} 
image).
+Later in the same command, we then stitch the top three crops horizontally 
(again, into a @mymath{300\times300} image)
+This leaves the the two @mymath{300\times300} images on the stack (see 
@ref{Reverse polish notation}).
+We finally stitch those two along the second dimension.
+This operator is therefore useful in scenarios like placing the CCD amplifiers 
into one image.
+
 @item collapse-sum
 Collapse the given dataset (second popped operand), by summing all elements 
along the first popped operand (a dimension in FITS standard: counting from 
one, from fastest dimension).
 The returned dataset has one dimension less compared to the input.
@@ -31235,12 +31296,20 @@ The latter two are the opposite.
 @deffnx Macro GAL_ARITHMETIC_OP_MEANVAL
 @deffnx Macro GAL_ARITHMETIC_OP_STDVAL
 @deffnx Macro GAL_ARITHMETIC_OP_MEDIANVAL
-Unary operand statistical operators that will return a single value for
-datasets of any size. These are just wrappers around similar functions in
-@ref{Statistical operations} and are included in @code{gal_arithmetic} only
-for completeness (to use easily in @ref{Arithmetic}). In your programs, it
-will probably be easier if you use those @code{gal_statistics_} functions
-directly.
+Unary operand statistical operators that will return a single value for 
datasets of any size.
+These are just wrappers around similar functions in @ref{Statistical 
operations} and are included in @code{gal_arithmetic} only for completeness (to 
use easily in @ref{Arithmetic}).
+In your programs, it will probably be easier if you use those 
@code{gal_statistics_} functions directly.
+@end deffn
+
+@deffn  Macro GAL_ARITHMETIC_OP_UNIQUE
+@deffnx Macro GAL_ARITHMETIC_OP_NOBLANK
+Unary operands that will remove some elements from the input dataset.
+The first will return the unique elements, and the second will return the 
non-blank elements.
+Due to the removal of elements, the dimensionality of the output will be lost.
+
+These are just wrappers over the @code{gal_statistics_unique} and 
@code{gal_blank_remove}.
+These are just wrappers around similar functions in @ref{Statistical 
operations} and are included in @code{gal_arithmetic} only for completeness (to 
use easily in @ref{Arithmetic}).
+In your programs, it will probably be easier if you use those 
@code{gal_statistics_} functions directly.
 @end deffn
 
 @deffn Macro GAL_ARITHMETIC_OP_ABS
@@ -31254,20 +31323,13 @@ Unary operand absolute-value operator.
 @deffnx Macro GAL_ARITHMETIC_OP_MEAN
 @deffnx Macro GAL_ARITHMETIC_OP_STD
 @deffnx Macro GAL_ARITHMETIC_OP_MEDIAN
-Multi-operand statistical operations. When @code{gal_arithmetic} is called
-with any of these operators, it will expect only a single operand that will
-be interpreted as a list of datasets (see @ref{List of gal_data_t}). These
-operators can work on multiple threads using the @code{numthreads}
-argument. See the discussion under the @code{min} operator in
-@ref{Arithmetic operators}.
+Multi-operand statistical operations.
+When @code{gal_arithmetic} is called with any of these operators, it will 
expect only a single operand that will be interpreted as a list of datasets 
(see @ref{List of gal_data_t}).
+These operators can work on multiple threads using the @code{numthreads} 
argument.
+See the discussion under the @code{min} operator in @ref{Arithmetic operators}.
 
-The output will be a single dataset with each of its elements replaced by
-the respective statistical operation on the whole list. The type of the
-output is determined from the operator (irrespective of the input type):
-for @code{GAL_ARITHMETIC_OP_MIN} and @code{GAL_ARITHMETIC_OP_MAX}, it will
-be the same type as the input, for @code{GAL_ARITHMETIC_OP_NUMBER}, the
-output will be @code{GAL_TYPE_UINT32} and for the rest, it will be
-@code{GAL_TYPE_FLOAT32}.
+The output will be a single dataset with each of its elements replaced by the 
respective statistical operation on the whole list.
+The type of the output is determined from the operator (irrespective of the 
input type): for @code{GAL_ARITHMETIC_OP_MIN} and @code{GAL_ARITHMETIC_OP_MAX}, 
it will be the same type as the input, for @code{GAL_ARITHMETIC_OP_NUMBER}, the 
output will be @code{GAL_TYPE_UINT32} and for the rest, it will be 
@code{GAL_TYPE_FLOAT32}.
 @end deffn
 
 @deffn Macro GAL_ARITHMETIC_OP_QUANTILE
@@ -31280,14 +31342,10 @@ The output type is the same as the inputs.
 @deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_MEAN
 @deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN
 @deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_NUMBER
-Similar to the operands above (including @code{GAL_ARITHMETIC_MIN}), except
-that when @code{gal_arithmetic} is called with these operators, it requires
-two arguments. The first is the list of datasets like before, and the
-second is the 2-element list of @mymath{\sigma}-clipping parameters. The first
-element in the parameters list is the multiple of sigma and the second is
-the termination criteria (see @ref{Sigma clipping}). The output type of
-@code{GAL_ARITHMETIC_OP_SIGCLIP_NUMBER} will be @code{GAL_TYPE_UINT32} and
-for the rest it will be @code{GAL_TYPE_FLOAT32}.
+Similar to the operands above (including @code{GAL_ARITHMETIC_MIN}), except 
that when @code{gal_arithmetic} is called with these operators, it requires two 
arguments.
+The first is the list of datasets like before, and the second is the 2-element 
list of @mymath{\sigma}-clipping parameters.
+The first element in the parameters list is the multiple of sigma and the 
second is the termination criteria (see @ref{Sigma clipping}).
+The output type of @code{GAL_ARITHMETIC_OP_SIGCLIP_NUMBER} will be 
@code{GAL_TYPE_UINT32} and for the rest it will be @code{GAL_TYPE_FLOAT32}.
 @end deffn
 
 @deffn  Macro GAL_ARITHMETIC_OP_MKNOISE_SIGMA
@@ -31311,6 +31369,11 @@ The first is the dataset, and the second is a single 
integer value.
 The output type is a single integer.
 @end deffn
 
+@deffn Macro GAL_ARITHMETIC_OP_STITCH
+Stitch a list of input datasets along the requested dimension.
+See the description of the @code{stitch} operator in Arithmetic 
(@ref{Dimensionality changing operators}).
+@end deffn
+
 @deffn Macro GAL_ARITHMETIC_OP_POW
 Binary operator to-power operator. When @code{gal_arithmetic} is called with 
any of these operators, it will expect two operands: raising the first by the 
second (returning a floating point, inputs can be integers).
 @end deffn
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 1b353853..032ff8df 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -705,6 +705,37 @@ arithmetic_from_statistics(int operator, int flags, 
gal_data_t *input)
 
 
 
+/* Call functions in the 'gnuastro/statistics' library. */
+static gal_data_t *
+arithmetic_to_oned(int operator, int flags, gal_data_t *input)
+{
+  gal_data_t *out=NULL;
+  int inplace=flags & GAL_ARITHMETIC_FLAG_FREE;
+
+  switch(operator)
+    {
+    case GAL_ARITHMETIC_OP_UNIQUE:
+      out=gal_statistics_unique(input, inplace);
+      break;
+    case GAL_ARITHMETIC_OP_NOBLANK:
+      out = inplace ? input : gal_data_copy(input);
+      gal_blank_remove(out);
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
+            "'%s' to fix the problem. The value '%d' to "
+            "'operator' is not recognized", __func__,
+            PACKAGE_BUGREPORT, operator);
+    }
+
+  /* Return the output dataset. */
+  return out;
+}
+
+
+
+
+
 
 
 
@@ -920,6 +951,150 @@ arithmetic_size(int operator, int flags, gal_data_t *in, 
gal_data_t *arg)
 
 
 
+static size_t
+arithmetic_stitch_sanity_check(gal_data_t *list, gal_data_t *fdim)
+{
+  float *fitsdim;
+  size_t c, dim;
+  gal_data_t *tmp;
+
+  /* Currently we only have the stitch operator for 2D datasets. */
+  if(list->ndim!=2)
+    error(EXIT_FAILURE, 0, "currently the 'stitch' operator only "
+          "works with 2D datasets (images) but you have given a "
+          "%zu dimensional dataset. Please get in touch with us "
+          "at '%s' to add this feature", list->ndim, PACKAGE_BUGREPORT);
+
+  /* Make sure that the dimention is an integer. */
+  fitsdim=fdim->array;
+  if( fitsdim[0] != ceil(fitsdim[0]) )
+    error(EXIT_FAILURE, 0, "the dimension operand (first popped "
+          "operand) to 'stitch' should be an integer, but it is "
+          "'%g'", fitsdim[0]);
+
+  /* Make sure that the requested dimension is not larger than the input's
+     number of dimensions. */
+  if(fitsdim[0]>list->ndim)
+    error(EXIT_FAILURE, 0, "the dimension operand (first popped "
+          "operand) to 'stitch' is %g, but your input datasets "
+          "only have %zu dimensions", fitsdim[0], list->ndim);
+
+  /* Convert the given dimension (in FITS standard) to C standard. */
+  dim = list->ndim - fitsdim[0];
+
+  /* Go through the list of datasets and make sure the dimensionality is
+     fine.*/
+  c=0;
+  for(tmp=list; tmp!=NULL; tmp=tmp->next)
+    {
+      /* Increment the counter. */
+      ++c;
+
+      /* Make sure they have the same numerical data type. */
+      if(tmp->type!=list->type)
+        error(EXIT_FAILURE, 0, "input dataset number %zu has a "
+              "numerical data type of '%s' while the first "
+              "input has a type of '%s'",
+              c, gal_type_name(tmp->type,1), gal_type_name(list->ndim,1));
+
+      /* Make sure they have the same number of dimensions. */
+      if(tmp->ndim!=list->ndim)
+        error(EXIT_FAILURE, 0, "input dataset number %zu has %zu "
+              "dimensions, while the first input has %zu dimensions",
+              c, tmp->ndim, list->ndim);
+
+      /* Make sure the length along the non-requested dimension in all the
+         inputs are the same. Note that since we currently only support 2D
+         arrays, and that 'dim' is in C standard (has a value of 0 or 1),
+         we can simply say '!dim' to point to the other dimension.*/
+      if( tmp->dsize[!dim]!=list->dsize[!dim])
+        error(EXIT_FAILURE, 0, "input dataset number %zu has %zu "
+              "elements along dimension number %d, while the first "
+              "has %zu elements along that dimension", c,
+              tmp->dsize[!dim], !dim==0 ? 2 : 1, list->dsize[!dim]);
+    }
+
+  /* Return the dimension number that must be used. */
+  return dim;
+}
+
+
+
+
+
+/* Stitch multiple operands along a given dimension. */
+static gal_data_t *
+arithmetic_stitch(int operator, int flags, gal_data_t *list,
+                  gal_data_t *fdim)
+{
+  void *oarr;
+  gal_data_t *tmp, *out;
+  uint8_t type=list->type;
+  size_t i, dim, sum, dsize[2];
+
+  /* In case we are dealing with a single-element list, just return it! */
+  if(list->next==NULL) return list;
+
+  /* Make sure everything is good! */
+  dim=arithmetic_stitch_sanity_check(list, fdim);
+
+  /* Find the size of the final output dataset. */
+  sum=0; for(tmp=list; tmp!=NULL; tmp=tmp->next) sum+=tmp->dsize[dim];
+  dsize[0] = dim==0 ? sum            : list->dsize[0];
+  dsize[1] = dim==0 ? list->dsize[1] : sum;
+
+  /* Allocate the output dataset. */
+  out=gal_data_alloc(NULL, list->type, list->ndim, dsize, list->wcs,
+                     0, list->minmapsize, list->quietmmap, "Stitched",
+                     list->unit, "Stitched dataset");
+
+  /* Write the individual datasets into the output. Note that 'dim' is the
+     dimension counter of the C standard. */
+  oarr=out->array;
+  for(tmp=list; tmp!=NULL; tmp=tmp->next)
+    {
+      switch(dim)
+        {
+
+        /* Vertical stitching (second FITS axis is the vertical axis). */
+        case 0:
+          memcpy(oarr,tmp->array, gal_type_sizeof(type)*tmp->size);
+          oarr += gal_type_sizeof(type)*tmp->size;
+          break;
+
+        /* Horizontal stitching (first FITS axis is the horizontal axis) */
+        case 1:
+
+          /* Copy row-by-row. */
+          for(i=0;i<dsize[0];++i)
+            memcpy(gal_pointer_increment(oarr,       i*dsize[1],      type),
+                   gal_pointer_increment(tmp->array, i*tmp->dsize[1], type),
+                   gal_type_sizeof(type)*tmp->dsize[1]);
+
+          /* Copying has finished, increment the start for the next
+             image. Note that in this scenario, the starting pixel for the
+             next image is on the first row, but tmp->dsize[1] pixels
+             away. */
+          oarr += gal_type_sizeof(type)*tmp->dsize[1];
+          break;
+
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
+                "to find and fix the problem. The value of 'dim' is "
+                "'%zu' is not understood", __func__, PACKAGE_BUGREPORT,
+                dim);
+        }
+    }
+
+  /* Clean up and return */
+  if(flags & GAL_ARITHMETIC_FLAG_FREE)
+    gal_list_data_free(list);
+  return out;
+}
+
+
+
+
 
 
 
@@ -2436,6 +2611,12 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
   else if (!strcmp(string, "sigclip-std"))
     { op=GAL_ARITHMETIC_OP_SIGCLIP_STD;       *num_operands=-1; }
 
+  /* To one-dimension (only based on values). */
+  else if (!strcmp(string, "unique"))
+    { op=GAL_ARITHMETIC_OP_UNIQUE;            *num_operands=1;  }
+  else if (!strcmp(string, "noblank"))
+    { op=GAL_ARITHMETIC_OP_NOBLANK;           *num_operands=1;  }
+
   /* Adding noise operators. */
   else if (!strcmp(string, "mknoise-sigma"))
     { op=GAL_ARITHMETIC_OP_MKNOISE_SIGMA;     *num_operands=2; }
@@ -2448,6 +2629,10 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
   else if (!strcmp(string, "size"))
     { op=GAL_ARITHMETIC_OP_SIZE;              *num_operands=2;  }
 
+  /* Stitching */
+  else if (!strcmp(string, "stitch"))
+    { op=GAL_ARITHMETIC_OP_STITCH;            *num_operands=-1; }
+
   /* Conditional operators. */
   else if (!strcmp(string, "lt" ))
     { op=GAL_ARITHMETIC_OP_LT;                *num_operands=2;  }
@@ -2598,6 +2783,9 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_STDVAL:          return "stdvalue";
     case GAL_ARITHMETIC_OP_MEDIANVAL:       return "medianvalue";
 
+    case GAL_ARITHMETIC_OP_UNIQUE:          return "unique";
+    case GAL_ARITHMETIC_OP_NOBLANK:         return "noblank";
+
     case GAL_ARITHMETIC_OP_MIN:             return "min";
     case GAL_ARITHMETIC_OP_MAX:             return "max";
     case GAL_ARITHMETIC_OP_NUMBER:          return "number";
@@ -2616,6 +2804,7 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM: return "mknoise-uniform";
 
     case GAL_ARITHMETIC_OP_SIZE:            return "size";
+    case GAL_ARITHMETIC_OP_STITCH:          return "stitch";
 
     case GAL_ARITHMETIC_OP_TO_UINT8:        return "uchar";
     case GAL_ARITHMETIC_OP_TO_INT8:         return "char";
@@ -2743,6 +2932,13 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       out=arithmetic_from_statistics(operator, flags, d1);
       break;
 
+    /* Return 1D array (only values). */
+    case GAL_ARITHMETIC_OP_UNIQUE:
+    case GAL_ARITHMETIC_OP_NOBLANK:
+      d1 = va_arg(va, gal_data_t *);
+      out=arithmetic_to_oned(operator, flags, d1);
+      break;
+
     /* Absolute operator. */
     case GAL_ARITHMETIC_OP_ABS:
       d1 = va_arg(va, gal_data_t *);
@@ -2799,6 +2995,13 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       out=arithmetic_size(operator, flags, d1, d2);
       break;
 
+    /* Stitch multiple datasets. */
+    case GAL_ARITHMETIC_OP_STITCH:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      out=arithmetic_stitch(operator, flags, d1, d2);
+      break;
+
     /* Conversion operators. */
     case GAL_ARITHMETIC_OP_TO_UINT8:
     case GAL_ARITHMETIC_OP_TO_INT8:
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index aafcd908..e21c2cc4 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -143,6 +143,8 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_MEANVAL,      /* Mean value of array.                  */
   GAL_ARITHMETIC_OP_STDVAL,       /* Standard deviation value of array.    */
   GAL_ARITHMETIC_OP_MEDIANVAL,    /* Median value of array.                */
+  GAL_ARITHMETIC_OP_UNIQUE,       /* Only return unique elements.          */
+  GAL_ARITHMETIC_OP_NOBLANK,      /* Only keep non-blank elements.         */
 
   GAL_ARITHMETIC_OP_MIN,          /* Minimum per pixel of multiple arrays. */
   GAL_ARITHMETIC_OP_MAX,          /* Maximum per pixel of multiple arrays. */
@@ -163,6 +165,8 @@ enum gal_arithmetic_operators
 
   GAL_ARITHMETIC_OP_SIZE,         /* Size of the dataset along an axis     */
 
+  GAL_ARITHMETIC_OP_STITCH,       /* Stitch multiple datasets together.    */
+
   GAL_ARITHMETIC_OP_TO_UINT8,     /* Convert to uint8_t.                   */
   GAL_ARITHMETIC_OP_TO_INT8,      /* Convert to int8_t.                    */
   GAL_ARITHMETIC_OP_TO_UINT16,    /* Convert to uint16_t.                  */