Surfaces
A surface (also referred to as a triangulation), contains points and triangular facets (or faces), and edges that form the boundary of each facet.
Some examples in this page use libraries not installed by default.
Libraries referenced:
Surface examples
Creating a surface
This example will demonstrates basic Surface construction by making a tetrahedron.
from mapteksdk.data import Surface
from mapteksdk.project import Project
project = Project() # Connect to default project
path = "/surfaces/tetra"
with project.new(path, Surface, overwrite=True) as tetra:
size = 3
tetra.points = [[-size, -size, -size], # [x, y, z]
[-size, size, size],
[size, -size, size],
[size, size, -size]]
tetra.facets = [[0, 1, 2], [0, 2, 3], [0, 1, 3], [1, 2, 3]] # [i1, i2, i3]
tetra.facet_colours = [[255, 0, 0, 255], # Red [R, G, B, A]
[0, 255, 0, 255], # Blue
[0, 0, 255, 255], # Green
[255, 255, 255, 255]] # White
Creating a surface from a 00t file
This example will import data to use for the rest of the examples, starting with a Vulcan 00t triangulation. For more examples on converting between Surface objects and 00t files, see Built-In Input/Output.
Note: The data used in this and the following examples can be downloaded here: pit_model.00t
from mapteksdk.project import Project
from mapteksdk.data import io
project = Project() # Connect to default project
# The file path you have exported the sample data to:
file_path = "F:\\Python SDK Help\\data\\pit_model.00t"
# Import the 00t:
imported_surface = io.import_00t(file_path)
save_as = "/surfaces/pit_model"
if imported_surface:
# Successful imported, now assign a path
# After add_object, set imported_surface as the new ObjectID:
imported_surface = project.add_object(save_as, imported_surface, overwrite=True)
else:
print("{} failed to import.".format(file_path))
Getting information about a surface
This example returns basic information about the Surface object including primitive counts and extents.
from mapteksdk.project import Project
from mapteksdk.data import Surface
project = Project() # Connect to default project
# Where you stored the surface:
path = "/surfaces/pit_model"
found = project.find_object(path)
if found and found.is_a(Surface):
with project.read(path) as surface:
print("Name: {}\nPath: {}\nParent: {}".format(
found.name, found.path, found.parent.path))
print("Number of points: {}\nNumber of edges: {}\nNumber of facets: {}".format(
surface.point_count, surface.edge_count, surface.facet_count
))
extent = surface.extent
print("Extent - Minimum: X: {}, Y: {}, Z: {}".format(
extent.minimum[0], extent.minimum[1], extent.minimum[2]
))
print("Extent - Maximum: X: {}, Y: {}, Z: {}".format(
extent.maximum[0], extent.maximum[1], extent.maximum[2]
))
print("Centre: X: {}, Y: {}, Z: {}".format(
extent.centre[0], extent.centre[1], extent.centre[2]
))
print("Maximum length between points: {}".format(
extent.length
))
else:
print("Please import Surface to '{}' first.".format(path))
# Expected output:
# >>> Name: pit_model
# >>> Path: /surfaces/pit_model
# >>> Parent: /surfaces
# >>> Number of points: 64268
# >>> Number of edges: 192775
# >>> Number of facets: 128508
# >>> Extent - Minimum: X: 939.0660536488951, Y: 1654.8985445524079, Z: 284.4358990415397
# >>> Extent - Maximum: X: 1425.4774194945908, Y: 2208.9500778135116, Z: 384.14591444660067
# >>> Centre: X: 1182.271736571743, Y: 1931.9243111829596, Z: 334.2909067440702
# >>> Maximum length between points: 554.0515332611037
Converting a surface into a point set
Converting a Surface object into a PointSet object involves copying the points array, and optionally the point_colours array.
from mapteksdk.project import Project
from mapteksdk.data import Surface, PointSet
project = Project() # Connect to default project
# Where you stored the surface:
path = "/surfaces/pit_model"
output = path + " (points)"
found = project.find_object(path)
if found and found.is_a(Surface):
with project.read(path) as surface:
with project.new(output, PointSet, overwrite=True) as points:
points.points, points.point_colours = surface.points, surface.point_colours
else:
print("Please import Surface to '{}' first.".format(path))
Converting a surface into a an edge network
A surface can be converted into a wireframe representation using an EdgeNetwork. This involves copying the points and edges arrays and optionally the point_colours array.
from mapteksdk.project import Project
from mapteksdk.data import Surface, EdgeNetwork
project = Project() # Connect to default project
# Where you stored the surface:
path = "/surfaces/pit_model"
output = path + " (wireframe)"
found = project.find_object(path)
if found and found.is_a(Surface):
with project.read(path) as surface:
with project.new(output, EdgeNetwork, overwrite=True) as edges:
edges.points = surface.points
edges.point_colours = surface.point_colours
edges.edges = surface.edges
else:
print("Please import Surface to '{}' first.".format(path))
Rotating, translating, or scaling a surface
Rotating, translating, or scaling a Surface can be performed exactly the same way as rotating, translating, or scaling a PointSet. See the following topics for examples:
Determining the elevation of a coordinate on a surface
This example demonstrates how to find the height on a Surface that a point will intersect with it vertically.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import numpy as np
def find_intersection(triangles, point):
""" From an array of triangles, find which one the point of interest intersects """
p1, p2, p3 = triangles[:, 0, :], triangles[:, 1, :], triangles[:, 2, :]
v1 = p3 - p1
v2 = p2 - p1
cp = np.cross(v1, v2)
a, b, c = cp[:, 0], cp[:, 1], cp[:, 2]
# Create an array of the same point for the length of the triangles array
# This is to allow vectorisation of the intersection calculations
point_as_array = np.empty([triangles.shape[0], 3], dtype=np.float)
point_as_array[:] = point
# Calculate intersection height at each facet
d = np.einsum("ij,ij->i", cp, p3)
intersect_height = (d - a * point_as_array[:,0] - b * point_as_array[:, 1]) / c
# Determine which triangles the point sits inside of
Area = 0.5 *(-p2[:,1]*p3[:,0] + p1[:,1]*(-p2[:,0] + p3[:,0]) + p1[:,0]*(p2[:,1] - p3[:,1]) + p2[:,0]*p3[:,1])
s = 1/(2*Area)*(p1[:,1]*p3[:,0] - p1[:,0]*p3[:,1] + (p3[:,1] - p1[:,1])*point[0] + (p1[:,0] - p3[:,0])*point[1])
t = 1/(2*Area)*(p1[:,0]*p2[:,1] - p1[:,1]*p2[:,0] + (p1[:,1] - p2[:,1])*point[0] + (p2[:,0] - p1[:,0])*point[1])
possible_intersections = np.where(((s>=0) & (s<=1)) & ((t>=0) & (t<=1)) & (1-s-t > 0))
if len(possible_intersections[0]) > 0:
return intersect_height[possible_intersections[0]][0]
return None
if __name__ == "__main__":
project = Project() # Connect to default project
# Where you stored the surface:
path = "/surfaces/pit_model"
# The horizontal coordinate to determine intersection height on the surface
intersect_at = np.array([1165, 1985, 0])
found = project.find_object(path)
if found and found.is_a(Surface):
with project.read(path) as surface:
# Generate array of 'triangles' [[[x1,y1,z1],[x2,y2,z2],[x3,y3,z3]],]
triangles = surface.points[surface.facets]
# Query the intersection height (if any)
result = find_intersection(triangles, intersect_at)
if result:
print("Intersected X: {}, Y: {}, with height: {}".format(intersect_at[0], intersect_at[1], result))
else:
print("Couldn't find any facets to intersect with.")
else:
print("Please import Surface to '{}' first.".format(path))
# Example output:
# Intersected X: 1165, Y: 1985, with height: 289.4469219564403
Converting a surface to an STL file
This example demonstrates converting a Surface object to an STL file.
Important: This makes use of the trimesh library, which is not included by default with Python. You may need to install it before trying the script.
Note: The following example imports a Maptek object file containing a single surface to apply to the marker. It can be downloaded here: PythonSDK_SampleData_Truck.maptekobj.
More comprehensive examples for importing Maptek object files (.maptekobj) are provided in Built-In Input/Output.
from mapteksdk.project import Project
from mapteksdk.data import io
import trimesh
project = Project() # Connect to default project
# Use the path you have saved the example maptekobj file:
maptekobj_path = "F:\\Python SDK Help\\data\\PythonSDK_SampleData_Truck.maptekobj"
export_path = "F:\\Python SDK Help\\data\\Truck.stl"
# Import the data
imported_data = io.import_maptekobj(maptekobj_path)
path = "surfaces/truck"
# Assign a path to the surface
project.rename_object(imported_data, path, overwrite=True)
with project.read(path) as surface:
# Convert the Surface to trimesh
output = trimesh.Trimesh(vertices=surface.points,
faces=surface.facets)
# Convert to binary representation of stl
output_bytes = trimesh.exchange.stl.export_stl(output)
# Save to stl file
with open(export_path, 'wb') as f:
f.write(output_bytes)
Creating a surface from an STL file
Using the STL file created in the previous example, the following example will import the STL data to create a Surface.
Important: This makes use of the trimesh library, which is not included by default with Python. You may need to install it before trying the script.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import trimesh
project = Project() # Connect to default project
# Path to your stl file:
import_path = "F:\\Python SDK Help\\data\\Truck.stl"
mesh = trimesh.load_mesh(import_path)
with project.new("surfaces/imported stl", Surface, overwrite=True) as surface:
# Copy the points and facets from the trimesh object
surface.points, surface.facets = mesh.vertices, mesh.faces
Conditionally colouring a facet selection
This example demonstrates iterating over selected facet primitives, calculating dip angle of each facet and highlighting those with a dip greater than 70 degrees by colouring them red.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import numpy as np
import math
project = Project() # Connect to default project
def calculate_dip(p1, p2, p3):
""" Calculate dip angle of a facet,
where p1, p2, p3 are the 3d points it is created from """
AC = p1 - p3
AB = p1 - p2
ui, uj, uk = np.cross(AC, AB)
return abs(math.degrees(math.asin(math.sqrt((ui**2 + uj**2)) / math.sqrt(ui**2 + uj**2 + uk**2))))
# For all surfaces in selection that have selected facets:
for item in project.get_selected():
if item.is_a(Surface):
with project.edit(item) as surface:
# Get list of indices of selected facets
selection = np.where(surface.facet_selection)
# Iterate the indices as i
for i in np.nditer(selection):
# If the dip is >70 degres, colour the facet red
# surface.facets[i] -> [i1,i2,i3]
# surface.points[surface.facets[i][0]] -> point[i1] -> p1 [x, y, z]
if calculate_dip(surface.points[surface.facets[i][0]],
surface.points[surface.facets[i][1]],
surface.points[surface.facets[i][2]]) > 70:
surface.facet_colours[i] = [255, 0, 0, 255]
Colouring a surface with two-sided colouring
The Surface.back_colour property allows surfaces to be coloured differently on the front and back. The front is coloured using the point or facet colours (depending on which are defined) and the back is coloured uniformly using the back colour. The following example script demonstrates creating a surface with two-sided colouring:
from mapteksdk.project import Project
from mapteksdk.data import Surface
from mapteksdk.operations import open_new_view
if __name__ == "__main__":
width = 5.0
height = 2.5
with Project() as project:
with project.new("surfaces/two_sided_open_box", Surface) as surface:
surface.points = [
[-width, -width, 0], [-width, width, 0],
[width, width, 0], [width, -width, 0],
[-width, -width, height], [-width, width, height],
[width, width, height], [width, -width, height],
]
surface.facets = [
[1, 2, 0], [2, 3, 0], # Bottom.
[7, 4, 0], [3, 7, 0], # Front.
[4, 1, 0], [4, 5, 1], # Left.
[5, 2, 1], [5, 6, 2], # Back.
[6, 3, 2], [7, 3, 6] # Right.
]
surface.facet_colours = [
[255, 0, 0, 255], [255, 0, 0, 255],
[0, 255, 0, 255], [0, 255, 0, 255],
[0, 0, 255, 255], [0, 0, 255, 255],
[255, 255, 0, 255], [255, 255, 0, 255],
[255, 0, 255, 255], [255, 0, 255, 255],
]
# Set the back of each facet to be grey.
surface.back_colour = [127, 127, 127]
open_new_view(surface.id)
The above script creates the following surface:
-
The back colour does not have an alpha component.
-
The back colour is the same for all facets, regardless of how the surface is coloured on the front.
-
The front and back of a surface is determined by the normals of the facet.
-
If a surface has inconsistent normals, then what is the front and what is the back of a facet will be inconsistent across the surface. This will cause inconsistent colouring on the surface.
-
The fix surface tool (e.g. fix_surface()) can be used to resolve this issue.
Deleting a facet selection
This example demonstrates efficient removal of facets from a Surface that have been selected in the view.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import numpy as np
project = Project() # Connect to default project
# For all surfaces in selection that have selected facets:
for item in project.get_selected():
if item.is_a(Surface):
with project.edit(item) as surface:
# Get list of indices of selected facets
selection = np.where(surface.facet_selection)
# Remove selected facet indices
surface.remove_facets(selection)
# Clear selection
surface.facet_selection = [False]
Creating a new surface from a facet selection
This example demonstrates how to copy selected primitives from an object, however will create a new object for each object selected if there are multiple. Other techniques can be used to merge multiple object primitive selections into a single new object.
from mapteksdk.project import Project
from mapteksdk.data import Surface
project = Project() # Connect to default project
# For all surfaces in selection that have selected facets:
for item in project.get_selected():
if item.is_a(Surface):
# Prepare path to store new object:
new_path = item.path + "_selection"
with project.edit(item) as surface:
with project.new(new_path, Surface, overwrite=True) as surface_new:
# Orphaned points will be removed automatically
surface_new.points = surface.points
surface_new.point_colours = surface.point_colours
# Copy only the selected facets and any colour information
surface_new.facets = surface.facets[surface.facet_selection]
surface_new.facet_colours = surface.facet_colours[surface.facet_selection]
Create facet primitive attributes
This example will demonstrate applying a "dip" attribute to each facet on the surface, which can be referenced by the Custom Colour tool in the software or used in conjunction with colour maps.
The Maptek Python SDK supports directly applying colour maps to points, but not edges or facets.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import numpy as np
project = Project() # Connect to default project
# Where you stored the surface:
path = "/surfaces/pit_model"
found = project.find_object(path)
if found and found.is_a(Surface):
with project.edit(path) as surface:
# Below is a vectorised calculation for dip angle for each facet
# Generate array of 'triangles' [[[x1,y1,z1],[x2,y2,z2],[x3,y3,z3]],]
triangles = surface.points[surface.facets]
p1, p2, p3 = triangles[:, 0, :], triangles[:, 1, :], triangles[:, 2, :]
AC = p1 - p3
AB = p1 - p2
cp = np.cross(AC, AB)
a, b, c = cp[:, 0], cp[:, 1], cp[:, 2]
dip_results = np.abs(np.degrees(np.arcsin(np.sqrt((a**2 + b**2)) / np.sqrt((a**2 + b**2 + c**2)))))
# Save the dip values in an attribute "dip" for the software to refer to
surface.save_facet_attribute("dip", dip_results)
else:
print("Please import Surface to '{}' first.".format(path))
Subdividing a surface
This example shows how to subdivide a Surface such that each facet is replaced with four smaller facets.
Important: This makes use of the trimesh library, which is not included by default with Python. You may need to install it before trying the script.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import trimesh.remesh
project = Project() # Connect to default project
path = "/surfaces/pit_model"
if project.find_object(path):
with project.read(path) as surface:
with project.new(path + " divided", Surface, overwrite=True) as new_surface:
new_surface.points, new_surface.facets = trimesh.remesh.subdivide(surface.points, surface.facets)
Creating a Delaunay surface from a CSV file
The following is a basic example that imports a CSV (format: X, Y, Z, R, G, B), triangulates it and creates a Surface object.
Note: The data used in this and the following examples can be downloaded here: pointset_3d.csv
Important: This makes use of the scipy.spatial.Delaunay library, which may not be included by default with Python. You may need to install it before trying the script.
from mapteksdk.project import Project
from mapteksdk.data import Surface
from numpy import genfromtxt
from scipy.spatial import Delaunay
import numpy as np
# The file path you have exported the sample data to:
file_path = "F:\\Python SDK Help\\data\\pointset_3d.csv"
# Read the CSV using numpy
raw_data = genfromtxt(file_path, delimiter=',') # [x,y,z,r,g,b]
# Extract the points and colours into separate arrays
points = raw_data[:,[0, 1, 2]] # x,y,z
colours = raw_data[:,[3, 4, 5]] # r,g,b
# Attach a column of same length of 255 to populate alpha colour
alpha = np.empty([colours.shape[0], 1])
alpha.fill(255)
colours = np.column_stack((colours, alpha)) # r,g,b,a
# Compute triangulation - Only x, y needed for delaunay triangulation
tri = Delaunay(points[:, :2]) # x, y
project = Project() # Connect to default project
with project.new("/surfaces/imported csv", Surface, overwrite=True) as surface:
surface.points = points
surface.point_colours = colours
surface.facets = tri.simplices
Creating a surface from visible points in selected objects
This example builds on examples provided in the example Creating a point set, creating a Delaunay triangulation using the scipy.spatial.Delaunay module.
Important: This makes use of the scipy.spatial.Delaunay library, which may not be included by default with Python. You may need to install it before trying the script.
from mapteksdk.project import Project
from mapteksdk.data import Surface
import numpy as np
from scipy.spatial import Delaunay
project = Project() # Connect to default project
save_as = "scrapbook/create_surface_example"
# Temporarily store points and colours in these arrays
all_points = np.empty([0,3], dtype=np.float64)
all_colours = np.empty([0,4], dtype=np.uint8)
for item in project.get_selected():
with project.read(item) as obj:
if not hasattr(obj, 'points'): continue
# Append visible points from obj to buffer object
# Use .point_selection to get indices for all True values from the selection array
all_points = np.vstack((all_points, obj.points[obj.point_visibility]))
all_colours = np.vstack((all_colours, obj.point_colours[obj.point_visibility]))
if len(all_points) > 3:
with project.new(save_as, Surface, overwrite=True) as surface:
surface.points = all_points
surface.point_colours = all_colours
print("Calculate triangulation of {} points".format(len(all_points)))
# for scipy.spatial.Delaunay we want [[x,y],[x,y],etc]
x_y_array = all_points[:, :2] # x, y (:2) columns of take all points (:)
tri = Delaunay(x_y_array)
# Store facet indices:
surface.facets = tri.simplices
else:
print("Not enough points in selected objects to triangulate.")
Creating a hyperbolic paraboloid
This example generates a hyperbolic paraboloid and outputs a representation of it as a PointSet , EdgeNetwork (wireframe) and Surface, using four input boundary points. Smoothness of the result is relative to the number of segments each line is divided by. The example also illustrates a way of subdividing a line into multiple segments.
The example also illustrates how the SDK will handle cases where a list type may be used for convenience instead of a NumPy array for providing coordinates and edge arrays.
from mapteksdk.project import Project
from mapteksdk.data import Surface, EdgeNetwork, PointSet
from scipy.spatial import Delaunay
import numpy as np
project = Project() # Connect to default project
def subdivide_line(p1, p2, segments):
""" Divide a line between p1 and p2 into a number of segments """
x_delta = (p2[0] - p1[0]) / float(segments)
y_delta = (p2[1] - p1[1]) / float(segments)
z_delta = (p2[2] - p1[2]) / float(segments)
points = []
for i in range(1, segments):
points.append([p1[0] + i * x_delta, p1[1] + i * y_delta, p1[2] + i * z_delta])
return [p1] + points + [p2]
path_prefix = "/surfaces/hyperbolic paraboloid"
NW = [0, 35, 17] # NW point
SW = [0, 0, 2] # SW point
SE = [15, 0, 13] # SE point
NE = [15, 35, 5] # NE point
segments = 40 # Number of segments to subdivide each line
# Create set of points alone the outline of the shape
NW_NE = subdivide_line(NW, NE, segments)
SW_SE = subdivide_line(SW, SE, segments)
NW_SW = subdivide_line(NW, SW, segments)
NE_SE = subdivide_line(NE, SE, segments)
# Temporary storage for points and edge joins
points = []
edges = []
# Temporary storage for all intersection points
intersect_points = []
for i in range(len(NW_NE)):
# NW -> NE joins to SW -> SE
points.append(NW_NE[i])
points.append(SW_SE[i])
edges.append([len(points)-1, len(points)-2])
# NW -> SW joins to NE -> SE
points.append(NW_SW[i])
points.append(NE_SE[i])
edges.append([len(points)-1, len(points)-2])
# Generate points for all intersections along these lines
intersect_points += subdivide_line(SW_SE[i], NW_NE[i], segments)
intersect_points += subdivide_line(NE_SE[i], NW_SW[i], segments)
# Store all the bounary and intersection points
with project.new(path_prefix + "_points", PointSet, overwrite=True) as pts:
intersect_points += points
pts.points = intersect_points
# Store wireframe of the shape
with project.new(path_prefix + "_edges", EdgeNetwork, overwrite=True) as lines:
lines.points = points
lines.edges = edges
# Generate a surface of the shape using the points
with project.new(path_prefix + "_surface", Surface, overwrite=True) as sfc:
# Convert the list() of points to a numpy array of [X, Y] for the Delaunay function
tri = Delaunay(np.array(intersect_points)[:, :2])
sfc.points = intersect_points
sfc.facets = tri.simplices
