[gnuastro-commits] master 1cb1e40: Library (dimension.h): new elliptical distance function

[Mohammad Akhlaghi]

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

    Library (dimension.h): new elliptical distance function
    
    Until now, Gnuastro's MakeProfiles program would internally calculate the
    elliptical distance of a point from the center of an ellipse. But it was
    too deeply intertwined with its internal structure.
    
    With this commit, a copy of that function has been brought into the
    library, enabling other programs to easily use this feature also (as well
    as users in their custom programs).
---
 NEWS                     |  1 +
 doc/gnuastro.texi        | 20 ++++++++++++++++++
 lib/dimension.c          | 54 ++++++++++++++++++++++++++++++++++++++++++++++++
 lib/gnuastro/dimension.h |  4 ++++
 4 files changed, 79 insertions(+)

diff --git a/NEWS b/NEWS
index 04a2954..0930039 100644
--- a/NEWS
+++ b/NEWS
@@ -112,6 +112,7 @@ See the end of the file for license conditions.
      without a square adjacency matrix (which can consume major RAM).
    - gal_blank_flag_remove: Remove all flagged elements in a dataset.
    - gal_blank_remove_rows: remove all rows that have at least one blank.
+   - gal_dimension_dist_elliptical: Elliptical dist. of a point from center.
    - gal_fits_hdu_datasum: calculate DATASUM of given HDU in given FITS file.
    - gal_fits_hdu_datasum_ptr: calculate DATASUM of opened FITS file pointer.
    - gal_pointer_allocate_ram_or_mmap: allocate space either in RAM or as a
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f047a3d..8220749 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -17551,6 +17551,8 @@ Everything else after that (when the pixel index and 
its radial profile have ent
 
 
 
+
+
 @node PSF, Stars, Defining an ellipse and ellipsoid, Modeling basics
 @subsubsection Point spread function
 
@@ -21628,6 +21630,8 @@ already be allocated before calling this function.
 @end deftypefun
 
 @deftypefun size_t gal_dimension_dist_manhattan (size_t @code{*a}, size_t 
@code{*b}, size_t @code{ndim})
+@cindex Manhattan distance
+@cindex Distance, Manhattan
 Return the manhattan distance (see
 @url{https://en.wikipedia.org/wiki/Taxicab_geometry, Wikipedia}) between
 the two coordinates @code{a} and @code{b} (each an array of @code{ndim}
@@ -21639,6 +21643,22 @@ Return the radial distance between the two coordinates 
@code{a} and
 @code{b} (each an array of @code{ndim} elements).
 @end deftypefun
 
+@deftypefun float gal_dimension_dist_elliptical (double @code{*center}, double 
@code{*pa_deg}, double @code{*q}, size_t @code{ndim}, double @code{*point})
+@cindex Ellipse
+@cindex Ellipoid
+@cindex Axis ratio
+@cindex Position angle
+@cindex Elliptical distance
+@cindex Ellipsoidal distance
+@cindex Distance, elliptical/ellipsoidal
+Return the elliptical/ellipsoidal distance of the single point @code{point} 
(containing @code{ndim} values: coordinates of the point in each dimension) 
from an ellipse that is defined by @code{center}, @code{pa_deg} and @code{q}.
+@code{center} is the coordinates of the ellipse center (also with @code{ndim} 
elements). @code{pa} is the position-angle in degrees (the angle of the 
semi-major axis from the first dimension in a 2D ellipse) and @code{q} is the 
axis ratio.
+
+In a 2D ellipse, @code{pa} and @code{q} are a single-element array.
+However, in a 3D ellipsoid, @code{pa} must have three elements, and @code{q} 
must have 2 elements.
+For more see @ref{Defining an ellipse and ellipsoid}.
+@end deftypefun
+
 @deftypefun {gal_data_t *} gal_dimension_collapse_sum (gal_data_t @code{*in}, 
size_t @code{c_dim}, gal_data_t @code{*weight})
 Collapse the input dataset (@code{in}) along the given dimension
 (@code{c_dim}, in C definition: starting from zero, from the slowest
diff --git a/lib/dimension.c b/lib/dimension.c
index 8187fb4..d0ec430 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -285,6 +285,60 @@ gal_dimension_dist_radial(size_t *a, size_t *b, size_t 
ndim)
 
 
 
+/* Elliptical distance of a point from a given center. */
+#define DEGREESTORADIANS   M_PI/180.0
+float
+gal_dimension_dist_elliptical(double *center, double *pa_deg, double *q,
+                              size_t ndim, double *point)
+{
+  double Xr, Yr, Zr;        /* Rotated x, y, z. */
+  double q1=q[0], q2;
+  double c1=cos( pa_deg[0] * DEGREESTORADIANS ), c2, c3;
+  double s1= sin( pa_deg[0] * DEGREESTORADIANS ), s2, s3;
+  double x=center[0]-point[0], y=center[1]-point[1], z;
+
+  /* Find the distance depending on the dimension. */
+  switch(ndim)
+    {
+    case 2:
+      /* The parenthesis aren't necessary, but help in readability and
+         avoiding human induced bugs. */
+      Xr = x * ( c1       )     +   y * ( s1 );
+      Yr = x * ( -1.0f*s1 )     +   y * ( c1 );
+      return sqrt( Xr*Xr + Yr*Yr/q1/q1 );
+      break;
+
+    case 3:
+      /* Define the necessary parameters. */
+      q2=q[1];
+      z=center[2]-point[2];
+      c2 = cos( pa_deg[1] * DEGREESTORADIANS );
+      s2 = sin( pa_deg[1] * DEGREESTORADIANS );
+      c3 = cos( pa_deg[2] * DEGREESTORADIANS );
+      s3 = sin( pa_deg[2] * DEGREESTORADIANS );
+
+      /* Calculate the distance. */
+      Xr = x*(  c3*c1   - s3*c2*s1 ) + y*( c3*s1   + s3*c2*c1) + z*( s3*s2 );
+      Yr = x*( -1*s3*c1 - c3*c2*s1 ) + y*(-1*s3*s1 + c3*c2*c1) + z*( c3*s2 );
+      Zr = x*(  s1*s2              ) + y*(-1*s2*c1           ) + z*( c2    );
+      return sqrt( Xr*Xr + Yr*Yr/q1/q1 + Zr*Zr/q2/q2 );
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%s: currently only 2 and 3 dimensional "
+            "distances are supported, you have asked for an %zu-dimensional "
+            "dataset", __func__, ndim);
+
+    }
+
+  /* Control should never reach here. */
+  error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address the "
+        "problem. Control should not reach the end of this function",
+        __func__, PACKAGE_BUGREPORT);
+  return NAN;
+}
+
+
 
 
 
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index c06b2e1..ea996db 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -95,6 +95,10 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t 
ndim);
 float
 gal_dimension_dist_radial(size_t *a, size_t *b, size_t ndim);
 
+float
+gal_dimension_dist_elliptical(double *center, double *pa_deg, double *q,
+                              size_t ndim, double *point);
+
 
 
 /************************************************************************/