Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update and fix connectivity docs #4872

Merged
merged 8 commits into from
Sep 15, 2023
Merged
Show file tree
Hide file tree
Changes from 7 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
5 changes: 3 additions & 2 deletions examples/01-filter/compute-volume.py
Original file line number Diff line number Diff line change
Expand Up @@ -77,10 +77,11 @@
###############################################################################
# Or better yet, you could simply extract the largest volume from your
# dataset directly by passing ``'largest'`` to the ``connectivity`` and
# specify the same scalar range used earlier for thresholding.
# specifying the scalar range of interest.

# Grab the largest connected volume within a scalar range
largest = threshed.connectivity('largest', scalar_range=[0.15, 0.50])
scalar_range = [0, 77] # Range corresponding to bottom 15% of values
largest = threshed.connectivity('largest', scalar_range=scalar_range)

# Get volume as numeric value
large_volume = largest.volume
Expand Down
158 changes: 143 additions & 15 deletions examples/01-filter/connectivity.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,39 +4,167 @@
Connectivity
~~~~~~~~~~~~

Use the connectivity filter to remove noisy isosurfaces.
This example highlights some applications of the :func:`~pyvista.DataSetFilters.connectivity`
filter.

This example is very similar to `this VTK example <https://kitware.github.io/vtk-examples/site/Python/VisualizationAlgorithms/PineRootConnectivity/>`__
"""

###############################################################################
# Remove Noisy Isosurfaces
# ~~~~~~~~~~~~~~~~~~~~~~~~
#
# Use connectivity to remove noisy isosurfaces.
#
# This section is similar to `this VTK example <https://kitware.github.io/vtk-examples/site/Python/VisualizationAlgorithms/PineRootConnectivity/>`__.

# sphinx_gallery_thumbnail_number = 2
import numpy as np

import pyvista as pv
from pyvista import examples

###############################################################################
# Load a dataset that has noisy isosurfaces
mesh = examples.download_pine_roots()

cpos = [(40.6018, -280.533, 47.0172), (40.6018, 37.2813, 50.1953), (0.0, 0.0, 1.0)]
# Load a dataset with noisy isosurfaces.
pine_roots = examples.download_pine_roots()

# Plot the raw data
cpos = [(40.6018, -280.533, 47.0172), (40.6018, 37.2813, 50.1953), (0.0, 0.0, 1.0)]
p = pv.Plotter()
p.add_mesh(mesh, color='#965434')
p.add_mesh(mesh.outline())
p.add_mesh(pine_roots, color='#965434')
p.add_mesh(pine_roots.outline())
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'`` to the ``connectivity``
# filter or by using the :func:`pyvista.DataSetFilters.extract_largest` filter
# (both are equivalent).
# The plotted mesh is very noisy. We can extract the largest connected
# isosurface using the ``'largest'`` ``extraction_mode`` of the
# :func:`~pyvista.DataSetFilters.connectivity` filter. Equivalently,
# :func:`~pyvista.DataSetFilters.extract_largest` may also be used.

# Grab the largest connected volume present
largest = mesh.connectivity('largest')
largest = pine_roots.connectivity('largest')
# or: largest = mesh.extract_largest()

p = pv.Plotter()
p.add_mesh(largest, color='#965434')
p.add_mesh(mesh.outline())
p.add_mesh(pine_roots.outline())
p.camera_position = cpos
p.show()


###############################################################################
# Extract Small Regions
# ~~~~~~~~~~~~~~~~~~~~~
#
# Use connectivity to extract the smaller 'noisy' regions that were
# removed in the remove noisy isosurfaces example above.
#
# First, get a list of all region ids.
all_regions = pine_roots.connectivity('all')
region_ids = np.unique(all_regions['RegionId'])

###############################################################################
# Since the region IDs are sorted in descending order (by cell count),
# we can extract all regions *except* for the largest one using the
# ``'specified'`` ``extraction_mode`` of the :func:`~pyvista.DataSetFilters.connectivity`
# filter.
noise_region_ids = region_ids[1::] # All region ids except '0'
noise = pine_roots.connectivity('specified', noise_region_ids)

###############################################################################
# Plot the noisy regions. For context, also plot the largest region.
p = pv.Plotter()
p.add_mesh(noise, cmap='glasbey', categories=True)
p.add_mesh(largest, color='lightgray', opacity=0.1)
p.add_mesh(pine_roots.outline())
p.camera_position = cpos
p.show()


###############################################################################
# Label Disconnected Regions
# ~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Use connectivity to label all disconnected regions.
#
# This section is similar to `this VTK example <https://examples.vtk.org/site/Cxx/PolyData/ColorDisconnectedRegionsDemo/>`__.
#
# First, load a dataset with disconnected regions.
mesh = examples.download_foot_bones()

###############################################################################
# Extract all regions.
conn = mesh.connectivity('all')

###############################################################################
# Plot the labeled regions.

# Format scalar bar text for integer values.
scalar_bar_args = dict(
fmt='%.f',
)

cpos = [(10.5, 12.2, 18.3), (0.0, 0.0, 0.0), (0.0, 1.0, 0.0)]

conn.plot(
categories=True,
cmap='glasbey',
scalar_bar_args=scalar_bar_args,
cpos=cpos,
)


###############################################################################
# Extract Regions From Seed Points
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Use connectivity to extract regions of interest using scalar data and
# seed points.
#
# First, create a dataset with salient features. Here, we create hills
# and use curvature to define its peaks and valleys.
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',
)

###############################################################################
# Extract a region of interest using the
# ``'point_seed'`` ``extraction_mode`` of the :func:`~pyvista.DataSetFilters.connectivity`
# filter. Let's extract the steepest peak using a seed point where the
# curvature is maximized.

# Get seed point
peak_point_id = np.argmax(mesh['Curvature'])

# Define a scalar range of the region to extract
data_min, data_max = mesh.get_data_range()
peak_range = [0.2, data_max] # Peak if curvature > 0.2

peak_mesh = mesh.connectivity('point_seed', peak_point_id, scalar_range=peak_range)

###############################################################################
# Let's also extract the closest valley to the steepest peak using the
# ``'closest'`` ``extraction_mode`` of the :func:`~pyvista.DataSetFilters.connectivity`
# filter.
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()