DM-32008: Wrap test_spoke in C++/pybind11

include/lsst/meas/astrom/pessimisticPatternMatcherUtils.h

@@ -0,0 +1,71 @@
+// -*- LSST-C++ -*-
+
+/*
+ * This file is part of meas_astrom.
+ *
+ * Developed for the LSST Data Management System.
+ * This product includes software developed by the LSST Project
+ * (https://www.lsst.org).
+ * See the COPYRIGHT file at the top-level directory of this distribution
+ * for details of code ownership.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <ndarray.h>
+
+namespace lsst {
+namespace meas {
+namespace astrom {
+
+/**
+ * Check the opening angle between the first spoke of our pattern for the
+ * source object against the reference object.
+ *
+ * This method makes heavy use of the small angle approximation to perform
+ * the comparison.
+ *
+ * @param[in] cos_theta_src Cosine of the angle between the current candidate
+ *     source spoke and the first spoke.
+ * @param[in] sin_theta_src Sine of the angle between the current candidate
+ *     source spoke and the first spoke.
+ * @param[in] ref_ctr 3 vector of the candidate reference center
+ * @param[in] ref_ctr_id id lookup of the ref_ctr into the master reference
+ *     array. uint16 type is to limit the amount of memory used and is set
+ *     in the pessimistic_pattern_matcher_3d python class with reference
+ *     catalogs trimmed by the matchPessimisticB runner class.
+ * @param[in] proj_ref_ctr_delta Plane projected first spoke in the reference
+ *     pattern using the pattern center as normal.
+ * @param[in] proj_ref_ctr_dist_sq Squared length of the projected vector.
+ * @param[in] ref_dist_idx_array Indices sorted by the delta distance between
+ *     the source spoke we are trying to test and the candidate reference
+ *     spokes.
+ * @param[in] ref_id_array Array of id lookups into the master reference array
+ *     that our center id object is paired with.
+ * @param[in] reference_array Array of three vectors representing the locations
+ *     of all reference objects.
+ * @param[in] src_sin_tol Sine of tolerance allowed between source and
+ *     reference spoke opening angles.
+ * @return ID of the candidate reference object successfully matched or -1 if
+ *     no match is found.
+ */
+int check_spoke(double cos_theta_src, double sin_theta_src, ndarray::Array<double, 1, 1> const& ref_ctr,
+                ndarray::Array<double, 1, 1> const& proj_ref_ctr_delta, double proj_ref_ctr_dist_sq,
+                ndarray::Array<long int, 1, 1> const& ref_dist_idx_array,
+                ndarray::Array<uint16_t, 1, 1> const& ref_id_array,
+                ndarray::Array<double, 2, 1> const& reference_array, double src_sin_tol);
+
+}  // namespace astrom
+}  // namespace meas
+}  // namespace lsst

python/lsst/meas/astrom/SConscript

@@ -6,4 +6,5 @@ scripts.BasicSConscript.pybind11([
     "polynomialTransform",
     "scaledPolynomialTransformFitter",
     "sipTransform",
+    "pessimisticPatternMatcherUtils",
 ], addUnderscore=False)

python/lsst/meas/astrom/pessimisticPatternMatcherUtils.cc

@@ -0,0 +1,46 @@
+// -*- LSST-C++ -*-
+
+/*
+ * This file is part of meas_astrom.
+ *
+ * Developed for the LSST Data Management System.
+ * This product includes software developed by the LSST Project
+ * (https://www.lsst.org).
+ * See the COPYRIGHT file at the top-level directory of this distribution
+ * for details of code ownership.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
+#include "numpy/arrayobject.h"
+#include "ndarray/pybind11.h"
+
+#include "lsst/meas/astrom/pessimisticPatternMatcherUtils.h"
+
+namespace py = pybind11;
+
+namespace lsst {
+namespace meas {
+namespace astrom {
+
+PYBIND11_MODULE(pessimisticPatternMatcherUtils, mod) {
+    mod.def("check_spoke", &check_spoke,
+            "Test the opening angle between the spokes of our source pattern against the reference.");
+}
+
+}  // namespace astrom
+}  // namespace meas
+}  // namespace lsst

python/lsst/meas/astrom/pessimistic_pattern_matcher_b_3D.py

@@ -4,6 +4,7 @@
 from scipy.spatial import cKDTree
 from scipy.stats import sigmaclip
 
+from .pessimisticPatternMatcherUtils import check_spoke
 import lsst.pipe.base as pipeBase
 
 
@@ -605,7 +606,7 @@ def _find_candidate_reference_pairs(self, src_dist, ref_dist_array,
 
         # If these are equal there are no candidates and we exit.
         if start_idx == end_idx:
-            return []
+            return np.array([], dtype="int")
 
         # Make sure the endpoints of the input array are respected.
         if start_idx < 0:
@@ -839,17 +840,17 @@ def _create_pattern_spokes(self, src_ctr, src_delta_array, src_dist_array,
 
             # Test the spokes and return the id of the reference object.
             # Return None if no match is found.
-            ref_id = self._test_spoke(
+            ref_id = check_spoke(
                 cos_theta_src,
                 sin_theta_src,
                 ref_ctr,
-                ref_ctr_id,
                 proj_ref_ctr_delta,
                 proj_ref_ctr_dist_sq,
                 ref_dist_idx_array,
                 ref_id_array,
+                self._reference_array,
                 src_sin_tol)
-            if ref_id is None:
+            if ref_id < 0:
                 n_fail += 1
                 continue
 
@@ -953,14 +954,13 @@ def _test_spoke(self, cos_theta_src, sin_theta_src, ref_ctr, ref_ctr_id,
             else:
                 sin_comparison = \
                     (sin_theta_src - sin_theta_ref) / cos_theta_ref
-
             if abs(sin_comparison) > src_sin_tol:
                 continue
 
             # Return the correct id of the candidate we found.
             return ref_id_array[ref_dist_idx]
 
-        return None
+        return -1
 
     def _create_shift_rot_matrix(self, cos_rot_sq, shift_matrix, src_delta,
                                  ref_ctr, ref_delta):

src/pessimisticPatternMatcherUtils.cc

@@ -0,0 +1,101 @@
+// -*- LSST-C++ -*-
+
+/*
+ * This file is part of meas_astrom.
+ *
+ * Developed for the LSST Data Management System.
+ * This product includes software developed by the LSST Project
+ * (https://www.lsst.org).
+ * See the COPYRIGHT file at the top-level directory of this distribution
+ * for details of code ownership.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <cmath>
+#include <Eigen/Dense>
+#include "ndarray/eigen.h"
+#include "lsst/meas/astrom/pessimisticPatternMatcherUtils.h"
+
+namespace lsst {
+namespace meas {
+namespace astrom {
+
+int check_spoke(double cos_theta_src, double sin_theta_src, ndarray::Array<double, 1, 1> const& ref_ctr,
+                ndarray::Array<double, 1, 1> const& proj_ref_ctr_delta, double proj_ref_ctr_dist_sq,
+                ndarray::Array<long int, 1, 1> const& ref_dist_idx_array,
+                ndarray::Array<uint16_t, 1, 1> const& ref_id_array,
+                ndarray::Array<double, 2, 1> const& reference_array, double src_sin_tol) {
+    // Loop over our candidate reference objects. ref_dist_idx_array is a view into ref_id_array
+    // and is not the same length as ref_id_array.
+    for (auto idx = ref_dist_idx_array.begin(); idx != ref_dist_idx_array.end(); idx++) {
+        // Compute the delta vector from the pattern center.
+        unsigned int ref_id = ref_id_array[*idx];
+        ndarray::Array<double, 1, 1> ref_delta = copy(reference_array[ref_id] - ref_ctr);
+
+        double ref_dot = ndarray::asEigenMatrix(ref_delta).dot(ndarray::asEigenMatrix(ref_ctr));
+        ndarray::Array<double, 1, 1> proj_ref_delta = copy(ref_delta - ref_dot * ref_ctr);
+        // Compute the cos between our "center" reference vector and the
+        // current reference candidate.
+        double proj_delta_dist_sq =
+                ndarray::asEigenMatrix(proj_ref_delta).dot(ndarray::asEigenMatrix(proj_ref_delta));
+        double geom_dist_ref = sqrt(proj_ref_ctr_dist_sq * proj_delta_dist_sq);
+        double cos_theta_ref =
+                ndarray::asEigenMatrix(proj_ref_delta).dot(ndarray::asEigenMatrix(proj_ref_ctr_delta)) /
+                geom_dist_ref;
+
+        // Make sure we can safely make the comparison in case
+        // our "center" and candidate vectors are mostly aligned.
+        double cos_sq_comparison;
+        if (cos_theta_ref * cos_theta_ref < (1 - src_sin_tol * src_sin_tol)) {
+            cos_sq_comparison = (cos_theta_src - cos_theta_ref) * (cos_theta_src - cos_theta_ref) /
+                                (1 - cos_theta_ref * cos_theta_ref);
+        } else {
+            cos_sq_comparison = (cos_theta_src - cos_theta_ref) * (cos_theta_src - cos_theta_ref) /
+                                (src_sin_tol * src_sin_tol);
+        }
+        // Test the difference of the cosine of the reference angle against
+        // the source angle. Assumes that the delta between the two is
+        // small.
+        if (cos_sq_comparison > src_sin_tol * src_sin_tol) {
+            continue;
+        }
+        // The cosine tests the magnitude of the angle but not
+        // its direction. To do that we need to know the sine as well.
+        // This cross product calculation does that.
+        Eigen::Vector3d cross_ref = ndarray::asEigenMatrix(proj_ref_delta)
+                                            .head<3>()
+                                            .cross(ndarray::asEigenMatrix(proj_ref_ctr_delta).head<3>()) /
+                                    geom_dist_ref;
+        double sin_theta_ref = cross_ref.dot(ndarray::asEigenMatrix(ref_ctr));
+        // Check the value of the cos again to make sure that it is not
+        // near zero.
+        double sin_comparison;
+        if (abs(cos_theta_src) < src_sin_tol) {
+            sin_comparison = (sin_theta_src - sin_theta_ref) / src_sin_tol;
+        } else {
+            sin_comparison = (sin_theta_src - sin_theta_ref) / cos_theta_ref;
+        }
+
+        // Return the correct id of the candidate we found.
+        if (abs(sin_comparison) < src_sin_tol) {
+            return ref_id;
+        }
+    }
+    return -1;
+}
+
+}  // namespace astrom
+}  // namespace meas
+}  // namespace lsst