Refactor cross product calculations to use numpy.linalg.det for impro…

…ved accuracy and to address deprecation of 2D cross product

romancal/skymatch/region.py

@@ -400,25 +400,37 @@ def intersection(self, edge):
         v = edge._stop - edge._start
         w = self._start - edge._start
 
-        # Make problem a proper 3D problem, since these are 2D vectors in the
-        # xy plane we expect the cross product to be in the z direction only.
-        u = np.insert(u, 2, 0)
-        v = np.insert(v, 2, 0)
-        w = np.insert(w, 2, 0)
-
-        # We only care about the z component of the cross product
-        D = np.cross(u, v)[-1]
+        # Find the determinant of the matrix formed by the vectors u and v
+        #    Note: Originally this was computed using a numpy "2D" cross product,
+        #          however, this functionality has been deprecated and slated for
+        #          removal.
+        #    Note: np.linalg.det([u, v]) can also be used however, it changes
+        #          the order of the arithmetic used to compute the value from
+        #          what the original 2D cross product was using
+        #          -> np.cross(u, v) = u[0]*v[1] - u[1]*v[0]
+        #          where as this would give
+        #          -> np.linalg.det([u, v]) = u[0]*v[1] - v[0]*u[1]
+        #          Note the second term is reversed. Due to the non-commutative
+        #          nature of floating point arithmetic, this can lead to
+        #          different results in the final value. This is why the transpose
+        #          is taken to ensure the original arithmetic order is preserved.
+        D = np.linalg.det(np.array([u, v]).transpose())
 
         if np.allclose(D, 0, rtol=0, atol=1e2 * np.finfo(float).eps):
             return np.array(self._start)
 
-        # We only care about the z component of the cross product
-        return np.cross(v, w)[-1] / D * u[:-1] + self._start
+        # See note above
+        return np.linalg.det(np.array([v, w]).transpose()) / D * u[:-1] + self._start
 
     def is_parallel(self, edge):
         u = self._stop - self._start
         v = edge._stop - edge._start
-        return np.allclose(np.cross(u, v), 0, rtol=0, atol=1e2 * np.finfo(float).eps)
+        return np.allclose(
+            np.linalg.det(np.array([u, v]).transpose()),
+            0,
+            rtol=0,
+            atol=1e2 * np.finfo(float).eps,
+        )
 
 
 def _round_vertex(v):