diff --git a/pylinalg/matrix.py b/pylinalg/matrix.py index 2cc79ab..1e95548 100644 --- a/pylinalg/matrix.py +++ b/pylinalg/matrix.py @@ -737,8 +737,7 @@ def _mat_inv(m) -> np.ndarray: det = n11 * t11 + n21 * t12 + n31 * t13 + n41 * t14 if det == 0: - out.fill(0) - return out # singular matrix + raise np.linalg.LinAlgError("Singular matrix") det_inv = 1 / det @@ -853,11 +852,11 @@ def _mat_inv(m) -> np.ndarray: if int(np.__version__.split(".")[0]) >= 2: _default_mat_inv_method = "numpy" else: - _default_mat_inv_method = "manual" + _default_mat_inv_method = "python" def mat_inverse( - matrix, /, *, method=_default_mat_inv_method, dtype=None, out=None + matrix, /, *, method=_default_mat_inv_method, raise_err=False, dtype=None, out=None ) -> np.ndarray: """ Compute the inverse of a matrix. @@ -868,7 +867,9 @@ def mat_inverse( The matrix to invert. method : str, optional The method to use for inversion. The default is "numpy" when - numpy version is 2.0.0 or newer, otherwise "manual". + numpy version is 2.0.0 or newer, otherwise "python". + raise_err : bool, optional + Raise a ValueError if the matrix is singular. Default is False. dtype : data-type, optional Overrides the data type of the result. out : ndarray, optional @@ -886,11 +887,11 @@ def mat_inverse( ----- The default method is "numpy" when numpy version >= 2.0.0, which uses the `numpy.linalg.inv` function. - The alternative method is "manual", which uses a manual implementation of - the inversion algorithm. The manual method is used to avoid a performance + The alternative method is "python", which uses a pure python implementation of + the inversion algorithm. The python method is used to avoid a performance issue with `numpy.linalg.inv` on some platforms when numpy version < 2.0.0. See: https://github.com/pygfx/pygfx/issues/763 - The manual method is slower than the numpy method, but it is guaranteed to work. + The python method is slower than the numpy method, but it is guaranteed to work. When the matrix is singular, it will return a matrix filled with zeros, This is a common behavior in real-time graphics applications. @@ -899,18 +900,18 @@ def mat_inverse( fn = { "numpy": np.linalg.inv, - "manual": _mat_inv, + "python": _mat_inv, }[method] matrix = np.asarray(matrix) try: inverse = fn(matrix) - except np.linalg.LinAlgError: + except np.linalg.LinAlgError as err: + if raise_err: + raise ValueError("The provided matrix is not invertible.") from err inverse = np.zeros_like(matrix, dtype=dtype) if out is None: - if dtype is not None: - return inverse.astype(dtype, copy=False) - return inverse + return np.asarray(inverse, dtype=dtype) else: out[:] = inverse return out diff --git a/pylinalg/misc.py b/pylinalg/misc.py index 9fc8409..01ef15b 100644 --- a/pylinalg/misc.py +++ b/pylinalg/misc.py @@ -131,7 +131,7 @@ def quat_to_axis_angle(quaternion, /, *, out=None, dtype=None) -> np.ndarray: quaternion = np.asarray(quaternion) if out is None: - quaternion = quaternion.astype(dtype) + quaternion = quaternion.astype(dtype, copy=False) out = ( quaternion[..., :3] / np.sqrt(1 - quaternion[..., 3] ** 2), 2 * np.arccos(quaternion[..., 3]), diff --git a/pylinalg/vector.py b/pylinalg/vector.py index d20cc26..a024b0a 100644 --- a/pylinalg/vector.py +++ b/pylinalg/vector.py @@ -115,7 +115,9 @@ def vec_transform(vectors, matrix, /, *, w=1, out=None, dtype=None) -> np.ndarra return out -def vec_unproject(vector, matrix, /, *, depth=0, out=None, dtype=None) -> np.ndarray: +def vec_unproject( + vector, matrix, /, *, matrix_is_inv=False, depth=0, out=None, dtype=None +) -> np.ndarray: """ Un-project a vector from 2D space to 3D space. @@ -130,6 +132,9 @@ def vec_unproject(vector, matrix, /, *, depth=0, out=None, dtype=None) -> np.nda The vector to be un-projected. matrix: ndarray, [4, 4] The camera's intrinsic matrix. + matrix_is_inv: bool, optional + Default is False. If True, the provided matrix is assumed to be the + inverse of the camera's intrinsic matrix. depth : number, optional The distance of the unprojected vector from the camera. out : ndarray, optional @@ -162,10 +167,12 @@ def vec_unproject(vector, matrix, /, *, depth=0, out=None, dtype=None) -> np.nda if out is None: out = np.empty((*result_shape, 3), dtype=dtype) - try: - inverse_projection = np.linalg.inv(matrix) - except np.linalg.LinAlgError as err: - raise ValueError("The provided matrix is not invertible.") from err + if matrix_is_inv: + inverse_projection = matrix + else: + from .matrix import mat_inverse + + inverse_projection = mat_inverse(matrix, raise_err=True) vector_hom = np.empty((*result_shape, 4), dtype=dtype) vector_hom[..., 2] = depth @@ -295,7 +302,7 @@ def vec_dist(vector_a, vector_b, /, *, out=None, dtype=None) -> np.ndarray: shape = vector_a.shape[:-1] if out is None: - out = np.linalg.norm(vector_a - vector_b, axis=-1).astype(dtype) + out = np.linalg.norm(vector_a - vector_b, axis=-1).astype(dtype, copy=False) elif len(shape) >= 0: out[:] = np.linalg.norm(vector_a - vector_b, axis=-1) else: @@ -346,7 +353,7 @@ def vec_angle(vector_a, vector_b, /, *, out=None, dtype=None) -> np.ndarray: ) if out is None: - out = np.arccos(the_cos).astype(dtype) + out = np.arccos(the_cos).astype(dtype, copy=False) elif len(shape) >= 0: out[:] = np.arccos(the_cos) else: diff --git a/pyproject.toml b/pyproject.toml index 6380216..c0420c4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -2,7 +2,7 @@ [project] name = "pylinalg" -version = "0.6.2" +version = "0.6.3" description = "Linear algebra utilities for Python" readme = "README.md" license = { file = "LICENSE" } diff --git a/tests/test_matrix.py b/tests/test_matrix.py index 04662cb..8415c19 100644 --- a/tests/test_matrix.py +++ b/tests/test_matrix.py @@ -385,15 +385,15 @@ def test_mat_euler_vs_scipy(): def test_mat_inverse(): a = la.mat_from_translation([1, 2, 3]) np_inv = la.mat_inverse(a, method="numpy") - manual_inv = la.mat_inverse(a, method="manual") - npt.assert_array_equal(np_inv, manual_inv) + python_inv = la.mat_inverse(a, method="python") + npt.assert_array_equal(np_inv, python_inv) # test for singular matrix b = la.mat_from_scale([0, 2, 3]) np_inv = la.mat_inverse(b, method="numpy") - manual_inv = la.mat_inverse(b, method="manual") + python_inv = la.mat_inverse(b, method="python") npt.assert_array_equal(np_inv, np.zeros((4, 4))) - npt.assert_array_equal(manual_inv, np.zeros((4, 4))) + npt.assert_array_equal(python_inv, np.zeros((4, 4))) def test_mat_has_shear(): diff --git a/tests/test_pylinalg.py b/tests/test_pylinalg.py index d67e854..1afa66c 100644 --- a/tests/test_pylinalg.py +++ b/tests/test_pylinalg.py @@ -58,9 +58,9 @@ def popattr(key): # confirm the signature of each callable sig = inspect.signature(getattr(la, key)) - assert ( - sig.return_annotation is not inspect.Signature.empty - ), f"Missing return annotation on {key}" + assert sig.return_annotation is not inspect.Signature.empty, ( + f"Missing return annotation on {key}" + ) if sig.return_annotation is bool: key_parts = key.split("_") assert key_parts[1] in ("is", "has") @@ -77,9 +77,9 @@ def popattr(key): assert param.KEYWORD_ONLY has_out = True assert has_out, f"Function {key} does not have 'out' keyword-only argument" - assert ( - has_dtype - ), f"Function {key} does not have 'dtype' keyword-only argument" + assert has_dtype, ( + f"Function {key} does not have 'dtype' keyword-only argument" + ) # assert that all callables are available in __all__ assert set(__all__) == set(callables) diff --git a/tests/test_vectors.py b/tests/test_vectors.py index e69dd68..b5e923b 100644 --- a/tests/test_vectors.py +++ b/tests/test_vectors.py @@ -303,6 +303,16 @@ def test_vec_unproject_exceptions(): la.vec_unproject(vector, matrix) +def test_vec_unproject_is_inverse(): + a = la.mat_perspective(-1, 1, -1, 1, 1, 100) + a_inv = la.mat_inverse(a) + vecs = np.array([[1, 2], [4, 5], [7, 8]]) + + expected = la.vec_unproject(vecs, a) + actual = la.vec_unproject(vecs, a_inv, matrix_is_inv=True) + npt.assert_array_equal(expected, actual) + + def test_vector_apply_rotation_ordered(): """Test that a positive pi/2 rotation about the z-axis and then the y-axis results in a different output then in standard rotation ordering."""