Repository: pyvista/pyvista

Issue:
Add `covariance` method to `fit_plane_to_points` and expose `Plane` arguments
### Describe the feature you would like to be added.

This is a follow-up to #4650.

### Links to VTK Documentation, Examples, or Class Definitions.

Currently, `fit_plane_to_points` uses SVD. However, calculating fit using covariance may be more "VTK-like" since:
[`vtkPlane`](https://vtk.org/doc/nightly/html/classvtkPlane.html#a7b88495d93a5aa5053b1044965df67bd) has a fitting method which uses covariance, and [`vtkOBBTree`](https://vtk.org/doc/nightly/html/classvtkOBBTree.html) also computes fit using covariance. I provided a graphical example of the different fit methods [here](https://github.com/pyvista/pyvista/issues/555#issuecomment-1646108382).

### Pseudocode or Screenshots

I propose adding `method='covariance'` as the default arg for `fit_plane_to_points` and having `method='svd'` as a non-default alternative.

In addition, sometimes it's useful to extract only the corners of the best-fit plane. In this case, we want to set `i_resolution=1` and `j_resolution=1` in the `Plane` constructor so that the plane mesh only has four points. However, we don't have direct access to these parameters in `fit_plane_to_points`. So I think we should expose these parameters directly.

This would look something like:
``` python 
def fit_plane_to_points(
    points,
    method='covariance',
    return_meta=False,
     i_resolution=10, 
    j_resolution=10):
```

Any thoughts on this or suggestions for an alternative API (e.g. use a `plane_parameters` dict instead)? Should the default resolution be 1 instead of 10?

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pyvista/pyvista/pull/6528?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 91.50%. Comparing base [(`99a71b5`)](https://app.codecov.io/gh/pyvista/pyvista/commit/99a71b5096d44299a78fba7ef428b0f9e14fa04b?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) to head [(`8f69d82`)](https://app.codecov.io/gh/pyvista/pyvista/commit/8f69d82aa2c351a601d2fb71a80e0099fd688745?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista).
> Report is 1 commits behind head on main.

> :exclamation:  There is a different number of reports uploaded between BASE (99a71b5) and HEAD (8f69d82). Click for more details.
> 
> <details><summary>HEAD has 5 uploads less than BASE</summary>
>
>| Flag | BASE (99a71b5) | HEAD (8f69d82) |
>|------|------|------|
>||6|1|
></details>

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


```diff
@@            Coverage Diff             @@
##             main    #6528      +/-   ##
==========================================
- Coverage   97.38%   91.50%   -5.89%     
==========================================
  Files         143      143              
  Lines       27386    27388       +2     
==========================================
- Hits        26671    25062    -1609     
- Misses        715     2326    +1611     
```

</details>
pyvista-bot: <!-- NETLIFY DEPLOY COMMENT GENERATED BY ACTIONS_NETLIFY - APP ID SHA256: d3574f413fe51079937885070831e027ca5b10f7facb9e609494e0a481a31f79 -->
🚀 Deployed on https://66c18db03f837a5464c50e35--pyvista-dev.netlify.app

Patch:
diff --git a/pyvista/core/utilities/points.py b/pyvista/core/utilities/points.py
index 57e4428e..a07daf29 100644
--- a/pyvista/core/utilities/points.py
+++ b/pyvista/core/utilities/points.py
@@ -185,7 +185,7 @@ def lines_from_points(points, close=False):
     return poly
 
 
-def fit_plane_to_points(points, return_meta=False):
+def fit_plane_to_points(points, return_meta=False, resolution=10):
     """Fit a plane to a set of points using the SVD algorithm.
 
     The plane is automatically sized and oriented to fit the extents of
@@ -208,6 +208,13 @@ def fit_plane_to_points(points, return_meta=False):
         If ``True``, also returns the center and normal of the
         generated plane.
 
+    resolution : int, default: 10
+        Number of points on the plane mesh along its edges. Specify two numbers to
+        set the resolution along the plane's long and short edge (respectively) or
+        a single number to set both edges to have the same resolution.
+
+        .. versionadded:: 0.45.0
+
     Returns
     -------
     pyvista.PolyData
@@ -225,22 +232,18 @@ def fit_plane_to_points(points, return_meta=False):
 
     >>> import pyvista as pv
     >>> import numpy as np
+    >>> from pyvista import examples
     >>>
-    >>> # Create point cloud
     >>> rng = np.random.default_rng(seed=0)
     >>> cloud = rng.random((10, 3))
     >>> cloud[:, 2] *= 0.1
     >>>
-    >>> # Fit plane
-    >>> plane, center, normal = pv.fit_plane_to_points(
-    ...     cloud, return_meta=True
-    ... )
-    >>>
-    >>> # Plot the fitted plane
+    >>> plane = pv.fit_plane_to_points(cloud)
+
+    Plot the point cloud and fitted plane.
+
     >>> pl = pv.Plotter()
-    >>> _ = pl.add_mesh(
-    ...     plane, color='lightblue', style='wireframe', line_width=4
-    ... )
+    >>> _ = pl.add_mesh(plane, style='wireframe', line_width=4)
     >>> _ = pl.add_points(
     ...     cloud,
     ...     render_points_as_spheres=True,
@@ -249,22 +252,18 @@ def fit_plane_to_points(points, return_meta=False):
     ... )
     >>> pl.show()
 
-    Fit a plane to a mesh.
+    Fit a plane to a mesh and return its meta-data. Set the plane resolution to 1
+    so that the plane has no internal points or edges.
 
-    >>> import pyvista as pv
-    >>> from pyvista import examples
-    >>>
-    >>> # Create mesh
     >>> mesh = examples.download_shark()
-    >>>
-    >>> # Fit plane
-    >>> plane = pv.fit_plane_to_points(mesh.points)
-    >>>
-    >>> # Plot the fitted plane
-    >>> pl = pv.Plotter()
-    >>> _ = pl.add_mesh(
-    ...     plane, show_edges=True, color='lightblue', opacity=0.25
+    >>> plane, center, normal = pv.fit_plane_to_points(
+    ...     mesh.points, return_meta=True, resolution=1
     ... )
+
+    Plot the mesh and fitted plane.
+
+    >>> pl = pv.Plotter()
+    >>> _ = pl.add_mesh(plane, show_edges=True, opacity=0.25)
     >>> _ = pl.add_mesh(mesh, color='gray')
     >>> pl.camera_position = [
     ...     (-117, 76, 235),
@@ -273,7 +272,34 @@ def fit_plane_to_points(points, return_meta=False):
     ... ]
     >>> pl.show()
 
+    Use the meta data with :meth:`pyvista.DataSetFilter.clip` to split the mesh into
+    two.
+
+    >>> first_half, second_half = mesh.clip(
+    ...     origin=center, normal=normal, return_clipped=True
+    ... )
+
+    Plot the two halves of the clipped mesh.
+
+    >>> pl = pv.Plotter()
+    >>> _ = pl.add_mesh(first_half, color='red')
+    >>> _ = pl.add_mesh(second_half, color='blue')
+    >>> pl.camera_position = [
+    ...     (-143, 43, 40),
+    ...     (-8.7, -11, -14),
+    ...     (0.25, 0.92, -0.29),
+    ... ]
+    >>> pl.show()
     """
+    valid_resolution = _validation.validate_array(
+        resolution,
+        must_have_shape=[(), (2,)],
+        must_be_integer=True,
+        broadcast_to=(2,),
+        dtype_out=int,
+    )
+    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)
@@ -304,7 +330,14 @@ def fit_plane_to_points(points, return_meta=False):
     center = rotate_transform_inv[:3, :3] @ projected.center + data_center
 
     # Initialize plane then move to final position
-    plane = pyvista.Plane(center=(0, 0, 0), direction=(0, 0, 1), i_size=i_size, j_size=j_size)
+    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
 
@@ -503,6 +536,9 @@ def principal_axes(points: MatrixLike[float], *, return_std: bool = False):
 
     See Also
     --------
+    fit_plane_to_points
+        Fit a plane to points using the first two principal axes.
+
     pyvista.DataSetFilters.align_xyz
         Filter which aligns principal axes to the x-y-z axes.
 


Test patch:
diff --git a/tests/core/test_utilities.py b/tests/core/test_utilities.py
index c2dbd292..6ce35dea 100644
--- a/tests/core/test_utilities.py
+++ b/tests/core/test_utilities.py
@@ -902,6 +902,16 @@ def test_has_module():
     assert not has_module('not_a_module')
 
 
+def test_fit_plane_to_points_resolution(airplane):
+    DEFAULT_RESOLUTION = 10
+    plane = fit_plane_to_points(airplane.points)
+    assert plane.n_points == (DEFAULT_RESOLUTION + 1) ** 2
+
+    resolution = (1.0, 2.0)  # Test with integer-valued floats
+    plane = fit_plane_to_points(airplane.points, resolution=resolution)
+    assert plane.n_points == (resolution[0] + 1) * (resolution[1] + 1)
+
+
 def test_fit_plane_to_points():
     points = ex.load_airplane().points
     plane, center, normal = fit_plane_to_points(points, return_meta=True)