diff --git a/resqpy/surface/_surface.py b/resqpy/surface/_surface.py index b0d9236d..531d6c9c 100644 --- a/resqpy/surface/_surface.py +++ b/resqpy/surface/_surface.py @@ -8,6 +8,7 @@ log = logging.getLogger(__name__) import numpy as np +import math as maths import resqpy.crs as rqc import resqpy.lines as rql @@ -47,7 +48,7 @@ def __init__(self, originator = None, extra_metadata = {}): """Create an empty Surface object (RESQML TriangulatedSetRepresentation). - + Optionally populates from xml, point set or mesh. arguments: @@ -706,6 +707,187 @@ def set_from_point_set(self, self.set_from_triangles_and_points(t, p_e) return flange_array + def extend_surface_with_flange(self, + convexity_parameter = 5.0, + reorient = False, + reorient_max_dip = None, + flange_point_count = 11, + flange_radial_factor = 10.0, + flange_radial_distance = None, + flange_inner_ring = False, + saucer_parameter = None, + make_clockwise = False, + retriangulate = False): + """Returns a new Surface object where the original surface has been extended with a flange with a Delaunay triangulation of points in a PointSet object. + + arguments: + 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 + flange_point_count (int, default 11): the number of points to generate in the flange ring; ignored if + retriangulate 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. Ignored if retriangulate + is False + 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 + retriangulate (bool, default False): if True, the surface will be generated with a retriangulation of + the existing points. If False, the surface will be generated by adding flange points and triangles directly + from the original surface edges, and will no retriangulate the input surface. If False the surface must not + contain tears + + returns: + a new surface, and a boolean array of length N, where N is the number of triangles on the surface. This boolean + array is False on original triangle points, and True for extended flange triangles + + notes: + 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 + """ + prev_t, prev_p = self.triangles_and_points() + point_set = rqs.PointSet(self.model, crs_uuid = self.crs_uuid, title = self.title, points_array = prev_p) + if retriangulate: + out_surf = Surface(self.model, crs_uuid = self.crs_uuid, title = self.title) + return out_surf, out_surf.set_from_point_set(point_set, convexity_parameter, reorient, reorient_max_dip, + True, flange_point_count, flange_radial_factor, + flange_radial_distance, flange_inner_ring, saucer_parameter, + make_clockwise) + else: + 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) + assert prev_p.ndim >= 2 + assert prev_p.shape[-1] == 3 + p = prev_p.reshape((-1, 3)) + 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, _, reorient_matrix = triangulate.reorient(unit_adjusted_p, max_dip = reorient_max_dip) + else: + p_xy = unit_adjusted_p + reorient_matrix = None + + centre_point = np.nanmean(p_xy.reshape((-1, 3)), axis = 0) # work out the radius for the flange points + p_radius_v = np.nanmax(np.abs(p.reshape((-1, 3)) - np.expand_dims(centre_point, axis = 0)), axis = 0)[:2] + p_radius = maths.sqrt(np.sum(p_radius_v * p_radius_v)) + radius = p_radius * flange_radial_factor + if flange_radial_distance is not None and flange_radial_distance > radius: + radius = flange_radial_distance + + de, dc = self.distinct_edges_and_counts() # find the distinct edges and counts + unique_edge = de[dc == 1] # find hull edges (edges on only a single triangle) + hull_points = p_xy[unique_edge] # find points defining the hull edges + hull_centres = np.mean(hull_points, axis = 1) # find the centre of each edge + + flange_points = np.empty( + shape = (hull_centres.shape), dtype = float + ) # loop over all the hull centres, generating a flange point and finding the azimuth from the centre to the hull centre point + az = np.empty(shape = len(hull_centres), dtype = float) + for i, c in enumerate(hull_centres): + v = [centre_point[0] - c[0], centre_point[1] - c[1], centre_point[2] - c[2]] + uv = -vec.unit_vector(v) + az[i] = vec.azimuth(uv) + flange_point = centre_point + radius * uv + if simple_saucer_angle is not None: + z_shift = radius * maths.tan(vec.radians_from_degrees(simple_saucer_angle)) + flange_point[2] -= z_shift + flange_points[i] = flange_point + + sort_az_ind = np.argsort(np.array(az)) # sort by azimuth, to run through the hull points + new_points = np.empty(shape = (len(flange_points), 3), dtype = float) + new_triangles = np.empty(shape = (len(flange_points) * 2, 3), dtype = int) + point_offset = len(p_xy) # the indices of the new triangles will begin after this + for i, ind in enumerate(sort_az_ind): # loop over each point in azimuth order + new_points[i] = flange_points[ind] + this_hull_edge = unique_edge[ind] + + def az_for_point(c): + v = [centre_point[0] - c[0], centre_point[1] - c[1], centre_point[2] - c[2]] + uv = -vec.unit_vector(v) + return vec.azimuth(uv) + + this_edge_az_sort = np.array( + [az_for_point(p_xy[this_hull_edge[0]]), + az_for_point(p_xy[this_hull_edge[1]])]) + if np.min(this_edge_az_sort) < az[ind] < np.max(this_edge_az_sort): + first, second = np.argsort(this_edge_az_sort) + else: + second, first = np.argsort(this_edge_az_sort) + if i != len(sort_az_ind) - 1: + new_triangles[2 * i] = np.array( + [this_hull_edge[first], this_hull_edge[second], + i + point_offset]) # add a triangle between the two hull points and the flange point + new_triangles[(2 * i) + 1] = np.array( + [this_hull_edge[second], i + point_offset, + i + point_offset + 1]) # for all but the last point, hookup triangle to the next flange point + else: + new_triangles[2 * i] = np.array( + [this_hull_edge[first], this_hull_edge[second], + i + point_offset]) # add a triangle between the two hull points and the first flange point + new_triangles[(2 * i) + 1] = np.array( + [this_hull_edge[second], point_offset, + i + point_offset]) # add in the final triangle between the first and last flange points + + all_points = np.concatenate((p_xy, new_points)) # concatenate triangle and points arrays + all_triangles = np.concatenate((prev_t, new_triangles)) + + flange_array = np.zeros(shape = all_triangles.shape[0], dtype = bool) + flange_array[ + len(prev_t):] = True # make a flange bool array, where all new triangles are flange and therefore True + + assert len(all_points) == ( + point_offset + len(flange_points)), "New point count should be old point count + flange point count" + assert len(all_triangles) == ( + len(prev_t) + + 2 * len(flange_points)), "New triangle count should be old triangle count + 2 x #flange points" + + if saucer_parameter is not None: + _adjust_flange_z(self.model, crs.uuid, all_points, len(all_points), all_triangles, flange_array, + saucer_parameter) # adjust the flange points if in saucer mode + if reorient: + all_points = vec.rotate_array(reorient_matrix.T, all_points) + if crs.xy_units != crs.z_units and reorient: + wam.convert_lengths(all_points[:, 2], crs.xy_units, crs.z_units) + + if make_clockwise: + triangulate.make_all_clockwise_xy(all_triangles, all_points) # modifies t in situ + + out_surf = Surface(self.model, crs_uuid = self.crs_uuid, title = self.title) + out_surf.set_from_triangles_and_points(all_triangles, all_points) # create the new surface + return out_surf, flange_array + def make_all_clockwise_xy(self, reorient = False): """Reorders cached triangles data such that all triangles are clockwise when viewed from -ve z axis. @@ -845,7 +1027,7 @@ def set_to_single_cell_faces_from_corner_points(self, cp, quad_triangles = True) 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 @@ -881,7 +1063,7 @@ def cell_axis_and_polarity_from_triangle_index(self, triangle_index): 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: @@ -985,7 +1167,7 @@ def set_from_rms_file(self, filename, quad_triangles = False): 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: @@ -1170,11 +1352,11 @@ def resampled_surface(self, title = None): return resampled def resample_surface_unique_edges(self): - """Returns a new surface, with the same model, title and crs as the original surface, but with additional refined points along original surface tears and edges. - + """Returns a new surface, with the same model, title and crs as the original surface, but with additional refined points along original surface tears and edges. + Each edge forming a tear or outer edge in the surface will have 3 additional points added, with 2 additional points on each edge of the original triangle. The output surface is re-triangulated using these new points (tears will be filled) - returns: + returns: resqpy.surface.Surface object with extra_metadata ('unique edges resampled from surface': uuid), where uuid is for the original surface uuid """ _, op = self.triangles_and_points() diff --git a/tests/unit_tests/surface/test_surface.py b/tests/unit_tests/surface/test_surface.py index caf52c43..8faba90d 100644 --- a/tests/unit_tests/surface/test_surface.py +++ b/tests/unit_tests/surface/test_surface.py @@ -1209,3 +1209,139 @@ def test_edge_lengths(example_model_and_crs): assert_array_almost_equal(lengths, np.array([(3.0, 4.0, 5.0), (5.0, 3.0, 4.0)], dtype = float)) lengths = surf.edge_lengths(required_uom = 'cm') assert_array_almost_equal(lengths, np.array([(300.0, 400.0, 500.0), (500.0, 300.0, 400.0)], dtype = float)) + + +def test_extended_surface_with_flange_extension(example_model_and_crs): + model, crs = example_model_and_crs + + # create a set of points (must have a convex hull) + x = np.array([-1, -1.5, -2, -1.5, -1, 0, 1, 1.5, 2, 1.5, 1, 0], dtype = float) + 5 + y = np.array([-1.1, -0.5, 0, 0.5, 1, 0.5, 1.1, 0.5, 0, -0.5, -1, 0], dtype = float) + 5 + z = np.zeros(shape = (12), dtype = float) + 5 + p = np.stack((x, y, z), axis = -1) + + # make a PointSet object + ps = rqs.PointSet(model, crs_uuid = crs.uuid, points_array = p, title = 'points') + ps.write_hdf5() + ps.create_xml() + + # try out flange extension + surf = rqs.Surface(model, crs_uuid = ps.crs_uuid, title = 'surface from ' + str(ps.title)) + assert surf is not None + surf.set_from_point_set(ps, reorient = False, reorient_max_dip = None, extend_with_flange = False) + surf.write_hdf5() + surf.create_xml() + orig_t, orig_p = surf.triangles_and_points() + + new_surf, flange_bool = surf.extend_surface_with_flange(reorient = False, + flange_radial_factor = 10.0, + flange_radial_distance = None, + saucer_parameter = None, + retriangulate = False) + + new_t, new_p = new_surf.triangles_and_points() + + assert np.all(flange_bool[len(orig_t):]) + assert not np.any(flange_bool[:len(orig_t)]) + assert np.all(np.isin(orig_t, new_t)) + + np.testing.assert_array_almost_equal( + np.array([[5, 28.07078471], [26.81071628, 12.43307607], [27.71578211, 1.25570298], [24.45543821, -7.28011087], + [5, -17.98745138], [-16.42279299, -3.40843501], [-17.86764418, 7.76400526], + [-15.03584362, 16.39531138]]), new_p[len(orig_p):, :2]) + + +def test_extended_surface_with_flange_extension_saucer(example_model_and_crs): + model, crs = example_model_and_crs + + # create a set of points (must have a convex hull) + x = np.array([-1, -1.5, -2, -1.5, -1, 0, 1, 1.5, 2, 1.5, 1, 0], dtype = float) + 5 + y = np.array([-1.1, -0.5, 0, 0.5, 1, 0.5, 1.1, 0.5, 0, -0.5, -1, 0], dtype = float) + 5 + z = np.zeros(shape = (12), dtype = float) + 5 + p = np.stack((x, y, z), axis = -1) + + # make a PointSet object + ps = rqs.PointSet(model, crs_uuid = crs.uuid, points_array = p, title = 'points') + ps.write_hdf5() + ps.create_xml() + + # try out flange extension + surf = rqs.Surface(model, crs_uuid = ps.crs_uuid, title = 'surface from ' + str(ps.title)) + assert surf is not None + surf.set_from_point_set(ps, reorient = False, reorient_max_dip = None, extend_with_flange = False) + surf.write_hdf5() + surf.create_xml() + orig_t, orig_p = surf.triangles_and_points() + + new_surf, flange_bool = surf.extend_surface_with_flange(reorient = False, + flange_radial_factor = 10.0, + flange_radial_distance = None, + saucer_parameter = 45.0, + retriangulate = False) + + new_t, new_p = new_surf.triangles_and_points() + + assert np.all(flange_bool[len(orig_t):]) + assert not np.any(flange_bool[:len(orig_t)]) + assert np.all(np.isin(orig_t, new_t)) + + np.testing.assert_array_almost_equal( + np.array([[5, 28.07078471, -18.029118], [26.81071628, 12.43307607, -18.029118], + [27.71578211, 1.25570298, -18.029118], [24.45543821, -7.28011087, -18.029118], + [5, -17.98745138, -18.029118], [-16.42279299, -3.40843501, -18.029118], + [-17.86764418, 7.76400526, -18.029118], [-15.03584362, 16.39531138, -18.029118]]), + new_p[len(orig_p):]) + + +def test_extended_surface_with_flange_extension_retriangulate(example_model_and_crs): + model, crs = example_model_and_crs + + # number of random points to use + n = 20 + + # create a set of random points + x = np.random.random(n) * 1000.0 - 500.0 + y = np.random.random(n) * 1000.0 - 500.0 + z = np.random.random(n) #  note: triangulation does not use z values + p = np.stack((x, y, z), axis = -1) + centre = np.mean(p, axis = 0) + + # make a PointSet object + ps = rqs.PointSet(model, crs_uuid = crs.uuid, points_array = p, title = 'random points in square') + ps.write_hdf5() + ps.create_xml() + + # try out flange extension + orig_surf = rqs.Surface(model, crs_uuid = ps.crs_uuid, title = 'surface from ' + str(ps.title)) + assert orig_surf is not None + orig_surf.set_from_point_set(ps, + reorient = True, + reorient_max_dip = None, + extend_with_flange = False, + make_clockwise = False) + orig_surf.write_hdf5() + orig_surf.create_xml() + + # flange extension with simple saucer + new_surf, flange_bool = orig_surf.extend_surface_with_flange(reorient = True, + reorient_max_dip = None, + flange_point_count = 12, + flange_radial_factor = 2.0, + flange_radial_distance = 2000.0, + flange_inner_ring = True, + saucer_parameter = -60.0, + make_clockwise = False, + retriangulate = True) + new_surf.write_hdf5() + new_surf.create_xml() + + assert new_surf.node_count() == 56 + _, p = new_surf.triangles_and_points() + min_p = np.min(p, axis = 0) + max_p = np.max(p, axis = 0) + assert -2010.0 <= min_p[0] - centre[0] < -1900.0 + assert -2010.0 <= min_p[1] - centre[1] < -1900.0 + assert 1900.0 < max_p[0] - centre[0] <= 2010.0 + assert 1900.0 < max_p[1] - centre[1] <= 2010.0 + assert 0.0 <= min_p[2] <= 1.0 + assert 3464.0 < max_p[2] < 3465.5