[gnuastro-commits] master a2cce53: MakeNoise: error printed with background is less than one

[Mohammad Akhlaghi]

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

    MakeNoise: error printed with background is less than one
    
    Until now there was no sanity check on the value to the given
    background. It may happen that the user mistakenly gives a value less than
    1. But since hte background is used in a Poisson distribution (which is a
    counting distribution), it must be larger than one.
    
    With this commit, a sanity check has been added to check for this. Also
    some cleaning up was done in the code and the default values to the
    '--background', '--zeropoint' and '--instrument' options has been removed:
    the user should specify the type of the background/noise on every
    run. Should it be Poisson? Is it in magnitudes?
---
 bin/mknoise/args.h           |  2 +-
 bin/mknoise/astmknoise.conf  |  3 --
 bin/mknoise/main.h           |  3 +-
 bin/mknoise/mknoise.c        | 70 ++++++++++++++++++++++++++------------------
 bin/mknoise/ui.c             | 56 ++++++++++++++++++++---------------
 tests/mknoise/addnoise-3d.sh |  3 +-
 tests/mknoise/addnoise.sh    |  5 ++--
 7 files changed, 81 insertions(+), 61 deletions(-)

diff --git a/bin/mknoise/args.h b/bin/mknoise/args.h
index eda7581..22c6ad6 100644
--- a/bin/mknoise/args.h
+++ b/bin/mknoise/args.h
@@ -64,7 +64,7 @@ struct argp_option program_options[] =
       0,
       "Fixed background magnitude for whole input.",
       GAL_OPTIONS_GROUP_INPUT,
-      &p->background_mag,
+      &p->background,
       GAL_TYPE_FLOAT64,
       GAL_OPTIONS_RANGE_ANY,
       GAL_OPTIONS_NOT_MANDATORY,
diff --git a/bin/mknoise/astmknoise.conf b/bin/mknoise/astmknoise.conf
index b029737..8295fd5 100644
--- a/bin/mknoise/astmknoise.conf
+++ b/bin/mknoise/astmknoise.conf
@@ -20,9 +20,6 @@
 # warranty.
 
 # Input:
- background     -10.00
- instrumental   0.000
- zeropoint      0.00
 
 # Output:
  type           float32
diff --git a/bin/mknoise/main.h b/bin/mknoise/main.h
index 3f357f9..024c662 100644
--- a/bin/mknoise/main.h
+++ b/bin/mknoise/main.h
@@ -47,13 +47,12 @@ struct mknoiseparams
   double           sigma;    /* Total noise sigma (ignoring others).     */
   double    instrumental;    /* Standard deviation constants.            */
   double       zeropoint;    /* Zeropoint magnitude of image.            */
-  double  background_mag;    /* Background in magnitudes.                */
+  double      background;    /* Background in magnitudes.                */
   uint8_t bgisbrightness;    /* Background is brightness, not magnitude. */
   uint8_t        envseed;    /* ==1, generate a random seed.             */
 
   /* Internal */
   gal_data_t      *input;    /* Input image data in double precision.    */
-  double      background;    /* Background in units of brightness.       */
   gsl_rng           *rng;    /* Main instance of random number generator.*/
   const char   *rng_name;    /* The type/name of the Random number gen.  */
   unsigned long rng_seed;    /* Seed of Random number generator.         */
diff --git a/bin/mknoise/mknoise.c b/bin/mknoise/mknoise.c
index 3ccbba2..7a16c9e 100644
--- a/bin/mknoise/mknoise.c
+++ b/bin/mknoise/mknoise.c
@@ -35,6 +35,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gsl/gsl_randist.h>     /* To make noise.        */
 
 #include <gnuastro-internal/timing.h>
+#include <gnuastro-internal/checkset.h>
 
 #include "main.h"
 
@@ -54,45 +55,56 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 void
 convertsaveoutput(struct mknoiseparams *p)
 {
-  char keyname1[FLEN_KEYWORD];
+  double tmp;
+  char *keyname;
   gal_fits_list_key_t *headers=NULL;
-  char keyname2[FLEN_KEYWORD], keyname3[FLEN_KEYWORD];
-  char keyname4[FLEN_KEYWORD], keyname5[FLEN_KEYWORD];
-
 
   /* Add the proper information to the header of the output: */
   gal_fits_key_write_filename("INF", p->inputname, &headers, 0);
-  if( !isnan(p->background_mag) )
+  if( !isnan(p->background) )
     {
-      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, 0);
-      strcpy(keyname2, "BZRPNT");
-      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname2, 0,
-                                &p->zeropoint, 0,
-                                "Zeropoint magnitude of image.", 0, NULL, 0);
-      strcpy(keyname3, "INSTRU");
-      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname3, 0,
-                                &p->instrumental, 0,
-                                "Instrumental noise in units of flux.",
+      gal_checkset_allocate_copy("BCKGRND", &keyname);
+      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname, 1,
+                                &p->background, 0,
+                                "Background value for Poisson noise.",
                                 0, NULL, 0);
+      if( !isnan(p->zeropoint) )
+        {
+          tmp=-2.5 * log10(p->background) + p->zeropoint;
+          gal_checkset_allocate_copy("BCKGMAG", &keyname);
+          gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname, 1,
+                                    &tmp, 0,
+                                    "Background value in magnitudes",
+                                    0, NULL, 0);
+          gal_checkset_allocate_copy("BCKGZP", &keyname);
+          gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname, 1,
+                                    &p->zeropoint, 0,
+                                    "Zeropoint for interpreting magnitudes.",
+                                    0, NULL, 0);
+        }
+      if( !isnan(p->instrumental) )
+        {
+          gal_checkset_allocate_copy("INSTRU", &keyname);
+          gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname, 1,
+                                    &p->instrumental, 0,
+                                    "Instrumental noise in units of flux.",
+                                    0, NULL, 0);
+        }
     }
   else
     {
-      strcpy(keyname1, "SIGMA");
-      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname1, 0,
+      gal_checkset_allocate_copy("SIGMA", &keyname);
+      gal_fits_key_list_add_end(&headers, GAL_TYPE_FLOAT64, keyname, 1,
                                 &p->sigma, 0, "Total noise sigma", 0,
                                 NULL, 0);
     }
-  strcpy(keyname4, "RNGTYPE");
-  gal_fits_key_list_add_end(&headers, GAL_TYPE_STRING, keyname4, 0,
+  gal_checkset_allocate_copy("RNGTYPE", &keyname);
+  gal_fits_key_list_add_end(&headers, GAL_TYPE_STRING, keyname, 1,
                             (void *)(p->rng_name), 0,
                             "Random number generator (by GSL) type.",
                             0, NULL, 0);
-  strcpy(keyname5, "RNGSEED");
-  gal_fits_key_list_add_end(&headers, GAL_TYPE_ULONG, keyname5, 0,
+  gal_checkset_allocate_copy("RNGSEED", &keyname);
+  gal_fits_key_list_add_end(&headers, GAL_TYPE_ULONG, keyname, 1,
                             &p->rng_seed, 0,
                             "Random number generator (by GSL) seed.",
                             0, NULL, 0);
@@ -116,8 +128,10 @@ convertsaveoutput(struct mknoiseparams *p)
 void
 mknoise(struct mknoiseparams *p)
 {
-  double *d, *df, background=p->background;
-  double instpowtwo = p->instrumental*p->instrumental;
+  double *d, *df, back=p->background;
+  double inst = ( isnan(p->instrumental)
+                  ? 0.0f
+                  : p->instrumental*p->instrumental );
 
   /* Add the noise: */
   df=(d=p->input->array)+p->input->size;
@@ -130,9 +144,7 @@ mknoise(struct mknoiseparams *p)
   else
     {
       do
-        *d += ( background
-                + gsl_ran_gaussian(p->rng,
-                                   sqrt( instpowtwo + background + *d )) );
+        *d += back + gsl_ran_gaussian(p->rng, sqrt( inst + back + *d ));
       while(++d<df);
     }
 
diff --git a/bin/mknoise/ui.c b/bin/mknoise/ui.c
index c299288..428d5d9 100644
--- a/bin/mknoise/ui.c
+++ b/bin/mknoise/ui.c
@@ -116,7 +116,8 @@ ui_initialize_options(struct mknoiseparams *p,
   /* Initialize options for this program. */
   p->sigma               = NAN;
   p->zeropoint           = NAN;
-  p->background_mag      = NAN;
+  p->background          = NAN;
+  p->instrumental        = NAN;
 
 
   /* Modify common options. */
@@ -222,20 +223,41 @@ 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( isnan(p->sigma) && isnan(p->background) )
+    error(EXIT_FAILURE, 0, "noise not defined: please define the noise "
+          "level with either '--sigma' (for a fixed noise STD for all "
+          "the pixels, irrespective of their value) or '--background' "
+          "(to use in a Poisson noise model, where the noise will differ "
+          "based on pixel value)");
 
 
   /* If a background magnitude is given ('--bgbrightness' hasn't been
      called), the zeropoint is necessary. */
-  if( !isnan(p->background_mag) && p->bgisbrightness==0 && 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. Alternatively, if your background value "
-          "is in units of brightness (which doesn't need a zeropoint), "
-          "please use '--bgisbrightness'");
+  if( !isnan(p->background) )
+    {
+      /* Make sure that the background can be interpretted properly. */
+      if( p->bgisbrightness==0 && isnan(p->zeropoint) )
+        error(EXIT_FAILURE, 0, "missing background information. When the "
+              "noise is identified by the background, a zeropoint magnitude "
+              "is mandatory. Please use the '--zeropoint' option to specify "
+              "a zeropoint magnitude. Alternatively, if your background value "
+              "is brightness (which is in linear scale and doesn't need a "
+              "zeropoint), please use '--bgisbrightness'");
+
+      /* If the background is in units of magnitudes, convert it to
+         brightness. */
+      if( p->bgisbrightness==0 )
+        p->background = pow(10, (p->zeropoint-p->background)/2.5f);
+
+      /* Make sure that the background is larger than 1 (where Poisson
+         noise is actually defined). */
+      if( p->background < 1 )
+        error(EXIT_FAILURE, 0, "backgroud value is smaller than 1. "
+              "Poisson noise is only defined on a positive distribution "
+              "with values larger than 1. You can use the '--sigma' "
+              "option to add a fixed noise level (with any positive value) "
+              "to any pixel.");
+    }
 }
 
 
@@ -307,18 +329,6 @@ ui_preparations(struct mknoiseparams *p)
     p->cp.output=gal_checkset_automatic_output(&p->cp, p->inputname,
                                                "_noised.fits");
 
-
-  /* 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).*/
-  if( !isnan(p->background_mag) )
-    p->background = ( p->bgisbrightness
-                      ? p->background_mag
-                      : pow(10, (p->zeropoint-p->background_mag)/2.5f) );
-
-
   /* Allocate the random number generator: */
   p->rng=gal_checkset_gsl_rng(p->envseed, &p->rng_name, &p->rng_seed);
 }
diff --git a/tests/mknoise/addnoise-3d.sh b/tests/mknoise/addnoise-3d.sh
index ec010f4..2a48601 100755
--- a/tests/mknoise/addnoise-3d.sh
+++ b/tests/mknoise/addnoise-3d.sh
@@ -56,4 +56,5 @@ if [ ! -f $img      ]; then echo "$img does not exist.";     
exit 77; fi
 # ==================
 export GSL_RNG_SEED=1
 export GSL_RNG_TYPE=ranlxs2
-$execname --envseed $img
+options="--background=-10 --zeropoint=0 --envseed"
+$check_with_program $execname $img $options
diff --git a/tests/mknoise/addnoise.sh b/tests/mknoise/addnoise.sh
index 969948d..0e86273 100755
--- a/tests/mknoise/addnoise.sh
+++ b/tests/mknoise/addnoise.sh
@@ -63,5 +63,6 @@ if [ ! -f $img2     ]; then echo "$img2 does not exist.";    
exit 77; fi
 # debugging when the developer doesn't have access to the user's system.
 export GSL_RNG_SEED=1
 export GSL_RNG_TYPE=ranlxs2
-$check_with_program $execname --envseed $img1
-$check_with_program $execname --envseed $img2
+options="--background=-10 --zeropoint=0 --envseed"
+$check_with_program $execname $img1 $options
+$check_with_program $execname $img2 $options