Repository: pyvista/pyvista

Issue:
Confusing behaviour of ParametricEllipsoid
### Describe the bug, what's wrong, and what you expected.

When creating a ParametricEllispoid using a direction of [0, 1, 0], the ellipsoid is rotated along the y axis.  
For example if setting the direction to [1e-5, 1, 0], which corresponds to approximately similar direction, the ellipsoid displays then the correct behaviour.

### Steps to reproduce the bug.

```python
import pyvista as pv

ellipsoid = pv.ParametricEllipsoid(300, 100, 10, direction=[0, 1, 0])
```

### System Information

```shell
Date: Wed Sep 06 14:07:38 2023 CEST

                OS : Linux
            CPU(s) : 8
           Machine : x86_64
      Architecture : 64bit
               RAM : 31.2 GiB
       Environment : Jupyter
       File system : ext4
        GPU Vendor : Intel
      GPU Renderer : Mesa Intel(R) UHD Graphics (CML GT2)
       GPU Version : 4.6 (Core Profile) Mesa 22.0.1
  MathText Support : False

  Python 3.8.13 (default, Apr 19 2022, 02:32:06)  [GCC 11.2.0]

           pyvista : 0.41.1
               vtk : 9.2.6
             numpy : 1.24.4
        matplotlib : 3.3.4
            scooby : 0.7.2
             pooch : v1.7.0
           imageio : 2.31.1
           IPython : 8.12.2
        ipywidgets : 8.0.7
             scipy : 1.10.1
              tqdm : 4.65.0
        jupyterlab : 3.6.5
             trame : 2.5.2
      trame_client : 2.10.0
      trame_server : 2.11.7
         trame_vtk : 2.5.8
      nest_asyncio : 1.5.6
```


### Screenshots

Here is the given ellipsoid
![confusing_ellipsoid](https://github.com/pyvista/pyvista/assets/57682091/f0e1b5f7-eca1-4224-a020-df44385ed68b)
Here is what is expected
![expected_ellipsoid](https://github.com/pyvista/pyvista/assets/57682091/d4f67ead-9928-4af3-9c3a-b6121180b780)


Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pyvista/pyvista/pull/4853?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) Report
> Merging [#4853](https://app.codecov.io/gh/pyvista/pyvista/pull/4853?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) (36158be) into [main](https://app.codecov.io/gh/pyvista/pyvista/commit/f0db627a9b3eecb0ae95a338fc06c5740ea868d6?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) (f0db627) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

```diff
@@           Coverage Diff           @@
##             main    #4853   +/-   ##
=======================================
  Coverage   95.78%   95.79%           
=======================================
  Files         130      130           
  Lines       21305    21305           
=======================================
+ Hits        20408    20410    +2     
+ Misses        897      895    -2     
```


Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index d202452d..6132cd55 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -6,7 +6,7 @@ ci:
 
 repos:
 - repo: https://github.com/psf/black
-  rev: 23.7.0
+  rev: 23.9.1
   hooks:
   - id: black
 
diff --git a/examples/01-filter/compute-volume.py b/examples/01-filter/compute-volume.py
index 3633107a..981f8b31 100644
--- a/examples/01-filter/compute-volume.py
+++ b/examples/01-filter/compute-volume.py
@@ -76,12 +76,11 @@ print(f"Original volume: {dataset.volume}")
 
 ###############################################################################
 # Or better yet, you could simply extract the largest volume from your
-# thresholded dataset by passing ``largest=True`` to the ``connectivity``
-# filter or by using ``extract_largest`` filter (both are equivalent).
+# dataset directly by passing ``'largest'`` to the ``connectivity`` and
+# specify the same scalar range used earlier for thresholding.
 
-# Grab the largest connected volume present
-largest = threshed.connectivity(largest=True)
-# or: largest = threshed.extract_largest()
+# Grab the largest connected volume within a scalar range
+largest = threshed.connectivity('largest', scalar_range=[0.15, 0.50])
 
 # Get volume as numeric value
 large_volume = largest.volume
diff --git a/examples/01-filter/connectivity.py b/examples/01-filter/connectivity.py
index 23c70007..8a3acbcf 100644
--- a/examples/01-filter/connectivity.py
+++ b/examples/01-filter/connectivity.py
@@ -27,12 +27,12 @@ p.show(cpos=cpos)
 ###############################################################################
 # The mesh plotted above is very noisy. We can extract the largest connected
 # isosurface in that mesh using the :func:`pyvista.DataSetFilters.connectivity`
-# filter and passing ``largest=True`` to the ``connectivity``
+# filter and passing ``'largest'`` to the ``connectivity``
 # filter or by using the :func:`pyvista.DataSetFilters.extract_largest` filter
 # (both are equivalent).
 
 # Grab the largest connected volume present
-largest = mesh.connectivity(largest=True)
+largest = mesh.connectivity('largest')
 # or: largest = mesh.extract_largest()
 
 p = pv.Plotter()
diff --git a/pyvista/core/_vtk_core.py b/pyvista/core/_vtk_core.py
index e17db5d8..dfa4d8b6 100644
--- a/pyvista/core/_vtk_core.py
+++ b/pyvista/core/_vtk_core.py
@@ -253,6 +253,11 @@ from vtkmodules.vtkFiltersExtraction import (
     vtkExtractSelection,
 )
 from vtkmodules.vtkFiltersFlowPaths import vtkEvenlySpacedStreamlines2D, vtkStreamTracer
+
+try:  # Introduced VTK v9.1.0
+    from vtkmodules.vtkFiltersGeneral import vtkRemovePolyData
+except ImportError:  # pragma: no cover
+    pass
 from vtkmodules.vtkFiltersGeneral import (
     vtkAxes,
     vtkBooleanOperationPolyDataFilter,
diff --git a/pyvista/core/filters/data_set.py b/pyvista/core/filters/data_set.py
index c1f9712d..c935cd4a 100644
--- a/pyvista/core/filters/data_set.py
+++ b/pyvista/core/filters/data_set.py
@@ -1,6 +1,6 @@
 """Filters module with a class of common filters that can be applied to any vtkDataSet."""
 import collections.abc
-from typing import Optional, Sequence, Union
+from typing import Literal, Optional, Sequence, Union
 import warnings
 
 import matplotlib.pyplot as plt
@@ -2373,58 +2373,483 @@ class DataSetFilters:
         _update_alg(alg, progress_bar, 'Computing Glyphs')
         return _get_output(alg)
 
-    def connectivity(self, largest=False, progress_bar=False):
-        """Find and label connected bodies/volumes.
+    def connectivity(
+        self,
+        extraction_mode: Literal[
+            'all', 'largest', 'specified', 'cell_seed', 'point_seed', 'closest'
+        ] = 'all',
+        variable_input=None,
+        scalar_range=None,
+        scalars=None,
+        label_regions=True,
+        region_ids=None,
+        point_ids=None,
+        cell_ids=None,
+        closest_point=None,
+        inplace=False,
+        progress_bar=False,
+        **kwargs,
+    ):
+        """Find and label connected regions.
+
+        This filter extracts cell regions based on a specified connectivity
+        criterion. The extraction criterion can be controlled with
+        ``extraction_mode`` to extract the largest region or the closest
+        region to a seed point, for example.
 
-        This adds an ID array to the point and cell data to
-        distinguish separate connected bodies. This applies a
-        ``vtkConnectivityFilter`` filter which extracts cells that
-        share common points and/or meet other connectivity criterion.
+        In general, cells are considered to be connected if they
+        share a point. However, if a ``scalar_range`` is provided, cells
+        must also have at least one point with scalar values in the
+        specified range to be considered connected.
 
-        Cells that share vertices and meet other connectivity
-        criterion such as scalar range are known as a region.
+        See :ref:`connectivity_example` and :ref:`volumetric_example` for
+        more examples using this filter.
+
+        .. versionadded:: 0.43.0
+           * New extraction modes: ``'specified'``, ``'cell_seed'``, ``'point_seed'``,
+             and ``'closest'``.
+           * Extracted regions are now sorted in descending order by
+             cell count.
+           * Region connectivity can be controlled using ``scalar_range``.
+
+        .. deprecated:: 0.43.0
+           Parameter ``largest`` is deprecated. Use ``'largest'`` or
+           ``extraction_mode='largest'`` instead.
 
         Parameters
         ----------
-        largest : bool, default: False
-            Extract the largest connected part of the mesh.
+        extraction_mode : str, default: "all"
+            * ``'all'``: Extract all connected regions.
+            * ``'largest'`` : Extract the largest connected region (by cell
+              count).
+            * ``'specified'``: Extract specific region IDs. Use ``region_ids``
+              to specify the region IDs to extract.
+            * ``'cell_seed'``: Extract all regions sharing the specified cell
+              ids. Use ``cell_ids`` to specify the cell ids.
+            * ``'point_seed'`` : Extract all regions sharing the specified
+              point ids. Use ``point_ids`` to specify the point ids.
+            * ``'closest'`` : Extract the region closest to the specified
+              point. Use ``closest_point`` to specify the point.
+
+        variable_input : float | sequence[float], optional
+            The convenience parameter used for specifying any required input values
+            for some values of ``extraction_mode``. Setting
+            ``extraction_input`` is equivalent to setting:
+            * ``'region_ids'`` if mode is ``'specified'``.
+            * ``'cell_ids'`` if mode is ``'cell_seed'``.
+            * ``'point_ids'`` if mode is ``'point_seed'``.
+            * ``'closest_point'`` if mode is ``'closest'``. It has no effect if the mode is ``'all'`` or ``'largest'``.
+
+        scalar_range : float | sequence[float], optional
+            Single value or ``(min, max)``. If set, the connectivity is
+            restricted to cells with at least one point with scalar values in
+            the specified range.
+
+        scalars : str, optional
+            Name of scalars to use if ``scalar_range`` is specified. Defaults
+            to currently active scalars.
+
+            .. note::
+               This filter requires point scalars to determine region
+               connectivity. If cell scalars are provided, they are first
+               converted to point scalars with :func:`cell_data_to_point_data`
+               before applying the filter. The converted point scalars are
+               removed from the output after applying the filter.
+
+        label_regions : bool, default: True
+            If ``True``, scalar point and cell arrays ``RegionId`` are stored.
+            Each region is assigned a unique ID. IDs are zero-indexed and are
+            assigned by region cell count in descending order (i.e. the largest
+            region has ID ``0``).
+
+        region_ids : sequence[int], optional
+            Region ids to extract. Only used if ``extraction_mode`` is
+            ``specified``.
+
+        point_ids : sequence[int], optional
+            Point ids to use as seeds. Only used if ``extraction_mode`` is
+            ``point_seed``.
+
+        cell_ids : sequence[int], optional
+            Cell ids to use as seeds. Only used if ``extraction_mode`` is
+            ``cell_seed``.
+
+        closest_point : sequence[int], optional
+            Point coordinates in ``(x, y, z)``. Only used if
+            ``extraction_mode`` is ``closest``.
+
+        inplace : bool, default: False
+            If ``True`` the mesh is updated in-place, otherwise a copy
+            is returned. A copy is always returned if the input type is
+            not ``pyvista.PolyData`` or ``pyvista.UnstructuredGrid``.
 
         progress_bar : bool, default: False
-            Display a progress bar to indicate progress.
+            Display a progress bar.
+
+        **kwargs : dict, optional
+            Used for handling deprecated parameters.
 
         Returns
         -------
         pyvista.DataSet
-            Dataset with labeled connected bodies.  Return type
-            matches input.
+            Dataset with labeled connected bodies. Return type is
+            ``pyvista.PolyData`` if input type is ``pyvista.PolyData`` and
+            ``pyvista.UnstructuredGrid`` otherwise.
+
+        See Also
+        --------
+            :func:`pyvista.DataSetFilter.connectivity`
 
         Examples
         --------
-        Join two meshes together and plot their connectivity.
+        Label all connected regions.
 
-        >>> import pyvista
-        >>> mesh = pyvista.Sphere() + pyvista.Sphere(center=(2, 0, 0))
-        >>> conn = mesh.connectivity(largest=False)
-        >>> conn.plot(cmap=['red', 'blue'])
+        Load data.
+
+        >>> import numpy as np
+        >>> from pyvista import examples
+        >>> import pyvista as pv
+        >>> mesh = examples.download_foot_bones()
+
+        Extract all disconnected regions.
+
+        >>> conn = mesh.connectivity()
+
+        Configure plotting parameters.
+
+        >>> categories = True
+        >>> cmap = 'glasbey'
+
+        Format scalar bar text for integer values.
+
+        >>> scalar_bar_args = dict(
+        ...     fmt='%.f',  # Do not show decimals
+        ... )
+
+        Set the plot's camera position.
+
+        >>> cpos = [(10.5, 12.2, 18.3), (0.0, 0.0, 0.0), (0.0, 1.0, 0.0)]
+
+        Plot the regions by ID.
+
+        Note that region IDs are sorted in descending order by
+        cell count.
+
+        >>> conn.plot(
+        ...     categories=categories,
+        ...     cmap=cmap,
+        ...     scalar_bar_args=scalar_bar_args,
+        ...     cpos=cpos,
+        ... )
+
+        Extract specific regions by size.
+        Calculate region sizes.
+
+        Get counts using the previous `connectivity('all')` results.
+
+        >>> regions, region_sizes = np.unique(
+        ...     conn['RegionId'], return_counts=True
+        ... )
+
+        Extract regions with fewer than 150 cells.
+
+        >>> region_ids = regions[np.flatnonzero(region_sizes < 150)]
+        >>> conn = mesh.connectivity('specified', region_ids)
+
+        Plot result with plotting parameters defined earlier.
+
+        >>> conn.plot(
+        ...     categories=categories,
+        ...     cmap=cmap,
+        ...     scalar_bar_args=scalar_bar_args,
+        ...     cpos=cpos,
+        ... )
+
+        Extract the largest region.
+
+        >>> conn = mesh.connectivity('largest', label_regions=False)
+
+        Plot the largest region and show the input mesh for reference.
+
+        >>> p = pv.Plotter()
+        >>> _ = p.add_mesh(conn)
+        >>> _ = p.add_mesh(mesh, style='wireframe')
+        >>> p.show(cpos=cpos)
+
+        Extract regions using seed points.
+
+        Create hills and use curvature to define their peaks and valleys.
+
+        >>> import pyvista as pv
+        >>> mesh = pv.ParametricRandomHills()
+        >>> mesh["Curvature"] = mesh.curvature()
+
+        Visualize the peaks and valleys.
+
+        Peaks have large positive curvature (i.e. are convex)
+        Valleys have large negative curvature (i.e. are concave)
+        Flat regions have curvature close to zero.
+
+        >>> mesh.plot(
+        ...     clim=[-0.5, 0.5],
+        ...     cmap='bwr',
+        ...     below_color='blue',
+        ...     above_color='red',
+        ... )
 
-        See :ref:`volumetric_example` for more examples using this filter.
+        Extract the steepest peak using a seed point.
+
+        >>> data_min, data_max = mesh.get_data_range()
+        >>> peak_range = [0.2, data_max]  # Peak if curvature > 0.2
+        >>> peak_point_id = np.argmax(mesh['Curvature'])  # Steepest point
+        >>> peak_mesh = mesh.connectivity(
+        ...     'point_seed', peak_point_id, scalar_range=peak_range
+        ... )
+
+        Extract the valley region closest to the steepest peak.
+
+        >>> valley_range = [data_min, -0.2]  # Valley if curvature < 0.2
+        >>> peak_point = mesh.points[peak_point_id]
+        >>> valley_mesh = mesh.connectivity(
+        ...     'closest', peak_point, scalar_range=valley_range
+        ... )
+
+        Plot extracted regions.
+
+        >>> p = pv.Plotter()
+        >>> _ = p.add_mesh(mesh, style='wireframe', color='lightgray')
+        >>> _ = p.add_mesh(peak_mesh, color='red', label='Steepest Peak')
+        >>> _ = p.add_mesh(
+        ...     valley_mesh, color='blue', label='Closest Valley'
+        ... )
+        >>> _ = p.add_legend()
+        >>> p.show()
 
         """
+        # Deprecated on v0.43.0
+        keep_largest = kwargs.pop('largest', False)
+        if keep_largest:  # pragma: no cover
+            warnings.warn(
+                "Use of `largest=True` is deprecated. Use 'largest' or "
+                "`extraction_mode='largest'` instead.",
+                PyVistaDeprecationWarning,
+            )
+            extraction_mode = 'largest'
+
+        def _unravel_and_validate_ids(ids):
+            ids = np.asarray(ids).ravel()
+            is_all_integers = np.issubdtype(ids.dtype, np.integer)
+            is_all_positive = not np.any(ids < 0)
+            if not (is_all_positive and is_all_integers):
+                raise ValueError('IDs must be positive integer values.')
+            return np.unique(ids)
+
+        def _extract_points(_input, points):
+            # This method is similar to `mesh.extract_points` but the
+            # output arrays are removed and the output type is PolyData
+            # if the input type is PolyData
+            _output = _input.extract_points(points, progress_bar=progress_bar)
+            if isinstance(_input, pyvista.PolyData):
+                # Output is UnstructuredGrid, so apply vtkRemovePolyData
+                # to input to make the output as PolyData type instead
+                all_ids = set(range(_input.n_cells))
+                ids_to_keep = set(_output['vtkOriginalCellIds'])
+                ids_to_remove = list(all_ids - ids_to_keep)
+                if len(ids_to_remove) == 0:
+                    _output = _input
+                else:
+                    if pyvista.vtk_version_info < (9, 1, 0):
+                        raise VTKVersionError(
+                            '`connectivity` with PolyData requires vtk>=9.1.0'
+                        )  # pragma: no cover
+                    remove = _vtk.vtkRemovePolyData()
+                    remove.SetInputData(_input)
+                    remove.SetCellIds(numpy_to_idarr(ids_to_remove))
+                    _update_alg(remove, progress_bar, "Removing Cells.")
+                    _output = _get_output(remove)
+                    _output.clean(
+                        point_merging=False, inplace=True, progress_bar=progress_bar
+                    )  # remove unused points
+            _output.point_data.remove('vtkOriginalPointIds')
+            _output.cell_data.remove('vtkOriginalCellIds')
+            return _output
+
+        # Store active scalars info to restore later if needed
+        active_field, active_name = self.active_scalars_info  # type: ignore
+
+        # Set scalars
+        if scalar_range is None:
+            input_mesh = self.copy(deep=False)  # type: ignore
+        else:
+            if scalar_range[0] > scalar_range[1]:
+                raise ValueError(
+                    f"Lower value of scalar range {scalar_range[0]} cannot be greater than the upper value {scalar_range[0]}"
+                )
+
+            # Input will be modified, so copy first
+            input_mesh = self.copy()  # type: ignore
+            if scalars is None:
+                set_default_active_scalars(input_mesh)
+            else:
+                input_mesh.set_active_scalars(scalars)
+            # Make sure we have point data (required by the filter)
+            field, name = input_mesh.active_scalars_info
+            if field == FieldAssociation.CELL:
+                # Convert to point data with a unique name
+                # The point array will be removed later
+                point_data = input_mesh.cell_data_to_point_data(progress_bar=progress_bar)[name]
+                input_mesh.point_data['__point_data'] = point_data
+                input_mesh.set_active_scalars('__point_data')
+
+            if extraction_mode in ['all', 'specified', 'closest']:
+                # Scalar connectivity has no effect if SetExtractionModeToAllRegions
+                # (which applies to 'all' and 'specified') and 'closest'
+                # can sometimes fail for some datasets/scalar values.
+                # So, we filter scalar values beforehand
+                if scalar_range is not None:
+                    # Use extract_points to ensure that cells with at least one
+                    # point within the range are kept (this is consistent
+                    # with how the filter operates for other modes)
+                    lower = input_mesh.active_scalars >= scalar_range[0]
+                    upper = input_mesh.active_scalars <= scalar_range[1]
+                    input_mesh = _extract_points(input_mesh, np.logical_and(upper, lower))
+
         alg = _vtk.vtkConnectivityFilter()
-        alg.SetInputData(self)
-        if largest:
+        alg.SetInputDataObject(input_mesh)
+
+        # Due to inconsistent/buggy output, always keep this on and
+        # remove scalars later as needed
+        alg.ColorRegionsOn()  # This will create 'RegionId' scalars
+
+        # Sort region ids
+        alg.SetRegionIdAssignmentMode(alg.CELL_COUNT_DESCENDING)
+
+        if scalar_range is not None:
+            alg.ScalarConnectivityOn()
+            alg.SetScalarRange(*scalar_range)
+
+        if extraction_mode == 'all':
+            alg.SetExtractionModeToAllRegions()
+
+        elif extraction_mode == 'largest':
             alg.SetExtractionModeToLargestRegion()
-        else:
+
+        elif extraction_mode == 'specified':
+            if region_ids is None:
+                if variable_input is None:
+                    raise ValueError(
+                        "`region_ids` must be specified when `extraction_mode='specified'`."
+                    )
+                else:
+                    region_ids = variable_input
+            # this mode returns scalar data with shape that may not match
+            # the number of cells/points, so we extract all and filter later
+            # alg.SetExtractionModeToSpecifiedRegions()
+            region_ids = _unravel_and_validate_ids(region_ids)
+            # [alg.AddSpecifiedRegion(i) for i in region_ids]
             alg.SetExtractionModeToAllRegions()
-        alg.SetColorRegions(True)
-        _update_alg(alg, progress_bar, 'Finding and Labeling Connected Bodies/Volumes.')
+
+        elif extraction_mode == 'cell_seed':
+            if cell_ids is None:
+                if variable_input is None:
+                    raise ValueError(
+                        "`cell_ids` must be specified when `extraction_mode='cell_seed'`."
+                    )
+                else:
+                    cell_ids = variable_input
+            alg.SetExtractionModeToCellSeededRegions()
+            alg.InitializeSeedList()
+            [alg.AddSeed(i) for i in _unravel_and_validate_ids(cell_ids)]
+
+        elif extraction_mode == 'point_seed':
+            if point_ids is None:
+                if variable_input is None:
+                    raise ValueError(
+                        "`point_ids` must be specified when `extraction_mode='point_seed'`."
+                    )
+                else:
+                    point_ids = variable_input
+            alg.SetExtractionModeToPointSeededRegions()
+            alg.InitializeSeedList()
+            [alg.AddSeed(i) for i in _unravel_and_validate_ids(point_ids)]
+
+        elif extraction_mode == 'closest':
+            if closest_point is None:
+                if variable_input is None:
+                    raise ValueError(
+                        "`closest_point` must be specified when `extraction_mode='closest'`."
+                    )
+                else:
+                    closest_point = variable_input
+            alg.SetExtractionModeToClosestPointRegion()
+            alg.SetClosestPoint(*closest_point)
+
+        else:
+            raise ValueError(
+                f"Invalid value for `extraction_mode` '{extraction_mode}'. Expected one of the following: 'all', 'largest', 'specified', 'cell_seed', 'point_seed', or 'closest'"
+            )
+
+        _update_alg(alg, progress_bar, 'Finding and Labeling Connected Regions.')
         output = _get_output(alg)
 
-        # remove regionID from the output when extraction mode is set to the
-        # largest to avoid the VTK warning:
-        # the Cell array RegionId ... has X tuples but there are only Y cells
-        if largest:
-            output.cell_data.pop('RegionId', None)
+        # Process output
+        output_needs_fixing = False  # initialize flag if output needs to be fixed
+        if extraction_mode == 'all':
+            pass  # Output is good
+        elif extraction_mode == 'specified':
+            # All regions were initially extracted, so extract only the
+            # specified regions
+            point_scalars = output.point_data['RegionId']
+            point_id_mask = np.zeros_like(point_scalars)
+            for i in region_ids:
+                point_id_mask[point_scalars == i] = True
+            output = _extract_points(output, np.nonzero(point_id_mask)[0])
+
+            if label_regions:
+                # Tresholded regions may not be contiguous and zero-based
+                # which will need to be fixed
+                output_needs_fixing = True
+
+        elif extraction_mode == 'largest' and isinstance(output, pyvista.PolyData):
+            # PolyData with 'largest' mode generates bad output with unreferenced points
+            output_needs_fixing = True
+
+        else:
+            # All other extraction modes / cases may generate incorrect scalar arrays
+            # e.g. 'largest' may output scalars with shape that does not match output mesh
+            # e.g. 'seed' method scalars may have one RegionId, yet may contain many
+            # disconnected regions. Therefore, check for correct scalars size
+            if label_regions:
+                invalid_cell_scalars = output.n_cells != output.cell_data['RegionId'].size
+                invalid_point_scalars = output.n_points != output.point_data['RegionId'].size
+                if invalid_cell_scalars or invalid_point_scalars:
+                    output_needs_fixing = True
+
+        if output_needs_fixing:
+            # Fix bad output recursively using 'all' mode which has known good output
+            output.point_data.remove('RegionId')
+            output.cell_data.remove('RegionId')
+            output = output.connectivity('all', label_regions=True, inplace=inplace)
+
+        # Remove temp point array
+        try:
+            output.point_data.remove('__point_data')
+        except KeyError:
+            pass
+
+        if not label_regions:
+            output.point_data.remove('RegionId')
+            output.cell_data.remove('RegionId')
+
+            # restore previously active scalars
+            output.set_active_scalars(active_name, preference=active_field)
+
+        if inplace:
+            try:
+                self.copy_from(output, deep=False)  # type: ignore
+                return self
+            except:
+                pass
         return output
 
     def extract_largest(self, inplace=False, progress_bar=False):
@@ -2454,18 +2879,18 @@ class DataSetFilters:
         >>> import pyvista
         >>> mesh = pyvista.Sphere() + pyvista.Cube()
         >>> largest = mesh.extract_largest()
-        >>> largest.point_data.clear()
-        >>> largest.cell_data.clear()
         >>> largest.plot()
 
-        See :ref:`volumetric_example` for more examples using this filter.
+        See :ref:`connectivity_example` and :ref:`volumetric_example` for
+        more examples using this filter.
+
+        .. seealso::
+            :func:`pyvista.DataSetFilter.connectivity`
 
         """
-        mesh = DataSetFilters.connectivity(self, largest=True, progress_bar=progress_bar)
-        if inplace:
-            self.copy_from(mesh, deep=False)
-            return self
-        return mesh
+        return DataSetFilters.connectivity(
+            self, 'largest', label_regions=False, inplace=inplace, progress_bar=progress_bar
+        )
 
     def split_bodies(self, label=False, progress_bar=False):
         """Find, label, and split connected bodies/volumes.
diff --git a/pyvista/core/utilities/geometric_objects.py b/pyvista/core/utilities/geometric_objects.py
index 3185a677..b4710bc8 100644
--- a/pyvista/core/utilities/geometric_objects.py
+++ b/pyvista/core/utilities/geometric_objects.py
@@ -62,7 +62,7 @@ def translate(surf, center=(0.0, 0.0, 0.0), direction=(1.0, 0.0, 0.0)):
     if np.array_equal(normz, (0.0, 0.0, 0.0)):
         # the assumed normy axis is parallel to normx, so shift its
         # axis and recalculate normz
-        norm_y_temp = np.roll(norm_y_temp, 1)
+        norm_y_temp = [-1.0, 0.0, 0.0]
         normz = np.cross(normx, norm_y_temp)
     normz /= np.linalg.norm(normz)
     normy = np.cross(normz, normx)


Test patch:
diff --git a/tests/core/test_dataset_filters.py b/tests/core/test_dataset_filters.py
index 95bd4f06..eab48cda 100644
--- a/tests/core/test_dataset_filters.py
+++ b/tests/core/test_dataset_filters.py
@@ -1003,7 +1003,343 @@ def test_glyph_orient_and_scale():
     assert glyph4.bounds[0] == geom.bounds[0] and glyph4.bounds[1] == geom.bounds[1]
 
 
-def test_split_and_connectivity():
+@pytest.fixture
+def connected_datasets():
+    # This is similar to the datasets fixture, but the PolyData is fully connected
+    return [
+        examples.load_uniform(),  # ImageData
+        examples.load_rectilinear(),  # RectilinearGrid
+        examples.load_hexbeam(),  # UnstructuredGrid
+        pv.Sphere(),  # PolyData
+        examples.load_structured(),  # StructuredGrid
+    ]
+
+
+@pytest.fixture
+def foot_bones():
+    return examples.download_foot_bones()
+
+
+@pytest.fixture
+def connected_datasets_single_disconnected_cell(connected_datasets):
+    # Create datasets of MultiBlocks with either 'point' or 'cell' scalar
+    # data, where a single cell (or its points) has a completely different
+    # scalar value to all other cells in the dataset
+    for i, dataset in enumerate(connected_datasets):
+        dataset.clear_data()
+        dataset_composite = pv.MultiBlock()
+        single_cell_id = dataset.n_cells - 1
+        single_cell_point_ids = dataset.get_cell(single_cell_id).point_ids
+        for association in ['point', 'cell']:
+            # Add copy as block
+            dataset_copy = dataset.copy()
+            dataset_composite[association] = dataset_copy
+
+            # Make scalar data such that values for a single cell or its points
+            # fall outside of the extraction range
+
+            if association == 'point':
+                num_scalars = dataset_copy.n_points
+                node_dataset = dataset_copy.point_data
+                ids = single_cell_point_ids
+            else:
+                num_scalars = dataset_copy.n_cells
+                node_dataset = dataset_copy.cell_data
+                ids = single_cell_id
+
+            # Assign non-zero floating point scalar data with range [-10, 10]
+            # which is intended to be extracted later
+            scalar_data = np.linspace(-10, 10, num=num_scalars)
+            node_dataset['data'] = scalar_data
+
+            # Assign a totally different scalar value to the single node
+            node_dataset['data'][ids] = -1000.0
+
+        connected_datasets[i] = dataset_composite
+    return connected_datasets
+
+
+@pytest.mark.parametrize('dataset_index', list(range(5)))
+@pytest.mark.parametrize(
+    'extraction_mode', ['all', 'largest', 'specified', 'cell_seed', 'point_seed', 'closest']
+)
+@pytest.mark.parametrize('label_regions', [True, False])
+@pytest.mark.parametrize('scalar_range', [True, False])
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_inplace_and_output_type(
+    datasets, dataset_index, extraction_mode, label_regions, scalar_range
+):
+    # paramaterize with label_regions and scalar_range as these parameters
+    # have branches which may modify input/input type
+    dataset = datasets[dataset_index]
+
+    # ensure we have scalars and set a restricted range
+    if scalar_range:
+        if len(dataset.array_names) == 0:
+            dataset.point_data['data'] = np.arange(0, dataset.n_points)
+        pv.set_default_active_scalars(dataset)
+        scalar_range = [np.mean(dataset.active_scalars), np.max(dataset.active_scalars)]
+    else:
+        scalar_range = None
+
+    common_args = dict(
+        extraction_mode=extraction_mode,
+        point_ids=0,
+        cell_ids=0,
+        region_ids=0,
+        closest_point=(0, 0, 0),
+        label_regions=label_regions,
+        scalar_range=scalar_range,
+    )
+    conn = dataset.connectivity(inplace=False, **common_args)
+    assert conn is not dataset
+
+    conn = dataset.connectivity(inplace=True, **common_args)
+    if isinstance(dataset, (pv.UnstructuredGrid, pv.PolyData)):
+        assert conn is dataset
+    else:
+        assert conn is not dataset
+
+    # test correct output type
+    if isinstance(dataset, pv.PolyData):
+        assert isinstance(conn, pv.PolyData)
+    else:
+        assert isinstance(conn, pv.UnstructuredGrid)
+
+
+@pytest.mark.parametrize('dataset_index', list(range(5)))
+@pytest.mark.parametrize(
+    'extraction_mode', ['all', 'largest', 'specified', 'cell_seed', 'point_seed', 'closest']
+)
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_label_regions(datasets, dataset_index, extraction_mode):
+    # the connectivity filter is known to output incorrectly sized scalars
+    # test all modes and datasets for correct scalar size
+    dataset = datasets[dataset_index]
+    common_args = dict(
+        extraction_mode=extraction_mode,
+        point_ids=0,
+        cell_ids=0,
+        region_ids=0,
+        closest_point=(0, 0, 0),
+    )
+    conn = dataset.connectivity(**common_args, label_regions=True)
+    assert 'RegionId' in conn.point_data.keys()
+    assert 'RegionId' in conn.cell_data.keys()
+
+    expected_cell_scalars_size = conn.n_cells
+    actual_cell_scalars_size = conn.cell_data['RegionId'].size
+    assert expected_cell_scalars_size == actual_cell_scalars_size
+
+    expected_point_scalars_size = conn.n_points
+    actual_point_scalars_size = conn.point_data['RegionId'].size
+    assert expected_point_scalars_size == actual_point_scalars_size
+
+    # test again but without labels
+    active_scalars_info = dataset.active_scalars_info
+    conn = dataset.connectivity(**common_args, label_regions=False)
+    assert 'RegionId' not in conn.point_data.keys()
+    assert 'RegionId' not in conn.cell_data.keys()
+
+    assert conn.n_cells == expected_cell_scalars_size
+    assert conn.n_points == expected_point_scalars_size
+
+    # test previously active scalars are restored
+    assert conn.active_scalars_info[0] == active_scalars_info[0]
+    assert conn.active_scalars_info[1] == active_scalars_info[1]
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_raises(
+    connected_datasets_single_disconnected_cell,
+):
+    dataset = connected_datasets_single_disconnected_cell[0]['point']
+    with pytest.raises(ValueError, match='Lower value'):
+        dataset.connectivity(scalar_range=[1, 0])
+
+    with pytest.raises(ValueError, match='Invalid value for `extraction_mode`'):
+        dataset.connectivity(extraction_mode='foo')
+
+    with pytest.raises(ValueError, match='`closest_point` must be specified'):
+        dataset.connectivity(extraction_mode='closest')
+
+    with pytest.raises(ValueError, match='`point_ids` must be specified'):
+        dataset.connectivity(extraction_mode='point_seed')
+
+    with pytest.raises(ValueError, match='`cell_ids` must be specified'):
+        dataset.connectivity(extraction_mode='cell_seed')
+
+    with pytest.raises(ValueError, match='`region_ids` must be specified'):
+        dataset.connectivity(extraction_mode='specified')
+
+    with pytest.raises(ValueError, match='positive integer values'):
+        dataset.connectivity(extraction_mode='cell_seed', cell_ids=[-1, 2])
+
+
+@pytest.mark.parametrize('dataset_index', list(range(5)))
+@pytest.mark.parametrize(
+    'extraction_mode', ['all', 'largest', 'specified', 'cell_seed', 'point_seed', 'closest']
+)
+@pytest.mark.parametrize('association', ['cell', 'point'])
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_scalar_range(
+    connected_datasets_single_disconnected_cell, dataset_index, extraction_mode, association
+):
+    dataset = connected_datasets_single_disconnected_cell[dataset_index][association]
+
+    common_args = dict(
+        extraction_mode=extraction_mode,
+        point_ids=dataset.get_cell(0).point_ids[0],
+        cell_ids=0,
+        region_ids=0,
+        closest_point=dataset.get_cell(0).points[0],
+        label_regions=True,
+    )
+
+    # test a single cell is removed
+    conn_no_range = dataset.connectivity(**common_args)
+    conn_with_range = dataset.connectivity(**common_args, scalar_range=[-10, 10])
+    assert conn_with_range.n_cells == conn_no_range.n_cells - 1
+
+    # test no cells are removed
+    conn_with_full_range = dataset.connectivity(
+        **common_args, scalar_range=dataset.get_data_range()
+    )
+    assert conn_with_full_range.n_cells == dataset.n_cells
+
+    # test input scalars are passed to output
+    assert len(conn_no_range.array_names) == 3  # ['data', 'RegionId', 'RegionId']
+    assert len(conn_with_range.array_names) == 3
+    assert len(conn_with_full_range.array_names) == 3
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_all(foot_bones):
+    conn = foot_bones.connectivity('all')
+    assert conn.n_cells == foot_bones.n_cells
+
+    # test correct labels
+    conn = foot_bones.connectivity('all', label_regions=True)
+    region_ids, counts = np.unique(conn.cell_data['RegionId'], return_counts=True)
+    assert np.array_equal(region_ids, list(range(26)))
+    assert np.array_equal(
+        counts,
+        [
+            598,
+            586,
+            392,
+            360,
+            228,
+            212,
+            154,
+            146,
+            146,
+            146,
+            134,
+            134,
+            134,
+            126,
+            124,
+            74,
+            66,
+            60,
+            60,
+            60,
+            48,
+            46,
+            46,
+            46,
+            46,
+            32,
+        ],
+    )
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_largest(foot_bones):
+    conn = foot_bones.connectivity('largest')
+    assert conn.n_cells == 598
+
+    # test correct labels
+    conn = foot_bones.connectivity('largest', label_regions=True)
+    region_ids, counts = np.unique(conn.cell_data['RegionId'], return_counts=True)
+    assert region_ids == [0]
+    assert counts == [598]
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_specified(foot_bones):
+    # test all regions
+    all_regions = list(range(26))
+    conn = foot_bones.connectivity('specified', region_ids=all_regions)
+    assert conn.n_cells == foot_bones.n_cells
+
+    # test irrelevant region IDs
+    test_regions = all_regions + [77, 99]
+    conn = foot_bones.connectivity('specified', test_regions)
+    assert conn.n_cells == foot_bones.n_cells
+
+    # test some regions
+    some_regions = [1, 2, 4, 5]
+    expected_n_cells = 586 + 392 + 228 + 212
+    conn = foot_bones.connectivity('specified', some_regions)
+    assert conn.n_cells == expected_n_cells
+
+    # test correct labels
+    conn = foot_bones.connectivity('specified', some_regions, label_regions=True)
+    region_ids, counts = np.unique(conn.cell_data['RegionId'], return_counts=True)
+    assert np.array_equal(region_ids, [0, 1, 2, 3])
+    assert np.array_equal(counts, [586, 392, 228, 212])
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_point_seed(foot_bones):
+    conn = foot_bones.connectivity('point_seed', point_ids=1598)
+    assert conn.n_cells == 598
+    conn = foot_bones.connectivity('point_seed', [1326, 1598])
+    assert conn.n_cells == 598 + 360
+    assert conn.n_points == 301 + 182
+
+    # test correct labels
+    conn = foot_bones.connectivity('point_seed', [1326, 1598], label_regions=True)
+    region_ids, counts = np.unique(conn.cell_data['RegionId'], return_counts=True)
+    assert np.array_equal(region_ids, [0, 1])
+    assert np.array_equal(counts, [598, 360])
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_cell_seed(foot_bones):
+    conn = foot_bones.connectivity('cell_seed', cell_ids=3122)
+    assert conn.n_cells == 598
+    assert conn.n_points == 301
+    conn = foot_bones.connectivity('cell_seed', [2588, 3122])
+    assert conn.n_cells == 598 + 360
+    assert conn.n_points == 301 + 182
+
+    # test correct labels
+    conn = foot_bones.connectivity('point_seed', [1326, 1598], label_regions=True)
+    region_ids, counts = np.unique(conn.cell_data['RegionId'], return_counts=True)
+    assert np.array_equal(region_ids, [0, 1])
+    assert np.array_equal(counts, [598, 360])
+
+
+@pytest.mark.needs_vtk_version(9, 1, 0)
+def test_connectivity_closest_point(foot_bones):
+    conn = foot_bones.connectivity('closest', closest_point=(-3.5, -0.5, -0.5))
+    assert conn.n_cells == 598
+    assert conn.n_points == 301
+    conn = foot_bones.connectivity('closest', (-1.5, -0.5, 0.05))
+    assert conn.n_cells == 360
+    assert conn.n_points == 182
+
+    # test correct labels
+    conn = foot_bones.connectivity('closest', (-3.5, -0.5, -0.5), label_regions=True)
+    region_ids, counts = np.unique(conn.cell_data['RegionId'], return_counts=True)
+    assert region_ids == [0]
+    assert counts == [598]
+
+
+def test_split_bodies():
     # Load a simple example mesh
     dataset = examples.load_uniform()
     dataset.set_active_scalars('Spatial Cell Data')
diff --git a/tests/core/test_parametric_geometry.py b/tests/core/test_parametric_geometry.py
index 925c90f8..ed219166 100644
--- a/tests/core/test_parametric_geometry.py
+++ b/tests/core/test_parametric_geometry.py
@@ -148,3 +148,15 @@ def test_ParametricSuperToroid():
 def test_ParametricTorus():
     geom = pv.ParametricTorus()
     assert geom.n_points
+
+
+def test_direction():
+    geom1 = pv.ParametricEllipsoid(300, 100, 10, direction=[1, 0, 0])
+    geom2 = pv.ParametricEllipsoid(300, 100, 10, direction=[0, 1, 0])
+    assert geom1.n_points
+    assert geom2.n_points
+    points1 = geom1.points
+    points2 = geom2.points
+    assert np.allclose(points1[:, 0], points2[:, 1])
+    assert np.allclose(points1[:, 1], -points2[:, 0])
+    assert np.allclose(points1[:, 2], points2[:, 2])