[gnuastro-commits] master b42b89b: Final noise sigma as input in MakeNoise

[Mohammad Akhlaghi]

branch: master
commit b42b89b5630989b35794a7b4c72ec67fe60d8eec
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Final noise sigma as input in MakeNoise
    
    Until now, it was necessary to calculate the background magnitude for the
    noise-level you wanted. This is indeed the most accurate way to estimate
    the noise (since noise depends on the background and the pixel
    values). However, when users are just interested in the low surface
    brightness regions and (where the noise levels don't change much), it is
    much more convenient to just input the noise sigma instead of going through
    all the trouble of converting it to a background magnitude with a certain
    zeropoint.
    
    Also with this commit, the old `--stdadd' (`-s') option is now called
    `--instrumental' (`-i') to be more clear.
---
 NEWS                        |  8 ++++++
 bin/mknoise/args.h          | 25 +++++++++++++-----
 bin/mknoise/astmknoise.conf |  2 +-
 bin/mknoise/main.h          |  3 ++-
 bin/mknoise/mknoise.c       | 63 ++++++++++++++++++++++++++++++---------------
 bin/mknoise/ui.c            | 40 ++++++++++++++++++++--------
 bin/mknoise/ui.h            | 11 ++++----
 doc/gnuastro.texi           | 39 +++++++++++++++++-----------
 8 files changed, 131 insertions(+), 60 deletions(-)

diff --git a/NEWS b/NEWS
index 8d23e44..76f1563 100644
--- a/NEWS
+++ b/NEWS
@@ -10,6 +10,11 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   Arithmetic: now has a new `--globalhdu' (`-g') option which can be used
   once for all the input images.
 
+  MakeNoise: with the new `--sigma' (`-s') option, it is now possible to
+  directly request the noise sigma or standard deviation. When this option
+  is called, the `--background', `--zeropoint' and other option values will
+  be ignored.
+
   MakeProfiles: the new `--kernel' option can make a kernel image without
   the need to define a catalog. With this option, a catalog (or
   accompanying background image) must not be given.
@@ -51,6 +56,9 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   multiple times and the order of its calling will be used for the column
   containing the center in the respective dimension (in FITS format).
 
+  MakeNoise: the old `--stdadd' (`-s') option has been renamed to
+  `--instrumental' (`-i') to be more clear.
+
   MakeProfiles: The new `--naxis' and `--shift' options can take multiple
   values for each dimension (separated by a comma). This replaces the old
   `--naxis1', `--naxis2' and `--xshift' and `--yshift' options.
diff --git a/bin/mknoise/args.h b/bin/mknoise/args.h
index ba19980..797952e 100644
--- a/bin/mknoise/args.h
+++ b/bin/mknoise/args.h
@@ -32,13 +32,26 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 struct argp_option program_options[] =
   {
     {
-      "stdadd",
-      UI_KEY_STDADD,
+      "sigma",
+      UI_KEY_SIGMA,
       "FLT",
       0,
-      "Standard deviation addition constant.",
+      "Total noise sigma, ignore other options.",
       GAL_OPTIONS_GROUP_INPUT,
-      &p->stdadd,
+      &p->sigma,
+      GAL_TYPE_FLOAT64,
+      GAL_OPTIONS_RANGE_GE_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "instrumental",
+      UI_KEY_INSTRUMENTAL,
+      "FLT",
+      0,
+      "Instrument noise level (in pixel value units).",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->instrumental,
       GAL_TYPE_FLOAT64,
       GAL_OPTIONS_RANGE_GE_0,
       GAL_OPTIONS_NOT_MANDATORY,
@@ -54,7 +67,7 @@ struct argp_option program_options[] =
       &p->background_mag,
       GAL_TYPE_FLOAT64,
       GAL_OPTIONS_RANGE_ANY,
-      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
     {
@@ -67,7 +80,7 @@ struct argp_option program_options[] =
       &p->zeropoint,
       GAL_TYPE_FLOAT64,
       GAL_OPTIONS_RANGE_GE_0,
-      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
 
diff --git a/bin/mknoise/astmknoise.conf b/bin/mknoise/astmknoise.conf
index a29f1b3..96e011c 100644
--- a/bin/mknoise/astmknoise.conf
+++ b/bin/mknoise/astmknoise.conf
@@ -19,7 +19,7 @@
 
 # Input:
  background     -10.00
- stdadd         0.000
+ instrumental   0.000
  zeropoint      0.00
 
 # Output:
diff --git a/bin/mknoise/main.h b/bin/mknoise/main.h
index d143cdd..fd4c386 100644
--- a/bin/mknoise/main.h
+++ b/bin/mknoise/main.h
@@ -44,7 +44,8 @@ struct mknoiseparams
   /* From command-line */
   struct gal_options_common_params cp;   /* Common parameters.           */
   char        *inputname;    /* Input filename.                          */
-  double          stdadd;    /* Standard deviation constants.            */
+  double           sigma;    /* Total noise sigma (ignoring others).     */
+  double    instrumental;    /* Standard deviation constants.            */
   double       zeropoint;    /* Zeropoint magnitude of image.            */
   double  background_mag;    /* Background in magnitudes.                */
   uint8_t        envseed;    /* ==1, generate a random seed.             */
diff --git a/bin/mknoise/mknoise.c b/bin/mknoise/mknoise.c
index e17b8ad..f91c265 100644
--- a/bin/mknoise/mknoise.c
+++ b/bin/mknoise/mknoise.c
@@ -62,32 +62,42 @@ convertsaveoutput(struct mknoiseparams *p)
 
   /* Add the proper information to the header of the output: */
   gal_fits_key_write_filename("INF", p->inputname, &headers);
-  strcpy(keyname1, "BCKGRND");
-  gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname1, 0,
-                            &p->background_mag, 0, "Background "
-                            "value (in magnitude) for noise.",
-                            0, NULL);
-  strcpy(keyname2, "BZRPNT");
-  gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname2, 0,
-                            &p->zeropoint, 0, "Zeropoint magnitude of image.",
-                            0, NULL);
-  strcpy(keyname3, "STDADD");
-  gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname3, 0,
-                            &p->stdadd, 0, "Instrumental noise in units of "
-                            "flux.", 0, NULL);
+  if( !isnan(p->background_mag) )
+    {
+      strcpy(keyname1, "BCKGRND");
+      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname1, 0,
+                                &p->background_mag, 0, "Background "
+                                "value (in magnitude) for noise.",
+                                0, NULL);
+      strcpy(keyname2, "BZRPNT");
+      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname2, 0,
+                                &p->zeropoint, 0,
+                                "Zeropoint magnitude of image.", 0, NULL);
+      strcpy(keyname3, "INSTRU");
+      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname3, 0,
+                                &p->instrumental, 0,
+                                "Instrumental noise in units of flux.",
+                                0, NULL);
+    }
+  else
+    {
+      strcpy(keyname1, "SIGMA");
+      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname1, 0,
+                                &p->sigma, 0, "Total noise sigma", 0, NULL);
+    }
   strcpy(keyname4, "RNGTYPE");
   gal_fits_key_list_add_end(&headers, GAL_TYPE_STRING, keyname4, 0,
                             p->rng_type, 0, "Random number generator (by "
                             "GSL) type.",  0, NULL);
   strcpy(keyname5, "RNGSEED");
-  gal_fits_key_list_add_end(&headers, GAL_TYPE_INT64, keyname5, 0, 
&p->rng_seed,
-                             0, "Random number generator (by GSL) seed.",
-                             0, NULL);
+  gal_fits_key_list_add_end(&headers, GAL_TYPE_INT64, keyname5, 0,
+                            &p->rng_seed, 0,
+                            "Random number generator (by GSL) seed.",
+                            0, NULL);
 
   /* Save the output: */
   p->input=gal_data_copy_to_new_type_free(p->input, p->cp.type);
   gal_fits_img_write(p->input, p->cp.output, headers, PROGRAM_STRING);
-
 }
 
 
@@ -97,13 +107,24 @@ convertsaveoutput(struct mknoiseparams *p)
 void
 mknoise(struct mknoiseparams *p)
 {
-  double *d, *df, background=p->background, stdadd=p->stdadd;
+  double *d, *df, background=p->background, instrumental=p->instrumental;
 
   /* Add the noise: */
   df=(d=p->input->array)+p->input->size;
-  do
-    *d += background+gsl_ran_gaussian(p->rng, sqrt(stdadd+background+*d));
-  while(++d<df);
+  if( !isnan(p->sigma) )
+    {
+      do
+        *d += gsl_ran_gaussian(p->rng, p->sigma);
+      while(++d<df);
+    }
+  else
+    {
+      do
+        *d += ( background
+                + gsl_ran_gaussian(p->rng,
+                                   sqrt( instrumental + background + *d )) );
+      while(++d<df);
+    }
 
   /* Convert and save the output in the proper format: */
   convertsaveoutput(p);
diff --git a/bin/mknoise/ui.c b/bin/mknoise/ui.c
index 14db897..79ab5a0 100644
--- a/bin/mknoise/ui.c
+++ b/bin/mknoise/ui.c
@@ -111,6 +111,12 @@ ui_initialize_options(struct mknoiseparams *p,
   cp->coptions           = gal_commonopts_options;
 
 
+  /* Initialize options for this program. */
+  p->sigma               = NAN;
+  p->zeropoint           = NAN;
+  p->background_mag      = NAN;
+
+
   /* Modify common options. */
   for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
     {
@@ -208,13 +214,25 @@ parse_opt(int key, char *arg, struct argp_state *state)
 /***************       Sanity Check         *******************/
 /**************************************************************/
 /* Read and check ONLY the options. When arguments are involved, do the
-   check in `ui_check_options_and_arguments'.
+   check in `ui_check_options_and_arguments'. */
 static void
 ui_read_check_only_options(struct mknoiseparams *p)
 {
-
+  /* At leaset one of `--sigma' or `--background' are necessary. */
+  if( isnan(p->sigma) && isnan(p->background_mag) )
+    error(EXIT_FAILURE, 0, "at least one of `--sigma' or `--background' "
+          "must be given to identify the noise level");
+
+
+  /* If a background magnitude is given (and the user hasn't given a
+     `--sigma'), the zeropoint is necessary. */
+  if( isnan(p->sigma) && !isnan(p->background_mag) && isnan(p->zeropoint) )
+    error(EXIT_FAILURE, 0, "no zeropoint magnitude given. When the noise is "
+          "identified by the background magnitude, a zeropoint magnitude "
+          "is mandatory. Please use the `--zeropoint' option to specify "
+          "a zeropoint magnitude");
 }
-*/
+
 
 
 
@@ -286,12 +304,12 @@ ui_preparations(struct mknoiseparams *p)
 
 
   /* Convert the background value from magnitudes to flux. Note that
-     magnitudes are actually calculated from the ratio of brightness,
-     not flux. But in the context of MakeNoise where everything is
-     done on pixels independently, brightness and flux are the same
-     (flux is multiplied by the area of one pixel (=1) to give
-     brightness).*/
-  p->background=pow(10, (p->zeropoint-p->background_mag)/2.5f);
+     magnitudes are actually calculated from the ratio of brightness, not
+     flux. But in the context of MakeNoise where everything is done on
+     pixels independently, brightness and flux are the same (flux is
+     multiplied by the area of one pixel (=1) to give brightness).*/
+  if( !isnan(p->background_mag) )
+    p->background=pow(10, (p->zeropoint-p->background_mag)/2.5f);
 
 
   /* Allocate the random number generator: */
@@ -359,9 +377,9 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
mknoiseparams *p)
 
 
   /* Read the options into the program's structure, and check them and
-     their relations prior to printing.
+     their relations prior to printing. */
   ui_read_check_only_options(p);
-  */
+
 
   /* Print the option values if asked. Note that this needs to be done
      after the option checks so un-sane values are not printed in the
diff --git a/bin/mknoise/ui.h b/bin/mknoise/ui.h
index 137b237..38c1d1f 100644
--- a/bin/mknoise/ui.h
+++ b/bin/mknoise/ui.h
@@ -29,16 +29,17 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 /* Available letters for short options:
 
-   a c d f g i j k l m n p r t u v w x y
+   a c d f g j k l m n p r t u v w x y
    A B C E G H J L O Q R W X Y
 */
 enum option_keys_enum
 {
   /* With short-option version. */
-  UI_KEY_STDADD      = 's',
-  UI_KEY_BACKGROUND  = 'b',
-  UI_KEY_ZEROPOINT   = 'z',
-  UI_KEY_ENVSEED     = 'e',
+  UI_KEY_SIGMA        = 's',
+  UI_KEY_INSTRUMENTAL = 'i',
+  UI_KEY_BACKGROUND   = 'b',
+  UI_KEY_ZEROPOINT    = 'z',
+  UI_KEY_ENVSEED      = 'e',
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index bfa68bc..e30d7a0 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15650,26 +15650,35 @@ $ astmknoise [OPTION ...] InputImage.fits
 One line examples:
 
 @example
-## Add noise to input image assuming background and instrumental noise:
-$ astmknoise --background=1000 --stdadd=20 mockimage.fits
+## Add noise with a standard deviation of 100 to image:
+$ astmknoise --sigma=100 image.fits
+
+## Add noise to input image assuming a background magnitude (with zeropoint
+## magnitude of 0) and a certain instrumental noise:
+$ astmknoise --background=-10 -z0 --instrumental=20 mockimage.fits
 @end example
 
 @noindent
 If actual processing is to be done, the input image is a mandatory
-argument. The full list of options common to all the programs in
-Gnuastro can be seen in @ref{Common options}. The output will have the
-same type as the input image, however the internal processing is done
-on a double precision floating point format. If the input values were
-integer types, then each floating point number will be rounded to the
-nearest integer away from zero. This might cause integer overflow if
-types with small ranges are used (for example images with a
address@hidden of @code{8} which can only keep 256 values). This can be
-disabled with the @option{doubletype} option.  The header of the
-output FITS file keeps all the parameters that were influential in
-making it. This is done for future reproducibility.
+argument. The full list of options common to all the programs in Gnuastro
+can be seen in @ref{Common options}. The type (see @ref{Numeric data
+types}) of the output can be specified with the @option{--type} option, see
address@hidden output options}. The header of the output FITS file keeps all
+the parameters that were influential in making it. This is done for future
+reproducibility.
 
 @table @option
 
address@hidden -s FLT
address@hidden --sigma=FLT
+The total noise sigma in the same units as the pixel values. With this
+option, the @option{--background}, @option{--zeropoint} and
address@hidden will be ignored. With this option, the noise will
+be independent of the pixel values (which is not realistic, see @ref{Photon
+counting noise}). Hence it is only useful if you are working on low surface
+brightness regions where the change in pixel value (and thus real noise) is
+insignificant.
+
 @item -b FLT
 @itemx --background=FLT
 The background pixel value for the image in units of magnitudes, see
@@ -15680,8 +15689,8 @@ The background pixel value for the image in units of 
magnitudes, see
 The zeropoint magnitude used to convert the value of @option{--background}
 (in units of magnitude) to flux, see @ref{Flux Brightness and magnitude}.
 
address@hidden -s FLT
address@hidden --stdadd=FLT
address@hidden -i FLT
address@hidden --instrumental=FLT
 The instrumental noise which is in units of flux, see @ref{Instrumental
 noise}.