import math
from typing import Optional, Tuple, Union
import numpy as np
import pyvista as pv
from pyvista import MultiBlock, PolyData, UnstructuredGrid
try:
from typing import Literal
except ImportError:
from typing_extensions import Literal
from ...logging import logger_manager as lm
from ..models import collect_models, multiblock2model
from .utils import _interactive_plotter
###############
# Create line #
###############
[docs]def find_plane_equation(point1: np.ndarray, point2: np.ndarray, point3: np.ndarray):
xo1, yo1, zo1 = point1
xo2, yo2, zo2 = point2
xo3, yo3, zo3 = point3
a = (yo2 - yo1) * (zo3 - zo1) - (zo2 - zo1) * (yo3 - yo1)
b = (xo3 - xo1) * (zo2 - zo1) - (xo2 - xo1) * (zo3 - zo1)
c = (xo2 - xo1) * (yo3 - yo1) - (xo3 - xo1) * (yo2 - yo1)
d = -(a * xo1 + b * yo1 + c * zo1)
equation_args = np.array([a, b, c, d])
return equation_args
[docs]def find_model_outline_planes(model) -> dict:
x1, x2, y1, y2, z1, z2 = model.bounds
vertices = np.asarray(
[
[x1, y1, z1],
[x1, y1, z2],
[x1, y2, z1],
[x1, y2, z2],
[x2, y1, z1],
[x2, y1, z2],
[x2, y2, z1],
[x2, y2, z2],
]
)
x_plane = find_plane_equation(point1=vertices[0, :], point2=vertices[1, :], point3=vertices[2, :])
x_plane_opposite = find_plane_equation(point1=vertices[4, :], point2=vertices[5, :], point3=vertices[6, :])
y_plane = find_plane_equation(point1=vertices[0, :], point2=vertices[1, :], point3=vertices[4, :])
y_plane_opposite = find_plane_equation(point1=vertices[2, :], point2=vertices[3, :], point3=vertices[6, :])
z_plane = find_plane_equation(point1=vertices[0, :], point2=vertices[2, :], point3=vertices[4, :])
z_plane_opposite = find_plane_equation(point1=vertices[1, :], point2=vertices[3, :], point3=vertices[5, :])
planes = {
"x": [x_plane, x_plane_opposite],
"y": [y_plane, y_plane_opposite],
"z": [z_plane, z_plane_opposite],
}
return planes
[docs]def find_intersection(model, vec, center, plane):
# Normalize the vector
normal = vec / np.linalg.norm(vec)
normal_x, normal_y, normal_z = normal
center = model.center if center is None else center
center_x, center_y, center_z = center
a, b, c, d = plane
t = (-a * center_x - b * center_y - c * center_z - d) / (a * normal_x + b * normal_y + c * normal_z)
intersection_x = normal_x * t + center_x
intersection_y = normal_y * t + center_y
intersection_z = normal_z * t + center_z
intersection = np.asarray([intersection_x, intersection_y, intersection_z])
return intersection
[docs]def euclidean_distance(instance1, instance2, dimension):
distance = 0
for i in range(dimension):
distance += (instance1[i] - instance2[i]) ** 2
return math.sqrt(distance)
[docs]def create_line(model, vec, center, n_points):
planes = find_model_outline_planes(model=model)
line_dict = {}
for key, value in planes.items():
center = model.center if center is None else center
intersection1 = find_intersection(model=model, vec=vec, center=center, plane=value[0])
intersection2 = find_intersection(model=model, vec=vec, center=center, plane=value[1])
ed = euclidean_distance(instance1=model.center, instance2=intersection2, dimension=3)
if not math.isnan(ed):
line_dict[ed] = [intersection1, intersection2]
if len(line_dict.keys()) != 0:
min_ed = np.min([i for i in line_dict.keys()])
intersection1 = line_dict[min_ed][0]
intersection2 = line_dict[min_ed][1]
# line_new = pv.Line(intersection1, intersection2, n_points - 1)
line = pv.Line(intersection1, intersection2, n_points + 1)
line_new = pv.Line(line.points[1], line.points[-2], n_points - 1)
return line_new, intersection1, intersection2
else:
raise ValueError("`vec` value is wrong. Please input another `vec` value.")
#########
# Slice #
#########
[docs]def three_d_slice(
model: Union[PolyData, UnstructuredGrid],
method: Literal["axis", "orthogonal", "line"] = "axis",
n_slices: int = 10,
axis: Literal["x", "y", "z"] = "x",
vec: Union[tuple, list] = (1, 0, 0),
center: Union[tuple, list] = None,
) -> Union[PolyData, Tuple[MultiBlock, MultiBlock, PolyData]]:
"""
Create many slices of the input dataset along a specified axis or
create three orthogonal slices through the dataset on the three cartesian planes or
slice a model along a vector direction perpendicularly.
Args:
model: Reconstructed 3D model.
method: The methods of slicing a model. Available `method` are:
* `'axis'`: Create many slices of the input dataset along a specified axis.
* `'orthogonal'`: Create three orthogonal slices through the dataset on the three cartesian planes.
This method is usually used interactively without entering a position which slices are taken.
* `'line'`: Slice a model along a vector direction perpendicularly.
n_slices: The number of slices to create along a specified axis. Only works when `method` is 'axis' or 'line'.
axis: The axis to generate the slices along. Only works when `method` is 'axis'.
vec: The vector direction. Only works when `method` is 'line'.
center: A 3-length sequence specifying the position which slices are taken. Defaults to the center of the model.
Returns:
If method is 'axis' or 'orthogonal', return a MultiBlock that contains all models you sliced; else return a
tuple that contains line model, all models you sliced and intersections of slices model and line model.
"""
# Check input model.
model = multiblock2model(model=model, message="slicing") if isinstance(model, MultiBlock) else model.copy()
center = model.center if center is None else center
if method == "axis":
# Create many slices of the input dataset along a specified axis.
return model.slice_along_axis(n=n_slices, axis=axis, center=center)
elif method == "orthogonal":
# Create three orthogonal slices through the dataset on the three cartesian planes.
return model.slice_orthogonal(x=center[0], y=center[1], z=center[2])
elif method == "line":
# Slice a model along a vector direction perpendicularly.
line, intersection1, intersection2 = create_line(model=model, vec=vec, center=center, n_points=n_slices)
slices = pv.MultiBlock()
line_points = pv.MultiBlock()
for point in line.points:
normal = vec / np.linalg.norm(vec)
point_slice = model.slice(normal=normal, origin=point)
if point_slice.n_points != 0:
slices.append(point_slice)
line_points.append(point)
lm.main_info(
f"Slice the model uniformly along the vector `vec` and generate {n_slices} slices. "
f"There are {n_slices-len(slices)} empty slices, {len(slices)} valid slices in all slices.",
indent_level=1,
)
return slices, line_points, line
else:
raise ValueError("`method` value is wrong. " "\nAvailable `method` are: `'axis'`, `'orthogonal'`, `'line'`.")
#####################
# Interactive slice #
#####################
[docs]def interactive_slice(
model: Union[PolyData, UnstructuredGrid, MultiBlock],
key: str = "groups",
method: Literal["axis", "orthogonal"] = "axis",
axis: Literal["x", "y", "z"] = "x",
) -> MultiBlock:
"""
Create a slice of the input dataset along a specified axis or
create three orthogonal slices through the dataset on the three cartesian planes.
Args:
model: Reconstructed 3D model.
key: The key under which are the labels.
method: The methods of slicing a model. Available `method` are:
* `'axis'`: Create a slice of the input dataset along a specified axis.
* `'orthogonal'`: Create three orthogonal slices through the dataset on the three cartesian planes.
axis: The axis to generate the slices along. Only works when `method` is 'axis'.
Returns:
sliced_model: A MultiBlock that contains all models you sliced.
"""
# Check input model.
model = multiblock2model(model=model, message="slicing") if isinstance(model, MultiBlock) else model.copy()
# Create an interactive window for using widgets.
p = _interactive_plotter()
p.add_mesh(model, opacity=0.2, scalars=f"{key}_rgba", rgba=True, style="surface")
# Slice a model using a slicing widget.
if method == "axis":
p.add_mesh_slice(
model,
normal=axis,
scalars=f"{key}_rgba",
rgba=True,
tubing=True,
widget_color="black",
color="black",
line_width=3.0,
)
else:
p.add_mesh_slice_orthogonal(
model,
scalars=f"{key}_rgba",
rgba=True,
tubing=True,
widget_color="black",
color="black",
line_width=3.0,
)
p.show(cpos="iso")
return collect_models(p.plane_sliced_meshes)