[gnuastro-commits] master 8df52c6 2/2: Operator to select blank pixels in Arithmetic

[Mohammad Akhlaghi]

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

    Operator to select blank pixels in Arithmetic
    
    With the equal operator in Arithmetic, we can choose certain pixels in the
    image except for blank pixels (NaN in float types and fixed constants for
    other types). But choosing blank pixels is necessary in many situations.
    
    The problem with blank values is that in floating points, the IEEE NaN
    standard defines this value as failing any conditional operator, so NaN is
    not equal to its self. The only reliable way we can check for it is the
    `isnan' function. For integer types the blank value is a constant (defined
    by the type in Gnuastro's `lib/gnuastro/fits.h'). But the user can't use
    that value manually because like all utilities, thanks to CFITSIO,
    Arithmetic will detect their blank-ness and since all internal operations
    are in double precision floating point type, they will be seen as NaN
    values. So even if the user reads the blank value for integer types from
    the FITS header, it is useless for them in Arithmetic.
    
    To address this problem, an `isblank' operator is now added to
    Arithmetic. Using C's `isnan' function, it will return 1 if the pixel is
    NaN and 0 if it isn't. Because of the internal conversion, the type of
    input (or value of blank pixels) is also irrelevant within Arithmetic.
    
    This finishes task #14146.
---
 doc/gnuastro.texi           |   18 ++++++++++++++----
 src/arithmetic/args.h       |    2 +-
 src/arithmetic/arithmetic.c |   35 +++++++++++++++++++++++++++++++++++
 3 files changed, 50 insertions(+), 5 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index b7794bb..db9644f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -6988,6 +6988,16 @@ Non-Equality: similar to @code{lt} (`less than' 
operator), but returning 1
 when the two popped operands are @emph{not} equal (to double precison
 floating point accuracy).
 
address@hidden Blank pixel
address@hidden isblank
+Test for a blank value (see @ref{Blank pixels}). In essence, this is very
+similar to the conditional operators: the output is either 1 or 0 (see the
+`less than' operator above). The difference is that it only needs one
+operand. Because of the definition of a blank pixel, a blank value is not
+even equal to itself, so you cannot use the equal operator above to select
+blank pixels. See the ``Blank pixels'' box below for more on Blank pixels
+in Arithmetic.
+
 @item where
 Change the pixel values `where' a certain condition holds. The conditional
 operators above can be used to define the condition. This operator requires
@@ -7028,9 +7038,9 @@ affected, they are only read.
 done in double precision floating points. So all blank pixels in the image
 (see @ref{Blank pixels}) will be stored as IEEE NaN values. One particular
 aspect of NaN values is that by definition they will fail on @emph{any}
-comparison, so a NaN is not equal to NaN! Currently, Gnuastro has no way to
-select only blank pixels but we have defined task 14146 to do it when we
-get the chance or someone offers to help.
+comparison. Hence both equal and not-equal operators will fail when both
+their operands are NaN! The only way to select blank pixels is through the
address@hidden operator explained above.
 
 One way you can exploit this property of the NaN value to your advantage is
 when you want a fully zero-valued image (even over the blank pixels) based
@@ -9839,7 +9849,7 @@ $ astsubtractsky --nch1=2 --nch2=2 image.fits
 $ astsubtractsky --mask=maskforimage.fits image.fits
 @end example
 
address@hidden Blank pixels
address@hidden Blank pixel
 The only required input to SubtractSky is the input data file that
 currently has to be only a 2D image. But in the future it might be
 useful to use it for 1D or 3D data too. Any pixels in the image with a
diff --git a/src/arithmetic/args.h b/src/arithmetic/args.h
index fcc8975..f46ce23 100644
--- a/src/arithmetic/args.h
+++ b/src/arithmetic/args.h
@@ -68,7 +68,7 @@ const char doc[] =
   "Please see the manual for more information. "
   "\n\nThe operators/functions recognized by "SPACK_NAME" are: +, -, *, /, "
   "abs, pow, sqrt, log, log10, minvalue, maxvalue, min, max, average, median, "
-  "lt, le, gt, ge, eq, neq. Please run `info gnuastro \"Arithmetic "
+  "lt, le, gt, ge, eq, neq, isblank. Please run `info gnuastro \"Arithmetic "
   "operators\"' for detailed information on each operator. Note that "
   "multiplication should be quoted (like \"*\", or '*') to avoid shell "
   "expansion.\n"
diff --git a/src/arithmetic/arithmetic.c b/src/arithmetic/arithmetic.c
index a99f7db..3909406 100644
--- a/src/arithmetic/arithmetic.c
+++ b/src/arithmetic/arithmetic.c
@@ -838,6 +838,40 @@ conditionals(struct imgarithparams *p, char *operator)
 
 
 
+/* In order to not conflict with the internal C `is...' functions, and in
+   particular the `isblank' function, we are calling this function
+   opisblank for operator-isblank. */
+void
+opisblank(struct imgarithparams *p)
+{
+  size_t size=p->s0*p->s1;
+  char *operator="isblank";
+  double *f, *ff, fnum, *farr;
+
+  /* Pop out the number of operands needed. */
+  pop_operand(p, &fnum, &farr, operator);
+
+  /* Do the operation: */
+  if(farr)                       /* Operand is array.        */
+    {
+      /* Do the operation, note that the output is stored in the first
+         input. Also note that since the linked list is
+         first-in-first-out, the second operand should be put first
+         here. */
+      ff=(f=farr)+size;
+      do *f = isnan(*f); while(++f<ff);
+
+      /* Push the output onto the stack. */
+      add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
+    }
+  else                          /* Operand is a number.      */
+    add_operand(p, NOOPTFILENAME, isnan(fnum), NOOPTARRAY);
+}
+
+
+
+
+
 /* Replace the pixels in the second popped element with the first. While
    choosing the pixels that are selected from the third The third popped
    element. The third is considered to be an array that can only be filled
@@ -964,6 +998,7 @@ reversepolish(struct imgarithparams *p)
                   || !strcmp(token->v, "ge")
                   || !strcmp(token->v, "eq")
                   || !strcmp(token->v, "neq")) conditionals(p, token->v);
+          else if(!strcmp(token->v, "isblank")) opisblank(p);
           else if(!strcmp(token->v, "where")) where(p);
           else
             error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "