Feat: Add errors for line and plane of best fit (#362)

Author: ajhynes7

.gitignore

@@ -4,6 +4,7 @@
 *.xml
 .*/
 .coverage
+.python-version
 build/
 dist/
 docs/build/

.python-version

@@ -542,7 +542,9 @@ def _best_fit(points_centered: Points, centroid: Point) -> Tuple[Vector, Point,
 
         def _compute_initial_direction(points: Points) -> np.ndarray:
             """Compute the initial direction as the best fit line."""
-            initial_direction = Line.best_fit(points).vector.unit()
+            line_best_fit = cast(Line, Line.best_fit(points))
+
+            initial_direction = line_best_fit.vector.unit()
             spherical_coordinates = _cartesian_to_spherical(*initial_direction)
 
             return np.array([spherical_coordinates.theta, spherical_coordinates.phi])

src/skspatial/objects/cylinder.py

@@ -657,23 +657,38 @@ def intersect_line(self, other: Line, check_coplanar: bool = True, **kwargs) ->
         return self.point + vector_a_scaled
 
     @classmethod
-    def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) -> Line:
+    def best_fit(
+        cls,
+        points: array_like,
+        tol: Optional[float] = None,
+        return_error: bool = False,
+        **kwargs,
+    ) -> Line | tuple[Line, float]:
         """
         Return the line of best fit for a set of points.
 
+        Also optionally return a value representing the error of the fit.
+        This is the sum of the squared singular values from SVD (excluding the first).
+
+        "The singular values reflect the amount of data variance captured by the bases.
+        The first basis (the one with largest singular value) lies in the direction of the greatest data variance.
+        The second basis captures the orthogonal direction with the second greatest variance, and so on." [1]_
+
         Parameters
         ----------
         points : array_like
              Input points.
         tol : float | None, optional
             Keyword passed to :meth:`Points.are_collinear` (default None).
+        return_error : bool, optional
+            If True, also return a value representing the error of the fit (default False).
         kwargs : dict, optional
             Additional keywords passed to :func:`numpy.linalg.svd`
 
         Returns
         -------
-        Line
-            The line of best fit.
+        Line | tuple[Line, float]
+            The line of best fit, and optionally the error of the fit.
 
         Raises
         ------
@@ -697,6 +712,10 @@ def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) ->
         >>> line.direction.round(3)
         Vector([0.707, 0.707])
 
+        References
+        ----------
+        .. [1] : "Singular Value Decomposition", Oracle, https://docs.oracle.com/en/database/oracle/machine-learning/oml4sql/23/dmcon/singular-value-decomposition.html#GUID-14AA4B45-3B36-4056-9B9A-BD9DC471F0AD
+
         """
         points_spatial = Points(points)
 
@@ -705,10 +724,17 @@ def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) ->
 
         points_centered, centroid = points_spatial.mean_center(return_centroid=True)
 
-        _, _, vh = np.linalg.svd(points_centered, **kwargs)
-        direction = vh[0, :]
+        _, S, Vh = np.linalg.svd(points_centered, **kwargs)
+
+        direction = Vh[0, :]
+        line_best_fit = cls(centroid, direction)
+
+        if return_error:
+            error = np.sum(S[1:] ** 2)
+
+            return line_best_fit, error
 
-        return cls(centroid, direction)
+        return line_best_fit
 
     def transform_points(self, points: array_like) -> np.ndarray:
         """

src/skspatial/objects/line.py

@@ -703,33 +703,44 @@ def intersect_plane(self, other: Plane, **kwargs) -> Line:
         return Line(point_line, direction_line)
 
     @classmethod
-    def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) -> Plane:
+    def best_fit(
+        cls,
+        points: array_like,
+        tol: Optional[float] = None,
+        return_error: bool = False,
+        **kwargs,
+    ) -> Plane | tuple[Plane, float]:
         """
         Return the plane of best fit for a set of 3D points.
 
+        Also optionally return a value representing the error of the fit.
+        This is the sum of the squared singular values from SVD (excluding the first two).
+
+        "The singular values reflect the amount of data variance captured by the bases.
+        The first basis (the one with largest singular value) lies in the direction of the greatest data variance.
+        The second basis captures the orthogonal direction with the second greatest variance, and so on." [1]_
+
         Parameters
         ----------
         points : array_like
              Input 3D points.
         tol : float | None, optional
             Keyword passed to :meth:`Points.are_collinear` (default None).
+        return_error : bool, optional
+            If True, also return a value representing the error of the fit (default False).
         kwargs : dict, optional
             Additional keywords passed to :func:`numpy.linalg.svd`
 
         Returns
         -------
-        Plane
-            The plane of best fit.
+        Plane | tuple[Plane, float]
+            The plane of best fit, and optionally the error of the fit.
 
         Raises
         ------
         ValueError
             If the points are collinear or are not 3D.
 
-        References
-        ----------
-        https://scicomp.stackexchange.com/a/6901
-
         Examples
         --------
         >>> from skspatial.objects import Plane
@@ -755,6 +766,11 @@ def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) ->
         >>> Plane.best_fit(points, full_matrices=False)
         Plane(point=Point([0.5, 0.5, 0. ]), normal=Vector([0., 0., 1.]))
 
+        References
+        ----------
+        .. [1] : "Singular Value Decomposition", Oracle, https://docs.oracle.com/en/database/oracle/machine-learning/oml4sql/23/dmcon/singular-value-decomposition.html#GUID-14AA4B45-3B36-4056-9B9A-BD9DC471F0AD
+        .. [2] : https://scicomp.stackexchange.com/a/6901
+
         """
         points = Points(points)
 
@@ -766,10 +782,16 @@ def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) ->
 
         points_centered, centroid = points.mean_center(return_centroid=True)
 
-        u, _, _ = np.linalg.svd(points_centered.T, **kwargs)
-        normal = Vector(u[:, 2])
+        U, S, _ = np.linalg.svd(points_centered.T, **kwargs)
+        normal = Vector(U[:, 2])
+
+        plane_fit = cls(centroid, normal)
+
+        if return_error:
+            error_fit = np.sum(S[2:] ** 2)
+            return plane_fit, error_fit
 
-        return cls(centroid, normal)
+        return plane_fit
 
     def to_mesh(
         self,

src/skspatial/objects/plane.py

@@ -323,6 +323,26 @@ def test_best_fit(points, line_expected):
     assert line_fit.point.is_close(line_expected.point)
 
 
+@pytest.mark.parametrize(
+    ("points", "line_expected", "error_expected"),
+    [
+        ([[0, 0], [1, 0]], Line([0.5, 0], [1, 0]), 0),
+        ([[1, 0], [0, 0]], Line([0.5, 0], [-1, 0]), 0),
+        ([[0, 0], [1, 1], [2, 2]], Line([1, 1], [1, 1]), 0),
+        ([[0, 0], [0, 1], [1, 0], [1, 1]], Line([0.5, 0.5], [1, 0]), 1),
+        ([[0, 0], [0, 2], [2, 0], [2, 2]], Line([1, 1], [1, 0]), 4),
+        ([[0, 0], [0, 3], [3, 0], [3, 3]], Line([1.5, 1.5], [1, 0]), 9),
+    ],
+)
+def test_best_fit_with_error(points, line_expected, error_expected):
+    line_fit, error_fit = Line.best_fit(np.array(points), return_error=True)
+
+    assert line_fit.is_close(line_expected)
+    assert line_fit.point.is_close(line_expected.point)
+
+    assert math.isclose(error_fit, error_expected, abs_tol=1e-9)
+
+
 @pytest.mark.parametrize(
     ("points", "message_expected"),
     [

tests/unit/objects/test_line.py

@@ -437,6 +437,32 @@ def test_best_fit(points, plane_expected):
     assert plane_fit.point.is_close(plane_expected.point)
 
 
+@pytest.mark.parametrize(
+    ("points", "plane_expected", "error_expected"),
+    [
+        ([[0, 0], [1, 1], [0, 2]], Plane([1 / 3, 1, 0], [0, 0, 1]), 0),
+        (
+            [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]],
+            Plane([0.25, 0.25, 0.25], [1, 1, 1]),
+            0.25,
+        ),
+        (
+            [[0, 0, 0], [2, 0, 0], [0, 2, 0], [0, 0, 2]],
+            Plane([0.5, 0.5, 0.5], [1, 1, 1]),
+            1,
+        ),
+    ],
+)
+def test_best_fit_with_error(points, plane_expected, error_expected):
+    points = Points(points).set_dimension(3)
+    plane_fit, error_fit = Plane.best_fit(points, return_error=True)
+
+    assert plane_fit.is_close(plane_expected)
+    assert plane_fit.point.is_close(plane_expected.point)
+
+    assert math.isclose(error_fit, error_expected, abs_tol=1e-9)
+
+
 @pytest.mark.parametrize(
     ("points", "message_expected"),
     [