[gnuastro-commits] master 5dc054e: Library (arithmetic.h): new operators to add noise

[Mohammad Akhlaghi]

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

    Library (arithmetic.h): new operators to add noise
    
    Until now, the only way to add noise was with the MakeNoise program. This
    had two limitations: 1) it didn't work on tables, 2) it wasn't easy to
    blend with other types of operations on the pixels.
    
    With this commit, the core operation of MakeNoise has now been moved into
    the arithmetic library as two operators ('mknoise-sigma' and
    'mknoise-poisson'). This enables both the features above.
    
    However, it doesn't yet support instrumental noise, but that can be added
    later, so its still too early to completely remove MakeNoise! Although it
    isn't too far out of sight now: possibly with a third operator, I just
    don't have enough time to implement it now.
---
 NEWS                        |   6 ++
 bin/arithmetic/args.h       |  13 ++++
 bin/arithmetic/arithmetic.c |  13 +++-
 bin/arithmetic/main.h       |   1 +
 bin/arithmetic/ui.h         |   1 +
 bin/table/arithmetic.c      |  12 ++++
 doc/gnuastro.texi           | 124 ++++++++++++++++++++++++--------------
 lib/arithmetic.c            | 141 ++++++++++++++++++++++++++++++++++++++++++++
 lib/gnuastro/arithmetic.h   |   3 +
 9 files changed, 269 insertions(+), 45 deletions(-)

diff --git a/NEWS b/NEWS
index ceb1926..a8ad81f 100644
--- a/NEWS
+++ b/NEWS
@@ -49,9 +49,15 @@ See the end of the file for license conditions.
      - asinh: Inverse of hyperbolic sine.
      - acosh: Inverse of hyperbolic cosine.
      - atanh: Inverse of hyperbolic tangent.
+     - mknoise-sigma: Add Gaussian noise with the fixed sigma.
+     - mknoise-poisson: Add Poisson noise with the given background.
      - counts-to-mag: Convert counts to magnitudes with given zero point.
      - counts-to-jy: Convert counts to Janskys through a zero point based
           on AB magnitudes.
+   --envseed: new option to get random number generator settings for the
+     new 'mknoise-sigma' and 'mknoise-poisson' operators from the
+     environment for reproducibility (see "Generating random numbers"
+     section in manual).
 
   ConvertType:
    --globalhdu: Use a single HDU identifier for all the input files
diff --git a/bin/arithmetic/args.h b/bin/arithmetic/args.h
index 8e765ee..355fe9f 100644
--- a/bin/arithmetic/args.h
+++ b/bin/arithmetic/args.h
@@ -69,6 +69,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "envseed",
+      UI_KEY_ENVSEED,
+      0,
+      0,
+      "Use GSL_RNG_SEED env. for '--rowrandom'.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->envseed,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 963633f..d8c2343 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1217,7 +1217,7 @@ static void
 arithmetic_operator_run(struct arithmeticparams *p, int operator,
                         char *operator_string, size_t num_operands)
 {
-  size_t i;
+  size_t i, one=1;
   unsigned int numop;
   gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
   int flags = ( GAL_ARITHMETIC_INPLACE | GAL_ARITHMETIC_FREE
@@ -1267,6 +1267,17 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
                 num_operands, operator_string);
         }
 
+      /* Save 'envseed' as third operand if necessary. */
+      switch(operator)
+        {
+        case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
+        case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
+          d3=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &one, NULL, 0, -1, 1,
+                            NULL, NULL, NULL);
+          ((uint8_t *)(d3->array))[0]=p->envseed;
+          break;
+        }
+
       /* Run the arithmetic operation. Note that 'gal_arithmetic'
          is a variable argument function (like printf). So the
          number of arguments it uses depend on the operator. So
diff --git a/bin/arithmetic/main.h b/bin/arithmetic/main.h
index f53c46a..bbe9b37 100644
--- a/bin/arithmetic/main.h
+++ b/bin/arithmetic/main.h
@@ -87,6 +87,7 @@ struct arithmeticparams
   int        wcs_collapsed;  /* If the internal WCS is already collapsed.*/
 
   /* Internal: */
+  uint8_t          envseed;  /* To setup the random number generator.   */
   struct operand *operands;  /* The operands linked list.               */
   time_t           rawtime;  /* Starting time of the program.           */
 };
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/ui.h
index 8a74d81..f9ab07e 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/ui.h
@@ -46,6 +46,7 @@ enum option_keys_enum
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
+  UI_KEY_ENVSEED         = 1000,
 };
 
 
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index b964535..d06baae 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -769,6 +769,7 @@ arithmetic_operator_run(struct tableparams *p,
                         struct gal_arithmetic_set_params *setprm,
                         gal_data_t **stack)
 {
+  size_t one=1;
   gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
   int flags = ( GAL_ARITHMETIC_INPLACE | GAL_ARITHMETIC_FREE
                 | GAL_ARITHMETIC_NUMOK );
@@ -812,6 +813,17 @@ arithmetic_operator_run(struct tableparams *p,
                 arithmetic_operator_name(token->operator));
         }
 
+      /* Save 'envseed' as third operand if necessary. */
+      switch(token->operator)
+        {
+        case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
+        case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
+          d3=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &one, NULL, 0, -1, 1,
+                            NULL, NULL, NULL);
+          ((uint8_t *)(d3->array))[0]=p->envseed;
+          break;
+        }
+
       /* Run the arithmetic operation. Note that 'gal_arithmetic' is a
          variable argument function (like printf). So the number of
          arguments it uses depend on the operator. In other words, when the
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f8241ad..55d6002 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12789,6 +12789,47 @@ The internal conversion of C will be used.
 Convert the type of the popped operand to 64-bit (double precision) floating 
point (see @ref{Numeric data types}).
 The internal conversion of C will be used.
 
+@item mknoise-sigma
+Add a fixed noise (Gaussian standard deviation) to each element of the input 
dataset.
+This operator takes two arguments: the top/first popped operand is the noise 
standard deviation, the next popped operand is the dataset that the noise 
should be added to.
+You can use the @option{--envseed} option to fix the random number generator 
seed (and thus get a reproducible result).
+For more on @option{--envseed}, see @ref{Generating random numbers}.
+
+For example with the first command below, @file{image.fits} will be degraded 
by a noise of standard deviation 3 units.
+@example
+$ astarithmetic image.fits 3 mknoise-sigma
+@end example
+
+Alternatively, you can use this operator within column arithmetic in the Table 
program, to generate a random number like below (centered on 0, with 
@mymath{\sigma=3}) like the first command below.
+With the second command, you can put it into a shell variable for later usage.
+
+@example
+$ echo 0 | asttable -c'arith $1 3 mknoise-sigma'
+$ value=$(echo 0 | asttable -c'arith $1 3 mknoise-sigma')
+$ echo $value
+@end example
+
+You can also use this operator in combination with AWK to easily generate an 
arbitrarily large table with random columns.
+In the example below, we'll create a two column table with 20 rows.
+The first column will be centered on 5 and @mymath{\sigma_1=2}, the second 
will be centered on 10 and @mymath{\sigma_2=3}:
+
+@example
+$ echo 5 10 \
+       | awk '@{for(i=0;i<20;++i) print $1, $2@}' \
+       | asttable -c'arith $1 2 mknoise-sigma' \
+                  -c'arith $2 3 mknoise-sigma'
+@end example
+
+By adding an extra @option{--output=random.fits}, the table will be saved into 
a file called @file{random.fits}, and you can change the @code{i<20} to 
@code{i<5000} to have 5000 rows instead.
+Of course, if your input table has different values in the desired column the 
noisy distribution will be centered on each input element, but all will have 
the same scatter/sigma.
+
+@item mknoise-poisson
+Add Poisson noise to each element of the input dataset (see @ref{Photon 
counting noise}).
+This operator takes two arguments: the top/first popped operand is the 
background, the next popped operand is the dataset that the noise should be 
added to.
+
+Except for the noise-model, this operator is very similar to 
@code{mknoise-sigma} and the examples there apply here too.
+The main difference with @code{mknoise-sigma} is that in a Poisson 
distribution the scatter/sigma will depend on each element's value.
+
 @item size
 Size of the dataset along a given FITS/Fortran dimension (counting from 1).
 The desired dimension should be the first popped operand and the dataset must 
be the second popped operand.
@@ -12986,6 +13027,10 @@ $ astarithmetic data.fits --wcsfile=other.fits 
-ofinal.fits
 HDU/extension to read the WCS within the file given to @option{--wcsfile}.
 For more, see the description of @option{--wcsfile}.
 
+@item --envseed
+Use the environment for the random number generator settings in operators that 
need them (for example @code{mknoise-sigma}).
+This is very important for obtaining reproducible results, for more see 
@ref{Generating random numbers}.
+
 @item -O
 @itemx --onedasimage
 When final dataset to write as output only has one dimension, write it as a 
FITS image/array.
@@ -19907,33 +19952,31 @@ The subsequent programs which use GSL's random number 
generators will hence fort
 In case you want to set fixed values for these parameters every time you use 
the GSL random number generator, you can add these two lines to your 
@file{.bashrc} startup script@footnote{Don't forget that if you are going to 
give your scripts (that use the GSL random number generator) to others you have 
to make sure you also tell them to set these environment variable separately.
 So for scripts, it is best to keep all such variable definitions within the 
script, even if they are within your @file{.bashrc}.}, see @ref{Installation 
directory}.
 
-@cartouche
-@noindent
-@strong{NOTE:} If the two environment variables @code{GSL_RNG_TYPE} and 
@code{GSL_RNG_SEED} are defined, GSL will report them by default, even if you 
don't use the @option{--envseed} option.
-For example you can see the top few lines of the output of MakeProfiles:
+@strong{IMPORTANT NOTE:} If the two environment variables @code{GSL_RNG_TYPE} 
and @code{GSL_RNG_SEED} are defined, GSL will report them by default, even if 
you don't use the @option{--envseed} option.
+For example, see this call to MakeProfiles:
 
 @example
-$ export GSL_RNG_TYPE="taus"
+$ export GSL_RNG_TYPE=taus
 $ export GSL_RNG_SEED=345
-$ astmkprof -s1 --kernel=gaussian,2,5 --envseed
+$ astmkprof -s1 --kernel=gaussian,2,5
 GSL_RNG_TYPE=taus
 GSL_RNG_SEED=345
-MakeProfiles A.B started on DDD MMM DD HH:MM:SS YYYY
+MakeProfiles V.VV started on DDD MMM DDD HH:MM:SS YYYY
   - Building one gaussian kernel
-  - Random number generator (RNG) type: ranlxs1
-  - RNG seed for all profiles: 345
+  - Random number generator (RNG) type: taus
+  - Basic RNG seed: 1618960836
   ---- ./kernel.fits created.
-MakeProfiles finished in 0.111271 seconds
+  -- Output: ./kernel.fits
+MakeProfiles finished in 0.068945 seconds
 @end example
 
 @noindent
 @cindex Seed, Random number generator
 @cindex Random number generator, Seed
-The first two output lines (showing the names of the environment variables) 
are printed by GSL before MakeProfiles actually starts generating random 
numbers.
-The Gnuastro programs will report the values they use independently, you 
should check them for the final values used.
-For example if @option{--envseed} is not given, @code{GSL_RNG_SEED} will not 
be used and the last line shown above will not be printed.
-In the case of MakeProfiles, each profile will get its own seed value.
-@end cartouche
+The first two output lines (showing the names and values of the GSL 
environment variables) are printed by GSL before MakeProfiles actually starts 
generating random numbers.
+Gnuastro's programs will report the actual values they use independently 
(after the name of the program), you should check them for the final values 
used, not GSL's printed values.
+In the example above, did you notice how the random number generator seed 
above is different between GSL and MakeProfiles?
+However, if @option{--envseed} was given, both printed seeds would be the same.
 
 
 @node Invoking astmknoise,  , Noise basics, MakeNoise
@@ -26120,35 +26163,19 @@ description of @code{gal_wcs_world_to_img} for more 
details.
 @node Arithmetic on datasets, Tessellation library, World Coordinate System, 
Gnuastro library
 @subsection Arithmetic on datasets (@file{arithmetic.h})
 
-When the dataset's type and other information are already known, any
-programming language (including C) provides some very good tools for
-various operations (including arithmetic operations like addition) on the
-dataset with a simple loop. However, as an author of a program, making
-assumptions about the type of data, its dimensions and other basic
-characteristics will come with a large processing burden.
-
-For example if you always read your data as double precision floating
-points for a simple operation like addition with an integer constant, you
-will be wasting a lot of CPU and memory when the input dataset is
-@code{int32} type for example (see @ref{Numeric data types}). This overhead
-may be small for small images, but as you scale your process up and work
-with hundred/thousands of files that can be very large, this overhead will
-take a significant portion of the processing power. The functions and
-macros in this section are designed precisely for this purpose: to allow
-you to do any of the defined operations on any dataset with no overhead (in
-the native type of the dataset).
-
-Gnuastro's Arithmetic program uses the functions and macros of this
-section, so please also have a look at the @ref{Arithmetic} program and in
-particular @ref{Arithmetic operators} for a better description of the
-operators discussed here.
-
-The main function of this library is @code{gal_arithmetic} that is
-described below. It can take an arbitrary number of arguments as operands
-(depending on the operator, similar to @code{printf}). Its first two
-arguments are integers specifying the flags and operator. So first we will
-review the constants for the recognized flags and operators and discuss
-them, then introduce the actual function.
+When the dataset's type and other information are already known, any 
programming language (including C) provides some very good tools for various 
operations (including arithmetic operations like addition) on the dataset with 
a simple loop.
+However, as an author of a program, making assumptions about the type of data, 
its dimensions and other basic characteristics will come with a large 
processing burden.
+
+For example if you always read your data as double precision floating points 
for a simple operation like addition with an integer constant, you will be 
wasting a lot of CPU and memory when the input dataset is @code{int32} type for 
example (see @ref{Numeric data types}).
+This overhead may be small for small images, but as you scale your process up 
and work with hundred/thousands of files that can be very large, this overhead 
will take a significant portion of the processing power.
+The functions and macros in this section are designed precisely for this 
purpose: to allow you to do any of the defined operations on any dataset with 
no overhead (in the native type of the dataset).
+
+Gnuastro's Arithmetic program uses the functions and macros of this section, 
so please also have a look at the @ref{Arithmetic} program and in particular 
@ref{Arithmetic operators} for a better description of the operators discussed 
here.
+
+The main function of this library is @code{gal_arithmetic} that is described 
below.
+It can take an arbitrary number of arguments as operands (depending on the 
operator, similar to @code{printf}).
+Its first two arguments are integers specifying the flags and operator.
+So first we will review the constants for the recognized flags and operators 
and discuss them, then introduce the actual function.
 
 @deffn  Macro GAL_ARITHMETIC_INPLACE
 @deffnx Macro GAL_ARITHMETIC_FREE
@@ -26354,6 +26381,15 @@ the termination criteria (see @ref{Sigma clipping}). 
The output type of
 for the rest it will be @code{GAL_TYPE_FLOAT32}.
 @end deffn
 
+@deffn  Macro GAL_ARITHMETIC_OP_MKNOISE_SIGMA
+@deffnx Macro GAL_ARITHMETIC_OP_MKNOISE_POISSON
+Add noise to the input dataset.
+Both operators take three arguments: the first is the input data set (can have 
any dimensionality or number of elements.
+The second argument is the noise specifier (a single element, of any type): 
for a fixed-sigma noise, it is the Gaussian standard deviation, for the Poisson 
noise, it is the background (see @ref{Photon counting noise}).
+The third argument should a dataset with only one element, a 
@code{GAL_TYPE_UINT8} type and only two acceptable values 0 (to signifiy that 
the random number generator, RNG, seed can be automatically set for each run) 
and 1 (to signify that the RNG seed should be read from the @code{GAL_RNG_SEED} 
environment variable for a reproducible result).
+For more on setting the RNG seed, see @ref{Generating random numbers}.
+@end deffn
+
 @deffn Macro GAL_ARITHMETIC_OP_SIZE
 Size operator that will return a single value for datasets of any kind. When 
@code{gal_arithmetic} is called with this operator, it requires two arguments.
 The first is the dataset, and the second is a single integer value.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 8f98575..cae5345 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -29,6 +29,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 #include <stdarg.h>
 
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
 #include <gnuastro/list.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/units.h>
@@ -39,6 +42,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/statistics.h>
 #include <gnuastro/arithmetic.h>
 
+#include <gnuastro-internal/checkset.h>
 #include <gnuastro-internal/arithmetic-internal.h>
 
 /* Headers for each binary operator. Since they heavily involve macros,
@@ -651,6 +655,125 @@ arithmetic_from_statistics(int operator, int flags, 
gal_data_t *input)
 
 
 /***********************************************************************/
+/***************               Adding noise               **************/
+/***********************************************************************/
+
+/* The size operator. Reports the size along a given dimension. */
+static gal_data_t *
+arithmetic_mknoise(int operator, int flags, gal_data_t *in, gal_data_t *arg,
+                   gal_data_t *envseed_in)
+{
+  gsl_rng *rng;
+  uint8_t *envseed;
+  const char *rng_name;
+  double *d, *df, arg_v;
+  unsigned long rng_seed;
+  gal_data_t *out, *targ;
+
+
+  /* Sanity checks. */
+  if(arg->size!=1)
+    error(EXIT_FAILURE, 0, "the first popped operand to the '%s' "
+          "operator should be a single number (specifying the fixed "
+          "sigma, or background level for Poisson noise), but it "
+          "has %zu elements, in %zu dimension(s)",
+          gal_arithmetic_operator_string(operator), arg->size,
+          arg->ndim);
+  if(envseed_in->size!=1)
+    error(EXIT_FAILURE, 0, "the third popped operand to the '%s' "
+          "operator should be a single number (with a value of 0 "
+          "or 1, specifying if the environment should be used for "
+          "the random number generator settings), but it has %zu "
+          "elements, in %zu dimension(s)",
+          gal_arithmetic_operator_string(operator), arg->size,
+          arg->ndim);
+  if(envseed_in->type!=GAL_TYPE_UINT8)
+    error(EXIT_FAILURE, 0, "the third popped operand to the '%s' "
+          "operator should have a type of 'uint8' (with "
+          "a value of 0 or 1, specifying if the environment should "
+          "be used for the random number generator settings), but "
+          "it has a type of '%s'",
+          gal_arithmetic_operator_string(operator),
+          gal_type_name(envseed_in->type, 1));
+  envseed=envseed_in->array;
+  if( *envseed>1 )
+    error(EXIT_FAILURE, 0, "the third popped operand to the '%s' "
+          "operator should only have a value of 0 or 1 (specifying "
+          "if the environment should be used for the random number "
+          "generator settings), but it has a value of %u",
+          gal_arithmetic_operator_string(operator), *envseed);
+
+
+  /* Convert the input and argument into 'double' (and immediately free it
+     if it is no longer necessary). */
+  if(in->type==GAL_TYPE_FLOAT64) out=in;
+  else
+    {
+      out=gal_data_copy_to_new_type(in, GAL_TYPE_FLOAT64);
+      if(flags & GAL_ARITHMETIC_FREE)
+        { gal_data_free(in); in=NULL; }
+    }
+  targ=gal_data_copy_to_new_type(arg, GAL_TYPE_FLOAT64);
+  arg_v=((double *)(targ->array))[0];
+  gal_data_free(targ);
+
+
+  /* Make sure the noise identifier is positive. */
+  if(arg_v<0)
+    error(EXIT_FAILURE, 0, "the noise identifier (sigma for "
+          "'mknoise-sigma' or background for 'mknoise-poisson') must "
+          "be positive (it is %g)", arg_v);
+
+
+  /* Setup the random number generator. */
+  rng=gal_checkset_gsl_rng(*envseed, &rng_name, &rng_seed);
+
+
+  /* Add the noise. */
+  df=(d=out->array)+out->size;
+  switch(operator)
+    {
+    case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
+      do *d += gsl_ran_gaussian(rng, arg_v); while(++d<df);
+      break;
+    case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
+      do
+        *d += arg_v + gsl_ran_gaussian(rng, sqrt( arg_v + *d ));
+      while(++d<df);
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+            "fix the problem. The operator code %d isn't recognized "
+            "in this function", __func__, PACKAGE_BUGREPORT, operator);
+    }
+
+  /* Clean up and return */
+  if(flags & GAL_ARITHMETIC_FREE)
+    {
+      gal_data_free(arg);
+      gal_data_free(envseed_in);
+    }
+  return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
 /***************                  Metadata                **************/
 /***********************************************************************/
 
@@ -2009,6 +2132,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; }
 
+  /* Adding noise operators. */
+  else if (!strcmp(string, "mknoise-sigma"))
+    { op=GAL_ARITHMETIC_OP_MKNOISE_SIGMA;     *num_operands=2; }
+  else if (!strcmp(string, "mknoise-poisson"))
+    { op=GAL_ARITHMETIC_OP_MKNOISE_POISSON;   *num_operands=2; }
+
   /* The size operator */
   else if (!strcmp(string, "size"))
     { op=GAL_ARITHMETIC_OP_SIZE;              *num_operands=2;  }
@@ -2165,6 +2294,9 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:    return "sigclip-mean";
     case GAL_ARITHMETIC_OP_SIGCLIP_STD:     return "sigclip-number";
 
+    case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:   return "mknoise-sigma";
+    case GAL_ARITHMETIC_OP_MKNOISE_POISSON: return "mknoise-poisson";
+
     case GAL_ARITHMETIC_OP_SIZE:            return "size";
 
     case GAL_ARITHMETIC_OP_TO_UINT8:        return "uchar";
@@ -2325,6 +2457,15 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       out=arithmetic_bitwise_not(flags, d1);
       break;
 
+    /* Adding noise. */
+    case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
+    case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      d3 = va_arg(va, gal_data_t *);
+      out=arithmetic_mknoise(operator, flags, d1, d2, d3);
+      break;
+
     /* Size operator */
     case GAL_ARITHMETIC_OP_SIZE:
       d1 = va_arg(va, gal_data_t *);
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 49d6930..ebe4e55 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -148,6 +148,9 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN,/* Sigma-clipped median of mult. arrays.*/
   GAL_ARITHMETIC_OP_SIGCLIP_STD,  /* Sigma-clipped STD of multiple arrays. */
 
+  GAL_ARITHMETIC_OP_MKNOISE_SIGMA,/* Fixed-sigma noise to every element.   */
+  GAL_ARITHMETIC_OP_MKNOISE_POISSON,/* Poission noise on every element.    */
+
   GAL_ARITHMETIC_OP_SIZE,         /* Size of the dataset along an axis     */
 
   GAL_ARITHMETIC_OP_TO_UINT8,     /* Convert to uint8_t.                   */