[gnuastro-commits] master 79acd68: MakeProfiles: new azimuthal angle profile

[Mohammad Akhlaghi]

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

    MakeProfiles: new azimuthal angle profile
    
    In order to make certain features (that only cover a certain angular
    range), it was necessary to have every pixel's azimuthal angle (along the
    ellipse of fixed radius that covers it). For example if it is necessary to
    mask a certain arc over an image, or to build tidal tails or tidal shocks
    over mock galaxies.
    
    With this commit, a new function has been added to MakeProfiles that is
    called 'azimuth' (or the code of '9'). It will return the azimuthal angle
    for every pixel (in degrees).
    
    While I was testing this new feature, I noticed that the function returning
    the size of a box (in pixels) containing an ellipse doesn't can miss one
    pixel along the length if its fractional value is more than 0.5 due to the
    conversion of floating point to integer. This has also been corrected with
    this commit, for both 2D ellipses and 3D ellipsoids.
    
    This feature was added after some very productive discussions with Fernando
    Buitrago and Matthias Kluge.
---
 NEWS                         |  9 ++++++++
 bin/mkprof/args.h            |  2 +-
 bin/mkprof/main.h            |  5 ++++-
 bin/mkprof/mkprof.h          |  2 +-
 bin/mkprof/oneprofile.c      | 10 +++++++++
 bin/mkprof/profiles.c        | 50 ++++++++++++++++++++++++++++++++++++++++++++
 bin/mkprof/profiles.h        |  3 +++
 bin/mkprof/ui.c              | 19 ++++++++++++++++-
 doc/announce-acknowledge.txt |  2 ++
 doc/gnuastro.texi            | 11 +++++++---
 lib/box.c                    | 20 +++++++++++++++---
 11 files changed, 123 insertions(+), 10 deletions(-)

diff --git a/NEWS b/NEWS
index 7db74de..6cab5ea 100644
--- a/NEWS
+++ b/NEWS
@@ -32,6 +32,15 @@ See the end of the file for license conditions.
      estimating the surface brightness error.
    --sberror: error in measuring the surface brightness (mag/arcsec^2).
 
+  MakeProfiles:
+   - New type of profile showing the azimuthal angle (in degrees, along the
+     elliptical circumference of fixed radius) of each pixel. In
+     combination with the radial distance profile, you can now create
+     complex features in polar coordinates, such as tidal tails or tidal
+     shocks (using the Arithmetic program to mix the radius and azimuthal
+     angle through a function to create your desired features). Implemented
+     after discussions with Fernando Buitrago and Matthias Kluge.
+
   Match:
    - When called with '--notmatched --outcols=AAA,BBB', Match will append
      non-matching rows of second table into first table's rows (for columns
diff --git a/bin/mkprof/args.h b/bin/mkprof/args.h
index d6a8887..78d849b 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -391,7 +391,7 @@ struct argp_option program_options[] =
       0,
       "sersic (1), moffat (2), gaussian (3), point (4), "
       "flat (5), circumference (6), distance (7), "
-      "radial-table (8).",
+      "radial-table (8), azimuth (9).",
       UI_GROUP_CATALOG,
       &p->fcol,
       GAL_TYPE_STRING,
diff --git a/bin/mkprof/main.h b/bin/mkprof/main.h
index d794946..ef37f06 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -41,6 +41,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* Some constants */
 #define EPSREL_FOR_INTEG   2
 #define DEGREESTORADIANS   M_PI/180.0
+#define RADIANSTODEGREES   180.0/M_PI
+
 
 
 /* Modes to interpret coordinates. */
@@ -66,7 +68,8 @@ enum profile_types
   PROFILE_FLAT,                 /* Flat profile.               */
   PROFILE_CIRCUMFERENCE,        /* Circumference profile.      */
   PROFILE_DISTANCE,             /* Elliptical radius of pixel. */
-  PROFILE_CUSTOM,          /* Radial prof. in file/table. */
+  PROFILE_CUSTOM,               /* Radial prof. in file/table. */
+  PROFILE_AZIMUTH,              /* Azimuthal angle at distance.*/
 
   PROFILE_MAXIMUM_CODE,         /* Just for a sanity check.    */
 };
diff --git a/bin/mkprof/mkprof.h b/bin/mkprof/mkprof.h
index 8411b10..3ebce60 100644
--- a/bin/mkprof/mkprof.h
+++ b/bin/mkprof/mkprof.h
@@ -31,7 +31,7 @@ struct mkonthread
 {
   /* General parameters: */
   double                r;   /* Elliptical radius at this point.      */
-  double         coord[3];   /* Pixel coordinate.                     */
+  double         coord[3];   /* Pixel coord. (from profile center).   */
   double         lower[3];   /* Coordinates of lower pixel position.  */
   double        higher[3];   /* Coordinates of higher pixel position. */
   double             c[3];   /* Cosine of position angle(s).          */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index caa3140..5da9e53 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -694,6 +694,16 @@ oneprofile_set_prof_params(struct mkonthread *mkp)
 
 
 
+
+
+    case PROFILE_AZIMUTH:
+      mkp->profile          = profiles_azimuth;
+      mkp->truncr           = tp ? p->t[id] : p->t[id]*p->r[id];
+      mkp->correction       = 0;
+      break;
+
+
+
     case PROFILE_CUSTOM:
       mkp->profile          = profiles_custom_table;
       mkp->truncr           = tp ? p->t[id] : p->t[id]*p->r[id];
diff --git a/bin/mkprof/profiles.c b/bin/mkprof/profiles.c
index d6a3984..ee40e6a 100644
--- a/bin/mkprof/profiles.c
+++ b/bin/mkprof/profiles.c
@@ -52,6 +52,56 @@ profiles_radial_distance(struct mkonthread *mkp)
 
 
 
+/* Azimuthal angle.
+
+   Assuming theta is the azimuthal angle (along a constant radius), then an
+   ellipse is defined by:
+
+      x=a*cos(theta)
+      y=b*sin(theta)
+
+   Now, let's take 'phi' to be the angle in the normal equi-distant
+   (non-elliptical coordinates). Therefore:
+
+      tan(phi)=y/x
+
+   So:
+                 b*sin(theta)
+      tan(phi) = ------------ = b/a * tan(theta)
+                 a*cos(theta)
+
+   Therefore we can find theta like this (where q=b/a):
+
+      theta = atan( a/b * tan(phi) )
+            = atan( a/b * y/x )
+            = atan( y / (x*q) )
+
+   However, the 'x' and 'y' above are only for the case where the ellipse
+   has a position angle of zero (its major axis is aligned with the
+   horizontal axis. When the ellipse is rotated, we first need to rotate
+   the 'mkp->coord' values, and then use 'x' and 'y' like above. */
+double
+profiles_azimuth(struct mkonthread *mkp)
+{
+  /* We are rotating by the inverse (multiplied by -1) position angle. */
+  double x =      mkp->coord[0]*mkp->c[0] + mkp->coord[1]*mkp->s[0];
+  double y = -1 * mkp->coord[0]*mkp->s[0] + mkp->coord[1]*mkp->c[0];
+
+  /* The normal 'atan' function will only return values between -90 to +90
+     degrees, or only two quadrants (because it only takes a single
+     value). However, with 'atan2' we can get the full -180 to +180 degrees
+     (all four quadrants). */
+  double d=atan2(y, x*mkp->q[0])*RADIANSTODEGREES;
+
+  /* 'atan2' returns values between -180 to 180 deg. But we want values
+     from 0 to 360, so we'll add the negative ones by 360. */
+  return d>0 ? d : d+360.0f;
+}
+
+
+
+
+
 /* Read the values based on the distance from a table. */
 double
 profiles_custom_table(struct mkonthread *mkp)
diff --git a/bin/mkprof/profiles.h b/bin/mkprof/profiles.h
index 72c7b9b..64e9afc 100644
--- a/bin/mkprof/profiles.h
+++ b/bin/mkprof/profiles.h
@@ -27,6 +27,9 @@ double
 profiles_radial_distance(struct mkonthread *mkp);
 
 double
+profiles_azimuth(struct mkonthread *mkp);
+
+double
 profiles_custom_table(struct mkonthread *mkp);
 
 double
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index acbc2d0..cc4b061 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -124,6 +124,9 @@ ui_profile_name_read(char *string, size_t row)
   else if ( !strcmp("distance", string) )
     return PROFILE_DISTANCE;
 
+  else if ( !strcmp("azimuth", string) )
+    return PROFILE_DISTANCE;
+
   else if ( !strcmp("custom", string) )
     return PROFILE_CUSTOM;
 
@@ -159,6 +162,7 @@ ui_profile_name_write(int profile_code)
     case PROFILE_FLAT:           return "flat";
     case PROFILE_CIRCUMFERENCE:  return "circum";
     case PROFILE_DISTANCE:       return "distance";
+    case PROFILE_AZIMUTH:        return "azimuth";
     default:
       error(EXIT_FAILURE, 0, "%s: %d not recognized as a profile code",
             __func__, profile_code);
@@ -428,6 +432,7 @@ ui_parse_kernel(struct argp_option *option, char *arg,
         case PROFILE_FLAT:          need = kernel->flag==2 ? 1 : 2;  break;
         case PROFILE_CIRCUMFERENCE: need = kernel->flag==2 ? 1 : 2;  break;
         case PROFILE_DISTANCE:      need = kernel->flag==2 ? 1 : 2;  break;
+        case PROFILE_AZIMUTH:       need = kernel->flag==2 ? 1 : 2;  break;
         default:
           error_at_line(EXIT_FAILURE, 0, filename, lineno, "%s: a bug! "
                         "Please contact us at %s to correct the issue. "
@@ -971,7 +976,9 @@ ui_read_cols_3d(struct mkprofparams *p)
               corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_UINT8);
               p->f=corrtype->array;
 
-              /* Check if they are in the correct range. */
+              /* Check if they are in the correct range. For profile names
+                 given as string, a non-matching string will result in an
+                 error, so there is no need for this in that scenario. */
               for(i=0;i<p->num;++i)
                 if(p->f[i]<=PROFILE_INVALID || p->f[i]>=PROFILE_MAXIMUM_CODE)
                   error(EXIT_FAILURE, 0, "%s: row %zu, the function "
@@ -985,6 +992,16 @@ ui_read_cols_3d(struct mkprofparams *p)
                         "    $ info %s\n", p->catname, i+1, p->f[i],
                         PROFILE_INVALID, PROFILE_MAXIMUM_CODE, PROGRAM_EXEC);
             }
+
+          /* General profile sanity checks. */
+          for(i=0;i<p->num;++i)
+            {
+              /* Azimuthal profile not yet supported for ellipsoids. */
+              if(p->f[i]==PROFILE_AZIMUTH)
+                error(EXIT_FAILURE, 0, "%s: row %zu: the azimuthal "
+                      "angle profile is not yet supported for 3D "
+                      "datasets", p->catname, i+1);
+            }
           break;
 
         case 5:
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 0a6e709..566334e 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,7 +1,9 @@
 Alphabetically ordered list to acknowledge in the next release.
 
+Fernando Buitrago
 Mark Calabretta
 Leslie Hunt
+Matthias Kluge
 Juan Miro
 Juan Molina Tobar
 Zahra Sharbaf
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index af397db..9bda624 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -20087,10 +20087,10 @@ See the description of @option{--config} in 
@ref{Operating mode options}.}.
 Please see @ref{Sufi simulates a detection} for a very complete tutorial 
explaining how one could use MakeProfiles in conjunction with other Gnuastro's 
programs to make a complete simulated image of a mock galaxy.
 
 @menu
-* MakeProfiles catalog::        Required catalog properties.
+* MakeProfiles catalog::           Required catalog properties.
 * MakeProfiles profile settings::  Configuration parameters for all profiles.
-* MakeProfiles output dataset::  The canvas/dataset to build profiles over.
-* MakeProfiles log file::       A description of the optional log file.
+* MakeProfiles output dataset::    The canvas/dataset to build profiles over.
+* MakeProfiles log file::          A description of the optional log file.
 @end menu
 
 @node MakeProfiles catalog, MakeProfiles profile settings, Invoking astmkprof, 
Invoking astmkprof
@@ -20170,6 +20170,11 @@ In the latter case, just note that the central values 
are going to be incorrect
 @item
 Custom profile with `@code{custom}' or `@code{8}'.
 The values to use for each radial interval should be in the table given to 
@option{--customtable}. For more, see @ref{MakeProfiles profile settings}.
+
+@item
+Azimuthal angle profile with `@code{azimuth}' or `@code{9}'.
+Every pixel within the truncation radius will be given its azimuthal angle (in 
degrees, from 0 to 360) from the major axis.
+In combination with the radial distance profile, you can now create complex 
features in polar coordinates, uch as tidal tails or tidal shocks (using the 
Arithmetic program to mix the radius and azimuthal angle through a function to 
create your desired features).
 @end itemize
 
 @item --rcol=STR/INT
diff --git a/lib/box.c b/lib/box.c
index 7f6353d..d5d4313 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -88,12 +88,18 @@ gal_box_bound_ellipse(double a, double b, double theta_deg, 
long *width)
   /* Find the extent of the ellipse. */
   gal_box_bound_ellipse_extent(a, b, theta_deg, extent);
 
+  /* The returned value is an integer. So if the fractional component of
+     the extent is larger than 0.5, add one to it so we still cover the
+     whole desired region. */
+  extent[0] += extent[0] - ((long)(extent[0])) > 0.5f ? 0.5f : 0.0f;
+  extent[1] += extent[1] - ((long)(extent[1])) > 0.5f ? 0.5f : 0.0f;
+
   /* max_x and max_y are calculated from the center of the ellipse. We
      want the final height and width of the box enclosing the
      ellipse. So we have to multiply them by two, then take one from
      them (for the center). */
-  width[0] = 2 * ( (long)extent[0] ) + 1;
-  width[1] = 2 * ( (long)extent[1] ) + 1;
+  width[0] = 2 * ( (long)(extent[0]) ) + 1;
+  width[1] = 2 * ( (long)(extent[1]) ) + 1;
 }
 
 
@@ -261,7 +267,15 @@ gal_box_bound_ellipsoid(double *semiaxes, double 
*euler_deg, long *width)
 
   /* Find the width along each dimension. */
   for(i=0;i<3;++i)
-    width[i] = 2 * ( (long)extent[i] ) + 1;
+    {
+      /* The returned width along each dimension is an integer. So if the
+         fractional component of the extent is larger than 0.5, add one to
+         it so we still cover the whole desired region. */
+      extent[i] += extent[i] - ((long)(extent[i])) > 0.5f ? 0.5f : 0.0f;
+
+      /* Set the width. */
+      width[i] = 2 * ( (long)extent[i] ) + 1;
+    }
 }