[gnuastro-commits] master 00d6f5f: Arithmetic: new operator to find bounding box of an ellipse

[Mohammad Akhlaghi]

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

    Arithmetic: new operator to find bounding box of an ellipse
    
    Until now, Gnuastro only had a library function to find the width and
    height of the bounding box of an ellipse (such that the box and ellipse
    share the same center point). This function is used heavily in several
    Gnuastro programs for low level operations, for example
    MakeProfiles. However, there was no easy command-line interface to it.
    
    With this commit, a wrapper for it has been written in the Arithmetic
    library, allowing easy access to the command-line. It can work on both
    table columns and images.
    
    It was the first image-operator that returns two datasets, so I had to add
    the new 'operands_add_list' function in Arithmetic's operands
    functions.
    
    During the testing, I also noticed a bug in the 'tofile-' and 'tofilefree-'
    operators of Arithmetic (they weren't deleting an existing file, just
    adding their dataset as new HDUs). So with this commit, that bug is also
    fixed.
    
    This fixes bug #60826.
---
 NEWS                        |  9 +++++
 bin/arithmetic/arithmetic.c | 10 +++---
 bin/arithmetic/operands.c   | 49 +++++++++++++++++++++++++
 doc/gnuastro.texi           | 36 ++++++++++++++++++-
 lib/arithmetic.c            | 87 +++++++++++++++++++++++++++++++++++++++++++++
 lib/gnuastro/arithmetic.h   |  2 ++
 6 files changed, 188 insertions(+), 5 deletions(-)

diff --git a/NEWS b/NEWS
index 015476c..c076db2 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,14 @@ See the end of the file for license conditions.
 
 ** New features
 
+  Arithmetic:
+   - New operands (also available in Table's column arithmetic):
+     - box-around-ellipse: width and height of the box covering an ellipse.
+
+  Library:
+   - Arithmetic macros:
+     - GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE
+
 ** Removed features
 
 ** Changed features
@@ -17,6 +25,7 @@ See the end of the file for license conditions.
               this bug was fixed by Zahra Sharbaf.
   bug #60778: Brightness error not NaN when all STD pixels are blank,
               this bug was reported by Zahra Sharbaf.
+  bug #60826: Arithmetic won't delete existing file with tofile operators.
 
 
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 7892d6b..6114062 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1019,6 +1019,10 @@ arithmetic_tofile(struct arithmeticparams *p, char 
*token, int freeflag)
                            ? OPERATOR_PREFIX_LENGTH_TOFILEFREE
                            : OPERATOR_PREFIX_LENGTH_TOFILE     ];
 
+  /* Make sure that if the file exists, it is deleted. */
+  gal_checkset_writable_remove(filename, p->cp.keep,
+                               p->cp.dontdelete);
+
   /* Save it to a file. */
   popped->wcs=p->refdata.wcs;
   if(popped->ndim==1 && p->onedasimage==0)
@@ -1031,10 +1035,8 @@ arithmetic_tofile(struct arithmeticparams *p, char 
*token, int freeflag)
 
   /* Reset the WCS to NULL and put it back on the stack. */
   popped->wcs=NULL;
-  if(freeflag)
-    gal_data_free(popped);
-  else
-    operands_add(p, NULL, popped);
+  if(freeflag) gal_data_free(popped);
+  else         operands_add(p, NULL, popped);
 }
 
 
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 7e29777..a947780 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -79,6 +79,47 @@ operands_num(struct arithmeticparams *p)
 /**********************************************************************/
 /************      Adding to and popping from stack     ***************/
 /**********************************************************************/
+/* When the operand to be added is a list, we don't have to worry about
+   filenames (external datasets) because lists can only be added */
+static void
+operands_add_list(struct arithmeticparams *p, gal_data_t *data)
+{
+  gal_data_t *tmp;
+  struct operand *newnode;
+
+  /* First, reverse the list so when we add them, they have the proper
+     order. */
+  gal_list_data_reverse(&data);
+
+  /* Go over the list and add them to the stack. */
+  while(data)
+    {
+      /* Shift the 'data' pointer to the next element. */
+      tmp=data;
+      data=data->next;
+
+      /* Allocate space for the new operand. */
+      errno=0;
+      newnode=malloc(sizeof *newnode);
+      if(newnode==NULL)
+        error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for 'newnode'",
+              __func__, sizeof *newnode);
+
+      /* Set the basic parameters. */
+      newnode->data=tmp;
+      newnode->hdu=NULL;
+      newnode->filename=NULL;
+
+      /* Add this dataset to the top of the stack. */
+      newnode->next=p->operands;
+      p->operands=newnode;
+    }
+}
+
+
+
+
+
 void
 operands_add(struct arithmeticparams *p, char *filename, gal_data_t *data)
 {
@@ -86,6 +127,14 @@ operands_add(struct arithmeticparams *p, char *filename, 
gal_data_t *data)
   size_t ndim, *dsize;
   struct operand *newnode;
 
+  /* In case 'data' is a list then there is no file name or HDU involved,
+     the given datasets were produced internally, so things are easier. */
+  if(data && data->next)
+    {
+      operands_add_list(p, data);
+      return;
+    }
+
   /* Some operators might not actually return any dataset (data=NULL), in
      such cases filename will also be NULL (since the operand was not added
      from the command-line). So, we shouldn't add anything to the stack. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 04de232..1cfa4c5 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -13072,7 +13072,34 @@ The output will be a single unsigned integer 
(dimensions cannot be negative).
 For example, the following command will produce the size of the first 
extension/HDU (the default HDU) of @file{a.fits} along the second FITS axis.
 
 @example
-astarithmetic a.fits 2 size
+$ astarithmetic a.fits 2 size
+@end example
+
+@item box-around-ellipse
+Return the width (along horizontal) and height (along vertical) of a box that 
encompases an ellipse with the same center point.
+The top-popped operand is assumed to be the position angle (angle from the 
horizontal axis) in @emph{degrees}.
+The second and third popped operands are the minor and major axis lengths 
respectively.
+This operator outputs two operands on the general stack.
+The first one is the width and the second (which will be the top one when this 
operator finishes) is the height.
+
+If the value to the second popped operand (minor axis) is larger than the 
third (major axis), a NaN value will be written for both the width and height 
of that element and a warning will be printed (the warning can be disabled with 
the @option{--quiet} option).
+
+As an example, if your ellipse has a major axis length of 10 units, a minor 
axis length of 4 units and a position angle of 20 degrees, you can estimate the 
bounding box with this command:
+
+@example
+$ echo "10 4 20" \
+       | asttable -c'arith $1 $2 $3 box-around-ellipse'
+@end example
+
+Alternatively if your three values are in separate FITS arrays/images, you can 
use the command below to have the width and height in similarly sized fits 
arrays.
+In this example @file{a.fits} and @file{b.fits} are respectively the major and 
minor axis lengths and @file{pa.fits} is the position angle (in degrees).
+Also, in all three, we assume the first extension is used.
+After its done, the height of the box will be put in @file{h.fits} and the 
width will be in @file{w.fits}.
+Just note that because this operator has two output datasets, you need to 
first write the height (top output operand) into a file and free it with the 
@code{tofilefree-} operator, then write the width in the file given to 
@option{--output}.
+
+@example
+$ astarithmetic a.fits b.fits pa.fits box-around-ellipse \
+                tofilefree-h.fits -ow.fits -g1
 @end example
 
 @item makenew
@@ -26779,6 +26806,13 @@ in @ref{Arithmetic}. So in your programs, it might be 
preferable to
 directly use those functions.
 @end deffn
 
+@deffn Macro GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE
+Return the width (along horizontal) and height (along vertical) of a box that 
encompases an ellipse with the same center point.
+For more on the three input operands to this operator see the description of 
@code{box-around-ellipse}.
+This function returns two datasets as a @code{gal_data_t} linked list.
+The top element of the list is the height and its next element is the width.
+@end deffn
+
 @deffn  Macro GAL_ARITHMETIC_OP_MAKENEW
 Create a new, zero-valued dataset with an unsigned 8-bit data type.
 The length along each dimension of the dataset should be given as a single 
list of @code{gal_data_t}s.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 4be99df..bce3c2f 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
 
+#include <gnuastro/box.h>
 #include <gnuastro/list.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/units.h>
@@ -1956,6 +1957,77 @@ arithmetic_function_binary_flt(int operator, int flags, 
gal_data_t *il,
 
 
 
+gal_data_t *
+arithmetic_box_around_ellipse(gal_data_t *d1, gal_data_t *d2,
+                              gal_data_t *d3, int operator, int flags)
+{
+  size_t i;
+  double *a, *b, *pa, out[2];
+  gal_data_t *a_data, *b_data, *pa_data;
+
+  /* Basic sanity check. */
+  if( d1->size != d2->size || d1->size != d3->size )
+    error(EXIT_FAILURE, 0, "%s: the three operands to this "
+          "function don't have the same number of elements",
+          __func__);
+
+  /* Convert the two inputs into double. Note that if the user doesn't want
+     to free the inputs, we should make a copy of 'a_data' and 'b_data'
+     because the output will also be written in them. */
+  a_data=( ( d1->type==GAL_TYPE_FLOAT64
+             || flags & GAL_ARITHMETIC_FLAG_FREE )
+           ? d1
+           : gal_data_copy_to_new_type(d1, GAL_TYPE_FLOAT64) );
+  b_data=( ( d2->type==GAL_TYPE_FLOAT64
+             || flags & GAL_ARITHMETIC_FLAG_FREE )
+           ? d2
+           : gal_data_copy_to_new_type(d2, GAL_TYPE_FLOAT64) );
+  pa_data=( d3->type==GAL_TYPE_FLOAT64
+            ? d3
+            : gal_data_copy_to_new_type(d3, GAL_TYPE_FLOAT64) );
+
+  /* Set the double pointers and do the calculation for each. */
+  a=a_data->array;
+  b=b_data->array;
+  pa=pa_data->array;
+  for(i=0;i<a_data->size;++i)
+    {
+      /* If the minor axis has a larger value, print a warning and set the
+         value to NaN. */
+      if(b[i]>a[i])
+        {
+          if( (flags & GAL_ARITHMETIC_FLAG_QUIET)==0 )
+            error(EXIT_SUCCESS, 0, "%s: the minor axis (%g) is larger "
+                  "than the major axis (%g), output for this element "
+                  "will be NaN", __func__, b[i], a[i]);
+          out[0]=out[1]=NAN;
+        }
+      else gal_box_bound_ellipse_extent(a[i], b[i], pa[i], out);
+
+      /* Write the output in the inputs (we don't need the inputs any
+         more). */
+      a[i]=out[0];
+      b[i]=out[1];
+    }
+
+  /* Clean up. */
+  gal_data_free(pa_data);
+  if(flags & GAL_ARITHMETIC_FLAG_FREE)
+    {
+      if(a_data !=d1) gal_data_free(d1);
+      if(b_data !=d2) gal_data_free(d2);
+      if(pa_data!=d3) gal_data_free(d3);
+    }
+
+  /* Make the output as a list and return it. */
+  b_data->next=a_data;
+  return b_data;
+}
+
+
+
+
+
 /* Make a new dataset. */
 gal_data_t *
 arithmetic_makenew(gal_data_t *sizes)
@@ -2214,6 +2286,10 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
   else if (!strcmp(string, "float64"))
     { op=GAL_ARITHMETIC_OP_TO_FLOAT64;        *num_operands=1;  }
 
+  /* Surrounding box. */
+  else if (!strcmp(string, "box-around-ellipse"))
+    { op=GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE;*num_operands=3;  }
+
   /* New dataset. */
   else if (!strcmp(string, "makenew"))
     { op=GAL_ARITHMETIC_OP_MAKENEW;           *num_operands=-1;  }
@@ -2324,6 +2400,8 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_TO_FLOAT32:      return "float32";
     case GAL_ARITHMETIC_OP_TO_FLOAT64:      return "float64";
 
+    case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE: return "box-around-ellipse";
+
     case GAL_ARITHMETIC_OP_MAKENEW:         return "makenew";
 
     default:                                return NULL;
@@ -2503,6 +2581,15 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
       out=arithmetic_change_type(d1, operator, flags);
       break;
 
+    /* Calculate the width and height of a box surrounding an ellipse with
+       a certain major axis, minor axis and position angle. */
+    case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      d3 = va_arg(va, gal_data_t *);
+      out=arithmetic_box_around_ellipse(d1, d2, d3, operator, flags);
+      break;
+
     /* Build dataset from scratch. */
     case GAL_ARITHMETIC_OP_MAKENEW:
       d1 = va_arg(va, gal_data_t *);
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 560601c..9f4c654 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -168,6 +168,8 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_TO_FLOAT32,   /* Convert to float32.                   */
   GAL_ARITHMETIC_OP_TO_FLOAT64,   /* Convert to float64.                   */
 
+  GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE, /* Width/Height of box over ellipse*/
+
   GAL_ARITHMETIC_OP_MAKENEW,      /* Build a new dataset, containing zeros.*/
 
   GAL_ARITHMETIC_OP_LAST_CODE,    /* Last code of the library operands.    */