Repository: pyvista/pyvista

Issue:
`fit_plane_to_points` crashes python with many points
### Describe the bug, what's wrong, and what you expected.

This code crashes python:
```python
import pyvista as pv
points = pv.examples.download_armadillo().points
len(points)  # 172974
pv.fit_plane_to_points(points)  # crashes
# Process finished with exit code 137 (interrupted by signal 9: SIGKILL)
```
I believe the issue stems from using `np.linalg.svd` to fit the points, which has a high memory footprint. Instead, PCA may be used for the computation. PCA is used in the `principal_axes` method from #6346, so that method can be used internally to replace the call to `np.linalg.svd`

### Steps to reproduce the bug.

See above

### System Information

```shell
OS : Darwin (macOS 14.4.1)
            CPU(s) : 8
           Machine : arm64
      Architecture : 64bit
       Environment : Python
        GPU Vendor : Apple
      GPU Renderer : Apple M2
       GPU Version : 4.1 Metal - 88
  MathText Support : True

  Python 3.9.13 (v3.9.13:6de2ca5339, May 17 2022, 11:37:23)  [Clang 13.0.0
  (clang-1300.0.29.30)]

           pyvista : 0.45.dev0
               vtk : 9.3.1
             numpy : 2.0.0
        matplotlib : 3.9.0
            scooby : 0.10.0
             pooch : 1.8.2
            pillow : 10.4.0
           IPython : 8.18.1
          colorcet : 3.1.0
             scipy : 1.13.1
              tqdm : 4.66.4
    pytest_pyvista : 0.1.9
```


### Screenshots

_No response_

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pyvista/pyvista/pull/6517?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 94.62%. Comparing base [(`a8bccf0`)](https://app.codecov.io/gh/pyvista/pyvista/commit/a8bccf007d2196ca233ebbebd1f0f1b5ecf2acc8?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) to head [(`60fb16d`)](https://app.codecov.io/gh/pyvista/pyvista/commit/60fb16d3abe053e52b4c1c330896c50842cb5a3f?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #6517      +/-   ##
==========================================
- Coverage   97.34%   94.62%   -2.72%     
==========================================
  Files         143      143              
  Lines       27339    27332       -7     
==========================================
- Hits        26613    25864     -749     
- Misses        726     1468     +742     
```

</details>
pyvista-bot: <!-- NETLIFY DEPLOY COMMENT GENERATED BY ACTIONS_NETLIFY - APP ID SHA256: d3574f413fe51079937885070831e027ca5b10f7facb9e609494e0a481a31f79 -->
🚀 Deployed on https://66bfde5b50dbc6c0d21e1d7e--pyvista-dev.netlify.app
user27182: I think I'll split this up to keep it maintenance only,  and add features in a separate PR
user27182: I think this may be a breaking change in the sense that the direction of the normal may be flipped between the old implementation and the new one. Any code fitting a plane, then relying on this plane for signed distance calcs or for clipping a dataset would now silently return a different result. 

Do we need a warning in the docs, or just a versionchanged directive?
tkoyama010: Please add a warning in the docs.
user27182: Keeping this as draft based on this comment: https://github.com/pyvista/pyvista/issues/6460#issuecomment-2299735699
I want to explore keeping SVD as an option for computing the principal axes
user27182: > I want to explore keeping SVD as an option for computing the principal axes

This is explored here: https://github.com/pyvista/pyvista/issues/6460#issuecomment-2302993313. We should not use SVD. This is ready for review as-is.

Patch:
diff --git a/pyvista/core/utilities/points.py b/pyvista/core/utilities/points.py
index a07daf29..8ebc433e 100644
--- a/pyvista/core/utilities/points.py
+++ b/pyvista/core/utilities/points.py
@@ -186,7 +186,7 @@ def lines_from_points(points, close=False):
 
 
 def fit_plane_to_points(points, return_meta=False, resolution=10):
-    """Fit a plane to a set of points using the SVD algorithm.
+    """Fit a plane to points using its :func:`principal_axes`.
 
     The plane is automatically sized and oriented to fit the extents of
     the points.
@@ -199,6 +199,16 @@ def fit_plane_to_points(points, return_meta=False, resolution=10):
         computed as the center of the generated plane mesh. In previous
         versions, the center of the input points was returned.
 
+    .. versionchanged:: 0.45.0
+        The internal method used for fitting the plane has changed. Previously, singular
+        value decomposition (SVD) was used, but eigenvectors are now used instead.
+        See warning below.
+
+    .. warning::
+        The normal direction of the plane prior to version 0.45 may differ
+        from the latest version of this function. This may impact methods which
+        rely on the plane's direction.
+
     Parameters
     ----------
     points : array_like[float]
@@ -215,6 +225,11 @@ def fit_plane_to_points(points, return_meta=False, resolution=10):
 
         .. versionadded:: 0.45.0
 
+    See Also
+    --------
+    principal_axes
+        Compute axes vectors which best fit a set of points.
+
     Returns
     -------
     pyvista.PolyData
@@ -300,48 +315,28 @@ def fit_plane_to_points(points, return_meta=False, resolution=10):
     )
     i_resolution, j_resolution = valid_resolution
 
-    # Apply SVD to get orthogonal basis vectors to define the plane
-    data = np.array(points)
-    data_center = data.mean(axis=0)
-    _, _, Vh = np.linalg.svd(data - data_center)
-    i_vector = Vh[0]
-    j_vector = Vh[1]
-    normal = np.cross(i_vector, j_vector)
-
-    # Create rotation matrix from basis vectors
-    rotate_transform = np.eye(4)
-    rotate_transform[:3, :3] = np.vstack((i_vector, j_vector, normal))
-    rotate_transform_inv = rotate_transform.T
-
-    # Project and transform points to align and center data to the XY plane
-    poly = pyvista.PolyData(points)
-    projected = poly.project_points_to_plane(origin=data_center, normal=normal)
-    projected.points -= data_center
-    projected.transform(rotate_transform)
-
-    # Compute size of the plane
-    i_size = projected.bounds.x_max - projected.bounds.x_min
-    j_size = projected.bounds.y_max - projected.bounds.y_min
-
-    # The center of the input data does not necessarily coincide with
-    # the center of the plane. The true center of the plane is the
-    # middle of the bounding box of the projected + transformed data
-    # relative to the input data's center
-    center = rotate_transform_inv[:3, :3] @ projected.center + data_center
-
-    # Initialize plane then move to final position
+    # Align points to the xyz-axes
+    aligned, matrix = pyvista.PolyData(points).align_xyz(return_matrix=True)
+
+    # Fit plane to xyz-aligned mesh
+    aligned_bnds = aligned.bounds
+    i_size = aligned_bnds.x_max - aligned_bnds.x_min
+    j_size = aligned_bnds.y_max - aligned_bnds.y_min
     plane = pyvista.Plane(
-        center=(0, 0, 0),
-        direction=(0, 0, 1),
         i_size=i_size,
         j_size=j_size,
         i_resolution=i_resolution,
         j_resolution=j_resolution,
     )
-    plane.transform(rotate_transform_inv)
-    plane.points += center
+
+    # Transform plane back to input points positioning
+    inverse_matrix = pyvista.Transform(matrix).inverse_matrix
+    plane.transform(inverse_matrix, inplace=True)
 
     if return_meta:
+        # Compute center and normal from the plane's points and normals
+        center = np.mean(plane.points, axis=0)
+        normal = np.mean(plane.point_normals, axis=0)
         return plane, center, normal
     return plane
 


Test patch:
diff --git a/tests/core/test_utilities.py b/tests/core/test_utilities.py
index 6ce35dea..fe1b5ad4 100644
--- a/tests/core/test_utilities.py
+++ b/tests/core/test_utilities.py
@@ -913,23 +913,27 @@ def test_fit_plane_to_points_resolution(airplane):
 
 
 def test_fit_plane_to_points():
-    points = ex.load_airplane().points
-    plane, center, normal = fit_plane_to_points(points, return_meta=True)
-
-    assert np.allclose(normal, [-2.5999512e-08, 0.121780515, -0.99255705])
-    assert np.allclose(center, [896.9954860028446, 686.6470205328502, 78.13187948615939])
-    assert np.allclose(
-        plane.bounds,
-        [
-            139.06036376953125,
-            1654.9306640625,
-            38.0776252746582,
-            1335.2164306640625,
-            -1.4434913396835327,
-            157.70724487304688,
-        ],
+    # Fit a plane to a plane's points
+    center = (1, 2, 3)
+    direction = np.array((4.0, 5.0, 6.0))
+    direction /= np.linalg.norm(direction)
+    expected_plane = pv.Plane(center=center, direction=direction, i_size=2, j_size=3)
+    fitted_plane, fitted_center, fitted_normal = fit_plane_to_points(
+        expected_plane.points, return_meta=True
     )
 
+    # Test bounds
+    assert np.allclose(fitted_plane.bounds, expected_plane.bounds, atol=1e-6)
+
+    # Test center
+    assert np.allclose(fitted_plane.center, center)
+    assert np.allclose(fitted_center, center)
+    assert np.allclose(fitted_plane.points.mean(axis=0), center)
+
+    # Test normal
+    assert np.allclose(fitted_normal, direction)
+    assert np.allclose(fitted_plane.point_normals.mean(axis=0), direction)
+
 
 # Default output from `np.linalg.eigh`
 DEFAULT_PRINCIPAL_AXES = [[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0]]
@@ -1021,11 +1025,17 @@ def one_million_points():
 
 
 def test_principal_axes_success_with_many_points(one_million_points):
-    # Use large mesh to verify no memory errors are raised
+    # Use many points to verify no memory errors are raised
     axes = pv.principal_axes(one_million_points)
     assert isinstance(axes, np.ndarray)
 
 
+def test_fit_plane_to_points_success_with_many_points(one_million_points):
+    # Use many points to verify no memory errors are raised
+    plane = pv.fit_plane_to_points(one_million_points)
+    assert isinstance(plane, pv.PolyData)
+
+
 @pytest.fixture()
 def no_new_attr_subclass():
     @no_new_attr