Fix crash in convex_hull when input includes repeated points (spacete…

…lescope#254)

* Fix crash in convex_hull when input includes repeated points

* add and update unit tests

CHANGES.rst

@@ -29,6 +29,13 @@ Release Notes
 - Replaced private ``numpy.core._umath_tests.inner1d()`` with an alternative
   implementation. [#253]
 
+- Fixed a bug in the ``math_util.angles()`` function that would result in crash
+  when one of the input vectors is a null vector. [#254]
+
+- Enhanced the code in the ``SphericalPolygon.convex_hull()`` to ignore points
+  that lead to null vectors in computations. [#254]
+
+
 1.2.23 (10-October-2022)
 ========================
 

spherical_geometry/polygon.py

@@ -348,7 +348,7 @@ def convex_hull(cls, points):
 
         duo = list(zip(pt, ang))
         duo = sorted(duo, key=lambda d: d[1])
-        points = np.asarray([d[0] for d in duo])
+        points = np.asarray([d[0] for d in duo if np.isfinite(d[1])])
 
         # Set the first point on the hull to the extreme point
 

spherical_geometry/tests/test_basic.py

@@ -477,12 +477,19 @@ def test_fast_area():
     assert carea > 2.0 * np.pi and carea < 4.0 * np.pi
 
 
-def test_convex_hull():
+@pytest.mark.parametrize(
+    'repeat_pts', [False, True]
+)
+def test_convex_hull(repeat_pts):
     lon = (0.02, 0.10, 0.05, 0.03, 0.04, 0.07, 0.00, 0.06, 0.08, 0.13,
            0.08, 0.14, 0.15, 0.12, 0.01, 0.11)
     lat = (0.06, 0.00, 0.05, 0.01, 0.12, 0.08, 0.03, 0.02, 0.04, 0.03,
            0.10, 0.11, 0.01, 0.13, 0.09, 0.07)
 
+    if repeat_pts:
+        lon = lon + lon[::-1]
+        lat = lat + lat[::-1]
+
     lon_lat = list(zip(lon, lat))
 
     points = []
@@ -518,8 +525,9 @@ def test_convex_hull():
 @pytest.mark.skipif(math_util is None, reason="math_util C-ext is missing")
 def test_math_util_angle_domain():
     # Before a fix, this would segfault
-    with pytest.raises(ValueError):
-        math_util.angle([[0, 0, 0]], [[0, 0, 0]], [[0, 0, 0]])
+    assert not np.isfinite(
+        math_util.angle([[0, 0, 0]], [[0, 0, 0]], [[0, 0, 0]])[0]
+    )
 
 
 @pytest.mark.skipif(math_util is None, reason="math_util C-ext is missing")

src/math_util.c

@@ -7,6 +7,8 @@
 #include "numpy/npy_math.h"
 
 #include "qd/c_qd.h"
+#include <math.h>
+#include <stdlib.h>
 
 
 /*
@@ -198,7 +200,7 @@ normalized_dot_qd(const qd *A, const qd *B, qd *dot_val) {
         /* return non-normalized value: */
         PyErr_SetString(PyExc_ValueError, "Null vector.");
         c_qd_copy(ab.x, dot_val->x);
-        return 1;
+        return 2;
     } else {
         c_qd_div(ab.x, norm, dot_val->x);
     }
@@ -707,7 +709,7 @@ char *angle_signature = "(i),(i),(i)->()";
 static void
 DOUBLE_angle(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
 {
-    int comp;
+    int comp, ret;
     qd A[3];
     qd B[3];
     qd C[3];
@@ -741,7 +743,18 @@ DOUBLE_angle(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func))
         cross_qd(C, B, BCX);
         cross_qd(ABX, BCX, X);
         dot_qd(B, X, &diff);
-        if (normalized_dot_qd(ABX, BCX, &inner)) return;
+        ret = normalized_dot_qd(ABX, BCX, &inner);
+        if (ret == 1)  {
+            return;
+        } else if (ret == 2) {
+            PyErr_Clear();
+#if defined(NAN)
+            *((double *)op) = NAN;
+#else
+            *((double *)op) = strtod("NaN", NULL);
+#endif
+            continue;
+        }
 
         c_qd_abs(inner.x, abs_inner);
         c_qd_comp(abs_inner, QD_ONE, &comp);