"""Surface class based on RESQML TriangulatedSetRepresentation class."""
# RMS and ROXAR are registered trademarks of Roxar Software Solutions AS, an Emerson company
# GOCAD is also a trademark of Emerson
import logging
log = logging.getLogger(__name__)
import numpy as np
import resqpy.crs as rqc
import resqpy.lines as rql
import resqpy.olio.intersection as meet
import resqpy.olio.triangulation as triangulate
import resqpy.olio.uuid as bu
import resqpy.olio.vector_utilities as vec
import resqpy.olio.write_hdf5 as rwh5
import resqpy.olio.xml_et as rqet
import resqpy.property as rqp
import resqpy.surface as rqs
import resqpy.surface._base_surface as rqsb
import resqpy.surface._triangulated_patch as rqstp
import resqpy.weights_and_measures as wam
from resqpy.olio.xml_namespaces import curly_namespace as ns
from resqpy.olio.zmap_reader import read_mesh
class Surface(rqsb.BaseSurface):
"""Class for RESQML triangulated set surfaces."""
resqml_type = 'TriangulatedSetRepresentation'
[docs] def __init__(self,
parent_model,
uuid = None,
point_set = None,
mesh = None,
mesh_file = None,
mesh_format = None,
tsurf_file = None,
quad_triangles = False,
title = None,
surface_role = 'map',
crs_uuid = None,
originator = None,
extra_metadata = {}):
"""Create an empty Surface object (RESQML TriangulatedSetRepresentation).
Optionally populates from xml, point set or mesh.
arguments:
parent_model (model.Model object): the model to which this surface belongs
uuid (uuid.UUID, optional): if present, the surface is initialised from an existing RESQML object with this uuid
point_set (PointSet object, optional): if present, the surface is initialised as a Delaunay
triangulation of the points in the point set; ignored if extracting from xml
mesh (Mesh object, optional): if present, the surface is initialised as a triangulation of
the mesh; ignored if extracting from xml or if point_set is present
mesh_file (string, optional): the path of an ascii file holding a mesh in RMS text or zmap+ format;
ignored if extracting from xml or point_set or mesh is present
mesh_format (string, optional): 'rms' or 'zmap'; required if initialising from mesh_file
tsurf_file (string, optional): the path of an ascii file holding details of triangles and points in GOCAD-Tsurf format;
ignored if extraction from xml or point_set or mesh is present
quad_triangles (boolean, default False): if initialising from mesh or mesh_file, each 'square'
is represented by 2 triangles if quad_triangles is False, 4 triangles if True
title (string, optional): used as the citation title for the new object, ignored if
extracting from xml
surface_role (string, default 'map'): 'map' or 'pick'; ignored if uuid is not None
crs_uuid (uuid.UUID, optional): if present and not extracting from xml, is set as the crs uuid
applicable to mesh etc. data
originator (str, optional): the name of the person creating the object; defaults to login id; ignored
when initialising from an existing RESQML object
extra_metadata (dict): items in this dictionary are added as extra metadata; ignored
when initialising from an existing RESQML object
returns:
a newly created surface object
notes:
there are 6 ways to initialise a surface object, in order of precendence:
1. extracting from xml
2. as a Delaunay triangulation of points in a PointSet
3. as a simple triangulation of a Mesh object
4. as a simple triangulation of a mesh in an ascii file
5. from a GOCAD-TSurf format file
5. as an empty surface
if an empty surface is created, set_from... methods are available to then set for one of:
- a horizontal plane
- a single triangle
- a 'sail' (a triangle wrapped onto a sphere)
- etc.
the quad_triangles option is only applied if initialising from a mesh or mesh_file that is fully
defined (ie. no NaN's)
:meta common:
"""
assert surface_role in ['map', 'pick']
self.surface_role = surface_role
self.patch_list = [] # ordered list of patches
self.crs_uuid = crs_uuid
self.triangles = None # composite triangles (all patches)
self.points = None # composite points (all patches)
self.boundaries = None # todo: read up on what this is for and look out for examples
self.represented_interpretation_root = None
self.normal_vector = None # a single derived vector that is roughly (true) normal to the surface
self.title = title
super().__init__(model = parent_model,
uuid = uuid,
title = title,
originator = originator,
extra_metadata = extra_metadata)
if self.root is not None:
pass
elif point_set is not None:
self.set_from_point_set(point_set)
elif mesh is not None:
self.set_from_mesh_object(mesh, quad_triangles = quad_triangles)
elif mesh_file and mesh_format:
self.set_from_mesh_file(mesh_file, mesh_format, quad_triangles = quad_triangles)
elif tsurf_file is not None:
self.set_from_tsurf_file(tsurf_file)
[docs] @classmethod
def from_tri_mesh(cls, tri_mesh, exclude_nans = False):
"""Create a Surface from a TriMesh.
arguments:
tri_mesh (TriMesh): the tri mesh for which an equivalent Surface is required
exclude_nans (bool, default False): if True, and tri mesh point involving a not-a-number is
excluded from the surface points, along with any triangle that has such a point as a vertex
returns:
a new Surface using the points and triangles of the tri mesh
note:
this method does not write hdf5 data nor create xml for the new Surface
"""
assert isinstance(tri_mesh, rqs.TriMesh)
surf = cls(tri_mesh.model,
crs_uuid = tri_mesh.crs_uuid,
title = tri_mesh.title,
surface_role = tri_mesh.surface_role,
extra_metadata = tri_mesh.extra_metadata if hasattr(tri_mesh, 'extra_metadata') else None)
t, p = tri_mesh.triangles_and_points()
if exclude_nans:
t, p = nan_removed_triangles_and_points(t, p)
surf.set_from_triangles_and_points(t, p)
surf.represented_interpretation_root = tri_mesh.represented_interpretation_root
return surf
[docs] @classmethod
def from_downsampling_surface(cls,
fine_surface,
point_count = None,
title = None,
target_model = None,
inherit_extra_metadata = True,
inherit_represented_interpretation = True,
convexity_parameter = 5.0):
"""Create a Surface by taking a random subset of points from an existing surface and re-triangulating.
arguments:
fine_surface (Surface): a finely triangulated existing surface
point_count (int, optional): the number of randomly selected points to keep; defaults to 1% of
the fine surface point count
title (str, optional): the title for the new (coarse) surface; defaults to the title of the
fine surface
target_model (Model, optional): the model to receive the new surface; defaults to that of the
fine surface
inherit_extra_metadata (bool, default True): if True, the coarse surface will inherit any extra
metadata that the fine surface has
inherit_represented_interpretation (bool, default True): if True and the fine surface has a
represented interpretation then it is inherited by the coarse surface
convexity_parameter (float, default 5.0): controls how likely the resulting triangulation is to be
convex; reduce to 1.0 to allow slightly more concavities; increase to 100.0 or more for very
little chance of even a slight concavity
returns:
new Surface, without its hdf5 having been written, nor xml created
note:
if the fine surface does not have more points than point count, then a copy of that surface
is returned without re-triangulating
"""
assert fine_surface is not None
if title is None:
title = fine_surface.title
inter_model = target_model is not None and target_model is not fine_surface.model
if target_model is None:
target_model = fine_surface.model
fine_t, fine_p = fine_surface.triangles_and_points()
if point_count is None:
point_count = len(fine_p) // 100
if point_count < 3:
point_count = 3
if inter_model:
target_model.copy_uuid_from_other_model(fine_surface.model, fine_surface.crs_uuid)
if inherit_represented_interpretation and fine_surface.represented_interpretation_root is not None:
ri_uuid = fine_surface.model.uuid_for_root(fine_surface.represented_interpretation_root)
target_model.copy_uuid_from_other_model(fine_surface.model, ri_uuid)
em = (fine_surface.extra_metadata
if inherit_extra_metadata and hasattr(fine_surface, 'extra_metadata') else None)
surf = cls(target_model,
crs_uuid = fine_surface.crs_uuid,
title = title,
surface_role = fine_surface.surface_role,
extra_metadata = em)
if point_count >= len(fine_p):
t, p = fine_t, fine_p
else:
p = fine_p.copy()
np.random.default_rng().shuffle(p, axis = 0)
p = p[:point_count]
t = triangulate.dt(p[:, :2], container_size_factor = convexity_parameter, algorithm = "scipy")
surf.set_from_triangles_and_points(t, p)
if inherit_represented_interpretation:
surf.represented_interpretation_root = fine_surface.represented_interpretation_root
return surf
def _load_from_xml(self):
root_node = self.root
assert root_node is not None
self.extract_patches(root_node)
ref_node = rqet.find_tag(root_node, 'RepresentedInterpretation')
if ref_node is not None:
interp_root = self.model.referenced_node(ref_node)
self.set_represented_interpretation_root(interp_root)
@property
def represented_interpretation_uuid(self):
"""Returns the uuid of the represented surface interpretation, or None."""
return rqet.uuid_for_part_root(self.represented_interpretation_root)
[docs] def set_represented_interpretation_root(self, interp_root):
"""Makes a note of the xml root of the represented interpretation."""
self.represented_interpretation_root = interp_root
[docs] def set_model(self, parent_model):
"""Associate the surface with a resqml model (does not create xml or write hdf5 data)."""
self.model = parent_model
[docs] def triangles_and_points(self):
"""Returns arrays representing combination of all the patches in the surface.
returns:
tuple (triangles, points):
triangles (int array of shape[:, 3]): integer indices into points array,
being the nodes of the corners of the triangles;
points (float array of shape[:, 3]): flat array of xyz points, indexed by triangles
:meta common:
"""
if self.triangles is not None:
return (self.triangles, self.points)
self.extract_patches(self.root)
points_offset = 0
for triangulated_patch in self.patch_list:
(t, p) = triangulated_patch.triangles_and_points()
if points_offset == 0:
self.triangles = t.copy()
self.points = p.copy()
else:
self.triangles = np.concatenate((self.triangles, t.copy() + points_offset))
self.points = np.concatenate((self.points, p.copy()))
points_offset += p.shape[0]
return (self.triangles, self.points)
[docs] def triangle_count(self):
"""Return the numner of triangles in this surface."""
self.extract_patches(self.root)
if not self.patch_list:
return 0
return np.sum([tp.triangle_count for tp in self.patch_list])
[docs] def node_count(self):
"""Return the number of nodes (points) used in this surface."""
self.extract_patches(self.root)
if not self.patch_list:
return 0
return np.sum([tp.node_count for tp in self.patch_list])
[docs] def change_crs(self, required_crs):
"""Changes the crs of the surface, also sets a new uuid if crs changed.
note:
this method is usually used to change the coordinate system for a temporary resqpy object;
to add as a new part, call write_hdf5() and create_xml() methods
"""
old_crs = rqc.Crs(self.model, uuid = self.crs_uuid)
self.crs_uuid = required_crs.uuid
if required_crs == old_crs or not self.patch_list:
log.debug(f'no crs change needed for {self.title}')
return
log.debug(f'crs change needed for {self.title} from {old_crs.title} to {required_crs.title}')
for patch in self.patch_list:
patch.triangles_and_points()
required_crs.convert_array_from(old_crs, patch.points)
patch.crs_uuid = self.crs_uuid
self.triangles = None # clear cached arrays for surface
self.points = None
self.uuid = bu.new_uuid() # hope this doesn't cause problems
assert self.root is None
[docs] def set_to_trimmed_surface(self, large_surface, xyz_box = None, xy_polygon = None):
"""Populate this (empty) surface with triangles and points which overlap with a trimming volume.
arguments:
large_surface (Surface): the larger surface, a copy of which is to be trimmed
xyz_box (numpy float array of shape (2, 3), optional): if present, a cuboid in xyz space
against which to trim the surface
xy_polygon (closed convex resqpy.lines.Polyline, optional): if present, an xy boundary
against which to trim
notes:
at least one of xyz_box or xy_polygon must be present; if both are present, a triangle
must have at least one point within both boundaries to survive the trimming;
xyz_box and xy_polygon must be in the same crs as the large surface
"""
assert xyz_box is not None or xy_polygon is not None
box = None
if xyz_box is not None:
assert xyz_box.shape == (2, 3)
# guard against reversed ranges in xyz_box
box = np.empty((2, 3), dtype = float)
box[0] = np.amin(xyz_box, axis = 0)
box[1] = np.amax(xyz_box, axis = 0)
log.debug(f'trimming surface {large_surface.title} from {large_surface.triangle_count()} triangles')
if not self.title:
self.title = str(large_surface.title) + ' trimmed'
self.crs_uuid = large_surface.crs_uuid
self.patch_list = []
for triangulated_patch in large_surface.patch_list:
trimmed_patch = rqstp.TriangulatedPatch(self.model,
patch_index = len(self.patch_list),
crs_uuid = self.crs_uuid)
trimmed_patch.set_to_trimmed_patch(triangulated_patch, xyz_box = box, xy_polygon = xy_polygon)
if trimmed_patch is not None and trimmed_patch.triangle_count > 0:
self.patch_list.append(trimmed_patch)
if len(self.patch_list):
log.debug(f'trimmed surface {self.title} has {self.triangle_count()} triangles')
else:
log.warning('surface does not intersect trimming volume')
self.uuid = bu.new_uuid()
[docs] def set_to_split_surface(self, large_surface, line, delta_xyz):
"""Populate this (empty) surface with a version of a larger surface split by an xy line.
arguments:
large_surface (Surface): the larger surface, a copy of which is to be split
line (numpy float array of shape (2, 3)): the end points of a line (z will be ignored)
delta_xyz (numpy float array of shape (2, 3)): an xyz offset to add to the points on the
right and left sides of the line repspectively
notes:
this method can be used to introduce a tear into a surface, typically to represent a horizon
being split by a planar fault with a throw
"""
assert line.shape == (2, 3)
assert delta_xyz.shape == (2, 3)
t, p = large_surface.triangles_and_points()
assert p.ndim == 2 and p.shape[1] == 3
pp = np.concatenate((p, line), axis = 0)
tp = np.empty(p.shape, dtype = int)
tp[:, 0] = len(p)
tp[:, 1] = len(p) + 1
tp[:, 2] = np.arange(len(p), dtype = int)
cw = vec.clockwise_triangles(pp, tp)
pai = np.where(cw >= 0.0, True, False) # bool mask over p
pbi = np.where(cw <= 0.0, True, False) # bool mask over p
tap = pai[t]
tbp = pbi[t]
ta = np.any(tap, axis = 1) # bool array over t
tb = np.any(tbp, axis = 1) # bool array over t
#Â here we stick the two halves together into a single patch
#Â todo: keep as two patches as required by RESQML business rules
p_combo = np.empty((0, 3))
t_combo = np.empty((0, 3), dtype = int)
for i, tab in enumerate((ta, tb)):
p_keep = np.unique(t[tab])
# note new point index for each old point that is being kept
p_map = np.full(len(p), -1, dtype = int)
p_map[p_keep] = np.arange(len(p_keep))
# copy those unique points into a trimmed points array
points_trimmed = p[p_keep].copy()
# copy selected triangles, replacing p indices with compressed indices
t_trim = t[tab]
triangles_trimmed = p_map[t_trim].copy()
assert np.all(triangles_trimmed >= 0)
assert np.all(triangles_trimmed < len(points_trimmed))
points_trimmed += np.expand_dims(delta_xyz[i], axis = 0)
t_shift = len(p_combo)
p_combo = np.concatenate((p_combo, points_trimmed), axis = 0)
triangles_trimmed += t_shift
t_combo = np.concatenate((t_combo, triangles_trimmed), axis = 0)
self.set_from_triangles_and_points(t_combo, p_combo)
[docs] def distinct_edges(self):
"""Returns a numpy int array of shape (N, 2) being the ordered node pairs of distinct edges of triangles."""
triangles, _ = self.triangles_and_points()
assert triangles is not None
unique_edges, _ = triangulate.edges(triangles)
return unique_edges
[docs] def distinct_edges_and_counts(self):
"""Returns unique edges as pairs of point indices, and a count of uses of each edge.
returns:
numpy int array of shape (N, 2), numpy int array of shape (N,)
where 2D array is list-like sorted points index pairs for unique edges
and 1D array contains corresponding edge usage count (usually 1 or 2)
notes:
first entry in each pair is always the lower of the two point indices;
for well formed surfaces, the count should everywhere be zero or one;
the function does not attempt to detect coincident points
"""
triangles, _ = self.triangles_and_points()
assert triangles is not None
return triangulate.edges(triangles)
[docs] def edge_lengths(self, required_uom = None):
"""Returns float array of shape (N, 3) being triangle edge lengths.
arguments:
required_uom (str, optional): the required length uom for the resulting edge lengths; default is crs xy units
returns:
numpy float array of shape (N, 3) where N is the number of triangles
"""
t, p = self.triangles_and_points()
crs = rqc.Crs(self.model, uuid = self.crs_uuid)
if required_uom is None:
required_uom = crs.xy_units
if crs.xy_units != required_uom or crs.z_units != required_uom:
p = p.copy()
wam.convert_lengths(p[:, :2], crs.xy_units, required_uom)
wam.convert_lengths(p[:, 2], crs.z_units, required_uom)
t_end = np.empty_like(t)
t_end[:, :2] = t[:, 1:]
t_end[:, 2] = t[:, 0]
edge_v = p[t_end, :] - p[t, :]
return vec.naive_lengths(edge_v)
[docs] def set_from_triangles_and_points(self, triangles, points):
"""Populate this (empty) Surface object from an array of triangle corner indices and an array of points."""
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_from_triangles_and_points(triangles, points)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
self.triangles = triangles.copy()
self.points = points.copy()
[docs] def set_from_point_set(self,
point_set,
convexity_parameter = 5.0,
reorient = False,
reorient_max_dip = None,
extend_with_flange = False,
flange_point_count = 11,
flange_radial_factor = 10.0,
flange_radial_distance = None,
flange_inner_ring = False,
saucer_parameter = None,
make_clockwise = False):
"""Populate this (empty) Surface object with a Delaunay triangulation of points in a PointSet object.
arguments:
point_set (PointSet): the set of points to be triangulated to form a surface
convexity_parameter (float, default 5.0): controls how likely the resulting triangulation is to be
convex; reduce to 1.0 to allow slightly more concavities; increase to 100.0 or more for very little
chance of even a slight concavity
reorient (bool, default False): if True, a copy of the points is made and reoriented to minimise the
z range (ie. z axis is approximate normal to plane of points), to enhace the triangulation
reorient_max_dip (float, optional): if present, the reorientation of perspective off vertical is
limited to this angle in degrees
extend_with_flange (bool, default False): if True, a ring of points is added around the outside of the
points before the triangulation, effectively extending the surface with a flange
flange_point_count (int, default 11): the number of points to generate in the flange ring; ignored if
extend_with_flange is False
flange_radial_factor (float, default 10.0): distance of flange points from centre of points, as a
factor of the maximum radial distance of the points themselves; ignored if extend_with_flange is False
flange_radial_distance (float, optional): if present, the minimum absolute distance of flange points from
centre of points; units are those of the crs
flange_inner_ring (bool, default False): if True, an inner ring of points, with double flange point counr,
is created at a radius just outside that of the furthest flung original point; this improves
triangulation of the extended point set when the original has a non-convex hull
saucer_parameter (float, optional): if present, and extend_with_flange is True, then a parameter
controlling the shift of flange points in a perpendicular direction away from the fault plane;
see notes for how this parameter is interpreted
make_clockwise (bool, default False): if True, the returned triangles will all be clockwise when
viewed in the direction -ve to +ve z axis; if reorient is also True, the clockwise aspect is
enforced in the reoriented space
returns:
if extend_with_flange is True, numpy bool array with a value per triangle indicating flange triangles;
if extent_with_flange is False, None
notes:
if extend_with_flange is True, then a boolean array is created for the surface, with a value per triangle,
set to False (zero) for non-flange triangles and True (one) for flange triangles; this array is
suitable for adding as a property for the surface, with indexable element 'faces';
when flange extension occurs, the radius is the greater of the values determined from the radial factor
and radial distance arguments;
the saucer_parameter is interpreted in one of two ways: (1) +ve fractoinal values between zero and one
are the fractional distance from the centre of the points to its rim at which to sample the surface for
extrapolation and thereby modify the recumbent z of flange points; 0 will usually give shallower and
smoother saucer; larger values (must be less than one) will lead to stronger and more erratic saucer
shape in flange; (2) other values between -90.0 and 90.0 are interpreted as an angle to apply out of
the plane of the original points, to give a simple (and less computationally demanding) saucer shape;
+ve angles result in the shift being in the direction of the -ve z hemisphere; -ve angles result in
the shift being in the +ve z hemisphere; in either case the direction of the shift is perpendicular
to the average plane of the original points
"""
simple_saucer_angle = None
if saucer_parameter is not None and (saucer_parameter > 1.0 or saucer_parameter < 0.0):
assert -90.0 < saucer_parameter < 90.0, f'simple saucer angle parameter must be less than 90 degrees; too big: {saucer_parameter}'
simple_saucer_angle = saucer_parameter
saucer_parameter = None
assert saucer_parameter is None or 0.0 <= saucer_parameter < 1.0
crs = rqc.Crs(self.model, uuid = point_set.crs_uuid)
p = point_set.full_array_ref()
assert p.ndim >= 2
assert p.shape[-1] == 3
p = p.reshape((-1, 3))
nan_mask = np.isnan(p)
if np.any(nan_mask):
row_mask = np.logical_not(np.any(nan_mask, axis = -1))
log.info(
f'removing {len(p) - np.count_nonzero(row_mask)} NaN points from point set {point_set.title} prior to surface triangulation'
)
p = p[row_mask, :]
if crs.xy_units == crs.z_units or not reorient:
unit_adjusted_p = p
else:
unit_adjusted_p = p.copy()
wam.convert_lengths(unit_adjusted_p[:, 2], crs.z_units, crs.xy_units)
if reorient:
p_xy, self.normal_vector, reorient_matrix = triangulate.reorient(unit_adjusted_p,
max_dip = reorient_max_dip)
else:
p_xy = unit_adjusted_p
if extend_with_flange:
if not reorient:
assert saucer_parameter is None and simple_saucer_angle is None, \
'flange saucer mode only available with reorientation active'
log.warning('extending point set with flange without reorientation')
flange_points = triangulate.surrounding_xy_ring(p_xy,
count = flange_point_count,
radial_factor = flange_radial_factor,
radial_distance = flange_radial_distance,
inner_ring = flange_inner_ring,
saucer_angle = simple_saucer_angle)
p_xy_e = np.concatenate((p_xy, flange_points), axis = 0)
if reorient:
# reorient back extenstion points into original p space
flange_points_reverse_oriented = vec.rotate_array(reorient_matrix.T, flange_points)
p_e = np.concatenate((unit_adjusted_p, flange_points_reverse_oriented), axis = 0)
else:
p_e = p_xy_e
else:
p_xy_e = p_xy
p_e = unit_adjusted_p
flange_array = None
log.debug('number of points going into dt: ' + str(len(p_xy_e)))
success = False
try:
t = triangulate.dt(p_xy_e[:, :2], container_size_factor = convexity_parameter, algorithm = "scipy")
success = True
except AssertionError:
pass
if not success:
log.warning('triangulation failed, trying again with tiny perturbation of points')
p_xy_e[:, :2] += (np.random.random((len(p_xy_e), 2)) - 0.5) * 0.001
t = triangulate.dt(p_xy_e[:, :2], container_size_factor = convexity_parameter * 1.1)
log.debug('number of triangles: ' + str(len(t)))
if make_clockwise:
triangulate.make_all_clockwise_xy(t, p_e) # modifies t in situ
if extend_with_flange:
flange_array = np.zeros(len(t), dtype = bool)
flange_array[:] = np.where(np.any(t >= len(p), axis = 1), True, False)
if saucer_parameter is not None:
_adjust_flange_z(self.model, self.crs_uuid, p_xy_e, len(p), t, flange_array, saucer_parameter)
p_e = vec.rotate_array(reorient_matrix.T, p_xy_e)
if crs.xy_units != crs.z_units and reorient:
wam.convert_lengths(p_e[:, 2], crs.xy_units, crs.z_units)
self.crs_uuid = point_set.crs_uuid
self.set_from_triangles_and_points(t, p_e)
return flange_array
[docs] def make_all_clockwise_xy(self, reorient = False):
"""Reorders cached triangles data such that all triangles are clockwise when viewed from -ve z axis.
note:
if reorient is set True, a copy of the points is first reoriented such that they lie roughly in
the xy plane, and the clockwise-ness of the triangles is effected in that space
"""
_, p = self.triangles_and_points()
if reorient:
unit_adjusted_p = self.unit_adjusted_points()
p_xy, self.normal_vector, _ = triangulate.reorient(unit_adjusted_p)
else:
p_xy = p
triangulate.make_all_clockwise_xy(self.triangles, p_xy) # modifies t in situ
# assert np.all(vec.clockwise_triangles(p_xy, self.triangles) >= 0.0)
[docs] def unit_adjusted_points(self):
"""Returns cached points or copy thereof with z values adjusted to xy units."""
_, p = self.triangles_and_points()
crs = rqc.Crs(self.model, uuid = self.crs_uuid)
if crs.xy_units == self.z_units:
return p
unit_adjusted_p = p.copy()
wam.convert_lengths(unit_adjusted_p[:, 2], crs.z_units, crs.xy_units)
return unit_adjusted_p
[docs] def normal(self):
"""Returns a vector that is roughly normal to the surface.
notes:
the result becomes more meaningless the less planar the surface is;
even for a parfectly planar surface, the result is approximate;
true normal vector is found when xy & z units differ
"""
if self.normal_vector is None:
p = self.unit_adjusted_points()
_, self.normal_vector, _ = triangulate.reorient(p)
return self.normal_vector
[docs] def set_from_irregular_mesh(self, mesh_xyz, quad_triangles = False):
"""Populate this (empty) Surface object from an untorn mesh array of shape (N, M, 3).
arguments:
mesh_xyz (numpy float array of shape (N, M, 3)): a 2D lattice of points in 3D space
quad_triangles: (boolean, optional, default False): if True, each quadrangle is represented by
4 triangles in the surface, with the mean of the 4 corner points used as a common centre node;
if False (the default), only 2 triangles are used for each quadrangle; note that the 2 triangle
mode gives a non-unique triangulated result
"""
mesh_shape = mesh_xyz.shape
assert len(mesh_shape) == 3 and mesh_shape[2] == 3
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_from_irregular_mesh(mesh_xyz, quad_triangles = quad_triangles)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_from_sparse_mesh(self, mesh_xyz):
"""Populate this (empty) Surface object from a mesh array of shape (N, M, 3) with NaNs.
arguments:
mesh_xyz (numpy float array of shape (N, M, 3)): a 2D lattice of points in 3D space, with NaNs in z
"""
mesh_shape = mesh_xyz.shape
assert len(mesh_shape) == 3 and mesh_shape[2] == 3
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_from_sparse_mesh(mesh_xyz)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_from_mesh_object(self, mesh, quad_triangles = False):
"""Populate the (empty) Surface object from a Mesh object."""
xyz = mesh.full_array_ref()
if np.any(np.isnan(xyz)):
self.set_from_sparse_mesh(xyz)
else:
self.set_from_irregular_mesh(xyz, quad_triangles = quad_triangles)
[docs] def set_from_torn_mesh(self, mesh_xyz, quad_triangles = False):
"""Populate this (empty) Surface object from a torn mesh array of shape (nj, ni, 2, 2, 3).
arguments:
mesh_xyz (numpy float array of shape (nj, ni, 2, 2, 3)): corner points of 2D faces in 3D space
quad_triangles: (boolean, optional, default False): if True, each quadrangle (face) is represented
by 4 triangles in the surface, with the mean of the 4 corner points used as a common centre node;
if False (the default), only 2 triangles are used for each quadrangle; note that the 2 triangle
mode gives a non-unique triangulated result
note:
this method uses a single patch to represent the torn surface, whereas strictly the RESQML standard
requires speparate patches where parts of a surface are completely disconnected
"""
mesh_shape = mesh_xyz.shape
assert len(mesh_shape) == 5 and mesh_shape[2:] == (2, 2, 3)
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_from_torn_mesh(mesh_xyz, quad_triangles = quad_triangles)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def column_from_triangle_index(self, triangle_index):
"""For surface freshly built from fully defined mesh, returns (j, i) for given triangle index.
argument:
triangle_index (int or numpy int array): the triangle index (or array of indices) for which column(s) is/are
being sought
returns:
pair of ints or pair of numpy int arrays: the (j0, i0) indices of the column(s) which the triangle(s) is/are
part of
notes:
this function will only work if the surface has been freshly constructed with data from a mesh without NaNs,
otherwise (None, None) will be returned;
the information needed to map from triangle to column is not persistently stored as part of a resqml surface;
if triangle_index is a numpy int array, a pair of similarly shaped numpy arrays is returned
:meta common:
"""
assert len(self.patch_list) == 1
return self.patch_list[0].column_from_triangle_index(triangle_index)
[docs] def set_to_single_cell_faces_from_corner_points(self, cp, quad_triangles = True):
"""Populates this (empty) surface to represent faces of a cell, from corner points of shape (2, 2, 2, 3)."""
assert cp.size == 24
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_to_cell_faces_from_corner_points(cp, quad_triangles = quad_triangles)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_to_multi_cell_faces_from_corner_points(self, cp, quad_triangles = True):
"""Populates this (empty) surface to represent faces of a set of cells.
From corner points of shape (N, 2, 2, 2, 3).
"""
assert cp.size % 24 == 0
cp = cp.reshape((-1, 2, 2, 2, 3))
self.patch_list = []
p_index = 0
for cell_cp in cp:
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = p_index, crs_uuid = self.crs_uuid)
tri_patch.set_to_cell_faces_from_corner_points(cell_cp, quad_triangles = quad_triangles)
self.patch_list.append(tri_patch)
p_index += 1
self.uuid = bu.new_uuid()
[docs] def cell_axis_and_polarity_from_triangle_index(self, triangle_index):
"""For surface freshly built for cell faces, returns (cell_number, face_axis, polarity).
arguments:
triangle_index (int or numpy int array): the triangle index (or array of indices) for which cell face
information is required
returns:
triple int: (cell_number, axis, polarity)
note:
if the surface was built for a single cell, the returned cell number will be zero
"""
triangles_per_face = 4 if self.patch_list[0].quad_triangles else 2
face_index = triangle_index // triangles_per_face
cell_number, remainder = divmod(face_index, 6)
axis, polarity = divmod(remainder, 2)
return cell_number, axis, polarity
[docs] def set_to_horizontal_plane(self, depth, box_xyz, border = 0.0, quad_triangles = False):
"""Populate this (empty) surface with a patch of two triangles.
Triangles define a flat, horizontal plane at a given depth.
arguments:
depth (float): z value to use in all points in the triangulated patch
box_xyz (float[2, 3]): the min, max values of x, y (&z) giving the area to be covered (z ignored)
border (float): an optional border width added around the x,y area defined by box_xyz
quad_triangles (bool, default False): if True, 4 triangles are used instead of 2
:meta common:
"""
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_to_horizontal_plane(depth, box_xyz, border = border, quad_triangles = quad_triangles)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_to_triangle(self, corners):
"""Populate this (empty) surface with a patch of one triangle."""
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_to_triangle(corners)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_to_triangle_pair(self, corners):
"""Populate this (empty) surface with a patch of two triangles.
arguments:
corners (numpy float array of shape [2, 2, 3] or [4, 3]): 4 corners in logical ordering
"""
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_to_triangle_pair(corners.reshape((4, 3)))
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_to_sail(self, n, centre, radius, azimuth, delta_theta):
"""Populate this (empty) surface with a patch representing a triangle wrapped on a sphere."""
tri_patch = rqstp.TriangulatedPatch(self.model, patch_index = 0, crs_uuid = self.crs_uuid)
tri_patch.set_to_sail(n, centre, radius, azimuth, delta_theta)
self.patch_list = [tri_patch]
self.uuid = bu.new_uuid()
[docs] def set_from_mesh_file(self, filename, format, quad_triangles = False):
"""Populate this (empty) surface from a zmap or RMS text mesh file."""
assert format in ['zmap', 'rms', 'roxar'] # 'roxar' is synonymous with 'rms'
x, y, z = read_mesh(filename, format = format)
assert x is not None and y is not None or z is not None, 'failed to read surface from zmap file'
assert x.size == y.size and x.size == z.size, 'non matching array sizes from zmap reader'
assert x.shape == y.shape and x.shape == z.shape, 'non matching array shapes from zmap reader'
xyz_mesh = np.stack((x, y, z), axis = -1)
if np.any(np.isnan(z)):
self.set_from_sparse_mesh(xyz_mesh)
else:
self.set_from_irregular_mesh(xyz_mesh, quad_triangles = quad_triangles)
[docs] def set_from_tsurf_file(self, filename):
"""Populate this (empty) surface from a GOCAD tsurf file."""
triangles, vertices = [], []
with open(filename, 'r') as fl:
lines = fl.readlines()
v_index = None
for line in lines:
if "VRTX" in line:
words = line.rstrip().split(" ")
v_i = int(words[1])
if v_index is None:
v_index = v_i
index_offset = v_index
else:
assert v_i == v_index + 1, 'Tsurf vertex indices out of sequence'
v_index = v_i
vertices.append(words[2:5])
elif "TRGL" in line:
triangles.append(line.rstrip().split(" ")[1:4])
assert len(vertices) >= 3, 'vertices missing'
assert len(triangles) > 0, 'triangles missing'
triangles = np.array(triangles, dtype = int) - index_offset
vertices = np.array(vertices, dtype = float)
assert np.all(triangles >= 0) and np.all(triangles < len(vertices)), 'triangle vertex indices out of range'
self.set_from_triangles_and_points(triangles = triangles, points = vertices)
[docs] def set_from_zmap_file(self, filename, quad_triangles = False):
"""Populate this (empty) surface from a zmap mesh file."""
self.set_from_mesh_file(filename, 'zmap', quad_triangles = quad_triangles)
[docs] def set_from_roxar_file(self, filename, quad_triangles = False):
"""Populate this (empty) surface from an RMS text mesh file."""
self.set_from_mesh_file(filename, 'rms', quad_triangles = quad_triangles)
[docs] def set_from_rms_file(self, filename, quad_triangles = False):
"""Populate this (empty) surface from an RMS text mesh file."""
self.set_from_mesh_file(filename, 'rms', quad_triangles = quad_triangles)
[docs] def vertical_rescale_points(self, ref_depth = None, scaling_factor = 1.0):
"""Modify the z values of points by rescaling.
Stretches the distance from reference depth by scaling factor.
"""
if scaling_factor == 1.0:
return
if ref_depth is None:
for patch in self.patch_list:
patch_min = np.min(patch.points[:, 2])
if ref_depth is None or patch_min < ref_depth:
ref_depth = patch_min
assert ref_depth is not None, 'no z values found for vertical rescaling of surface'
self.triangles = None # invalidate any cached triangles & points in surface object
self.points = None
for patch in self.patch_list:
patch.vertical_rescale_points(ref_depth, scaling_factor)
[docs] def line_intersection(self, line_p, line_v, line_segment = False):
"""Returns x,y,z of an intersection point of straight line with the surface, or None if no intersection found."""
t, p = self.triangles_and_points()
tp = p[t]
intersects = meet.line_triangles_intersects(line_p, line_v, tp, line_segment = line_segment)
indices = meet.intersects_indices(intersects)
if not indices or len(indices) == 0:
return None
return intersects[indices[0]]
[docs] def sample_z_at_xy_points(self, points, multiple_handling = 'any'):
"""Returns interpolated z values for an array of xy values.
arguments:
points (numpy float array of shape (..., 2 or 3)): xy points to sample surface at (z values ignored)
multiple_handling (str, default 'any'): one of 'any', 'minimum', 'maximum', 'exception'
returns:
numpy float array of shape points.shape[:-1] being z values interpolated from the surface z values
notes:
points must be in the same crs as the surface;
NaN will be set for any points that do not intersect with the surface in the xy projection;
multiple_handling argument controls behaviour when one sample point intersects surface more than
once: 'any' a random one of the intersection z values is returned; 'minimum' or 'maximum': the
numerical min or max of the z values is returned; 'exception': a ValueError is raised
"""
assert points.ndim > 1 and 2 <= points.shape[-1] <= 3
assert multiple_handling in ['any', 'minimum', 'maximum', 'exception'], \
f'invalid multiple handling mode: {multiple_handling}'
if points.shape[-1] == 3:
sample_xy = points.reshape((-1, 3))
else:
sample_xy = np.zeros((points.size // 2, 3), dtype = float)
sample_xy[:, :2] = points.reshape((-1, 2))
t, p = self.triangles_and_points()
p_list = vec.points_in_triangles_njit(sample_xy, p[t], 1)
vertical = np.array((0.0, 0.0, 1.0), dtype = float)
z = np.full(sample_xy.shape[0], np.NaN, dtype = float)
for triangle_index, p_index, _ in p_list:
# todo: replace following with cheaper trilinear interpolation, given vertical intersection line
xyz = meet.line_triangle_intersect_numba(sample_xy[p_index], vertical, p[t[triangle_index]], t_tol = 0.05)
if np.isnan(z[p_index]) or multiple_handling == 'any':
z[p_index] = xyz[2]
elif multiple_handling == 'minimum':
if xyz[2] < z[p_index]:
z[p_index] = xyz[2]
elif multiple_handling == 'maximum':
if xyz[2] > z[p_index]:
z[p_index] = xyz[2]
else:
raise ValueError(f'multiple {self.title} surface intersections at xy: {sample_xy[p_index]}')
return z.reshape(points.shape[:-1])
[docs] def normal_vectors(self, add_as_property: bool = False) -> np.ndarray:
"""Returns the normal vectors for each triangle in the surface.
arguments:
add_as_property (bool): if True, face_surface_normal_vectors_array is added as a property to the model.
returns:
normal_vectors_array (np.ndarray): the normal vectors corresponding to each triangle in the surface.
"""
crs = rqc.Crs(self.model, uuid = self.crs_uuid)
triangles, points = self.triangles_and_points()
if crs.xy_units != crs.z_units:
points = points.copy()
wam.convert_lengths(points[:, 2], crs.z_units, crs.xy_units)
n_triangles = len(triangles)
normal_vectors_array = np.empty((n_triangles, 3))
for triangle_num in range(n_triangles):
normal_vector = vec.triangle_normal_vector_numba(points[triangles[triangle_num]])
if (normal_vector[2] > 0.0) == crs.z_inc_down:
normal_vector *= -1.0
normal_vectors_array[triangle_num] = normal_vector
if add_as_property:
pc = rqp.PropertyCollection()
pc.set_support(support = self)
crs = rqc.Crs(self.model, uuid = self.crs_uuid)
pc.add_cached_array_to_imported_list(
normal_vectors_array,
"computed from surface",
"normal vector",
uom = None, # uom not used for points properties
property_kind = "normal vector",
indexable_element = "faces",
points = True)
pc.write_hdf5_for_imported_list()
pc.create_xml_for_imported_list_and_add_parts_to_model(find_local_property_kinds = True)
return normal_vectors_array
[docs] def axial_edge_crossings(self, axis, value = 0.0):
"""Returns list-like array of points interpolated from triangle edges that cross value in axis.
arguments:
axis (int): xyz axis; 0 for crossings in x axis, 1 for y, 2 for z
value (float, default 0.0): the value in axis at which crossing points are sought
returns:
numpy float array of shape (N, 3) being interpolated triangle edge points with the crossing
value in axis; or None if no crossings found
"""
t, p = self.triangles_and_points()
zero = p[:, axis] - value # +ve one side of value, -ve the other side
abs_zero = np.expand_dims(np.abs(zero), axis = -1) # used to weight ends of crossing edges to interpolate
e, _ = triangulate.edges(t) # list-like pairs of p indices for ends of distinct edges
# find crossing edges
crossing = np.where(zero[e[:, 0]] * zero[e[:, 1]] < 0.0) # zero crossing edge indices in e
if not len(crossing[0]):
return None
ec = e[crossing]
s = abs_zero[ec[:, 0]] + abs_zero[ec[:, 1]]
# TODO: ignore divide by zero
ep = (p[ec[:, 0]] * abs_zero[ec[:, 1]] + p[ec[:, 1]] * abs_zero[ec[:, 0]]) / s
# TODO: unignore, and use midpoint of those edges?
return ep
[docs] def resampled_surface(self, title = None):
"""Creates a new triangulated set which is a resampled version of the current triangulated set. Each existing triangle in the tset is divided equally into 4 new triangles.
arguments:
title (str): a new title for the output triangulated set, if None the title will have the same title as the input triangulated set
returns:
resqpy.surface.Surface object, with extra_metadata ('resampled from surface': <uuid>), where uuid is the origin surface uuid
"""
rt, rp = self.triangles_and_points()
edge1 = np.mean(rp[rt[:]][:, ::2, :], axis = 1)
edge2 = np.mean(rp[rt[:]][:, 1:, :], axis = 1)
edge3 = np.mean(rp[rt[:]][:, :2, :], axis = 1)
allpoints = np.concatenate((rp, edge1, edge2, edge3), axis = 0)
count1 = len(rp)
count2 = count1 + len(edge1)
count3 = count2 + len(edge2)
tris = []
for i in range(len(rt)):
tris.extend([[rt[i][0], count1 + i, count3 + i], [rt[i][1], count2 + i, count3 + i],
[rt[i][2], count1 + i, count2 + i], [count1 + i, count2 + i, count3 + i]])
# TODO: implement alternate solution using edge functions in olio triangulation to optimise
points_unique, inverse = np.unique(allpoints, axis = 0, return_inverse = True)
tris = np.array(tris)
tris_unique = np.empty(shape = tris.shape, dtype = int)
tris_unique[:, 0] = inverse[tris[:, 0]]
tris_unique[:, 1] = inverse[tris[:, 1]]
tris_unique[:, 2] = inverse[tris[:, 2]]
if title is None:
title = self.citation_title
resampled = rqs.Surface(self.model,
title = title,
crs_uuid = self.crs_uuid,
extra_metadata = {'resampled from surface': str(self.uuid)})
resampled.set_from_triangles_and_points(tris_unique, points_unique)
return resampled
[docs] def write_hdf5(self, file_name = None, mode = 'a'):
"""Create or append to an hdf5 file, writing datasets for the triangulated patches after caching arrays.
:meta common:
"""
if self.uuid is None:
self.uuid = bu.new_uuid()
# NB: patch arrays must all have been set up prior to calling this function
h5_reg = rwh5.H5Register(self.model)
# todo: sort patches by patch index and check sequence
for triangulated_patch in self.patch_list:
(t, p) = triangulated_patch.triangles_and_points()
h5_reg.register_dataset(self.uuid, 'points_patch{}'.format(triangulated_patch.patch_index), p)
h5_reg.register_dataset(self.uuid, 'triangles_patch{}'.format(triangulated_patch.patch_index), t)
h5_reg.write(file_name, mode = mode)
[docs] def create_xml(self,
ext_uuid = None,
add_as_part = True,
add_relationships = True,
crs_uuid = None,
title = None,
originator = None):
"""Creates a triangulated surface xml node from this surface object and optionally adds as part of model.
arguments:
ext_uuid (uuid.UUID): the uuid of the hdf5 external part holding the surface arrays
add_as_part (boolean, default True): if True, the newly created xml node is added as a part
in the model
add_relationships (boolean, default True): if True, a relationship xml part is created relating the
new triangulated representation part to the crs part (and optional interpretation part)
crs_uuid (optional): the uuid of the coordinate reference system applicable to the surface points data;
if None, the main crs for the model is assumed to apply
title (string): used as the citation Title text; should be meaningful to a human
originator (string, optional): the name of the human being who created the triangulated representation part;
default is to use the login name
returns:
the newly created triangulated representation (surface) xml node
:meta common:
"""
if ext_uuid is None:
ext_uuid = self.model.h5_uuid()
if not self.title:
self.title = 'surface'
tri_rep = super().create_xml(add_as_part = False, title = title, originator = originator)
# todo: if crs_uuid is None, attempt to set to surface patch crs uuid (within patch loop, below)
if crs_uuid is not None:
self.crs_uuid = crs_uuid
if self.crs_uuid is None:
self.crs_uuid = self.model.crs_uuid # maverick use of model's default crs
if self.crs_uuid is None:
crs = rqc.Crs(self.model)
crs.create_xml()
self.crs_uuid = crs.uuid
assert self.crs_uuid is not None
if self.represented_interpretation_root is not None:
interp_root = self.represented_interpretation_root
interp_uuid = bu.uuid_from_string(interp_root.attrib['uuid'])
interp_part = self.model.part_for_uuid(interp_uuid)
interp_title = rqet.find_nested_tags_text(interp_root, ['Citation', 'Title'])
self.model.create_ref_node('RepresentedInterpretation',
interp_title,
interp_uuid,
content_type = self.model.type_of_part(interp_part),
root = tri_rep)
if interp_title and not title:
title = interp_title
# if not title: title = 'surface'
# self.model.create_citation(root = tri_rep, title = title, originator = originator)
role_node = rqet.SubElement(tri_rep, ns['resqml2'] + 'SurfaceRole')
role_node.set(ns['xsi'] + 'type', ns['resqml2'] + 'SurfaceRole')
role_node.text = self.surface_role
for patch in self.patch_list:
p_node = rqet.SubElement(tri_rep, ns['resqml2'] + 'TrianglePatch')
p_node.set(ns['xsi'] + 'type', ns['resqml2'] + 'TrianglePatch')
p_node.text = '\n'
pi_node = rqet.SubElement(p_node, ns['resqml2'] + 'PatchIndex')
pi_node.set(ns['xsi'] + 'type', ns['xsd'] + 'nonNegativeInteger')
pi_node.text = str(patch.patch_index)
ct_node = rqet.SubElement(p_node, ns['resqml2'] + 'Count')
ct_node.set(ns['xsi'] + 'type', ns['xsd'] + 'positiveInteger')
ct_node.text = str(patch.triangle_count)
cn_node = rqet.SubElement(p_node, ns['resqml2'] + 'NodeCount')
cn_node.set(ns['xsi'] + 'type', ns['xsd'] + 'nonNegativeInteger')
cn_node.text = str(patch.node_count)
triangles_node = rqet.SubElement(p_node, ns['resqml2'] + 'Triangles')
triangles_node.set(ns['xsi'] + 'type', ns['resqml2'] + 'IntegerHdf5Array')
triangles_node.text = '\n'
# not sure if null value node is needed, not actually used in data
triangles_null = rqet.SubElement(triangles_node, ns['resqml2'] + 'NullValue')
triangles_null.set(ns['xsi'] + 'type', ns['xsd'] + 'integer')
triangles_null.text = '-1' # or set to number of points in surface coords?
triangles_values = rqet.SubElement(triangles_node, ns['resqml2'] + 'Values')
triangles_values.set(ns['xsi'] + 'type', ns['eml'] + 'Hdf5Dataset')
triangles_values.text = '\n'
self.model.create_hdf5_dataset_ref(ext_uuid,
self.uuid,
'triangles_patch{}'.format(patch.patch_index),
root = triangles_values)
geom = rqet.SubElement(p_node, ns['resqml2'] + 'Geometry')
geom.set(ns['xsi'] + 'type', ns['resqml2'] + 'PointGeometry')
geom.text = '\n'
self.model.create_crs_reference(crs_uuid = self.crs_uuid, root = geom)
points_node = rqet.SubElement(geom, ns['resqml2'] + 'Points')
points_node.set(ns['xsi'] + 'type', ns['resqml2'] + 'Point3dHdf5Array')
points_node.text = '\n'
coords = rqet.SubElement(points_node, ns['resqml2'] + 'Coordinates')
coords.set(ns['xsi'] + 'type', ns['eml'] + 'Hdf5Dataset')
coords.text = '\n'
self.model.create_hdf5_dataset_ref(ext_uuid,
self.uuid,
'points_patch{}'.format(patch.patch_index),
root = coords)
patch.node = p_node
if add_as_part:
self.model.add_part('obj_TriangulatedSetRepresentation', self.uuid, tri_rep)
if add_relationships:
# todo: add multiple crs'es (one per patch)?
crs_root = self.model.root_for_uuid(self.crs_uuid)
self.model.create_reciprocal_relationship(tri_rep, 'destinationObject', crs_root, 'sourceObject')
if self.represented_interpretation_root is not None:
self.model.create_reciprocal_relationship(tri_rep, 'destinationObject',
self.represented_interpretation_root, 'sourceObject')
ext_part = rqet.part_name_for_object('obj_EpcExternalPartReference', ext_uuid, prefixed = False)
ext_node = self.model.root_for_part(ext_part)
self.model.create_reciprocal_relationship(tri_rep, 'mlToExternalPartProxy', ext_node,
'externalPartProxyToMl')
return tri_rep
def distill_triangle_points(t, p):
"""Returns a (triangles, points) pair with points distilled as only those used from p."""
assert np.all(t < len(p))
# find unique points used by triangles
p_keep = np.unique(t)
# note new point index for each old point that is being kept
p_map = np.full(len(p), -1, dtype = int)
p_map[p_keep] = np.arange(len(p_keep))
# copy those unique points into a trimmed points array
points_distilled = p[p_keep]
# copy triangles, replacing p indices with compressed indices
triangles_mapped = p_map[t]
assert triangles_mapped.shape == t.shape
assert np.all(triangles_mapped >= 0)
assert np.all(triangles_mapped < len(points_distilled))
return triangles_mapped, points_distilled
def nan_removed_triangles_and_points(t, p):
"""Returns a (triangles, points) pair which excludes any point involving a NaN and related triangles."""
assert p.ndim == 2 and p.shape[1] == 3
assert t.ndim == 2 and t.shape[1] == 3
p_nan_mask = np.any(np.isnan(p), axis = 1)
p_non_nan_mask = np.logical_not(p_nan_mask)
expanded_mask = np.empty(p.shape, dtype = bool)
expanded_mask[:] = np.expand_dims(p_non_nan_mask, axis = -1)
p_filtered = p[expanded_mask].reshape((-1, 3))
t_nan_mask = np.any(p_nan_mask[t], axis = 1)
expanded_mask = np.empty(t.shape, dtype = bool)
expanded_mask[:] = np.expand_dims(np.logical_not(t_nan_mask), axis = -1)
t_filtered = t[expanded_mask].reshape((-1, 3))
# modified the filtered t values to adjust for the compression of filtered p
p_map = np.full(len(p), -1, dtype = int)
p_map[p_non_nan_mask] = np.arange(len(p_filtered), dtype = int)
t_filtered = p_map[t_filtered]
assert t_filtered.ndim == 2 and t_filtered.shape[1] == 3
assert not np.any(t_filtered < 0) and not np.any(t_filtered >= len(p_filtered))
return (t_filtered, p_filtered)
def _adjust_flange_z(model, crs_uuid, p_xy_e, flange_start_index, t, flange_array, saucer_parameter):
"""Adjust the flange point z values (in recumbent space) by extrapolation of pair of points on original."""
# reconstruct the hull (could be concave) of original points
all_edges, edge_use_count = triangulate.edges(t)
inner_edges = triangulate.internal_edges(all_edges, edge_use_count)
t_for_inner_edges = triangulate.triangles_using_edges(t, inner_edges)
assert np.all(t_for_inner_edges >= 0)
flange_pairs = flange_array[t_for_inner_edges]
rim_edges = inner_edges[np.where(flange_pairs[:, 0] != flange_pairs[:, 1])]
assert rim_edges.ndim == 2 and rim_edges.shape[1] == 2 and len(rim_edges) > 0
rim_edge_index_list, rim_point_index_list = triangulate.rims(rim_edges)
assert len(rim_edge_index_list) == 1 and len(rim_point_index_list) == 1
rim_edge_indices = rim_edge_index_list[0]
rim_point_indices = rim_point_index_list[0] # ordered list of points on original hull (could be concave)
rim_pl = rql.Polyline(model,
set_coord = p_xy_e[rim_point_indices],
set_crs = crs_uuid,
is_closed = True,
title = 'rim')
centre = np.mean(p_xy_e[:flange_start_index], axis = 0)
# for each flange point, intersect a line from centre with the rim, and sample surface at saucer parameter
for flange_pi in range(flange_start_index, len(p_xy_e)):
f_xyz = p_xy_e[flange_pi]
pl_seg, rim_x, rim_y = rim_pl.first_line_intersection(centre[0], centre[1], f_xyz[0], f_xyz[1])
assert pl_seg is not None
rim_xyz = rim_pl.segment_xyz_from_xy(pl_seg, rim_x, rim_y)
sample_p = (1.0 - saucer_parameter) * centre + saucer_parameter * rim_xyz
p_list = vec.points_in_triangles_njit(np.expand_dims(sample_p, axis = 0), p_xy_e[t], 1)
vertical = np.array((0.0, 0.0, 1.0), dtype = float)
assert len(p_list) > 0
triangle_index, p_index, _ = p_list[0]
start_xyz = meet.line_triangle_intersect_numba(sample_p, vertical, p_xy_e[t[triangle_index]], t_tol = 0.05)
v_to_rim = rim_xyz - start_xyz
v_to_flange_p = f_xyz - start_xyz
if abs(v_to_rim[0]) > abs(v_to_rim[1]):
f = (v_to_rim[0]) / (v_to_flange_p[0])
else:
f = (v_to_rim[1]) / (v_to_flange_p[1])
assert 0.0 < f < 1.0
z = (rim_xyz[2] - start_xyz[2]) / f + start_xyz[2]
p_xy_e[flange_pi, 2] = z