[gnuastro-commits] master 325d717 045/113: 3D matching now in Match program and library

[Mohammad Akhlaghi]

branch: master
commit 325d717c2ee231ed1a83eb5f2eb299175890ab36
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    3D matching now in Match program and library
    
    The Match program and library can now do a 3D match. To do this the new
    `gal_box_bound_ellipsoid_extent' has been defined in the libraries to make
    it easy to find the aperture extent.
---
 bin/match/args.h   |   2 +-
 bin/match/ui.c     | 141 +++++++++++++++++++++++++++++++++--
 doc/gnuastro.texi  |  95 ++++++++++++++++++------
 lib/box.c          |  49 ++++++++----
 lib/gnuastro/box.h |   4 +
 lib/match.c        | 214 ++++++++++++++++++++++++++++++++++++++++-------------
 6 files changed, 409 insertions(+), 96 deletions(-)

diff --git a/bin/match/args.h b/bin/match/args.h
index ab95b57..f2b46a1 100644
--- a/bin/match/args.h
+++ b/bin/match/args.h
@@ -104,7 +104,7 @@ struct argp_option program_options[] =
     {
       "aperture",
       UI_KEY_APERTURE,
-      "FLT[,FLT[,FLT]]",
+      "FLT[,...]",
       0,
       "Acceptable aperture for matching.",
       UI_GROUP_CATALOGMATCH,
diff --git a/bin/match/ui.c b/bin/match/ui.c
index b625545..eebb680 100644
--- a/bin/match/ui.c
+++ b/bin/match/ui.c
@@ -327,7 +327,7 @@ ui_read_columns_aperture_2d(struct matchparams *p)
 {
   size_t apersize=3;
   gal_data_t *newaper=NULL;
-  double *naper, *oaper=p->aperture->array;
+  double *naper=NULL, *oaper=p->aperture->array;
 
   /* A general sanity check: the first two elements of aperture cannot be
      zero or negative. */
@@ -339,7 +339,7 @@ ui_read_columns_aperture_2d(struct matchparams *p)
           "zero or negative");
 
   /* Will be needed in more than one case. */
-  if(p->aperture->size!=3)
+  if(p->aperture->size!=apersize)
     {
       newaper=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &apersize, NULL,
                              0, -1, NULL, NULL, NULL);
@@ -386,6 +386,136 @@ ui_read_columns_aperture_2d(struct matchparams *p)
 
 
 
+/* The 3D aperture must finally have these elements.
+
+      p->aperture[0] = Major axis length.
+      p->aperture[1] = Second axis ratio.
+      p->aperture[2] = Third axis ratio.
+      p->aperture[3] = First ZXZ Euler angle rotation.
+      p->aperture[4] = Second ZXZ Euler angle rotation.
+      p->aperture[5] = Third ZXZ Euler angle rotation.   */
+static void
+ui_read_columns_aperture_3d(struct matchparams *p)
+{
+  size_t apersize=6;
+  gal_data_t *newaper=NULL;
+  double *naper=NULL, *oaper=p->aperture->array;
+
+  /* A general sanity check: the first two elements of aperture cannot be
+     zero or negative. */
+  if( oaper[0]<=0 )
+    error(EXIT_FAILURE, 0, "the first value of `--aperture' cannot be "
+          "zero or negative");
+  if( p->aperture->size>2 && (oaper[1]<=0 || oaper[2]<=0) )
+    error(EXIT_FAILURE, 0, "the second and third values of `--aperture' "
+          "cannot be zero or negative");
+
+  /* Will be needed in more than one case. */
+  if(p->aperture->size!=apersize)
+    {
+      newaper=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &apersize, NULL,
+                             0, -1, NULL, NULL, NULL);
+      naper=newaper->array;
+    }
+
+  /* Different based on  */
+  switch(p->aperture->size)
+    {
+    case 1:
+      naper[0] = oaper[0];
+      naper[1] = naper[2] = 1;
+      naper[3] = naper[4] = naper[5] = 0;
+      break;
+
+    case 3:
+      /* Major axis is along the first dimension, no rotation necessary. */
+      if(oaper[0]>=oaper[1] && oaper[0]>=oaper[2])
+        {
+          naper[0] = oaper[0];
+          naper[1] = oaper[1] / oaper[0];
+          naper[2] = oaper[2] / oaper[0];
+          naper[3] = naper[4] = naper[5] = 0;
+        }
+
+      /* Major axis is along the second dimension. So we want `X' to be in
+         the direction of `y'. Therefore, just the first Eurler ZXZ
+         rotation is necessary by 90 degrees. Here is how the rotated
+         coordinates (X,Y,Z) look like (after one rotation about Z).
+
+                  |z (before)                   Z|  /Y
+                  |                              | /
+                  |                  ==>         |/
+                  .-------- y                    .------X
+                 /
+              x /
+
+         You see how the major axis (X) now lies in the second original
+         direction (y). The length along `x' is now along `Y' and the third
+         hasn't changed. Note that we are talking about the semi-axis
+         lengths (a scalar, not a vector), so +- is irrelevant. */
+      else if(oaper[1]>=oaper[0] && oaper[1]>=oaper[2])
+        {
+          naper[0] = oaper[1];
+          naper[1] = oaper[0] / oaper[1];
+          naper[2] = oaper[2] / oaper[1];
+          naper[3] = 90;
+          naper[4] = naper[5] = 0;
+        }
+
+      /* The major axis is along the third dimension. So we want `X' to
+         point in the direction of `z'. To get to that point, we need 90
+         degree rotations in all three Euler ZXZ rotations.
+
+              |z (before)     z'|  /y'            |y''                  |X
+              |                 | /               |                     |
+              |            ==>  |/        ==>     |         ==>         |
+              .------ y         .------x'         .------x'      Y------.
+             /                                   /                     /
+          x /                                z''/                    Z/
+
+         We thus see that the length along the second original direction
+         (y) doesn't change stays in the second direction (Y). But the
+         original first axis direction (x) is now represented by (Z). Note
+         that we are talking about the semi-axis lengths (a scalar, not a
+         vector), so +- is irrelevant. */
+      else
+        {
+          naper[0] = oaper[2];
+          naper[1] = oaper[1] / oaper[2];
+          naper[2] = oaper[0] / oaper[2];
+          naper[3] = naper[4] = naper[5] = 90;
+        }
+
+      break;
+
+    case 6:
+      if(oaper[1]>1 || oaper[2]>1)
+        error(EXIT_FAILURE, 0, "atleast one of the second or third values "
+              "to `--aperture' is larger than one. When size numbers are "
+              "given to this option (in 3D), the second and third are the "
+              "axis ratios (which must always be less than 1).");
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%zu values given to `--aperture'. In 3D, this "
+            "option can only take 1, 3, or 6 values. See the description of "
+            "this option in the book for more with this command:\n\n"
+            "    $ info astmatch\n", p->aperture->size);
+    }
+
+  /* If a new aperture was defined, then replace it with the exitsting
+     one. */
+  if(newaper)
+    {
+      gal_data_free(p->aperture);
+      p->aperture=newaper;
+    }
+}
+
+
+
+
+
 /* Read catalog columns */
 static void
 ui_read_columns(struct matchparams *p)
@@ -418,12 +548,9 @@ ui_read_columns(struct matchparams *p)
                 p->aperture->size);
         break;
 
-      case 2:
-        ui_read_columns_aperture_2d(p);
-        break;
-
+      case 2: ui_read_columns_aperture_2d(p); break;
+      case 3: ui_read_columns_aperture_3d(p); break;
       default:
-
         error(EXIT_FAILURE, 0, "%zu dimensional matches are not currently "
               "supported (maximum is 2 dimensions). The number of "
               "dimensions is deduced from the number of values given to "
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8fe764b..05c50b1 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15798,6 +15798,10 @@ $ astmatch --aperture=1/3600,2/3600 in1.fits in2.txt
 ## and DEC_D columns of the second within a 0.5 arcseconds aperture.
 $ astmatch --ccol1=RA,DEC --ccol2=RA_D,DEC_D --aperture0.5/3600  \
            in1.fits in2.fits
+
+## Match in 3D (RA, Dec and Wavelength).
+$ astmatch --ccol1=2,3,4 --ccol2=2,3,4 -a0.5/3600,0.5/3600,5e-10 \
+           in1.fits in2.txt
 @end example
 
 Two inputs are necessary for Match to start processing. The inputs can be
@@ -15866,8 +15870,8 @@ columns}. See the one-line examples above for some 
usages of this option.
 The coordinate columns of the second input. See the example in
 @option{--ccol1} for more.
 
-@item -a FLT[,FLT[,FLT]]
-@itemx --aperture=FLT[,FLT[,FLT]]
+@item -a FLT[,...]
+@itemx --aperture=FLT[,...]
 Parameters of the aperture for matching. The values given to this option
 can be fractions, for example when the position columns are in units of
 degrees, @option{1/3600} can be used to ask for one arcsecond. The
@@ -15875,15 +15879,23 @@ interpretation of the values depends on the requested 
dimensionality
 (determined from @option{--ccol1} and @code{--ccol2}) and how many values
 are given to this option.
 
+When multiple objects are found within the aperture, the match is defined
+as the nearest one. In a multi-dimensional dataset, when the aperture is a
+general ellipse or ellipsoid (and not a circle or sphere), the distance is
+calculated in the elliptical space along the major axis. For the defintion
+of this distance, see @mymath{r_{el}} in @ref{Defining an ellipse and
+ellipsoid}.
+
 @table @asis
-@item 1D match
+@item 1D
 The aperture/interval can only take one value: half of the interval around
 each point (maximum distance from each point).
 
-@item 2D match
-In a 2D match, the aperture can be a circle, an ellipse aligned in the axes
-or an ellipse with a rotated major axis. To simply the usage, you can
-determine the shape based on the number of free parameters for each.
+@item 2D
+The aperture can be a circle, an ellipse aligned in the axes or an ellipse
+with a rotated major axis. To simplify the usage, the shape can be
+determined based on the number of values given to this option.
+
 @table @asis
 @item 1 number
 For example @option{--aperture=2}. The aperture will be a circle of the
@@ -15891,24 +15903,48 @@ given radius. The value will be in the same units as 
the columns in
 @option{--ccol1} and @option{--ccol2}).
 
 @item 2 numbers
-For example @option{--aperture=3,4e-10}. The aperture will be an ellipse
-(if the two numbers are different) with the respective value along each
-dimension. The numbers are in units of the first and second axis. In the
-example above, the semi-axis value along the first axis will be 3 (in units
-of the first coordinate) and along the second axis will be
-@mymath{4\times10^{-10}} (in units of the second coordinate). Such values
-can happen if you are comparing catalogs of a spectra for example. If more
-than one object exists in the aperture, the nearest will be found along the
-major axis as described in @ref{Defining an ellipse and ellipsoid}.
+For example @option{--aperture=3,4e-10}. The aperture will be a general
+ellipse with the respective value along each dimension. The numbers are in
+units of the first and second axis. In the example above, the semi-axis
+value along the first axis will be 3 (in units of the first coordinate) and
+along the second axis will be @mymath{4\times10^{-10}} (in units of the
+second coordinate). Such values can happen if you are comparing catalogs of
+a spectra for example.
 
 @item 3 numbers
 For example @option{--aperture=2,0.6,30}. The aperture will be an ellipse
 (if the second value is not 1). The first number is the semi-major axis,
 the second is the axis ratio and the third is the position angle (in
-degrees). If multiple matches are found within the ellipse, the distance
-(to find the nearest) is calculated along the major axis in the elliptical
-space, see @ref{Defining an ellipse and ellipsoid}.
+degrees).
+@end table
+
+@item 3D
+The aperture (matching volume) can be a sphere, an ellipsoid aligned on the
+three axises or a genenral ellipsoid rotated in any direction. To simplify
+the usage, the shape can be determined based on the number of values given
+to this option.
+
+@table @asis
+@item 1 number
+For example @option{--aperture=3}. The matching volume will be a sphere of
+the given radius. The value is in the same units as the input coordinates.
+
+@item 3 numbers
+For example @option{--aperture=4,5,6e-10}. The aperture will be a general
+ellipsoid with the respective extent along each dimension. The numbers must
+be in the same units as each axis. This is very similar to the two number
+case of 2D inputs. See there for more.
+
+@item 6 numbers
+For example @option{--aperture=4,0.5,0.6,10,20,30}. The numbers represent
+the full general ellipsoid definition (in any orientation). For the
+definition of a general ellipsoid, see @ref{Defining an ellipse and
+ellipsoid}. The first number is the semi-major axis. The second and third
+are the two axis ratios. The last three are the three Euler angles in units
+of degrees in the ZXZ order as fully described in @ref{Defining an ellipse
+and ellipsoid}.
 @end table
+
 @end table
 @end table
 
@@ -23177,9 +23213,9 @@ must already be allocated before this function.
 
 @deftypefun void gal_box_bound_ellipse (double @code{a}, double @code{b}, 
double @code{theta_deg}, long @code{*width})
 Any ellipse can be enclosed into a rectangular box. The purpose of this
-function is to give the height and width of that box assuming the center of
-the ellipse is in the box center. @code{a} is the ellipse major axis,
-@code{b} is the minor axis, @code{theta_deg} is the position angle in
+function is to give the integer height and width of that box assuming the
+center of the ellipse is in the box center. @code{a} is the ellipse major
+axis, @code{b} is the minor axis, @code{theta_deg} is the position angle in
 degrees. The @code{width} array will contain the output size in long
 integer type. @code{width[0]}, and @code{width[1]} are the number of pixels
 along the first and second FITS axis. Since the ellipse center is assumed
@@ -23187,10 +23223,21 @@ to be in the center of the box, all the values in 
@code{width} will be an
 odd integer.
 @end deftypefun
 
+@deftypefun void gal_box_bound_ellipsoid_extent (double @code{*semiaxes}, 
double @code{*euler_deg}, double @code{*extent})
+Return the maximum extent along each dimension of the given ellipsoid from
+its center. Therefore this is half the extent of the box in each
+dimension. The semi-axis lengths of the ellipsoid must be present in the 3
+element @code{semiaxis} array. The @code{euler_deg} array contains the
+three ellipsoid Euler angles in degrees. For a description of the Euler
+angles, see description of @code{gal_box_bound_ellipsoid} below. The extent
+in each dimension is in floating point format and stored in @code{extent}
+which must already be allocated before this function.
+@end deftypefun
+
 @deftypefun void gal_box_bound_ellipsoid (double @code{*semiaxes}, double 
@code{*euler_deg}, long @code{*width})
 Any ellipsoid can be enclosed into a rectangular volume/box. The purpose of
-this function is to give the size/width of that box. The semi-axes lengths
-of the ellipse must be in the @code{semiaxes} array (with three
+this function is to give the integer size/width of that box. The semi-axes
+lengths of the ellipse must be in the @code{semiaxes} array (with three
 elements). The major axis length must be the first element of
 @code{semiaxes}. The only other condition is that the next two semi-axes
 must both be smaller than the first. The orientation of the major axis is
diff --git a/lib/box.c b/lib/box.c
index 6c80d3c..07080b2 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -92,8 +92,8 @@ gal_box_bound_ellipse(double a, double b, double theta_deg, 
long *width)
      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 * ( (size_t)extent[0] + 1 ) + 1;
-  width[1] = 2 * ( (size_t)extent[1] + 1 ) + 1;
+  width[0] = 2 * ( (long)extent[0] + 1 ) + 1;
+  width[1] = 2 * ( (long)extent[1] + 1 ) + 1;
 }
 
 
@@ -198,17 +198,12 @@ gal_box_bound_ellipse(double a, double b, double 
theta_deg, long *width)
       x = M[1,4] \pm sqrt( M[1,1]^2 + M[1,2]^2 + M[1,3]^2 )
       y = M[2,4] \pm sqrt( M[2,1]^2 + M[2,2]^2 + M[2,3]^2 )
       z = M[3,4] \pm sqrt( M[3,1]^2 + M[3,2]^2 + M[3,3]^2 )        */
-#define PRINT3BY3(C, A){                                                \
-    printf("%s: | %-15g%-15g%-15g |\n"                                  \
-           "   | %-15g%-15g%-15g |\n"                                   \
-           "   | %-15g%-15g%-15g |\n\n", (C), (A)[0], (A)[1], (A)[2],   \
-           (A)[3], (A)[4], (A)[5], (A)[6], (A)[7], (A)[8]);             \
-  }
 void
-gal_box_bound_ellipsoid(double *semiaxes, double *euler_deg, long *width)
+gal_box_bound_ellipsoid_extent(double *semiaxes, double *euler_deg,
+                               double *extent)
 {
   size_t i, j;
-  double sumsq, a=semiaxes[0], b=semiaxes[1], c=semiaxes[2];
+  double a=semiaxes[0], b=semiaxes[1], c=semiaxes[2];
   double c1=cos(euler_deg[0]*M_PI/180), s1=sin(euler_deg[0]*M_PI/180);
   double c2=cos(euler_deg[1]*M_PI/180), s2=sin(euler_deg[1]*M_PI/180);
   double c3=cos(euler_deg[2]*M_PI/180), s3=sin(euler_deg[2]*M_PI/180);
@@ -225,8 +220,10 @@ gal_box_bound_ellipsoid(double *semiaxes, double 
*euler_deg, long *width)
   /* Find the width along each dimension. */
   for(i=0;i<3;++i)
     {
-      sumsq=0.0f; for(j=0;j<3;++j) sumsq += R[i*3+j] * R[i*3+j];
-      width[i] = 2 * ((long)sqrt(sumsq)+1) + 1;
+      extent[i]=0.0;
+      for(j=0;j<3;++j)
+        extent[i] += R[i*3+j] * R[i*3+j];
+      extent[i] = sqrt(extent[i]);
     }
 
   /* For a check:
@@ -236,7 +233,7 @@ gal_box_bound_ellipsoid(double *semiaxes, double 
*euler_deg, long *width)
   printf("c1: %f, s1: %f\nc2: %f, s2: %f\nc3: %f, s3: %f\n",
          c1, s1, c2, s2, c3, s3);
   PRINT3BY3("R", R);
-  printf("width: %ld, %ld, %ld\n", width[0], width[1], width[2]);
+  printf("extent: %ld, %ld, %ld\n", extent[0], extent[1], extent[2]);
   exit(0);
   */
 }
@@ -245,6 +242,32 @@ gal_box_bound_ellipsoid(double *semiaxes, double 
*euler_deg, long *width)
 
 
 
+/* Using `gal_box_bound_ellipsoid_extent', find the integer width of a box
+   that contains the ellipsoid. */
+#define PRINT3BY3(C, A){                                                \
+    printf("%s: | %-15g%-15g%-15g |\n"                                  \
+           "   | %-15g%-15g%-15g |\n"                                   \
+           "   | %-15g%-15g%-15g |\n\n", (C), (A)[0], (A)[1], (A)[2],   \
+           (A)[3], (A)[4], (A)[5], (A)[6], (A)[7], (A)[8]);             \
+  }
+void
+gal_box_bound_ellipsoid(double *semiaxes, double *euler_deg, long *width)
+{
+  size_t i;
+  double extent[3];
+
+  /* Find the extent of the ellipsoid in each axis. */
+  gal_box_bound_ellipsoid_extent(semiaxes, euler_deg, extent);
+
+  /* Find the width along each dimension. */
+  for(i=0;i<3;++i)
+    width[i] = 2 * ((long)extent[i] + 1) + 1;
+}
+
+
+
+
+
 /* We have the central pixel and box size of the crop box, find the
    starting (first) and ending (last) pixels: */
 void
diff --git a/lib/gnuastro/box.h b/lib/gnuastro/box.h
index 78d9710..72f48a6 100644
--- a/lib/gnuastro/box.h
+++ b/lib/gnuastro/box.h
@@ -59,6 +59,10 @@ void
 gal_box_bound_ellipse(double a, double b, double theta_deg, long *width);
 
 void
+gal_box_bound_ellipsoid_extent(double *semiaxes, double *euler_deg,
+                              double *extent);
+
+void
 gal_box_bound_ellipsoid(double *semiaxes, double *euler_deg, long *width);
 
 void
diff --git a/lib/match.c b/lib/match.c
index 17de927..2c7467b 100644
--- a/lib/match.c
+++ b/lib/match.c
@@ -158,10 +158,10 @@ match_coordinaes_sanity_check(gal_data_t *coord1, 
gal_data_t *coord2,
           gal_list_data_number(coord2));
 
 
-  /* This function currently only works for 2 dimensions. */
-  if(ncoord1>2)
+  /* This function currently only works for less than 4 dimensions. */
+  if(ncoord1>3)
     error(EXIT_FAILURE, 0, "%s: %zu dimension matching requested, this "
-          "function currently only matches datasets with a maximum of 2 "
+          "function currently only matches datasets with a maximum of 3 "
           "dimensions", __func__, ncoord1);
 
   /* Check the column properties. */
@@ -172,11 +172,32 @@ match_coordinaes_sanity_check(gal_data_t *coord1, 
gal_data_t *coord2,
   if(aperture[0]<=0)
     error(EXIT_FAILURE, 0, "%s: the first value in the aperture (%g) cannot "
           "be zero or negative", __func__, aperture[0]);
-  if(ncoord1==2)
-    if(aperture[1]<=0 || aperture[1]>1)
-      error(EXIT_FAILURE, 0, "%s: the second value in the aperture (%g) "
-            "is the axis ratio, so it must be larger than zero and less "
-            "than 1", __func__, aperture[1]);
+  switch(ncoord1)
+    {
+    case 1:  /* We don't need any checks in a 1D match. */
+      break;
+
+    case 2:
+      if( (aperture[1]<=0 || aperture[1]>1))
+        error(EXIT_FAILURE, 0, "%s: the second value in the aperture (%g) "
+              "is the axis ratio, so it must be larger than zero and less "
+              "than 1", __func__, aperture[1]);
+      break;
+
+    case 3:
+      if(aperture[1]<=0 || aperture[1]>1 || aperture[2]<=0 || aperture[2]>1)
+        error(EXIT_FAILURE, 0, "%s: at least one of the second or third "
+              "values in the aperture (%g and %g respectively) is smaller "
+              "than zero or larger than one. In a 3D match, these are the "
+              "axis ratios, so they must be larger than zero and less than "
+              "1", __func__, aperture[1], aperture[2]);
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+            "the issue. The value %zu not recognized for `ndim'", __func__,
+            PACKAGE_BUGREPORT, ncoord1);
+    }
 }
 
 
@@ -288,6 +309,89 @@ match_coordinates_prepare(gal_data_t *coord1, gal_data_t 
*coord2,
 /********************************************************************/
 /*************            Coordinate matching           *************/
 /********************************************************************/
+/* Preparations for `match_coordinates_second_in_first'. */
+static void
+match_coordinates_sif_prepare(gal_data_t *A, gal_data_t *B,
+                              double *aperture, size_t ndim, double **a,
+                              double **b, double *dist, double *c,
+                              double *s, int *iscircle)
+{
+  double semiaxes[3];
+
+  /* These two are common for all dimensions. */
+  a[0]=A->array;
+  b[0]=B->array;
+
+  /* See if the aperture is a circle or not. */
+  switch(ndim)
+    {
+    case 1:
+      *iscircle = 0; /* Irrelevant for 1D. */
+      dist[0]=aperture[0];
+      break;
+
+    case 2:
+      /* Set the main coordinate arrays. */
+      a[1]=A->next->array;
+      b[1]=B->next->array;
+
+      /* See if the aperture is circular. */
+      if( (*iscircle=(aperture[1]==1)?1:0)==0 )
+        {
+          /* Using the box that encloses the aperture, calculate the
+             distance along each axis. */
+          gal_box_bound_ellipse_extent(aperture[0], aperture[0]*aperture[1],
+                                       aperture[2], dist);
+
+          /* Calculate the sin and cos of the given ellipse if necessary
+             for ease of processing later. */
+          c[0] = cos( aperture[2] * M_PI/180.0 );
+          s[0] = sin( aperture[2] * M_PI/180.0 );
+        }
+      else
+        dist[0]=dist[1]=aperture[0];
+      break;
+
+    case 3:
+      /* Set the main coordinate arrays. */
+      a[1]=A->next->array;
+      b[1]=B->next->array;
+      a[2]=A->next->next->array;
+      b[2]=B->next->next->array;
+
+      if( (*iscircle=(aperture[1]==1 && aperture[2]==1)?1:0)==0 )
+        {
+          /* Using the box that encloses the aperture, calculate the
+             distance along each axis. */
+          semiaxes[0]=aperture[0];
+          semiaxes[1]=aperture[1]*aperture[0];
+          semiaxes[2]=aperture[2]*aperture[0];
+          gal_box_bound_ellipsoid_extent(semiaxes, &aperture[3], dist);
+
+          /* Calculate the sin and cos of the given ellipse if necessary
+             for ease of processing later. */
+          c[0] = cos( aperture[3] * M_PI/180.0 );
+          s[0] = sin( aperture[3] * M_PI/180.0 );
+          c[1] = cos( aperture[4] * M_PI/180.0 );
+          s[1] = sin( aperture[4] * M_PI/180.0 );
+          c[2] = cos( aperture[5] * M_PI/180.0 );
+          s[2] = sin( aperture[5] * M_PI/180.0 );
+        }
+      else
+        dist[0]=dist[1]=dist[2]=aperture[0];
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+            "the problem. The value %zu is not recognized for ndim",
+            __func__, PACKAGE_BUGREPORT, ndim);
+    }
+}
+
+
+
+
+
 static double
 match_coordinates_elliptical_r_2d(double d1, double d2, double *ellipse,
                                   double c, double s)
@@ -302,22 +406,50 @@ match_coordinates_elliptical_r_2d(double d1, double d2, 
double *ellipse,
 
 
 static double
-match_coordinates_distance(double *delta, int iscircle, size_t ndim,
-                           double *aperture, double c, double s)
+match_coordinates_elliptical_r_3d(double *delta, double *ellipsoid,
+                                  double *c, double *s)
 {
-  /* In a one dimensional case, the distance is just the absolute value of
-     delta[0]. */
-  if(ndim==1) return fabs(delta[0]);
+  double Xr, Yr, Zr;
+  double c1=c[0], s1=s[0];
+  double c2=c[1], s2=s[1];
+  double c3=c[2], s3=s[2];
+  double q1=ellipsoid[1], q2=ellipsoid[2];
+  double x=delta[0], y=delta[1], z=delta[2];
+
+  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 );
+}
+
 
+
+
+
+static double
+match_coordinates_distance(double *delta, int iscircle, size_t ndim,
+                           double *aperture, double *c, double *s)
+{
   /* For more than one dimension, we'll need to calculate the distance from
      the deltas (differences) in each dimension. */
   switch(ndim)
     {
+    case 1:
+      return fabs(delta[0]);
+
     case 2:
       return ( iscircle
                ? sqrt( delta[0]*delta[0] + delta[1]*delta[1] )
                : match_coordinates_elliptical_r_2d(delta[0], delta[1],
-                                                   aperture, c, s) );
+                                                   aperture, c[0], s[0]) );
+
+    case 3:
+      return ( iscircle
+               ? sqrt( delta[0]*delta[0]
+                       + delta[1]*delta[1]
+                       + delta[2]*delta[2] )
+               : match_coordinates_elliptical_r_3d(delta, aperture, c, s) );
+
     default:
       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
             "the problem. The value %zu is not recognized for ndim",
@@ -348,41 +480,16 @@ match_coordinates_second_in_first(gal_data_t *A, 
gal_data_t *B,
      redundant variables (those that equal a previous value) are only
      defined to make it easy to read the code.*/
   int iscircle=0;
-  double delta[2];
-  double r, c=NAN, s=NAN, dist[2];
   size_t ai, bi, blow, prevblow=0;
   size_t i, ar=A->size, br=B->size;
   size_t ndim=gal_list_data_number(A);
-  double *a[2]={A->array, A->next ? A->next->array : NULL};
-  double *b[2]={B->array, B->next ? B->next->array : NULL};
-
-  /* See if the aperture is a circle or not. */
-  switch(ndim)
-    {
-    case 1: iscircle = 0; /* Irrelevant for 1D. */  break;
-    case 2: iscircle = aperture[1]==1 ? 1 : 0;      break;
-    default:
-      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
-            "the problem. The value %zu is not recognized for ndim",
-            __func__, PACKAGE_BUGREPORT, ndim);
-    }
-
-  /* Preparations for the shape of the aperture. */
-  if(ndim==1 || iscircle)
-    dist[0]=dist[1]=aperture[0];
-  else
-    {
-      /* Using the box that encloses the aperture, calculate the distance
-         along each axis. */
-      gal_box_bound_ellipse_extent(aperture[0], aperture[0]*aperture[1],
-                                   aperture[2], dist);
-
-      /* Calculate the sin and cos of the given ellipse if necessary for
-         ease of processing later. */
-      c = cos( aperture[2] * M_PI/180.0 );
-      s = sin( aperture[2] * M_PI/180.0 );
-    }
+  double r, c[3]={NAN, NAN, NAN}, s[3]={NAN, NAN, NAN};
+  double dist[3]={NAN, NAN, NAN}, delta[3]={NAN, NAN, NAN};
+  double *a[3]={NULL, NULL, NULL}, *b[3]={NULL, NULL, NULL};
 
+  /* Necessary preperations. */
+  match_coordinates_sif_prepare(A, B, aperture,  ndim, a, b, dist, c, s,
+                                &iscircle);
 
   /* For each row/record of catalog `a', make a list of the nearest records
      in catalog b within the maximum distance. Note that both catalogs are
@@ -430,7 +537,7 @@ match_coordinates_second_in_first(gal_data_t *A, gal_data_t 
*B,
             axis coordinates have EXACTLY the same floating point
             value as each other to easily define an independent
             sorting in the second axis. */
-         if( ndim==1
+         if( ndim<2
               || ( b[1][bi] >= a[1][ai]-dist[1]
                    && b[1][bi] <= a[1][ai]+dist[1] ) )
            {
@@ -462,11 +569,16 @@ match_coordinates_second_in_first(gal_data_t *A, 
gal_data_t *B,
 
                 The next two problems will be solved with the list after
                 parsing of the whole catalog is complete.*/
-              for(i=0;i<ndim;++i) delta[i]=b[i][bi]-a[i][ai];
-              r=match_coordinates_distance(delta, iscircle, ndim, aperture,
-                                           c, s);
-             if(r<aperture[0])
-               match_coordinate_add_to_sfll(&bina[ai], bi, r);
+              if( ndim<3
+                  || ( b[2][bi] >= a[2][ai]-dist[2]
+                       && b[2][bi] <= a[2][ai]+dist[2] ) )
+                {
+                  for(i=0;i<ndim;++i) delta[i]=b[i][bi]-a[i][ai];
+                  r=match_coordinates_distance(delta, iscircle, ndim,
+                                               aperture, c, s);
+                  if(r<aperture[0])
+                    match_coordinate_add_to_sfll(&bina[ai], bi, r);
+                }
            }
        }