[gnuastro-commits] master 94e6d5b7: Library (arithmetic.h): new operator to rotate a point

[Mohammad Akhlaghi]

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

    Library (arithmetic.h): new operator to rotate a point
    
    Until now, if you wanted to rotate a set of points within a table, you had
    to do the mathematics yourself in a complex column arithmetic operation!
    
    With this commit, to simplify the job, a new 'rotate-coord' has been added
    to the general Arithmetic operators (so it can be used in both the
    Arithmetic and Table programs).
---
 NEWS                        |   3 +
 bin/arithmetic/arithmetic.c |  12 +++-
 bin/table/arithmetic.c      |  12 +++-
 doc/gnuastro.texi           |  54 ++++++++++++++---
 lib/arithmetic.c            | 142 ++++++++++++++++++++++++++++++++++++++++++--
 lib/gnuastro/arithmetic.h   |   1 +
 6 files changed, 208 insertions(+), 16 deletions(-)

diff --git a/NEWS b/NEWS
index b1728c15..e56fdc25 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,9 @@ See the end of the file for license conditions.
 ** New features
 *** Arithmetic
   - New operators:
+    - rotate-coord: given a 2D point's coordinates, return the coordinates
+      after it has been rotated by your requested angle around your
+      requested center (see documentation for example).
     - mad: Median Absolute Deviation (MAD) stacking.
     - madclip-mad: MAD after MAD-clipping stacking.
     - madclip-std: Standard deviation after MAD-clipping stacking.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 5f82fbd5..89c45c5c 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1601,7 +1601,7 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
   size_t i;
   unsigned int numop;
   int flags = GAL_ARITHMETIC_FLAGS_BASIC;
-  gal_data_t *d1=NULL, *d2=NULL, *d3=NULL, *d4=NULL;
+  gal_data_t *d1=NULL, *d2=NULL, *d3=NULL, *d4=NULL, *d5=NULL;
 
   /* Set the operating-mode flags if necessary. */
   if(p->cp.quiet) flags |= GAL_ARITHMETIC_FLAG_QUIET;
@@ -1642,6 +1642,14 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
           d1=operands_pop(p, operator_string);
           break;
 
+        case 5:
+          d5=operands_pop(p, operator_string);
+          d4=operands_pop(p, operator_string);
+          d3=operands_pop(p, operator_string);
+          d2=operands_pop(p, operator_string);
+          d1=operands_pop(p, operator_string);
+          break;
+
         case -1:
           /* This case is when the number of operands is itself an
              operand. So except for operators that have high-level
@@ -1715,7 +1723,7 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
          when the operator doesn't need three operands, the extra
          arguments will be ignored. */
       operands_add(p, NULL, gal_arithmetic(operator, p->cp.numthreads,
-                                           flags, d1, d2, d3, d4));
+                                           flags, d1, d2, d3, d4, d5));
 
       /* Operators with special attention afterwards. */
       switch(operator)
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index a44e4d46..bef8c6f9 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -1052,7 +1052,7 @@ arithmetic_operator_run(struct tableparams *p,
                         gal_data_t **stack)
 {
   int flags=GAL_ARITHMETIC_FLAGS_BASIC;
-  gal_data_t *d1=NULL, *d2=NULL, *d3=NULL, *d4=NULL;
+  gal_data_t *d1=NULL, *d2=NULL, *d3=NULL, *d4=NULL, *d5=NULL;
 
   /* Set the operating-mode flags if necessary. */
   if(p->cp.quiet) flags |= GAL_ARITHMETIC_FLAG_QUIET;
@@ -1093,6 +1093,14 @@ arithmetic_operator_run(struct tableparams *p,
           d1=arithmetic_stack_pop(stack, token->operator, NULL);
           break;
 
+        case 5:
+          d5=arithmetic_stack_pop(stack, token->operator, NULL);
+          d4=arithmetic_stack_pop(stack, token->operator, NULL);
+          d3=arithmetic_stack_pop(stack, token->operator, NULL);
+          d2=arithmetic_stack_pop(stack, token->operator, NULL);
+          d1=arithmetic_stack_pop(stack, token->operator, NULL);
+          break;
+
         case -1:
           error(EXIT_FAILURE, 0, "operators with a variable number of "
                 "operands are not yet implemented. Please contact us at "
@@ -1114,7 +1122,7 @@ arithmetic_operator_run(struct tableparams *p,
          ignored. */
       gal_list_data_add(stack, gal_arithmetic(token->operator,
                                               p->cp.numthreads,
-                                              flags, d1, d2, d3, d4) );
+                                              flags, d1, d2, d3, d4, d5) );
 
       /* Reset the meta-data for the element that was just put on the
          stack. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index e50d48fa..24b3e23f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -572,7 +572,7 @@ Arithmetic operators
 * Bitwise operators::           Work on bits within one pixel.
 * Numerical type conversion operators::  Convert the numeric datatype of a 
dataset.
 * Random number generators::    Random numbers can be used to add noise for 
example.
-* Box shape operators::         Return edges of 2D boxes.
+* Coordinate and border operators::         Return edges of 2D boxes.
 * Loading external columns::    Read a column from a table into the stack.
 * Size and position operators::  Extracting image size and pixel positions.
 * Building new dataset and stack management::  How to construct an empty 
dataset from scratch.
@@ -20602,7 +20602,7 @@ Reading NaN as a floating point number in Gnuastro is 
not case-sensitive.
 * Bitwise operators::           Work on bits within one pixel.
 * Numerical type conversion operators::  Convert the numeric datatype of a 
dataset.
 * Random number generators::    Random numbers can be used to add noise for 
example.
-* Box shape operators::         Return edges of 2D boxes.
+* Coordinate and border operators::         Return edges of 2D boxes.
 * Loading external columns::    Read a column from a table into the stack.
 * Size and position operators::  Extracting image size and pixel positions.
 * Building new dataset and stack management::  How to construct an empty 
dataset from scratch.
@@ -22380,7 +22380,7 @@ Convert the type of the popped operand to 64-bit 
(double precision) floating poi
 The internal conversion of C will be used.
 @end table
 
-@node Random number generators, Box shape operators, Numerical type conversion 
operators, Arithmetic operators
+@node Random number generators, Coordinate and border operators, Numerical 
type conversion operators, Arithmetic operators
 @subsubsection Random number generators
 
 When you simulate data (for example, see @ref{Sufi simulates a detection}), 
everything is ideal and there is no noise!
@@ -22789,12 +22789,50 @@ $ echo 1000 \
 These columns can easily be placed in the format for @ref{MakeProfiles} to be 
inserted into an image automatically.
 @end table
 
-@node Box shape operators, Loading external columns, Random number generators, 
Arithmetic operators
-@subsubsection Box shape operators
+@node Coordinate and border operators, Loading external columns, Random number 
generators, Arithmetic operators
+@subsubsection Coordinate and border operators
 
-The operators here help you in defining or using coordinates that form a 
``box'' (a rectangular region).
+The operators here help you in defining or manipulating coordinates.
+For examples to define the ``box'' (a rectangular region) that surrounds an 
ellipse or to rotate a point around a reference point.
 
 @table @command
+@item rotate-coord
+Rotate the given point (horizontal and vertical coordinates given in 5th and 
4th popped operands) around a center/reference point (coordinates given in the 
3rd and 2nd popped operands) by a given angle (first popped operand).
+
+For example, if you want to trace the outer edge of a circle centered on 
(1.23,45.6) with a radius of 0.78, you can use this operator like below.
+The logic is that we assume a single point that is located on 0.78 units after 
the center on the horizontal axis (the point's vertical axis position is the 
same as the center).
+We then rotate this point in each row by one degree to build the circle's 
circumference.
+
+@verbatim
+$ cx=1.23
+$ cy=45.6
+$ rad=0.78
+$ seq 0 360 \
+      | awk '{print '$rad'+'$cx', '$cy', $1}' \
+      | asttable -c'arith $1 $2 '$cx' '$cy' $3 rotate-coord' \
+                 --output=circle.fits
+
+##  Within TOPCAT, after opening "Plane Plot", within "Axes" select
+## "Aspect lock" so the steps in both axis is the same.
+$ astscript-fits-view circle.fits
+@end verbatim
+
+If you want the points to create a circle on the celestial sphere, you can use 
the @code{eq-j2000-from-flat} operator after this one (see @ref{Column 
arithmetic}):
+
+@verbatim
+$ seq 0 360 \
+      | awk '{print '$rad'+'$cx', '$cy', $1}' \
+      | asttable -c'arith $1 $2 '$cx' '$cy' $3 rotate-coord \
+                          '$cx' '$cy' TAN eq-j2000-from-flat' \
+                          --output=circle-on-sky.fits
+@end verbatim
+
+When you open TOPCAT, if you open the ``Plane Plot'', you will see an ellipse.
+However, if you open ``Sky Plot'' (from the ``Graphics'' menu), and select the 
first and second columns respectively, you will see a circle.
+
+The center coordinates and angle can be fixed for all the rows (as in the 
example above) or be different for every row.
+Recall that if you want these to change on every row, you should give the 
column name (or number followed by @code{$}) for these operands instead of the 
constant number above.
+
 @item box-around-ellipse
 Return the width (along horizontal) and height (along vertical) of a box that 
encompasses 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}.
@@ -22911,7 +22949,7 @@ It is smaller by a factor of @mymath{\cos({\delta})}.
 Therefore, an angular change (let's call it @mymath{\Delta_{lon}}) along the 
small circle defined by the fixed declination of @mymath{\delta} corresponds to 
@mymath{\Delta_{lon}/\cos(\delta)} on the equator.
 @end table
 
-@node Loading external columns, Size and position operators, Box shape 
operators, Arithmetic operators
+@node Loading external columns, Size and position operators, Coordinate and 
border operators, Arithmetic operators
 @subsubsection Loading external columns
 
 In the Arithmetic program, you can always load new dataset by simply giving 
their name.
@@ -39541,7 +39579,7 @@ For the declination, a simple addition/subtraction is 
enough.
 Also, on the equator (where the RA is defined), a simple addition/subtraction 
along the RA is fine.
 However, at other declinations, the new RA after a shift needs special 
treatment, such that close to the poles, a shift of 1 degree can correspond to 
a new RA that is much more distant than the original RA.
 Assuming a point at Right Ascension (RA) and Declination of @mymath{\alpha} 
and @mymath{\delta}, a shift of @mymath{R} degrees along the positive RA 
direction corresponds to a right ascension of 
@mymath{\alpha+\frac{R}{\cos(\delta)}}.
-For more, see the description of @code{box-vertices-on-sphere} in @ref{Box 
shape operators}.
+For more, see the description of @code{box-vertices-on-sphere} in 
@ref{Coordinate and border operators}.
 
 The 8 coordinates of the 4 vertices of the box are written in the order below.
 Where ``bottom'' corresponds to a lower declination and ``top'' to higher 
declination, ``left'' corresponds to a larger RA and ``right'' corresponds to 
lower RA.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 25168406..b8e2f891 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -1483,9 +1483,9 @@ arithmetic_stitch(int flags, gal_data_t *list, gal_data_t 
*fdim)
                        gal_type_sizeof(type)*tmp->dsize[1]);
 
               /* Copying has finished, increment the start for the next
-                 image. Note that in this scenario, the starting pixel for the
-                 next image is on the first row, but tmp->dsize[1] pixels
-                 away. */
+                 image. Note that in this scenario, the starting pixel for
+                 the next image is on the first row, but tmp->dsize[1]
+                 pixels away. */
               oarr += gal_type_sizeof(type)*tmp->dsize[1];
               break;
 
@@ -3275,6 +3275,127 @@ arithmetic_box(gal_data_t *d1, gal_data_t *d2, 
gal_data_t *d3,
 
 
 
+static gal_data_t *
+arithmetic_rotate(int operator, int flags, gal_data_t *d1, gal_data_t *d2,
+                  gal_data_t *d3, gal_data_t *d4, gal_data_t *d5)
+{
+  size_t i;
+  gal_data_t *out=NULL;
+  double x, y, rx, ry, rot, onerot;
+  double     *d1a=NULL, *d2a=NULL, *d3a=NULL, *d4a=NULL, *d5a=NULL;
+  gal_data_t *d1d=NULL, *d2d=NULL, *d3d=NULL, *d4d=NULL, *d5d=NULL;
+
+  /* Basic sanity check. */
+  if(d1->size != d2->size)
+    error(EXIT_FAILURE, 0, "%s: the input coordinate operands (5th and "
+          "4th popped operands) to this function don't have the same "
+          "number of elements", __func__);
+  if(d3->size != d4->size)
+    error(EXIT_FAILURE, 0, "%s: the reference coordinate operands (3rd "
+          "and 2nd popped operands) to this function don't have the same "
+          "number of elements", __func__);
+  if(d3->size>1
+     && (d3->size != d1->size || d4->size!= d1->size) )
+    error(EXIT_FAILURE, 0, "%s: the reference coordinate operands (3rd "
+          "and 2th popped operands) have %zu elements, while the input "
+          "coordinate operands (5th and 4th operands) have %zu elements! "
+          "The reference coordinates should either have a single element "
+          "(to be used for all inputs) or have the same number of "
+          "elements as input elements", __func__, d3->size, d1->size);
+  if(d5->size>1 && d5->size!=d1->size)
+    error(EXIT_FAILURE, 0, "%s: the rotation angle operand (1st popped "
+          "operand) has %zu elements, while the input coordinate operands "
+          "(5th and 4th operands) have %zu elements! The rotation angle "
+          "should either have a single element (to be used for all inputs) "
+          "or have the same number of elements as input elements",
+          __func__, d5->size, d1->size);
+
+
+  /* The datasets may be empty. In this case the output should also be
+     empty (we can have tables and images with 0 rows or pixels!). */
+  if(    d1->size==0 || d1->array==NULL
+      || d2->size==0 || d2->array==NULL
+      || d3->size==0 || d3->array==NULL
+      || d4->size==0 || d4->array==NULL
+      || d5->size==0 || d5->array==NULL )
+    {
+      if(flags & GAL_ARITHMETIC_FLAG_FREE)
+        {
+          gal_data_free(d2); gal_data_free(d3);
+          gal_data_free(d4); gal_data_free(d5);
+        }
+      if(d1->array) {free(d1->array); d1->array=NULL;}
+      if(d1->dsize) for(i=0;i<d1->ndim;++i) d1->dsize[i]=0;
+      d1->size=0; return d1;
+    }
+
+ /* Convert the 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. */
+  d1d=( d1->type==GAL_TYPE_FLOAT64
+        ? d1
+        : gal_data_copy_to_new_type(d1, GAL_TYPE_FLOAT64) );
+  d2d=( d2->type==GAL_TYPE_FLOAT64
+        ? d2
+        : gal_data_copy_to_new_type(d2, GAL_TYPE_FLOAT64) );
+  d3d=( d3->type==GAL_TYPE_FLOAT64
+        ? d3
+        : gal_data_copy_to_new_type(d3, GAL_TYPE_FLOAT64) );
+  d4d=( d4->type==GAL_TYPE_FLOAT64
+        ? d4
+        : gal_data_copy_to_new_type(d4, GAL_TYPE_FLOAT64) );
+  d5d=( d5->type==GAL_TYPE_FLOAT64
+        ? d5
+        : gal_data_copy_to_new_type(d5, GAL_TYPE_FLOAT64) );
+  d1a=d1d->array;
+  d2a=d2d->array;
+  d3a=d3d->array;
+  d4a=d4d->array;
+  d5a=d5d->array;
+
+  /* Do the operation. We will do it in the same space  */
+  onerot = d5d->size==1 ? (d5a[0] * M_PI/180.0f) : NAN;
+  for(i=0; i<d1d->size; ++i)
+    {
+      /* To simplify the readability. */
+      x   =                           d1a[i];
+      y   =                           d2a[i];
+      rx  = d3d->size==1 ? d3a[0] :   d3a[i];
+      ry  = d4d->size==1 ? d4a[0] :   d4a[i];
+      rot = d5d->size==1 ? onerot : ( d5a[i] * M_PI/180.0f );
+
+      /* Write the outputs in the inputs (note that above, we copied the
+         inputs in the short variables). But first, subtract the reference
+         X and Y so we rotate around that (we will add the reference
+         afterwards). */
+      x = x-rx;
+      y = y-ry;
+      d1a[i] = x*cos(rot) - y*sin(rot) + rx;
+      d2a[i] = x*sin(rot) + y*cos(rot) + ry;
+    }
+
+  /* Set the output. */
+  out=d2d;
+  out->next=d1d;
+
+  /* Clean up. */
+  if(flags & GAL_ARITHMETIC_FLAG_FREE)
+    {
+      if(d1d!=d1)                         gal_data_free(d1);
+      if(d2d!=d2)                         gal_data_free(d2);
+      if(d3d!=d3) { gal_data_free(d3d); } gal_data_free(d3);
+      if(d4d!=d4) { gal_data_free(d4d); } gal_data_free(d4);
+      if(d5d!=d5) { gal_data_free(d5d); } gal_data_free(d5);
+    }
+
+  /* Return the output. */
+  return out;
+}
+
+
+
+
+
 static gal_data_t *
 arithmetic_constants_standard(int operator)
 {
@@ -4052,6 +4173,8 @@ gal_arithmetic_set_operator(char *string, size_t 
*num_operands)
     { op=GAL_ARITHMETIC_OP_FINESTRUCTURE;     *num_operands=0;  }
 
   /* Surrounding box. */
+  else if (!strcmp(string, "rotate-coord"))
+    { op=GAL_ARITHMETIC_OP_ROTATE_COORD;      *num_operands=5;  }
   else if (!strcmp(string, "box-around-ellipse"))
     { op=GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE;*num_operands=3;  }
   else if (!strcmp(string, "box-vertices-on-sphere"))
@@ -4254,6 +4377,7 @@ gal_arithmetic_operator_string(int operator)
     case GAL_ARITHMETIC_OP_AVOGADRO:        return "avogadro";
     case GAL_ARITHMETIC_OP_FINESTRUCTURE:   return "fine-structure";
 
+    case GAL_ARITHMETIC_OP_ROTATE_COORD:    return "rotate-coord";
     case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE: return "box-around-ellipse";
     case GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE: return "vertices-on-sphere";
 
@@ -4347,7 +4471,7 @@ gal_data_t *
 gal_arithmetic(int operator, size_t numthreads, int flags, ...)
 {
   va_list va;
-  gal_data_t *d1, *d2, *d3, *d4, *out=NULL;
+  gal_data_t *d1, *d2, *d3, *d4, *d5, *out=NULL;
 
   /* Prepare the variable arguments (starting after the flags argument). */
   va_start(va, flags);
@@ -4642,6 +4766,16 @@ gal_arithmetic(int operator, size_t numthreads, int 
flags, ...)
         }
       break;
 
+    /* Rotate coordinates about a reference. */
+    case GAL_ARITHMETIC_OP_ROTATE_COORD:
+      d1 = va_arg(va, gal_data_t *);
+      d2 = va_arg(va, gal_data_t *);
+      d3 = va_arg(va, gal_data_t *);
+      d4 = va_arg(va, gal_data_t *);
+      d5 = va_arg(va, gal_data_t *);
+      out=arithmetic_rotate(operator, flags, d1, d2, d3, d4, d5);
+      break;
+
     /* Size and position operators. */
     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 66933e22..a2bf91f5 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -228,6 +228,7 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_TO_FLOAT32,   /* Convert to float32.                   */
   GAL_ARITHMETIC_OP_TO_FLOAT64,   /* Convert to float64.                   */
 
+  GAL_ARITHMETIC_OP_ROTATE_COORD, /* Return input coords after rotation.   */
   GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE, /* Width/Height of box over ellipse*/
   GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE, /* Vert. from center and width */