diff --git a/mesh_inset/__init__.py b/mesh_inset/__init__.py
index 108ff3657ba071ddbc0a1583f0f072350128a19f..0a441d1b7e7665746389c9c5abb4ef638352b77b 100644
--- a/mesh_inset/__init__.py
+++ b/mesh_inset/__init__.py
@@ -16,26 +16,30 @@
#
# ##### END GPL LICENSE BLOCK #####
+# <pep8 compliant>
+
bl_info = {
- "name": "Inset Polygon",
- "author": "Howard Trickey",
- "version": (0, 2),
- "blender": (2, 5, 7),
- "api": 36147,
- "location": "View3D > Tools",
- "description": "Make an inset polygon inside selection.",
- "warning": "",
- "wiki_url": "http://wiki.blender.org/index.php/Extensions:2.5/Py/Scripts/Modeling/Inset-Polygon",
- "tracker_url": "http://projects.blender.org/tracker/index.php?func=detail&aid=27290&group_id=153&atid=468",
- "category": "Mesh"}
+ "name": "Inset Polygon",
+ "author": "Howard Trickey",
+ "version": (0, 2),
+ "blender": (2, 5, 7),
+ "api": 36147,
+ "location": "View3D > Tools",
+ "description": "Make an inset polygon inside selection.",
+ "warning": "",
+ "wiki_url": \
+ "http://wiki.blender.org/index.php/Extensions:2.5/Py/Scripts/Modeling/Inset-Polygon",
+ "tracker_url": \
+ "http://projects.blender.org/tracker/index.php?func=detail&aid=27290&group_id=153&atid=468",
+ "category": "Mesh"}
if "bpy" in locals():
- import imp
+ import imp
else:
- from . import geom
- from . import model
- from . import offset
- from . import triquad
+ from . import geom
+ from . import model
+ from . import offset
+ from . import triquad
import math
import bpy
@@ -44,134 +48,156 @@ from bpy.props import *
class Inset(bpy.types.Operator):
- bl_idname = "mesh.inset"
- bl_label = "Inset"
- bl_description = "Make an inset polygon inside selection"
- bl_options = {'REGISTER', 'UNDO'}
-
- inset_amount = FloatProperty(name="Amount",
- description="Distance of inset edges from outer ones",
- default = 0.05,
- min = 0.0,
- max = 1000.0,
- soft_min = 0.0,
- soft_max = 1.0,
- unit = 'LENGTH')
- inset_height = FloatProperty(name="Height",
- description="Distance to raise inset faces",
- default = 0.0,
- min = -1000.0,
- max = 1000.0,
- soft_min = -1.0,
- soft_max = 1.0,
- unit = 'LENGTH')
- region = BoolProperty(name="Region",
- description="Inset selection as one region?",
- default = True)
-
- @classmethod
- def poll(cls, context):
- obj = context.active_object
- return (obj and obj.type == 'MESH' and context.mode == 'EDIT_MESH')
-
- def draw(self, context):
- layout= self.layout
- box = layout.box()
- box.label("Inset Options")
- box.prop(self, "inset_amount")
- box.prop(self, "inset_height")
- box.prop(self, "region")
-
- def invoke(self, context, event):
- self.action(context)
- return {'FINISHED'}
-
- def execute(self, context):
- self.action(context)
- return {'FINISHED'}
-
- def action(self, context):
- save_global_undo = bpy.context.user_preferences.edit.use_global_undo
- bpy.context.user_preferences.edit.use_global_undo = False
+ bl_idname = "mesh.inset"
+ bl_label = "Inset"
+ bl_description = "Make an inset polygon inside selection"
+ bl_options = {'REGISTER', 'UNDO'}
+
+ inset_amount = FloatProperty(name="Amount",
+ description="Amount to move inset edges",
+ default=5.0,
+ min=0.0,
+ max=1000.0,
+ soft_min=0.0,
+ soft_max=100.0,
+ unit='LENGTH')
+ inset_height = FloatProperty(name="Height",
+ description="Amount to raise inset faces",
+ default=0.0,
+ min=-10000.0,
+ max=10000.0,
+ soft_min=-500.0,
+ soft_max=500.0,
+ unit='LENGTH')
+ region = BoolProperty(name="Region",
+ description="Inset selection as one region?",
+ default=True)
+ scale = EnumProperty(name="Scale",
+ description="Scale for amount",
+ items=[
+ ('PERCENT', "Percent",
+ "Percentage of maximum inset amount"),
+ ('ABSOLUTE', "Absolute",
+ "Length in blender units")
+ ],
+ default='PERCENT')
+
+ @classmethod
+ def poll(cls, context):
+ obj = context.active_object
+ return (obj and obj.type == 'MESH' and context.mode == 'EDIT_MESH')
+
+ def draw(self, context):
+ layout = self.layout
+ box = layout.box()
+ box.label("Inset Options")
+ box.prop(self, "scale")
+ box.prop(self, "inset_amount")
+ box.prop(self, "inset_height")
+ box.prop(self, "region")
+
+ def invoke(self, context, event):
+ self.action(context)
+ return {'FINISHED'}
+
+ def execute(self, context):
+ self.action(context)
+ return {'FINISHED'}
+
+ def action(self, context):
+ save_global_undo = bpy.context.user_preferences.edit.use_global_undo
+ bpy.context.user_preferences.edit.use_global_undo = False
+ bpy.ops.object.mode_set(mode='OBJECT')
+ obj = bpy.context.active_object
+ mesh = obj.data
+ do_inset(mesh, self.inset_amount, self.inset_height, self.region,
+ self.scale == 'PERCENT')
+ bpy.ops.object.mode_set(mode='EDIT')
+ bpy.context.user_preferences.edit.use_global_undo = save_global_undo
+
+
+def do_inset(mesh, amount, height, region, as_percent):
+ if amount <= 0.0:
+ return
+ pitch = math.atan(height / amount)
+ selfaces = []
+ selface_indices = []
+ for face in mesh.faces:
+ if face.select and not face.hide:
+ selfaces.append(face)
+ selface_indices.append(face.index)
+ m = geom.Model()
+ # if add all mesh.vertices, coord indices will line up
+ # Note: not using Points.AddPoint which does dup elim
+ # because then would have to map vertices in and out
+ m.points.pos = [v.co.to_tuple() for v in mesh.vertices]
+ for f in selfaces:
+ m.faces.append(list(f.vertices))
+ m.face_data.append(f.index)
+ orig_numv = len(m.points.pos)
+ orig_numf = len(m.faces)
+ model.BevelSelectionInModel(m, amount, pitch, True, region, as_percent)
+ if len(m.faces) == orig_numf:
+ # something went wrong with Bevel - just treat as no-op
+ return
+ # blender_faces: newfaces but all 4-tuples and no 0
+ # in 4th position if a 4-sided poly
+ blender_faces = []
+ blender_old_face_index = []
+ for i in range(orig_numf, len(m.faces)):
+ f = m.faces[i]
+ if len(f) == 3:
+ blender_faces.append(list(f) + [0])
+ blender_old_face_index.append(m.face_data[i])
+ elif len(f) == 4:
+ if f[3] == 0:
+ blender_faces.append([f[3], f[0], f[1], f[2]])
+ else:
+ blender_faces.append(f)
+ blender_old_face_index.append(m.face_data[i])
+ num_new_vertices = len(m.points.pos) - orig_numv
+ mesh.vertices.add(num_new_vertices)
+ for i in range(orig_numv, len(m.points.pos)):
+ mesh.vertices[i].co = mathutils.Vector(m.points.pos[i])
+ start_faces = len(mesh.faces)
+ mesh.faces.add(len(blender_faces))
+ for i, newf in enumerate(blender_faces):
+ mesh.faces[start_faces + i].vertices_raw = newf
+ # copy face attributes from old face that it was derived from
+ bfi = blender_old_face_index[i]
+ if bfi and 0 <= bfi < start_faces:
+ bfacenew = mesh.faces[start_faces + i]
+ bface = mesh.faces[bfi]
+ bfacenew.material_index = bface.material_index
+ bfacenew.use_smooth = bface.use_smooth
+ mesh.update(calc_edges=True)
+ # remove original faces
+ bpy.ops.object.mode_set(mode='EDIT')
+ save_select_mode = bpy.context.tool_settings.mesh_select_mode
+ bpy.context.tool_settings.mesh_select_mode = [False, False, True]
+ bpy.ops.mesh.select_all(action='DESELECT')
bpy.ops.object.mode_set(mode='OBJECT')
- obj = bpy.context.active_object
- mesh = obj.data
- do_inset(mesh, self.inset_amount, self.inset_height, self.region)
+ for fi in selface_indices:
+ mesh.faces[fi].select = True
bpy.ops.object.mode_set(mode='EDIT')
- bpy.context.user_preferences.edit.use_global_undo = save_global_undo
-
-
-def do_inset(mesh, amount, height, region):
- if amount <= 0.0:
- return
- pitch = math.atan(height / amount)
- selfaces = []
- selface_indices = []
- for face in mesh.faces:
- if face.select and not face.hide:
- selfaces.append(face)
- selface_indices.append(face.index)
- m = geom.Model()
- # if add all mesh.vertices, coord indices will line up
- # Note: not using Points.AddPoint which does dup elim
- # because then would have to map vertices in and out
- m.points.pos = [ v.co.to_tuple() for v in mesh.vertices ]
- for f in selfaces:
- m.faces.append(list(f.vertices))
- orig_numv = len(m.points.pos)
- orig_numf = len(m.faces)
- model.BevelSelectionInModel(m, m.faces, amount, pitch, True, region)
- if len(m.faces) == orig_numf:
- # something went wrong with Bevel - just treat as no-op
- return
- # blender_faces: newfaces but all 4-tuples and no 0
- # in 4th position if a 4-sided poly
- blender_faces = []
- for i in range(orig_numf, len(m.faces)):
- f = m.faces[i]
- if len(f) == 3:
- blender_faces.append(list(f) + [0])
- elif len(f) == 4:
- if f[3] == 0:
- blender_faces.append([f[3], f[0], f[1], f[2]])
- else:
- blender_faces.append(f)
- num_new_vertices = len(m.points.pos) - orig_numv
- mesh.vertices.add(num_new_vertices)
- for i in range(orig_numv, len(m.points.pos)):
- mesh.vertices[i].co = mathutils.Vector(m.points.pos[i])
- start_faces = len(mesh.faces)
- mesh.faces.add(len(blender_faces))
- for i, newf in enumerate(blender_faces):
- mesh.faces[start_faces + i].vertices_raw = newf
- mesh.update(calc_edges = True)
- # remove original faces
- bpy.ops.object.mode_set(mode='EDIT')
- save_select_mode = bpy.context.tool_settings.mesh_select_mode
- bpy.context.tool_settings.mesh_select_mode = [False, False, True]
- bpy.ops.mesh.select_all(action = 'DESELECT')
- bpy.ops.object.mode_set(mode = 'OBJECT')
- for fi in selface_indices:
- mesh.faces[fi].select = True
- bpy.ops.object.mode_set(mode = 'EDIT')
- bpy.ops.mesh.delete(type = 'FACE')
- bpy.context.tool_settings.mesh_select_mode = save_select_mode
+ bpy.ops.mesh.delete(type='FACE')
+ bpy.context.tool_settings.mesh_select_mode = save_select_mode
def panel_func(self, context):
- self.layout.label(text="Inset:")
- self.layout.operator("mesh.inset", text="Inset")
+ self.layout.label(text="Inset:")
+ self.layout.operator("mesh.inset", text="Inset")
def register():
- bpy.utils.register_class(Inset)
- bpy.types.VIEW3D_PT_tools_meshedit.append(panel_func)
+ bpy.utils.register_class(Inset)
+ bpy.types.VIEW3D_PT_tools_meshedit.append(panel_func)
def unregister():
- bpy.utils.unregister_class(Inset)
- bpy.types.VIEW3D_PT_tools_meshedit.remove(panel_func)
+ bpy.utils.unregister_class(Inset)
+ bpy.types.VIEW3D_PT_tools_meshedit.remove(panel_func)
if __name__ == "__main__":
- register()
+ register()
diff --git a/mesh_inset/geom.py b/mesh_inset/geom.py
index 7c21015efe80ff88a73f284617c520ae7eef1427..98059dcf29d40ef11d77b24157106d61aaf0db43 100644
--- a/mesh_inset/geom.py
+++ b/mesh_inset/geom.py
@@ -16,6 +16,8 @@
#
# ##### END GPL LICENSE BLOCK #####
+# <pep8 compliant>
+
"""Geometry classes and operations.
Also, vector file representation (Art).
"""
@@ -31,666 +33,687 @@ INVDISTTOL = 1e3
class Points(object):
- """Container of points without duplication, each mapped to an int.
+ """Container of points without duplication, each mapped to an int.
- Points are either have dimension at least 2, maybe more.
+ Points are either have dimension at least 2, maybe more.
- Implementation:
- In order to efficiently find duplicates, we quantize the points
- to triples of ints and map from quantized triples to vertex
- index.
+ Implementation:
+ In order to efficiently find duplicates, we quantize the points
+ to triples of ints and map from quantized triples to vertex
+ index.
- Attributes:
- pos: list of tuple of float - coordinates indexed by
- vertex number
- invmap: dict of (int, int, int) to int - quantized coordinates
- to vertex number map
- """
+ Attributes:
+ pos: list of tuple of float - coordinates indexed by
+ vertex number
+ invmap: dict of (int, int, int) to int - quantized coordinates
+ to vertex number map
+ """
- def __init__(self, initlist = []):
- self.pos = []
- self.invmap = dict()
- for p in initlist:
- self.AddPoint(p)
+ def __init__(self, initlist=[]):
+ self.pos = []
+ self.invmap = dict()
+ for p in initlist:
+ self.AddPoint(p)
- @staticmethod
- def Quantize(p):
- """Quantize the float tuple into an int tuple.
+ @staticmethod
+ def Quantize(p):
+ """Quantize the float tuple into an int tuple.
- Args:
- p: tuple of float
- Returns:
- tuple of int - scaled by INVDISTTOL and rounded p
- """
+ Args:
+ p: tuple of float
+ Returns:
+ tuple of int - scaled by INVDISTTOL and rounded p
+ """
- return tuple([int(round(v*INVDISTTOL)) for v in p])
+ return tuple([int(round(v * INVDISTTOL)) for v in p])
- def AddPoint(self, p):
- """Add point p to the Points set and return vertex number.
+ def AddPoint(self, p):
+ """Add point p to the Points set and return vertex number.
- If there is an existing point which quantizes the same,,
- don't add a new one but instead return existing index.
+ If there is an existing point which quantizes the same,,
+ don't add a new one but instead return existing index.
- Args:
- p: tuple of float - coordinates (2-tuple or 3-tuple)
- Returns:
- int - the vertex number of added (or existing) point
- """
+ Args:
+ p: tuple of float - coordinates (2-tuple or 3-tuple)
+ Returns:
+ int - the vertex number of added (or existing) point
+ """
- qp = Points.Quantize(p)
- if qp in self.invmap:
- return self.invmap[qp]
- else:
- self.invmap[qp] = len(self.pos)
- self.pos.append(p)
- return len(self.pos)-1
+ qp = Points.Quantize(p)
+ if qp in self.invmap:
+ return self.invmap[qp]
+ else:
+ self.invmap[qp] = len(self.pos)
+ self.pos.append(p)
+ return len(self.pos) - 1
- def AddPoints(self, points):
- """Add another set of points to this set.
+ def AddPoints(self, points):
+ """Add another set of points to this set.
- We need to return a mapping from indices
- in the argument points space into indices
- in this point space.
+ We need to return a mapping from indices
+ in the argument points space into indices
+ in this point space.
- Args:
- points: Points - to union into this set
- Returns:
- list of int: maps added indices to new ones
- """
+ Args:
+ points: Points - to union into this set
+ Returns:
+ list of int: maps added indices to new ones
+ """
- vmap = [ 0 ] * len(points.pos)
- for i in range(len(points.pos)):
- vmap[i] = self.AddPoint(points.pos[i])
- return vmap
+ vmap = [0] * len(points.pos)
+ for i in range(len(points.pos)):
+ vmap[i] = self.AddPoint(points.pos[i])
+ return vmap
- def AddZCoord(self, z):
- """Change this in place to have a z coordinate, with value z.
+ def AddZCoord(self, z):
+ """Change this in place to have a z coordinate, with value z.
- Assumes the coordinates are currently 2d.
+ Assumes the coordinates are currently 2d.
- Args:
- z: the value of the z coordinate to add
- Side Effect:
- self now has a z-coordinate added
- """
+ Args:
+ z: the value of the z coordinate to add
+ Side Effect:
+ self now has a z-coordinate added
+ """
- assert(len(self.pos) == 0 or len(self.pos[0]) == 2)
- newinvmap = dict()
- for i, (x,y) in enumerate(self.pos):
- newp = (x, y, z)
- self.pos[i] = newp
- newinvmap[self.Quantize(newp)] = i
- self.invmap = newinvmap
+ assert(len(self.pos) == 0 or len(self.pos[0]) == 2)
+ newinvmap = dict()
+ for i, (x, y) in enumerate(self.pos):
+ newp = (x, y, z)
+ self.pos[i] = newp
+ newinvmap[self.Quantize(newp)] = i
+ self.invmap = newinvmap
- def AddToZCoord(self, i, delta):
- """Change the z-coordinate of point with index i to add delta.
+ def AddToZCoord(self, i, delta):
+ """Change the z-coordinate of point with index i to add delta.
- Assumes the coordinates are currently 3d.
+ Assumes the coordinates are currently 3d.
- Args:
- i: int - index of a point
- delta: float - value to add to z-coord
- """
+ Args:
+ i: int - index of a point
+ delta: float - value to add to z-coord
+ """
- (x, y, z) = self.pos[i]
- self.pos[i] = (x, y, z + delta)
+ (x, y, z) = self.pos[i]
+ self.pos[i] = (x, y, z + delta)
class PolyArea(object):
- """Contains a Polygonal Area (polygon with possible holes).
-
- A polygon is a list of vertex ids, each an index given by
- a Points object. The list represents a CCW-oriented
- outer boundary (implicitly closed).
- If there are holes, they are lists of CW-oriented vertices
- that should be contained in the outer boundary.
- (So the left face of both the poly and the holes is
- the filled part.)
-
- Attributes:
- points: Points
- poly: list of vertex ids
- holes: list of lists of vertex ids (each a hole in poly)
- color: (float, float, float)- rgb color used to fill
- """
-
- def __init__(self, points = None, poly = None, holes = None):
- self.points = points if points else Points()
- self.poly = poly if poly else []
- self.holes = holes if holes else []
- self.color = (0.0, 0.0, 0.0)
-
- def AddHole(self, holepa):
- """Add a PolyArea's poly as a hole of self.
-
- Need to reverse the contour and
- adjust the the point indexes and self.points.
-
- Args:
- holepa: PolyArea
+ """Contains a Polygonal Area (polygon with possible holes).
+
+ A polygon is a list of vertex ids, each an index given by
+ a Points object. The list represents a CCW-oriented
+ outer boundary (implicitly closed).
+ If there are holes, they are lists of CW-oriented vertices
+ that should be contained in the outer boundary.
+ (So the left face of both the poly and the holes is
+ the filled part.)
+
+ Attributes:
+ points: Points
+ poly: list of vertex ids
+ holes: list of lists of vertex ids (each a hole in poly)
+ data: any - application data (can hold color, e.g.)
"""
- vmap = self.points.AddPoints(holepa.points)
- holepoly = [ vmap[i] for i in holepa.poly ]
- holepoly.reverse()
- self.holes.append(holepoly)
+ def __init__(self, points=None, poly=None, holes=None, data=None):
+ self.points = points if points else Points()
+ self.poly = poly if poly else []
+ self.holes = holes if holes else []
+ self.data = data
- def ContainsPoly(self, poly, points):
- """Tests if poly is contained within self.poly.
+ def AddHole(self, holepa):
+ """Add a PolyArea's poly as a hole of self.
- Args:
- poly: list of int - indices into points
- points: Points - maps to coords
- Returns:
- bool - True if poly is fully contained within self.poly
- """
+ Need to reverse the contour and
+ adjust the the point indexes and self.points.
- for v in poly:
- if PointInside(points.pos[v], self.poly, self.points) == -1:
- return False
- return True
+ Args:
+ holepa: PolyArea
+ """
- def Normal(self):
- """Returns the normal of the polyarea's main poly."""
+ vmap = self.points.AddPoints(holepa.points)
+ holepoly = [vmap[i] for i in holepa.poly]
+ holepoly.reverse()
+ self.holes.append(holepoly)
- pos = self.points.pos
- poly = self.poly
- if len(pos) == 0 or len(pos[0]) == 2 or len(poly) == 0:
- print("whoops, not enough info to calculate normal")
- return (0.0, 0.0, 1.0)
- return Newell(poly, self.points)
+ def ContainsPoly(self, poly, points):
+ """Tests if poly is contained within self.poly.
+ Args:
+ poly: list of int - indices into points
+ points: Points - maps to coords
+ Returns:
+ bool - True if poly is fully contained within self.poly
+ """
-class PolyAreas(object):
- """Contains a list of PolyAreas and a shared Points.
-
- Attributes:
- polyareas: list of PolyArea
- points: Points
- """
-
- def __init__(self):
- self.polyareas = []
- self.points = Points()
-
- def scale_and_center(self, scaled_side_target):
- """Adjust the coordinates of the polyareas so that
- it is centered at the origin and has its longest
- dimension scaled to be scaled_side_target."""
-
- if len(self.points.pos) == 0:
- return
- (minv, maxv) = self.bounds()
- maxside = max([ maxv[i]-minv[i] for i in range(2) ])
- if maxside > 0.0:
- scale = scaled_side_target / maxside
- else:
- scale = 1.0
- translate = [ -0.5*(maxv[i]+minv[i]) for i in range(2) ]
- dim = len(self.points.pos[0])
- if dim == 3:
- translate.append([0.0])
- for v in range(len(self.points.pos)):
- self.points.pos[v] = tuple([ scale*(self.points.pos[v][i] + translate[i]) for i in range(dim) ])
+ for v in poly:
+ if PointInside(points.pos[v], self.poly, self.points) == -1:
+ return False
+ return True
- def bounds(self):
- """Find bounding box of polyareas in xy.
+ def Normal(self):
+ """Returns the normal of the polyarea's main poly."""
- Returns:
- ([minx,miny],[maxx,maxy]) - all floats
- """
+ pos = self.points.pos
+ poly = self.poly
+ if len(pos) == 0 or len(pos[0]) == 2 or len(poly) == 0:
+ print("whoops, not enough info to calculate normal")
+ return (0.0, 0.0, 1.0)
+ return Newell(poly, self.points)
- huge = 1e100
- minv = [huge, huge]
- maxv = [-huge, -huge]
- for pa in self.polyareas:
- for face in [ pa.poly ] + pa.holes:
- for v in face:
- vcoords = self.points.pos[v]
- for i in range(2):
- if vcoords[i] < minv[i]:
- minv[i] = vcoords[i]
- if vcoords[i] > maxv[i]:
- maxv[i] = vcoords[i]
- if minv[0] == huge:
- minv = [0.0, 0.0]
- if maxv[0] == huge:
- maxv = [0.0, 0.0]
- return (minv, maxv)
+class PolyAreas(object):
+ """Contains a list of PolyAreas and a shared Points.
-class Model(object):
- """Contains a generic 3d model.
+ Attributes:
+ polyareas: list of PolyArea
+ points: Points
+ """
- A generic 3d model has vertices with 3d coordinates.
- Each vertex gets a 'vertex id', which is an index that
- can be used to refer to the vertex and can be used
- to retrieve the 3d coordinates of the point.
+ def __init__(self):
+ self.polyareas = []
+ self.points = Points()
+
+ def scale_and_center(self, scaled_side_target):
+ """Adjust the coordinates of the polyareas so that
+ it is centered at the origin and has its longest
+ dimension scaled to be scaled_side_target."""
+
+ if len(self.points.pos) == 0:
+ return
+ (minv, maxv) = self.bounds()
+ maxside = max([maxv[i] - minv[i] for i in range(2)])
+ if maxside > 0.0:
+ scale = scaled_side_target / maxside
+ else:
+ scale = 1.0
+ translate = [-0.5 * (maxv[i] + minv[i]) for i in range(2)]
+ dim = len(self.points.pos[0])
+ if dim == 3:
+ translate.append([0.0])
+ for v in range(len(self.points.pos)):
+ self.points.pos[v] = tuple([scale * (self.points.pos[v][i] + \
+ translate[i]) for i in range(dim)])
+
+ def bounds(self):
+ """Find bounding box of polyareas in xy.
+
+ Returns:
+ ([minx,miny],[maxx,maxy]) - all floats
+ """
+
+ huge = 1e100
+ minv = [huge, huge]
+ maxv = [-huge, -huge]
+ for pa in self.polyareas:
+ for face in [pa.poly] + pa.holes:
+ for v in face:
+ vcoords = self.points.pos[v]
+ for i in range(2):
+ if vcoords[i] < minv[i]:
+ minv[i] = vcoords[i]
+ if vcoords[i] > maxv[i]:
+ maxv[i] = vcoords[i]
+ if minv[0] == huge:
+ minv = [0.0, 0.0]
+ if maxv[0] == huge:
+ maxv = [0.0, 0.0]
+ return (minv, maxv)
- The actual visible part of the geometry are the faces,
- which are n-gons (n>2), specified by a vector of the
- n corner vertices.
- Faces may also have colors.
- Attributes:
- points: geom.Points - the 3d vertices
- faces: list of list of indices (each a CCW traversal of a face)
- colors: list of (float, float, float) - if present, is parallel to
- faces list and gives rgb colors of faces
- """
+class Model(object):
+ """Contains a generic 3d model.
+
+ A generic 3d model has vertices with 3d coordinates.
+ Each vertex gets a 'vertex id', which is an index that
+ can be used to refer to the vertex and can be used
+ to retrieve the 3d coordinates of the point.
+
+ The actual visible part of the geometry are the faces,
+ which are n-gons (n>2), specified by a vector of the
+ n corner vertices.
+ Faces may also have data associated with them,
+ and the data will be copied into newly created faces
+ from the most likely neighbor faces..
+
+ Attributes:
+ points: geom.Points - the 3d vertices
+ faces: list of list of indices (each a CCW traversal of a face)
+ face_data: list of any - if present, is parallel to
+ faces list and holds arbitrary data
+ """
- def __init__(self):
- self.points = Points()
- self.faces = []
- self.colors = []
+ def __init__(self):
+ self.points = Points()
+ self.faces = []
+ self.face_data = []
class Art(object):
- """Contains a vector art diagram.
+ """Contains a vector art diagram.
- Attributes:
- paths: list of Path objects
- """
+ Attributes:
+ paths: list of Path objects
+ """
- def __init__(self):
- self.paths = []
+ def __init__(self):
+ self.paths = []
class Paint(object):
- """A color or pattern to fill or stroke with.
+ """A color or pattern to fill or stroke with.
- For now, just do colors, but could later do
- patterns or images too.
+ For now, just do colors, but could later do
+ patterns or images too.
- Attributes:
- color: (r,g,b) triple of floats, 0.0=no color, 1.0=max color
- """
+ Attributes:
+ color: (r,g,b) triple of floats, 0.0=no color, 1.0=max color
+ """
- def __init__(self, r=0.0, g=0.0, b=0.0):
- self.color = (r, g, b)
+ def __init__(self, r=0.0, g=0.0, b=0.0):
+ self.color = (r, g, b)
- @staticmethod
- def CMYK(c, m, y, k):
- """Return Paint specified in CMYK model.
+ @staticmethod
+ def CMYK(c, m, y, k):
+ """Return Paint specified in CMYK model.
- Uses formula from 6.2.4 of PDF Reference.
+ Uses formula from 6.2.4 of PDF Reference.
- Args:
- c, m, y, k: float - in range [0, 1]
- Returns:
- Paint - with components in rgb form now
- """
+ Args:
+ c, m, y, k: float - in range [0, 1]
+ Returns:
+ Paint - with components in rgb form now
+ """
- return Paint(1.0 - min(1.0, c+k),
- 1.0 - min(1.0, m+k), 1.0 - min(1.0, y+k))
+ return Paint(1.0 - min(1.0, c + k),
+ 1.0 - min(1.0, m + k), 1.0 - min(1.0, y + k))
black_paint = Paint()
white_paint = Paint(1.0, 1.0, 1.0)
ColorDict = {
- 'aqua' : Paint(0.0, 1.0, 1.0),
- 'black' : Paint(0.0, 0.0, 0.0),
- 'blue' : Paint(0.0, 0.0, 1.0),
- 'fuchsia' : Paint(1.0, 0.0, 1.0),
- 'gray' : Paint(0.5, 0.5, 0.5),
- 'green' : Paint(0.0, 0.5, 0.0),
- 'lime' : Paint(0.0, 1.0, 0.0),
- 'maroon' : Paint(0.5, 0.0, 0.0),
- 'navy' : Paint(0.0, 0.0, 0.5),
- 'olive' : Paint(0.5, 0.5, 0.0),
- 'purple' : Paint(0.5, 0.0, 0.5),
- 'red' : Paint(1.0, 0.0, 0.0),
- 'silver' : Paint(0.75, 0.75, 0.75),
- 'teal' : Paint(0.0, 0.5, 0.5),
- 'white' : Paint(1.0, 1.0,1.0),
- 'yellow' : Paint(1.0, 1.0, 0.0)
+ 'aqua': Paint(0.0, 1.0, 1.0),
+ 'black': Paint(0.0, 0.0, 0.0),
+ 'blue': Paint(0.0, 0.0, 1.0),
+ 'fuchsia': Paint(1.0, 0.0, 1.0),
+ 'gray': Paint(0.5, 0.5, 0.5),
+ 'green': Paint(0.0, 0.5, 0.0),
+ 'lime': Paint(0.0, 1.0, 0.0),
+ 'maroon': Paint(0.5, 0.0, 0.0),
+ 'navy': Paint(0.0, 0.0, 0.5),
+ 'olive': Paint(0.5, 0.5, 0.0),
+ 'purple': Paint(0.5, 0.0, 0.5),
+ 'red': Paint(1.0, 0.0, 0.0),
+ 'silver': Paint(0.75, 0.75, 0.75),
+ 'teal': Paint(0.0, 0.5, 0.5),
+ 'white': Paint(1.0, 1.0, 1.0),
+ 'yellow': Paint(1.0, 1.0, 0.0)
}
class Path(object):
- """Represents a path in the PDF sense, with painting instructions.
-
- Attributes:
- subpaths: list of Subpath objects
- filled: True if path is to be filled
- fillevenodd: True if use even-odd rule to fill (else non-zero winding)
- stroked: True if path is to be stroked
- fillpaint: Paint to fill with
- strokepaint: Paint to stroke with
- """
+ """Represents a path in the PDF sense, with painting instructions.
+
+ Attributes:
+ subpaths: list of Subpath objects
+ filled: True if path is to be filled
+ fillevenodd: True if use even-odd rule to fill (else non-zero winding)
+ stroked: True if path is to be stroked
+ fillpaint: Paint to fill with
+ strokepaint: Paint to stroke with
+ """
- def __init__(self):
- self.subpaths = []
- self.filled = False
- self.fillevenodd = False
- self.stroked = False
- self.fillpaint = black_paint
- self.strokepaint = black_paint
+ def __init__(self):
+ self.subpaths = []
+ self.filled = False
+ self.fillevenodd = False
+ self.stroked = False
+ self.fillpaint = black_paint
+ self.strokepaint = black_paint
- def AddSubpath(self, subpath):
- """"Add a subpath."""
+ def AddSubpath(self, subpath):
+ """"Add a subpath."""
- self.subpaths.append(subpath)
+ self.subpaths.append(subpath)
- def Empty(self):
- """Returns True if this Path as no subpaths."""
+ def Empty(self):
+ """Returns True if this Path as no subpaths."""
- return not self.subpaths
+ return not self.subpaths
class Subpath(object):
- """Represents a subpath in PDF sense, either open or closed.
-
- We'll represent lines, bezier pieces, circular arc pieces
- as tuples with letters giving segment type in first position
- and coordinates (2-tuples of floats) in the other positions.
-
- Segment types:
- ('L', a, b) - line from a to b
- ('B', a, b, c, d) - cubic bezier from a to b, with control points c,d
- ('Q', a, b, c) - quadratic bezier from a to b, with 1 control point c
- ('A', a, b, rad, xrot, large-arc, ccw) - elliptical arc from a to b,
- with rad=(rx, ry) as radii, xrot is x-axis rotation in degrees,
- large-arc is True if arc should be >= 180 degrees,
- ccw is True if start->end follows counter-clockwise direction
- (see SVG spec); note that after rad,
- the rest are floats or bools, not coordinate pairs
- Note that s[1] and s[2] are the start and end points for any segment s.
-
- Attributes:
- segments: list of segment tuples (see above)
- closed: True if closed
- """
-
- def __init__(self):
- self.segments = []
- self.closed = False
-
- def Empty(self):
- """Returns True if this subpath as no segments."""
-
- return not self.segments
-
- def AddSegment(self, seg):
- """Add a segment."""
-
- self.segments.append(seg)
-
- @staticmethod
- def SegStart(s):
- """Return start point for segment.
-
- Args:
- s: a segment tuple
- Returns:
- (float, float): the coordinates of the segment's start point
+ """Represents a subpath in PDF sense, either open or closed.
+
+ We'll represent lines, bezier pieces, circular arc pieces
+ as tuples with letters giving segment type in first position
+ and coordinates (2-tuples of floats) in the other positions.
+
+ Segment types:
+ ('L', a, b) - line from a to b
+ ('B', a, b, c, d) - cubic bezier from a to b, with control points c,d
+ ('Q', a, b, c) - quadratic bezier from a to b, with 1 control point c
+ ('A', a, b, rad, xrot, large-arc, ccw) - elliptical arc from a to b,
+ with rad=(rx, ry) as radii, xrot is x-axis rotation in degrees,
+ large-arc is True if arc should be >= 180 degrees,
+ ccw is True if start->end follows counter-clockwise direction
+ (see SVG spec); note that after rad,
+ the rest are floats or bools, not coordinate pairs
+ Note that s[1] and s[2] are the start and end points for any segment s.
+
+ Attributes:
+ segments: list of segment tuples (see above)
+ closed: True if closed
"""
- return s[1]
+ def __init__(self):
+ self.segments = []
+ self.closed = False
- @staticmethod
- def SegEnd(s):
- """Return end point for segment.
+ def Empty(self):
+ """Returns True if this subpath as no segments."""
- Args:
- s: a segment tuple
- Returns:
- (float, float): the coordinates of the segment's end point
- """
+ return not self.segments
- return s[2]
+ def AddSegment(self, seg):
+ """Add a segment."""
+ self.segments.append(seg)
-class TransformMatrix(object):
- """Transformation matrix for 2d coordinates.
+ @staticmethod
+ def SegStart(s):
+ """Return start point for segment.
- The transform matrix is:
- [ a b 0 ]
- [ c d 0 ]
- [ e f 1 ]
- and coordinate tranformation is defined by:
- [x' y' 1] = [x y 1] x TransformMatrix
+ Args:
+ s: a segment tuple
+ Returns:
+ (float, float): the coordinates of the segment's start point
+ """
- Attributes:
- a, b, c, d, e, f: floats
- """
+ return s[1]
- def __init__(self, a=1.0, b=0.0, c=0.0, d=1.0, e=0.0, f=0.0):
- self.a = a
- self.b = b
- self.c = c
- self.d =d
- self.e = e
- self.f = f
+ @staticmethod
+ def SegEnd(s):
+ """Return end point for segment.
- def __str__(self):
- return str([self.a, self.b, self.c, self.d, self.e, self.f])
+ Args:
+ s: a segment tuple
+ Returns:
+ (float, float): the coordinates of the segment's end point
+ """
- def Copy(self):
- """Return a copy of this matrix."""
+ return s[2]
- return TransformMatrix(self.a, self.b, self.c, self.d, self.e, self.f)
- def ComposeTransform(self, a, b, c, d, e, f):
- """Apply the transform given the the arguments on top of this one.
+class TransformMatrix(object):
+ """Transformation matrix for 2d coordinates.
- This is accomplished by returning t x sel
- where t is the transform matrix that would be formed from the args.
+ The transform matrix is:
+ [ a b 0 ]
+ [ c d 0 ]
+ [ e f 1 ]
+ and coordinate tranformation is defined by:
+ [x' y' 1] = [x y 1] x TransformMatrix
- Arguments:
- a, b, c, d, e, f: float - defines a composing TransformMatrix
+ Attributes:
+ a, b, c, d, e, f: floats
"""
- newa = a * self.a + b*self.c
- newb = a*self.b + b*self.d
- newc = c*self.a + d*self.c
- newd = c*self.b + d*self.d
- newe = e*self.a + f*self.c + self.e
- newf = e*self.b + f*self.d + self.f
- self.a = newa
- self.b = newb
- self.c = newc
- self.d = newd
- self.e = newe
- self.f = newf
-
- def Apply(self, pt):
- """Return the result of applying this tranform to pt = (x,y).
-
- Arguments:
- (x, y) : (float, float)
- Returns:
- (x', y'): 2-tuple of floats, the result of [x y 1] x self
- """
+ def __init__(self, a=1.0, b=0.0, c=0.0, d=1.0, e=0.0, f=0.0):
+ self.a = a
+ self.b = b
+ self.c = c
+ self.d = d
+ self.e = e
+ self.f = f
- (x, y) = pt
- return (self.a*x + self.c*y + self.e, self.b*x + self.d*y + self.f)
+ def __str__(self):
+ return str([self.a, self.b, self.c, self.d, self.e, self.f])
+
+ def Copy(self):
+ """Return a copy of this matrix."""
+
+ return TransformMatrix(self.a, self.b, self.c, self.d, self.e, self.f)
+
+ def ComposeTransform(self, a, b, c, d, e, f):
+ """Apply the transform given the the arguments on top of this one.
+
+ This is accomplished by returning t x sel
+ where t is the transform matrix that would be formed from the args.
+
+ Arguments:
+ a, b, c, d, e, f: float - defines a composing TransformMatrix
+ """
+
+ newa = a * self.a + b * self.c
+ newb = a * self.b + b * self.d
+ newc = c * self.a + d * self.c
+ newd = c * self.b + d * self.d
+ newe = e * self.a + f * self.c + self.e
+ newf = e * self.b + f * self.d + self.f
+ self.a = newa
+ self.b = newb
+ self.c = newc
+ self.d = newd
+ self.e = newe
+ self.f = newf
+
+ def Apply(self, pt):
+ """Return the result of applying this tranform to pt = (x,y).
+
+ Arguments:
+ (x, y) : (float, float)
+ Returns:
+ (x', y'): 2-tuple of floats, the result of [x y 1] x self
+ """
+
+ (x, y) = pt
+ return (self.a * x + self.c * y + self.e, \
+ self.b * x + self.d * y + self.f)
def ApproxEqualPoints(p, q):
- """Return True if p and q are approximately the same points.
+ """Return True if p and q are approximately the same points.
- Args:
- p: n-tuple of float
- q: n-tuple of float
- Returns:
- bool - True if the 1-norm <= DISTTOL
- """
+ Args:
+ p: n-tuple of float
+ q: n-tuple of float
+ Returns:
+ bool - True if the 1-norm <= DISTTOL
+ """
- for i in range(len(p)):
- if abs(p[i] - q[i]) > DISTTOL:
- return False
- return True
+ for i in range(len(p)):
+ if abs(p[i] - q[i]) > DISTTOL:
+ return False
+ return True
def PointInside(v, a, points):
- """Return 1, 0, or -1 as v is inside, on, or outside polygon.
-
- Cf. Eric Haines ptinpoly in Graphics Gems IV.
-
- Args:
- v : (float, float) or (float, float, float) - coordinates of a point
- a : list of vertex indices defining polygon (assumed CCW)
- points: Points - to get coordinates for polygon
- Returns:
- 1, 0, -1: as v is inside, on, or outside polygon a
- """
-
- (xv, yv) = (v[0], v[1])
- vlast = points.pos[a[-1]]
- (x0, y0) = (vlast[0], vlast[1])
- if x0 == xv and y0 == yv:
- return 0
- yflag0 = y0 > yv
- inside = False
- n = len(a)
- for i in range(0, n):
- vi = points.pos[a[i]]
- (x1, y1) = (vi[0], vi[1])
- if x1 == xv and y1 == yv:
- return 0
- yflag1 = y1 > yv
- if yflag0 != yflag1:
- xflag0 = x0 > xv
- xflag1 = x1 > xv
- if xflag0 == xflag1:
- if xflag0:
- inside = not inside
- else:
- z = x1 - (y1-yv)*(x0-x1)/(y0-y1)
- if z >= xv:
- inside = not inside
- x0 = x1
- y0 = y1
- yflag0 = yflag1
- if inside:
- return 1
- else:
- return -1
+ """Return 1, 0, or -1 as v is inside, on, or outside polygon.
+
+ Cf. Eric Haines ptinpoly in Graphics Gems IV.
+
+ Args:
+ v : (float, float) or (float, float, float) - coordinates of a point
+ a : list of vertex indices defining polygon (assumed CCW)
+ points: Points - to get coordinates for polygon
+ Returns:
+ 1, 0, -1: as v is inside, on, or outside polygon a
+ """
+
+ (xv, yv) = (v[0], v[1])
+ vlast = points.pos[a[-1]]
+ (x0, y0) = (vlast[0], vlast[1])
+ if x0 == xv and y0 == yv:
+ return 0
+ yflag0 = y0 > yv
+ inside = False
+ n = len(a)
+ for i in range(0, n):
+ vi = points.pos[a[i]]
+ (x1, y1) = (vi[0], vi[1])
+ if x1 == xv and y1 == yv:
+ return 0
+ yflag1 = y1 > yv
+ if yflag0 != yflag1:
+ xflag0 = x0 > xv
+ xflag1 = x1 > xv
+ if xflag0 == xflag1:
+ if xflag0:
+ inside = not inside
+ else:
+ z = x1 - (y1 - yv) * (x0 - x1) / (y0 - y1)
+ if z >= xv:
+ inside = not inside
+ x0 = x1
+ y0 = y1
+ yflag0 = yflag1
+ if inside:
+ return 1
+ else:
+ return -1
def SignedArea(polygon, points):
- """Return the area of the polgon, positive if CCW, negative if CW.
+ """Return the area of the polgon, positive if CCW, negative if CW.
- Args:
- polygon: list of vertex indices
- points: Points
- Returns:
- float - area of polygon, positive if it was CCW, else negative
- """
+ Args:
+ polygon: list of vertex indices
+ points: Points
+ Returns:
+ float - area of polygon, positive if it was CCW, else negative
+ """
- a = 0.0
- n = len(polygon)
- for i in range(0, n):
- u = points.pos[polygon[i]]
- v = points.pos[polygon[(i+1) % n]]
- a += u[0]*v[1] - u[1]*v[0]
- return 0.5*a
+ a = 0.0
+ n = len(polygon)
+ for i in range(0, n):
+ u = points.pos[polygon[i]]
+ v = points.pos[polygon[(i + 1) % n]]
+ a += u[0] * v[1] - u[1] * v[0]
+ return 0.5 * a
def VecAdd(a, b):
- """Return vector a-b.
+ """Return vector a-b.
- Args:
- a: n-tuple of floats
- b: n-tuple of floats
- Returns:
- n-tuple of floats - pairwise addition a+b
- """
+ Args:
+ a: n-tuple of floats
+ b: n-tuple of floats
+ Returns:
+ n-tuple of floats - pairwise addition a+b
+ """
- n = len(a)
- assert(n == len(b))
- return tuple([ a[i]+b[i] for i in range(n)])
+ n = len(a)
+ assert(n == len(b))
+ return tuple([a[i] + b[i] for i in range(n)])
def VecSub(a, b):
- """Return vector a-b.
+ """Return vector a-b.
+
+ Args:
+ a: n-tuple of floats
+ b: n-tuple of floats
+ Returns:
+ n-tuple of floats - pairwise subtraction a-b
+ """
- Args:
- a: n-tuple of floats
- b: n-tuple of floats
- Returns:
- n-tuple of floats - pairwise subtraction a-b
- """
+ n = len(a)
+ assert(n == len(b))
+ return tuple([a[i] - b[i] for i in range(n)])
+
+
+def VecDot(a, b):
+ """Return the dot product of two vectors.
+
+ Args:
+ a: n-tuple of floats
+ b: n-tuple of floats
+ Returns:
+ n-tuple of floats - dot product of a and b
+ """
- n = len(a)
- assert(n == len(b))
- return tuple([ a[i]-b[i] for i in range(n)])
+ n = len(a)
+ assert(n == len(b))
+ sum = 0.0
+ for i in range(n):
+ sum += a[i] * b[i]
+ return sum
def VecLen(a):
- """Return the Euclidean lenght of the argument vector.
+ """Return the Euclidean lenght of the argument vector.
- Args:
- a: n-tuple of floats
- Returns:
- float: the 2-norm of a
- """
+ Args:
+ a: n-tuple of floats
+ Returns:
+ float: the 2-norm of a
+ """
- s = 0.0
- for v in a:
- s += v*v
- return math.sqrt(s)
+ s = 0.0
+ for v in a:
+ s += v * v
+ return math.sqrt(s)
def Newell(poly, points):
- """Use Newell method to find polygon normal.
-
- Assume poly has length at least 3 and points are 3d.
-
- Args:
- poly: list of int - indices into points.pos
- points: Points - assumed 3d
- Returns:
- (float, float, float) - the average normal
- """
-
- sumx = 0.0
- sumy = 0.0
- sumz = 0.0
- n = len(poly)
- pos = points.pos
- for i, ai in enumerate(poly):
- bi = poly[(i+1) % n]
- a = pos[ai]
- b = pos[bi]
- sumx += (a[1]-b[1]) * (a[2]+b[2])
- sumy += (a[2]-b[2]) * (a[0]+b[0])
- sumz += (a[0]-b[0]) * (a[1]+b[1])
- return Norm3(sumx, sumy, sumz)
+ """Use Newell method to find polygon normal.
+
+ Assume poly has length at least 3 and points are 3d.
+
+ Args:
+ poly: list of int - indices into points.pos
+ points: Points - assumed 3d
+ Returns:
+ (float, float, float) - the average normal
+ """
+
+ sumx = 0.0
+ sumy = 0.0
+ sumz = 0.0
+ n = len(poly)
+ pos = points.pos
+ for i, ai in enumerate(poly):
+ bi = poly[(i + 1) % n]
+ a = pos[ai]
+ b = pos[bi]
+ sumx += (a[1] - b[1]) * (a[2] + b[2])
+ sumy += (a[2] - b[2]) * (a[0] + b[0])
+ sumz += (a[0] - b[0]) * (a[1] + b[1])
+ return Norm3(sumx, sumy, sumz)
def Norm3(x, y, z):
- """Return vector (x,y,z) normalized by dividing by squared length.
- Return (0.0, 0.0, 1.0) if the result is undefined."""
- sqrlen = x * x + y * y + z * z
- if sqrlen < 1e-100:
- return (0.0, 0.0, 1.0)
- else:
- try:
- d = math.sqrt(sqrlen)
- return (x / d, y / d, z / d)
- except:
- return (0.0, 0.0, 1.0)
+ """Return vector (x,y,z) normalized by dividing by squared length.
+ Return (0.0, 0.0, 1.0) if the result is undefined."""
+ sqrlen = x * x + y * y + z * z
+ if sqrlen < 1e-100:
+ return (0.0, 0.0, 1.0)
+ else:
+ try:
+ d = math.sqrt(sqrlen)
+ return (x / d, y / d, z / d)
+ except:
+ return (0.0, 0.0, 1.0)
# We're using right-hand coord system, where
# forefinger=x, middle=y, thumb=z on right hand.
# Then, e.g., (1,0,0) x (0,1,0) = (0,0,1)
def Cross3(a, b):
- """Return the cross product of two vectors, a x b."""
+ """Return the cross product of two vectors, a x b."""
- (ax, ay, az) = a
- (bx, by, bz) = b
- return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
+ (ax, ay, az) = a
+ (bx, by, bz) = b
+ return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
def MulPoint3(p, m):
- """Return matrix multiplication of p times m
- where m is a 4x3 matrix and p is a 3d point, extended with 1."""
-
- (x, y, z) = p
- return (x * m[0] + y * m[3] + z * m[6] + m[9],
- x * m[1] + y * m[4] + z * m[7] + m[10],
- x * m[2] + y * m[5] + z * m[8] + m[11])
+ """Return matrix multiplication of p times m
+ where m is a 4x3 matrix and p is a 3d point, extended with 1."""
+ (x, y, z) = p
+ return (x * m[0] + y * m[3] + z * m[6] + m[9],
+ x * m[1] + y * m[4] + z * m[7] + m[10],
+ x * m[2] + y * m[5] + z * m[8] + m[11])
diff --git a/mesh_inset/model.py b/mesh_inset/model.py
index 4e3b0f664211205db385fb238bb936e9eaa5097a..a3eb2aacc5d68141459a608a145e44612f4c6a2c 100644
--- a/mesh_inset/model.py
+++ b/mesh_inset/model.py
@@ -16,6 +16,8 @@
#
# ##### END GPL LICENSE BLOCK #####
+# <pep8 compliant>
+
"""Manipulations of Models.
"""
@@ -28,480 +30,546 @@ import math
def PolyAreasToModel(polyareas, bevel_amount, bevel_pitch, quadrangulate):
- """Convert a PolyAreas into a Model object.
+ """Convert a PolyAreas into a Model object.
- Assumes polyareas are in xy plane.
+ Assumes polyareas are in xy plane.
- Args:
- polyareas: geom.PolyAreas
- bevel_amount: float - if > 0, amount of bevel
- bevel_pitch: float - if > 0, angle in radians of bevel
- quadrangulate: bool - should n-gons be quadrangulated?
- Returns:
- geom.Model
- """
+ Args:
+ polyareas: geom.PolyAreas
+ bevel_amount: float - if > 0, amount of bevel
+ bevel_pitch: float - if > 0, angle in radians of bevel
+ quadrangulate: bool - should n-gons be quadrangulated?
+ Returns:
+ geom.Model
+ """
- m = geom.Model()
- if not polyareas:
+ m = geom.Model()
+ if not polyareas:
+ return m
+ polyareas.points.AddZCoord(0.0)
+ m.points = polyareas.points
+ for pa in polyareas.polyareas:
+ PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate)
return m
- polyareas.points.AddZCoord(0.0)
- m.points = polyareas.points
- for pa in polyareas.polyareas:
- PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate)
- return m
def PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate):
- if bevel_amount > 0.0:
- BevelPolyAreaInModel(m, pa, bevel_amount, bevel_pitch, quadrangulate)
- elif quadrangulate:
- if len(pa.poly) == 0:
- return
- qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points)
- m.faces.extend(qpa)
- m.colors.extend([ pa.color ] * len(qpa))
- else:
- m.faces.append(pa.poly)
- # TODO: just the first part of QuadrangulateFaceWithHoles, to join
- # holes to outer poly
- m.colors.append(pa.color)
+ if bevel_amount > 0.0:
+ BevelPolyAreaInModel(m, pa, bevel_amount, bevel_pitch, quadrangulate,
+ False)
+ elif quadrangulate:
+ if len(pa.poly) == 0:
+ return
+ qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points)
+ m.faces.extend(qpa)
+ m.face_data.extend([pa.data] * len(qpa))
+ else:
+ m.faces.append(pa.poly)
+ # TODO: just the first part of QuadrangulateFaceWithHoles, to join
+ # holes to outer poly
+ m.face_data.append(pa.data)
def ExtrudePolyAreasInModel(mdl, polyareas, depth, cap_back):
- """Extrude the boundaries given by polyareas by -depth in z.
-
- Assumes polyareas are in xy plane.
-
- Arguments:
- mdl: geom.Model - where to do extrusion
- polyareas: geom.Polyareas
- depth: float
- cap_back: bool - if True, cap off the back
- Side Effects:
- For all edges in polys in polyareas, make quads in Model
- extending those edges by depth in the negative z direction.
- The color will be the color of the face that the edge is part of.
- """
-
- for pa in polyareas.polyareas:
- back_poly = _ExtrudePoly(mdl, pa.poly, depth, pa.color, True)
- back_holes = []
- for p in pa.holes:
- back_holes.append(_ExtrudePoly(mdl, p, depth, pa.color, False))
- if cap_back:
- qpa = triquad.QuadrangulateFaceWithHoles(back_poly, back_holes,
- polyareas.points)
- # need to reverse each poly to get normals pointing down
- for i, p in enumerate(qpa):
- t = list(p)
- t.reverse()
- qpa[i] = tuple(t)
- model.faces.extend(qpa)
- model.colors.extend([pa.color] * len(qpa))
-
-
-def _ExtrudePoly(mdl, poly, depth, color, isccw):
- """Extrude the poly by -depth in z
-
- Arguments:
- mdl: geom.Model - where to do extrusion
- poly: list of vertex indices
- depth: float
- color: tuple of three floats
- isccw: True if counter-clockwise
- Side Effects
- For all edges in poly, make quads in Model
- extending those edges by depth in the negative z direction.
- The color will be the color of the face that the edge is part of.
- Returns:
- list of int - vertices for extruded poly
- """
-
- if len(poly) < 2:
- return
- extruded_poly = []
- points = mdl.points
- if isccw:
- incr = 1
- else:
- incr = -1
- for i, v in enumerate(poly):
- vnext = poly[(i+incr) % len(poly)]
- (x0,y0,z0) = points.pos[v]
- (x1,y1,z1) = points.pos[vnext]
- vextrude = points.AddPoint((x0,y0,z0-depth))
- vnextextrude = points.AddPoint((x1,y1,z1-depth))
+ """Extrude the boundaries given by polyareas by -depth in z.
+
+ Assumes polyareas are in xy plane.
+
+ Arguments:
+ mdl: geom.Model - where to do extrusion
+ polyareas: geom.Polyareas
+ depth: float
+ cap_back: bool - if True, cap off the back
+ Side Effects:
+ For all edges in polys in polyareas, make quads in Model
+ extending those edges by depth in the negative z direction.
+ The application data will be the data of the face that the edge
+ is part of.
+ """
+
+ for pa in polyareas.polyareas:
+ back_poly = _ExtrudePoly(mdl, pa.poly, depth, pa.data, True)
+ back_holes = []
+ for p in pa.holes:
+ back_holes.append(_ExtrudePoly(mdl, p, depth, pa.data, False))
+ if cap_back:
+ qpa = triquad.QuadrangulateFaceWithHoles(back_poly, back_holes,
+ polyareas.points)
+ # need to reverse each poly to get normals pointing down
+ for i, p in enumerate(qpa):
+ t = list(p)
+ t.reverse()
+ qpa[i] = tuple(t)
+ mdl.faces.extend(qpa)
+ mdl.face_data.extend([pa.data] * len(qpa))
+
+
+def _ExtrudePoly(mdl, poly, depth, data, isccw):
+ """Extrude the poly by -depth in z
+
+ Arguments:
+ mdl: geom.Model - where to do extrusion
+ poly: list of vertex indices
+ depth: float
+ data: application data
+ isccw: True if counter-clockwise
+ Side Effects
+ For all edges in poly, make quads in Model
+ extending those edges by depth in the negative z direction.
+ The application data will be the data of the face that the edge
+ is part of.
+ Returns:
+ list of int - vertices for extruded poly
+ """
+
+ if len(poly) < 2:
+ return
+ extruded_poly = []
+ points = mdl.points
if isccw:
- sideface = [v, vextrude, vnextextrude, vnext]
+ incr = 1
else:
- sideface = [v, vnext, vnextextrude, vextrude]
- mdl.faces.append(sideface)
- mdl.colors.append(color)
- extruded_poly.append(vextrude)
- return extruded_poly
+ incr = -1
+ for i, v in enumerate(poly):
+ vnext = poly[(i + incr) % len(poly)]
+ (x0, y0, z0) = points.pos[v]
+ (x1, y1, z1) = points.pos[vnext]
+ vextrude = points.AddPoint((x0, y0, z0 - depth))
+ vnextextrude = points.AddPoint((x1, y1, z1 - depth))
+ if isccw:
+ sideface = [v, vextrude, vnextextrude, vnext]
+ else:
+ sideface = [v, vnext, vnextextrude, vextrude]
+ mdl.faces.append(sideface)
+ mdl.face_data.append(data)
+ extruded_poly.append(vextrude)
+ return extruded_poly
def BevelPolyAreaInModel(mdl, polyarea,
- bevel_amount, bevel_pitch, quadrangulate):
- """Bevel the interior of polyarea in model.
-
- This does smart beveling: advancing edges are merged
- rather than doing an 'overlap'. Advancing edges that
- hit an opposite edge result in a split into two beveled areas.
-
- If the polyarea is not in the xy plane, do the work in a
- transformed model, and then transfer the changes back.
-
- Arguments:
- mdl: geom.Model - where to do bevel
- polyarea geom.PolyArea - area to bevel into
- bevel_amount: float - if > 0, amount of bevel
- bevel_pitch: float - if > 0, angle in radians of bevel
- quadrangulate: bool - should n-gons be quadrangulated?
- Side Effects:
- Faces and points are added to model to model the
- bevel and the interior of the polyareas.
- """
-
- pa_norm = polyarea.Normal()
- if pa_norm == (0.0, 0.0, 1.0):
- m = mdl
- pa_rot = polyarea
- else:
- (pa_rot, inv_rot, inv_map) = _RotatedPolyAreaToXY(polyarea, pa_norm)
- # don't actually have to add the original faces into model, just their points.
- m = geom.Model()
- m.points = pa_rot.points
- vspeed = math.tan(bevel_pitch)
- off = offset.Offset(pa_rot, 0.0, vspeed)
- off.Build(bevel_amount)
- inner_pas = AddOffsetFacesToModel(m, off, polyarea.color)
- for pa in inner_pas.polyareas:
- if quadrangulate:
- if len(pa.poly) == 0:
- continue
- qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points)
- m.faces.extend(qpa)
- m.colors.extend([ pa.color ] * len(qpa))
+ bevel_amount, bevel_pitch, quadrangulate, as_percent):
+ """Bevel the interior of polyarea in model.
+
+ This does smart beveling: advancing edges are merged
+ rather than doing an 'overlap'. Advancing edges that
+ hit an opposite edge result in a split into two beveled areas.
+
+ If the polyarea is not in the xy plane, do the work in a
+ transformed model, and then transfer the changes back.
+
+ Arguments:
+ mdl: geom.Model - where to do bevel
+ polyarea geom.PolyArea - area to bevel into
+ bevel_amount: float - if > 0, amount of bevel
+ bevel_pitch: float - if > 0, angle in radians of bevel
+ quadrangulate: bool - should n-gons be quadrangulated?
+ as_percent: bool - if True, interpret amount as percent of max
+ Side Effects:
+ Faces and points are added to model to model the
+ bevel and the interior of the polyareas.
+ """
+
+ pa_norm = polyarea.Normal()
+ if pa_norm == (0.0, 0.0, 1.0):
+ m = mdl
+ pa_rot = polyarea
else:
- m.faces.append(pa.poly)
- m.colors.append(pa.color)
- if m != mdl:
- _AddTransformedPolysToModel(mdl, m.faces, m.points, inv_rot, inv_map)
-
-
-def AddOffsetFacesToModel(mdl, off, color = (0.0, 0.0, 0.0)):
- """Add the faces due to an offset into model.
-
- Returns the remaining interiors of the offset as a PolyAreas.
-
- Args:
- mdl: geom.Model - model to add offset faces into
- off: offset.Offset
- color: (float, float, float) - color to make the faces
- Returns:
- geom.PolyAreas
- """
-
- mdl.points = off.polyarea.points
- assert(len(mdl.points.pos) == 0 or len(mdl.points.pos[0]) == 3)
- o = off
- ostack = [ ]
- while o:
- if o.endtime != 0.0:
- for face in o.facespokes:
- n = len(face)
- for i, spoke in enumerate(face):
- nextspoke = face[(i+1) % n]
- v0 = spoke.origin
- v1 = nextspoke.origin
- v2 = nextspoke.dest
- v3 = spoke.dest
- if v2 == v3:
- mface = [v0, v1, v2]
- else:
- mface = [v0, v1, v2, v3]
- mdl.faces.append(mface)
- mdl.colors.append(color)
- ostack.extend(o.inneroffsets)
- if ostack:
- o = ostack.pop()
+ (pa_rot, inv_rot, inv_map) = _RotatedPolyAreaToXY(polyarea, pa_norm)
+ # don't have to add the original faces into model, just their points.
+ m = geom.Model()
+ m.points = pa_rot.points
+ vspeed = math.tan(bevel_pitch)
+ off = offset.Offset(pa_rot, 0.0, vspeed)
+ if as_percent:
+ bevel_amount = bevel_amount * off.MaxAmount() / 100.0
+ off.Build(bevel_amount)
+ inner_pas = AddOffsetFacesToModel(m, off, polyarea.data)
+ for pa in inner_pas.polyareas:
+ if quadrangulate:
+ if len(pa.poly) == 0:
+ continue
+ qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes,
+ pa.points)
+ m.faces.extend(qpa)
+ m.face_data.extend([pa.data] * len(qpa))
+ else:
+ m.faces.append(pa.poly)
+ m.face_data.append(pa.data)
+ if m != mdl:
+ _AddTransformedPolysToModel(mdl, m.faces, m.points, m.face_data,
+ inv_rot, inv_map)
+
+
+def AddOffsetFacesToModel(mdl, off, data=None):
+ """Add the faces due to an offset into model.
+
+ Returns the remaining interiors of the offset as a PolyAreas.
+
+ Args:
+ mdl: geom.Model - model to add offset faces into
+ off: offset.Offset
+ data: any - application data to be copied to the faces
+ Returns:
+ geom.PolyAreas
+ """
+
+ mdl.points = off.polyarea.points
+ assert(len(mdl.points.pos) == 0 or len(mdl.points.pos[0]) == 3)
+ o = off
+ ostack = []
+ while o:
+ if o.endtime != 0.0:
+ for face in o.facespokes:
+ n = len(face)
+ for i, spoke in enumerate(face):
+ nextspoke = face[(i + 1) % n]
+ v0 = spoke.origin
+ v1 = nextspoke.origin
+ v2 = nextspoke.dest
+ v3 = spoke.dest
+ if v2 == v3:
+ mface = [v0, v1, v2]
+ else:
+ mface = [v0, v1, v2, v3]
+ mdl.faces.append(mface)
+ mdl.face_data.append(data)
+ ostack.extend(o.inneroffsets)
+ if ostack:
+ o = ostack.pop()
+ else:
+ o = None
+ return off.InnerPolyAreas()
+
+
+def BevelSelectionInModel(mdl, bevel_amount, bevel_pitch, quadrangulate,
+ as_region, as_percent):
+ """Bevel all the faces in the model, perhaps as one region.
+
+ If as_region is False, each face is beveled individually,
+ otherwise regions of contiguous faces are merged into
+ PolyAreas and beveled as a whole.
+
+ TODO: something if extracted PolyAreas are not approximately
+ planar.
+
+ Args:
+ mdl: geom.Model
+ bevel_amount: float - amount to inset
+ bevel_pitch: float - angle of bevel side
+ quadrangulate: bool - should insides be quadrangulated?
+ as_region: bool - should faces be merged into regions?
+ as_percent: bool - should amount be interpreted as a percent
+ of the maximum amount (if True) or an absolute amount?
+ Side effect:
+ Beveling faces will be added to the model
+ """
+
+ pas = []
+ if as_region:
+ pas = RegionToPolyAreas(mdl.faces, mdl.points, mdl.face_data)
else:
- o = None
- return off.InnerPolyAreas()
-
-
-def BevelSelectionInModel(mdl, selected_faces,
- bevel_amount, bevel_pitch, quadrangulate, as_region):
- """Bevel the selected faces in the model.
-
- If as_region is False, each face is beveled individually,
- otherwise regions of contiguous faces are merged into
- PolyAreas and beveled as a whole.
-
- TODO: something if extracted PolyAreas are not approximately
- planar.
-
- Args:
- mdl: geom.Model
- selected_faces: list of list of int
- bevel_amount: float - amount to inset
- bevel_pitch: float - angle of bevel side
- quadrangulate: bool - should insides be quadrangulated?
- as_region: bool - should faces be merged into regions?
- Side effect:
- Beveling faces will be added to the model
- """
-
- pas = []
- if as_region:
- pas = RegionToPolyAreas(selected_faces, mdl.points)
- else:
- for face in selected_faces:
- pas.append(geom.PolyArea(mdl.points, face))
- for pa in pas:
- BevelPolyAreaInModel(mdl, pa,
- bevel_amount, bevel_pitch, quadrangulate)
-
-
-def RegionToPolyAreas(faces, points):
- """Find polygonal outlines induced by union of faces.
-
- Finds the polygons formed by boundary edges (those not
- sharing an edge with another face in region_faces), and
- turns those into PolyAreas.
- In the general case, there will be holes inside.
-
- Args:
- faces: list of list of int - each sublist is a face (indices into points)
- points: geom.Points - gives coordinates for vertices
- Returns:
- list of geom.PolyArea
- """
-
- ans = []
- (edges, vtoe) = _GetEdgeData(faces)
- (face_adj, is_interior_edge) = _GetFaceGraph(faces, edges, vtoe, points)
- (components, ftoc) = _FindFaceGraphComponents(faces, face_adj)
- for c in range(len(components)):
- boundary_edges = set()
- vstobe = dict()
- for e, ((vs, ve), f) in enumerate(edges):
- if ftoc[f] != c or is_interior_edge[e]:
- continue
- boundary_edges.add(e)
- vstobe[vs] = e
- polys = []
- while boundary_edges:
- e = boundary_edges.pop()
- ((vstart, ve), _) = edges[e]
- poly = [ vstart, ve ]
- while ve != vstart:
- if ve not in vstobe:
- print("whoops, couldn't close boundary")
- break
- nexte = vstobe[ve]
- ((_, ve), _) = edges[nexte]
- boundary_edges.remove(nexte)
- if ve != vstart:
- poly.append(ve)
- polys.append(poly)
- if len(polys) == 0:
- # can happen if an entire closed polytope is given
- # at least until we do an edge check
- return []
- elif len(polys) == 1:
- ans.append(geom.PolyArea(points, polys[0]))
- else:
- outerf = _FindOuterPoly(polys, points)
- pa = geom.PolyArea(points, polys[outerf])
- pa.holes = [ polys[i] for i in range(len(polys)) if i != outerf ]
- ans.append(pa)
- return ans
+ for f, face in enumerate(mdl.faces):
+ pas.append(geom.PolyArea(mdl.points, face, [],
+ mdl.face_data[f]))
+ for pa in pas:
+ BevelPolyAreaInModel(mdl, pa,
+ bevel_amount, bevel_pitch, quadrangulate, as_percent)
+
+
+def RegionToPolyAreas(faces, points, data):
+ """Find polygonal outlines induced by union of faces.
+
+ Finds the polygons formed by boundary edges (those not
+ sharing an edge with another face in region_faces), and
+ turns those into PolyAreas.
+ In the general case, there will be holes inside.
+ We want to associate data with the region PolyAreas.
+ Just choose a representative element of data[] when
+ more than one face is combined into a PolyArea.
+
+ Args:
+ faces: list of list of int - each sublist is a face (indices into points)
+ points: geom.Points - gives coordinates for vertices
+ data: list of any - parallel to faces, app data to put in PolyAreas
+ Returns:
+ list of geom.PolyArea
+ """
+
+ ans = []
+ (edges, vtoe) = _GetEdgeData(faces)
+ (face_adj, is_interior_edge) = _GetFaceGraph(faces, edges, vtoe, points)
+ (components, ftoc) = _FindFaceGraphComponents(faces, face_adj)
+ for c in range(len(components)):
+ boundary_edges = set()
+ betodata = dict()
+ vstobe = dict()
+ for e, ((vs, ve), f) in enumerate(edges):
+ if ftoc[f] != c or is_interior_edge[e]:
+ continue
+ boundary_edges.add(e)
+ # vstobe[v] is set of edges leaving v
+ # (could be more than one if boundary touches itself at a vertex)
+ if vs in vstobe:
+ vstobe[vs].append(e)
+ else:
+ vstobe[vs] = [e]
+ betodata[(vs, ve)] = data[f]
+ polys = []
+ poly_data = []
+ while boundary_edges:
+ e = boundary_edges.pop()
+ ((vstart, ve), face_i) = edges[e]
+ poly = [vstart, ve]
+ datum = betodata[(vstart, ve)]
+ while ve != vstart:
+ if ve not in vstobe:
+ print("whoops, couldn't close boundary")
+ break
+ nextes = vstobe[ve]
+ if len(nextes) == 1:
+ nexte = nextes[0]
+ else:
+ # find a next edge with face index face_i
+ # TODO: this is not guaranteed to work,
+ # as continuation edge may have been for a different
+ # face that is now combined with face_i by erasing
+ # interior edges. Find a better algorithm here.
+ nexte = -1
+ for ne_cand in nextes:
+ if edges[ne_cand][1] == face_i:
+ nexte = ne_cand
+ break
+ if nexte == -1:
+ # case mentioned in TODO may have happened;
+ # just choose any nexte - may mess things up
+ nexte = nextes[0]
+ ((_, ve), face_i) = edges[nexte]
+ if nexte not in boundary_edges:
+ print("whoops, nexte not a boundary edge", nexte)
+ break
+ boundary_edges.remove(nexte)
+ if ve != vstart:
+ poly.append(ve)
+ polys.append(poly)
+ poly_data.append(datum)
+ if len(polys) == 0:
+ # can happen if an entire closed polytope is given
+ # at least until we do an edge check
+ return []
+ elif len(polys) == 1:
+ ans.append(geom.PolyArea(points, polys[0], [], poly_data[0]))
+ else:
+ outerf = _FindOuterPoly(polys, points, faces)
+ pa = geom.PolyArea(points, polys[outerf], [], poly_data[outerf])
+ pa.holes = [polys[i] for i in range(len(polys)) if i != outerf]
+ ans.append(pa)
+ return ans
def _GetEdgeData(faces):
- """Find edges from faces, and some lookup dictionaries.
-
- Args:
- faces: list of list of int - each a closed CCW polygon of vertex indices
- Returns:
- (list of ((int, int), int), dict{ int->list of int}) -
- list elements are ((startv, endv), face index)
- dict maps vertices to edge indices
- """
-
- edges = []
- vtoe = dict()
- for findex, f in enumerate(faces):
- nf = len(f)
- for i, v in enumerate(f):
- endv = f[(i+1) % nf]
- edges.append(((v, endv), findex))
- eindex = len(edges)-1
- if v in vtoe:
- vtoe[v].append(eindex)
- else:
- vtoe[v] = [ eindex ]
- return (edges, vtoe)
+ """Find edges from faces, and some lookup dictionaries.
+
+ Args:
+ faces: list of list of int - each a closed CCW polygon of vertex indices
+ Returns:
+ (list of ((int, int), int), dict{ int->list of int}) -
+ list elements are ((startv, endv), face index)
+ dict maps vertices to edge indices
+ """
+
+ edges = []
+ vtoe = dict()
+ for findex, f in enumerate(faces):
+ nf = len(f)
+ for i, v in enumerate(f):
+ endv = f[(i + 1) % nf]
+ edges.append(((v, endv), findex))
+ eindex = len(edges) - 1
+ if v in vtoe:
+ vtoe[v].append(eindex)
+ else:
+ vtoe[v] = [eindex]
+ return (edges, vtoe)
def _GetFaceGraph(faces, edges, vtoe, points):
- """Find the face adjacency graph.
-
- Faces are adjacent if they share an edge,
- and the shared edge goes in the reverse direction,
- and if the angle between them isn't too large.
-
- Args:
- faces: list of list of int
- edges: list of ((int, int), int) - see _GetEdgeData
- vtoe: dict{ int->list of int } - see _GetEdgeData
- points: geom.Points
- Returns:
- (list of list of int, list of bool) -
- first list: each sublist is adjacent face indices for each face
- second list: maps edge index to True if it separates adjacent faces
- """
-
- face_adj = [ [] for i in range(len(faces)) ]
- is_interior_edge = [ False ] * len(edges)
- for e, ((vs, ve), f) in enumerate(edges):
- for othere in vtoe[ve]:
- ((_, we), g) = edges[othere]
- if we == vs:
- # face g is adjacent to face f
- # TODO: angle check
- if g not in face_adj[f]:
- face_adj[f].append(g)
- is_interior_edge[e] = True
- # Don't bother with mirror relations, will catch later
- return (face_adj, is_interior_edge)
+ """Find the face adjacency graph.
+
+ Faces are adjacent if they share an edge,
+ and the shared edge goes in the reverse direction,
+ and if the angle between them isn't too large.
+
+ Args:
+ faces: list of list of int
+ edges: list of ((int, int), int) - see _GetEdgeData
+ vtoe: dict{ int->list of int } - see _GetEdgeData
+ points: geom.Points
+ Returns:
+ (list of list of int, list of bool) -
+ first list: each sublist is adjacent face indices for each face
+ second list: maps edge index to True if it separates adjacent faces
+ """
+
+ face_adj = [[] for i in range(len(faces))]
+ is_interior_edge = [False] * len(edges)
+ for e, ((vs, ve), f) in enumerate(edges):
+ for othere in vtoe[ve]:
+ ((_, we), g) = edges[othere]
+ if we == vs:
+ # face g is adjacent to face f
+ # TODO: angle check
+ if g not in face_adj[f]:
+ face_adj[f].append(g)
+ is_interior_edge[e] = True
+ # Don't bother with mirror relations, will catch later
+ return (face_adj, is_interior_edge)
def _FindFaceGraphComponents(faces, face_adj):
- """Partition faces into connected components.
-
- Args:
- faces: list of list of int
- face_adj: list of list of int - see _GetFaceGraph
- Returns:
- (list of list of int, list of int) -
- first list partitions face indices into separate lists, each a component
- second list maps face indices into their component index
- """
-
- if not faces:
- return ([], [])
- components = []
- ftoc = [ -1 ] * len(faces)
- for i in range(len(faces)):
- if ftoc[i] == -1:
- compi = len(components)
- comp = []
- _FFGCSearch(i, faces, face_adj, ftoc, compi, comp)
- components.append(comp)
- return (components, ftoc)
+ """Partition faces into connected components.
+
+ Args:
+ faces: list of list of int
+ face_adj: list of list of int - see _GetFaceGraph
+ Returns:
+ (list of list of int, list of int) -
+ first list partitions face indices into separate lists,
+ each a component
+ second list maps face indices into their component index
+ """
+
+ if not faces:
+ return ([], [])
+ components = []
+ ftoc = [-1] * len(faces)
+ for i in range(len(faces)):
+ if ftoc[i] == -1:
+ compi = len(components)
+ comp = []
+ _FFGCSearch(i, faces, face_adj, ftoc, compi, comp)
+ components.append(comp)
+ return (components, ftoc)
def _FFGCSearch(findex, faces, face_adj, ftoc, compi, comp):
- """Depth first search helper function for _FindFaceGraphComponents
-
- Searches recursively through all faces connected to findex, adding
- each face found to comp and setting ftoc for that face to compi.
- """
-
- comp.append(findex)
- ftoc[findex] = compi
- for otherf in face_adj[findex]:
- if ftoc[otherf] == -1:
- _FFGCSearch(otherf, faces, face_adj, ftoc, compi, comp)
-
-def _FindOuterPoly(polys, points):
- """Assuming polys has one that contains the rest, find that one.
-
- Args:
- polys: list of list of int - list of polys given by vertex indices
- points: geom.Points
- Returns:
- int - the index in polys of the outermost one
- """
-
- if len(polys) < 2:
+ """Depth first search helper function for _FindFaceGraphComponents
+
+ Searches recursively through all faces connected to findex, adding
+ each face found to comp and setting ftoc for that face to compi.
+ """
+
+ comp.append(findex)
+ ftoc[findex] = compi
+ for otherf in face_adj[findex]:
+ if ftoc[otherf] == -1:
+ _FFGCSearch(otherf, faces, face_adj, ftoc, compi, comp)
+
+
+def _FindOuterPoly(polys, points, faces):
+ """Assuming polys has one CCW-oriented face when looking
+ down average normal of faces, return that one.
+
+ Only one of the faces should have a normal whose dot product
+ with the average normal of faces is positive.
+
+ Args:
+ polys: list of list of int - list of polys given by vertex indices
+ points: geom.Points
+ faces: list of list of int - original selected region, used to find
+ average normal
+ Returns:
+ int - the index in polys of the outermost one
+ """
+
+ if len(polys) < 2:
+ return 0
+ fnorm = (0.0, 0.0, 0.0)
+ for face in faces:
+ if len(face) > 2:
+ fnorm = geom.VecAdd(fnorm, geom.Newell(face, points))
+ if fnorm == (0.0, 0.0, 0.0):
+ return 0
+ # fnorm is really a multiple of the normal, but fine for test below
+ for i, poly in enumerate(polys):
+ if len(poly) > 2:
+ pnorm = geom.Newell(poly, points)
+ if geom.VecDot(fnorm, pnorm) > 0:
+ return i
+ print("whoops, couldn't find an outermost poly")
return 0
- for i, poly in enumerate(polys):
- otherpoly = polys[(i+1) % len(polys)]
- if geom.PointInside(points.pos[otherpoly[0]], poly, points) == 1:
- return i
- print("whoops, couldn't find an outermost poly")
- return 0
def _RotatedPolyAreaToXY(polyarea, norm):
- """Return a PolyArea rotated to xy plane.
-
- Only the points in polyarea will be transferred.
-
- Args:
- polyarea: geom.PolyArea
- norm: the normal for polyarea
- Returns:
- (geom.PolyArea, (float, ..., float), dict{ int -> int }) - new PolyArea,
- 4x3 inverse transform, dict mapping new verts to old ones
- """
-
- # find rotation matrix that takes norm to (0,0,1)
- (nx, ny, nz) = norm
- if abs(nx) < abs(ny) and abs(nx) < abs(nz):
- v = (vx, vy, vz) = geom.Norm3(0.0, nz, - ny)
- elif abs(ny) < abs(nz):
- v = (vx, vy, vz) = geom.Norm3(nz, 0.0, - nx)
- else:
- v = (vx, vy, vz) = geom.Norm3(ny, - nx, 0.0)
- (ux, uy, uz) = geom.Cross3(v, norm)
- rotmat = [ux, vx, nx, uy, vy, ny, uz, vz, nz, 0.0, 0.0, 0.0]
- # rotation matrices are orthogonal, so inverse is transpose
- invrotmat = [ux, uy, uz, vx, vy, vz, nx, ny, nz, 0.0, 0.0, 0.0]
- pointmap = dict()
- invpointmap = dict()
- newpoints = geom.Points()
- for poly in [polyarea.poly] + polyarea.holes:
- for v in poly:
- vcoords = polyarea.points.pos[v]
- newvcoords = geom.MulPoint3(vcoords, rotmat)
- newv = newpoints.AddPoint(newvcoords)
- pointmap[v] = newv
- invpointmap[newv] = v
- pa = geom.PolyArea(newpoints)
- pa.poly = [ pointmap[v] for v in polyarea.poly ]
- pa.holes = [ [ pointmap[v] for v in hole ] for hole in polyarea.holes ]
- return (pa, invrotmat, invpointmap)
-
-
-def _AddTransformedPolysToModel(mdl, polys, points, transform, pointmap):
- """Add (transformed) the points and faces to a model.
-
- Add polys to mdl. The polys have coordinates given by indices into points.pos;
- those need to be transformed by multiplying by the transform matrix.
- The vertices may already exist in mdl. Rather than relying on AddPoint to detect
- the duplicate (transform rounding error makes that dicey), the pointmap dictionary
- is used to map vertex indices in polys into those in mdl - if they exist already.
-
- Args:
- mdl: geom.Model - where to put new vertices, faces
- polys: list of list of int - each sublist a poly
- points: geom.Points - coords for vertices in polys
- transform: (float, ..., float) - 12-tuple, a 4x3 transform matrix
- pointmap: dict { int -> int } - maps new vertex indices to old ones
- Side Effects:
- The model gets new faces and vertices, based on those in polys.
- We are allowed to modify pointmap, as it will be discarded after call.
- """
-
- for i, coords in enumerate(points.pos):
- if i not in pointmap:
- p = geom.MulPoint3(coords, transform)
- pointmap[i] = mdl.points.AddPoint(p)
- for poly in polys:
- mpoly = [ pointmap[v] for v in poly ]
- mdl.faces.append(mpoly)
+ """Return a PolyArea rotated to xy plane.
+
+ Only the points in polyarea will be transferred.
+
+ Args:
+ polyarea: geom.PolyArea
+ norm: the normal for polyarea
+ Returns:
+ (geom.PolyArea, (float, ..., float), dict{ int -> int }) - new PolyArea,
+ 4x3 inverse transform, dict mapping new verts to old ones
+ """
+
+ # find rotation matrix that takes norm to (0,0,1)
+ (nx, ny, nz) = norm
+ if abs(nx) < abs(ny) and abs(nx) < abs(nz):
+ v = (vx, vy, vz) = geom.Norm3(0.0, nz, - ny)
+ elif abs(ny) < abs(nz):
+ v = (vx, vy, vz) = geom.Norm3(nz, 0.0, - nx)
+ else:
+ v = (vx, vy, vz) = geom.Norm3(ny, - nx, 0.0)
+ (ux, uy, uz) = geom.Cross3(v, norm)
+ rotmat = [ux, vx, nx, uy, vy, ny, uz, vz, nz, 0.0, 0.0, 0.0]
+ # rotation matrices are orthogonal, so inverse is transpose
+ invrotmat = [ux, uy, uz, vx, vy, vz, nx, ny, nz, 0.0, 0.0, 0.0]
+ pointmap = dict()
+ invpointmap = dict()
+ newpoints = geom.Points()
+ for poly in [polyarea.poly] + polyarea.holes:
+ for v in poly:
+ vcoords = polyarea.points.pos[v]
+ newvcoords = geom.MulPoint3(vcoords, rotmat)
+ newv = newpoints.AddPoint(newvcoords)
+ pointmap[v] = newv
+ invpointmap[newv] = v
+ pa = geom.PolyArea(newpoints)
+ pa.poly = [pointmap[v] for v in polyarea.poly]
+ pa.holes = [[pointmap[v] for v in hole] for hole in polyarea.holes]
+ pa.data = polyarea.data
+ return (pa, invrotmat, invpointmap)
+
+
+def _AddTransformedPolysToModel(mdl, polys, points, poly_data,
+ transform, pointmap):
+ """Add (transformed) the points and faces to a model.
+
+ Add polys to mdl. The polys have coordinates given by indices
+ into points.pos; those need to be transformed by multiplying by
+ the transform matrix.
+ The vertices may already exist in mdl. Rather than relying on
+ AddPoint to detect the duplicate (transform rounding error makes
+ that dicey), the pointmap dictionar is used to map vertex indices
+ in polys into those in mdl - if they exist already.
+
+ Args:
+ mdl: geom.Model - where to put new vertices, faces
+ polys: list of list of int - each sublist a poly
+ points: geom.Points - coords for vertices in polys
+ poly_data: list of any - parallel to polys
+ transform: (float, ..., float) - 12-tuple, a 4x3 transform matrix
+ pointmap: dict { int -> int } - maps new vertex indices to old ones
+ Side Effects:
+ The model gets new faces and vertices, based on those in polys.
+ We are allowed to modify pointmap, as it will be discarded after call.
+ """
+
+ for i, coords in enumerate(points.pos):
+ if i not in pointmap:
+ p = geom.MulPoint3(coords, transform)
+ pointmap[i] = mdl.points.AddPoint(p)
+ for i, poly in enumerate(polys):
+ mpoly = [pointmap[v] for v in poly]
+ mdl.faces.append(mpoly)
+ mdl.face_data.append(poly_data[i])
diff --git a/mesh_inset/offset.py b/mesh_inset/offset.py
index 7c1994696dde39ed0b4f80d5c8fee818dfc550d2..3a44b95b7aa1b7711ee5c72e9f52bdf0b20f2bc7 100644
--- a/mesh_inset/offset.py
+++ b/mesh_inset/offset.py
@@ -16,6 +16,8 @@
#
# ##### END GPL LICENSE BLOCK #####
+# <pep8 compliant>
+
"""Creating offset polygons inside faces."""
__author__ = "howard.trickey@gmail.com"
@@ -29,693 +31,725 @@ from .geom import Points
AREATOL = 1e-4
-class Spoke(object):
- """A Spoke is a line growing from an outer vertex to an inner one.
-
- A Spoke is contained in an Offset (see below).
-
- Attributes:
- origin: int - index of origin point in a Points
- dest: int - index of dest point
- is_reflex: bool - True if spoke grows from a reflex angle
- dir: (float, float, float) - direction vector (normalized)
- speed: float - at time t, other end of spoke is
- origin + t*dir. Speed is such that the wavefront
- from the face edges moves at speed 1.
- face: int - index of face containing this Spoke, in Offset
- index: int - index of this Spoke in its face
- destindex: int - index of Spoke dest in its face
- """
-
- def __init__(self, v, prev, next, face, index, points):
- """Set attribute of spoke from points making up initial angle.
-
- The spoke grows from an angle inside a face along the bisector
- of that angle. Its speed is 1/sin(.5a), where a is the angle
- formed by (prev, v, next). That speed means that the perpendicular
- from the end of the spoke to either of the prev->v or v->prev
- edges will grow at speed 1.
-
- Args:
- v: int - index of point spoke grows from
- prev: int - index of point before v on boundary (in CCW order)
- next: int - index of point after v on boundary (in CCW order)
- face: int - index of face containing this spoke, in containing offset
- index: int - index of this spoke in its face
- points: geom.Points - maps vertex indices to 3d coords
- """
- self.origin = v
- self.dest = v
- self.face = face
- self.index = index
- self.destindex = -1
- vmap = points.pos
- vp = vmap[v]
- prevp = vmap[prev]
- nextp = vmap[next]
- uin = Normalized2(Sub2(vp, prevp))
- uout = Normalized2(Sub2(nextp, vp))
- uavg = Normalized2((0.5*(uin[0]+uout[0]), 0.5*(uin[1]+uout[1])))
- if abs(Length2(uavg)) < TOL:
- # in and out vectors are reverse of each other
- self.dir = (uout[0], uout[1], 0.0)
- self.is_reflex = False
- self.speed = 1e7
- else:
- # bisector direction is 90 degree CCW rotation of average incoming/outgoing
- self.dir = (-uavg[1], uavg[0], 0.0)
- self.is_reflex = Ccw(next, v, prev, points)
- ang = Angle(prev, v, next, points) # in range [0, 180)
- sin_half_ang = math.sin(math.pi*ang / 360.0)
- if abs(sin_half_ang) < TOL:
- self.speed = 1e7
- else:
- self.speed = 1.0 / sin_half_ang
-
- def __repr__(self):
- """Printing representation of a Spoke."""
-
- return "@%d+%gt%s <%d,%d>" % (self.origin, \
- self.speed, str(self.dir), \
- self.face, self.index)
-
- def EndPoint(self, t, points, vspeed):
- """Return the coordinates of the non-origin point at time t.
-
- Args:
- t: float - time to end of spoke
- points: geom.Points - coordinate map
- vspeed: float - speed in z direction
- Returns:
- (float, float, float) - coords of spoke's endpoint at time t
+class Spoke(object):
+ """A Spoke is a line growing from an outer vertex to an inner one.
+
+ A Spoke is contained in an Offset (see below).
+
+ Attributes:
+ origin: int - index of origin point in a Points
+ dest: int - index of dest point
+ is_reflex: bool - True if spoke grows from a reflex angle
+ dir: (float, float, float) - direction vector (normalized)
+ speed: float - at time t, other end of spoke is
+ origin + t*dir. Speed is such that the wavefront
+ from the face edges moves at speed 1.
+ face: int - index of face containing this Spoke, in Offset
+ index: int - index of this Spoke in its face
+ destindex: int - index of Spoke dest in its face
"""
- p = points.pos[self.origin]
- d = self.dir
- v = self.speed
- return (p[0] + v*t*d[0], p[1] + v*t*d[1], p[2] + vspeed*t)
-
+ def __init__(self, v, prev, next, face, index, points):
+ """Set attribute of spoke from points making up initial angle.
+
+ The spoke grows from an angle inside a face along the bisector
+ of that angle. Its speed is 1/sin(.5a), where a is the angle
+ formed by (prev, v, next). That speed means that the perpendicular
+ from the end of the spoke to either of the prev->v or v->prev
+ edges will grow at speed 1.
+
+ Args:
+ v: int - index of point spoke grows from
+ prev: int - index of point before v on boundary (in CCW order)
+ next: int - index of point after v on boundary (in CCW order)
+ face: int - index of face containing this spoke, in containing offset
+ index: int - index of this spoke in its face
+ points: geom.Points - maps vertex indices to 3d coords
+ """
+
+ self.origin = v
+ self.dest = v
+ self.face = face
+ self.index = index
+ self.destindex = -1
+ vmap = points.pos
+ vp = vmap[v]
+ prevp = vmap[prev]
+ nextp = vmap[next]
+ uin = Normalized2(Sub2(vp, prevp))
+ uout = Normalized2(Sub2(nextp, vp))
+ uavg = Normalized2((0.5 * (uin[0] + uout[0]), \
+ 0.5 * (uin[1] + uout[1])))
+ if abs(Length2(uavg)) < TOL:
+ # in and out vectors are reverse of each other
+ self.dir = (uout[0], uout[1], 0.0)
+ self.is_reflex = False
+ self.speed = 1e7
+ else:
+ # bisector direction is 90 degree CCW rotation of
+ # average incoming/outgoing
+ self.dir = (-uavg[1], uavg[0], 0.0)
+ self.is_reflex = Ccw(next, v, prev, points)
+ ang = Angle(prev, v, next, points) # in range [0, 180)
+ sin_half_ang = math.sin(math.pi * ang / 360.0)
+ if abs(sin_half_ang) < TOL:
+ self.speed = 1e7
+ else:
+ self.speed = 1.0 / sin_half_ang
+
+ def __repr__(self):
+ """Printing representation of a Spoke."""
+
+ return "@%d+%gt%s <%d,%d>" % (self.origin, \
+ self.speed, str(self.dir), \
+ self.face, self.index)
+
+ def EndPoint(self, t, points, vspeed):
+ """Return the coordinates of the non-origin point at time t.
+
+ Args:
+ t: float - time to end of spoke
+ points: geom.Points - coordinate map
+ vspeed: float - speed in z direction
+ Returns:
+ (float, float, float) - coords of spoke's endpoint at time t
+ """
+
+ p = points.pos[self.origin]
+ d = self.dir
+ v = self.speed
+ return (p[0] + v * t * d[0], p[1] + v * t * d[1], p[2] + vspeed * t)
+
+ def VertexEvent(self, other, points):
+ """Intersect self with other spoke, and return the OffsetEvent, if any.
+
+ A vertex event is with one advancing spoke intersects an adjacent
+ adavancing spoke, forming a new vertex.
+
+ Args:
+ other: Spoke - other spoke to intersect with
+ points: Geom.points
+ Returns:
+ None or OffsetEvent - if there's an intersection in the growing
+ directions of the spokes, will return the OffsetEvent for
+ the intersection;
+ if lines are collinear or parallel, return None
+ """
+
+ vmap = points.pos
+ a = vmap[self.origin]
+ b = Add2(a, self.dir)
+ c = vmap[other.origin]
+ d = Add2(c, other.dir)
+ # find intersection of line ab with line cd
+ u = Sub2(b, a)
+ v = Sub2(d, c)
+ w = Sub2(a, c)
+ pp = Perp2(u, v)
+ if abs(pp) > TOL:
+ # lines or neither parallel nor collinear
+ si = Perp2(v, w) / pp
+ ti = Perp2(u, w) / pp
+ if si >= 0 and ti >= 0:
+ p = LinInterp2(a, b, si)
+ dist_ab = si * Length2(u)
+ dist_cd = ti * Length2(v)
+ time_ab = dist_ab / self.speed
+ time_cd = dist_cd / other.speed
+ time = max(time_ab, time_cd)
+ return OffsetEvent(True, time, p, self, other)
+ return None
- def VertexEvent(self, other, points):
- """Intersect self with other spoke, and return the OffsetEvent, if any.
+ def EdgeEvent(self, other, offset):
+ """Intersect self with advancing edge and return OffsetEvent, if any.
+
+ An edge event is when one advancing spoke intersects an advancing
+ edge. Advancing edges start out as face edges and move perpendicular
+ to them, at a rate of 1. The endpoints of the edge are the advancing
+ spokes on either end of the edge (so the edge shrinks or grows as
+ it advances). At some time, the edge may shrink to nothing and there
+ will be no EdgeEvent after that time.
+
+ We represent an advancing edge by the first spoke (in CCW order
+ of face) of the pair of defining spokes.
+
+ At time t, end of this spoke is at
+ o + d*s*t
+ where o=self.origin, d=self.dir, s= self.speed.
+ The advancing edge line has this equation:
+ oo + od*os*t + p*a
+ where oo, od, os are o, d, s for other spoke, and p is direction
+ vector parallel to advancing edge, and a is a real parameter.
+ Equating x and y of intersection point:
+
+ o.x + d.x*s*t = oo.x + od.x*os*t + p.x*w
+ o.y + d.y*s*t = oo.y + od.y*os*t + p.y*w
+
+ which can be rearranged into the form
+
+ a = bt + cw
+ d = et + fw
+
+ and solved for t, w.
+
+ Args:
+ other: Spoke - the edge out of this spoke's origin is the advancing
+ edge to be checked for intersection
+ offset: Offset - the containing Offset
+ Returns:
+ None or OffsetEvent - with data about the intersection, if any
+ """
+
+ vmap = offset.polyarea.points.pos
+ o = vmap[self.origin]
+ oo = vmap[other.origin]
+ otherface = offset.facespokes[other.face]
+ othernext = otherface[(other.index + 1) % len(otherface)]
+ oonext = vmap[othernext.origin]
+ p = Normalized2(Sub2(oonext, oo))
+ a = o[0] - oo[0]
+ d = o[1] - oo[1]
+ b = other.dir[0] * other.speed - self.dir[0] * self.speed
+ e = other.dir[1] * other.speed - self.dir[1] * self.speed
+ c = p[0]
+ f = p[1]
+ if abs(c) > TOL:
+ dem = e - f * b / c
+ if abs(dem) > TOL:
+ t = (d - f * a / c) / dem
+ w = (a - b * t) / c
+ else:
+ return None
+ elif abs(f) > TOL:
+ dem = b - c * e / f
+ if abs(dem) > TOL:
+ t = (a - c * d / f) / dem
+ w = (d - e * t) / f
+ else:
+ return None
+ else:
+ return None
+ if t < 0.0:
+ # intersection is in backward direction along self spoke
+ return None
+ if w < 0.0:
+ # intersection on wrong side of first end of advancing line segment
+ return None
+ # calculate the equivalent of w for the other end
+ aa = o[0] - oonext[0]
+ dd = o[1] - oonext[1]
+ bb = othernext.dir[0] * othernext.speed - self.dir[0] * self.speed
+ ee = othernext.dir[1] * othernext.speed - self.dir[1] * self.speed
+ cc = -p[0]
+ ff = -p[1]
+ if abs(cc) > TOL:
+ ww = (aa - bb * t) / cc
+ elif abs(ff) > TOL:
+ ww = (dd - ee * t) / ff
+ else:
+ return None
+ if ww < 0.0:
+ return None
+ evertex = (o[0] + self.dir[0] * self.speed * t, \
+ o[1] + self.dir[1] * self.speed * t)
+ return OffsetEvent(False, t, evertex, self, other)
- A vertex event is with one advancing spoke intersects an adjacent
- adavancing spoke, forming a new vertex.
- Args:
- other: Spoke - other spoke to intersect with
- points: Geom.points
- Returns:
- None or OffsetEvent - if there's an intersection in the growing
- directions of the spokes, will return the OffsetEvent for
- the intersection;
- if lines are collinear or parallel, return None
- """
-
- vmap = points.pos
- a = vmap[self.origin]
- b = Add2(a, self.dir)
- c = vmap[other.origin]
- d = Add2(c, other.dir)
- # find intersection of line ab with line cd
- u = Sub2(b, a)
- v = Sub2(d, c)
- w = Sub2(a, c)
- pp = Perp2(u, v)
- if abs(pp) > TOL:
- # lines or neither parallel nor collinear
- si = Perp2(v, w) / pp
- ti = Perp2(u, w) / pp
- if si >= 0 and ti >= 0:
- p = LinInterp2(a, b, si)
- dist_ab = si*Length2(u)
- dist_cd = ti*Length2(v)
- time_ab = dist_ab / self.speed
- time_cd = dist_cd / other.speed
- time = max(time_ab, time_cd)
- return OffsetEvent(True, time, p, self, other)
- return None
-
- def EdgeEvent(self, other, offset):
- """Intersect self with advancing edge and return OffsetEvent, if any.
-
- An edge event is when one advancing spoke intersects an advancing
- edge. Advancing edges start out as face edges and move perpendicular
- to them, at a rate of 1. The endpoints of the edge are the advancing
- spokes on either end of the edge (so the edge shrinks or grows as
- it advances). At some time, the edge may shrink to nothing and there
- will be no EdgeEvent after that time.
-
- We represent an advancing edge by the first spoke (in CCW order
- of face) of the pair of defining spokes.
-
- At time t, end of this spoke is at
- o + d*s*t
- where o=self.origin, d=self.dir, s= self.speed.
- The advancing edge line has this equation:
- oo + od*os*t + p*a
- where oo, od, os are o, d, s for other spoke, and p is direction
- vector parallel to advancing edge, and a is a real parameter.
- Equating x and y of intersection point:
-
- o.x + d.x*s*t = oo.x + od.x*os*t + p.x*w
- o.y + d.y*s*t = oo.y + od.y*os*t + p.y*w
-
- which can be rearranged into the form
-
- a = bt + cw
- d = et + fw
-
- and solved for t, w.
-
- Args:
- other: Spoke - the edge out of this spoke's origin is the advancing
- edge to be checked for intersection
- offset: Offset - the containing Offset
- Returns:
- None or OffsetEvent - with data about the intersection, if any
+class OffsetEvent(object):
+ """An event involving a spoke during offset computation.
+
+ The events kinds are:
+ vertex event: the spoke intersects an adjacent spoke and makes a new
+ vertex
+ edge event: the spoke hits an advancing edge and splits it
+
+ Attributes:
+ is_vertex_event: True if this is a vertex event (else it is edge event)
+ time: float - time at which it happens (edges advance at speed 1)
+ event_vertex: (float, float) - intersection point of event
+ spoke: Spoke - the spoke that this event is for
+ other: Spoke - other spoke involved in event; if vertex event, this will
+ be an adjacent spoke that intersects; if an edge event, this is the
+ spoke whose origin's outgoing edge grows to hit this event's spoke
"""
- vmap = offset.polyarea.points.pos
- o = vmap[self.origin]
- oo = vmap[other.origin]
- otherface = offset.facespokes[other.face]
- othernext = otherface[(other.index+1) % len(otherface)]
- oonext = vmap[othernext.origin]
- p = Normalized2(Sub2(oonext, oo))
- a = o[0] - oo[0]
- d = o[1] - oo[1]
- b = other.dir[0]*other.speed - self.dir[0]*self.speed
- e = other.dir[1]*other.speed - self.dir[1]*self.speed
- c = p[0]
- f = p[1]
- if abs(c) > TOL:
- dem = e - f*b/c
- if abs(dem) > TOL:
- t = (d - f*a/c) / dem
- w = (a - b*t) / c
- else:
- return None
- elif abs(f) > TOL:
- dem = b - c*e/f
- if abs(dem) > TOL:
- t = (a - c*d/f) / dem
- w = (d - e*t) / f
- else:
- return None
- else:
- return None
- if t < 0.0:
- # intersection is in backward direction along self spoke
- return None
- if w < 0.0:
- # intersection is on wrong side of first end of advancing line segment
- return None
- # calculate the equivalent of w for the other end
- aa = o[0] - oonext[0]
- dd = o[1] - oonext[1]
- bb = othernext.dir[0]*othernext.speed - self.dir[0]*self.speed
- ee = othernext.dir[1]*othernext.speed - self.dir[1]*self.speed
- cc = -p[0]
- ff = -p[1]
- if abs(cc) > TOL:
- ww = (aa - bb*t) / cc
- elif abs(ff) > TOL:
- ww = (dd - ee*t) / ff
- else:
- return None
- if ww < 0.0:
- return None
- evertex = (o[0] + self.dir[0]*self.speed*t, \
- o[1] + self.dir[1]*self.speed*t)
- return OffsetEvent(False, t, evertex, self, other)
-
+ def __init__(self, isv, time, evertex, spoke, other):
+ """Creates and initializes attributes of an OffsetEvent."""
-class OffsetEvent(object):
- """An event involving a spoke during offset computation.
-
- The events kinds are:
- vertex event: the spoke intersects an adjacent spoke and makes a new vertex
- edge event: the spoke hits an advancing edge and splits it
-
- Attributes:
- is_vertex_event: True if this is a vertex event (else it is edge event)
- time: float - time at which it happens (edges advance at speed 1)
- event_vertex: (float, float) - intersection point of event
- spoke: Spoke - the spoke that this event is for
- other: Spoke - other spoke involved in event; if vertex event, this will
- be an adjacent spoke that intersects; if an edge event, this is the
- spoke whose origin's outgoing edge grows to hit this event's spoke
- """
-
- def __init__(self, isv, time, evertex, spoke, other):
- """Creates and initializes attributes of an OffsetEvent."""
-
- self.is_vertex_event = isv
- self.time = time
- self.event_vertex = evertex
- self.spoke = spoke
- self.other = other
-
- def __repr__(self):
- """Printing representation of an event."""
-
- if self.is_vertex_event:
- c = "V"
- else:
- c = "E"
- return "%s t=%5f %s %s %s" % (c, self.time, str(self.event_vertex), \
- repr(self.spoke), repr(self.other))
+ self.is_vertex_event = isv
+ self.time = time
+ self.event_vertex = evertex
+ self.spoke = spoke
+ self.other = other
+ def __repr__(self):
+ """Printing representation of an event."""
-class Offset(object):
- """Represents an offset polygonal area, and used to construct one.
-
- Currently, the polygonal area must lie approximately in the XY plane.
- As well as growing inwards in that plane, the advancing lines also
- move in the Z direction at the rate of vspeed.
-
- Attributes:
- polyarea: geom.PolyArea - the area we are offsetting from.
- We share the polyarea.points, and add to it as points in
- the offset polygonal area are computed.
- facespokes: list of list of Spoke - each sublist is a closed face
- (oriented CCW); the faces may mutually interfere.
- These lists are spokes for polyarea.poly + polyarea.holes.
- endtime: float - time when this offset hits its first
- event (relative to beginning of this offset), or the amount
- that takes this offset to the end of the total Build time
- timesofar: float - sum of times taken by all containing Offsets
- vspeed: float - speed that edges move perpendicular to offset plane
- inneroffsets: list of Offset - the offsets that take over after this (inside it)
- """
-
- def __init__(self, polyarea, time, vspeed):
- """Set up initial state of Offset from a polyarea.
-
- Args:
- polyarea: geom.PolyArea
- time: float - time so far
- """
+ if self.is_vertex_event:
+ c = "V"
+ else:
+ c = "E"
+ return "%s t=%5f %s %s %s" % (c, self.time, str(self.event_vertex), \
+ repr(self.spoke), repr(self.other))
- self.polyarea = polyarea
- self.facespokes = []
- self.endtime = 0.0
- self.timesofar = time
- self.vspeed = vspeed
- self.inneroffsets = []
- self.InitFaceSpokes(polyarea.poly)
- for f in polyarea.holes:
- self.InitFaceSpokes(f)
-
- def __repr__(self):
- ans = ["Offset: endtime=%g" % self.endtime]
- for i, face in enumerate(self.facespokes):
- ans.append(("<%d>" % i) + str([ str(spoke) for spoke in face ]))
- return '\n'.join(ans)
-
- def PrintNest(self, indent_level=0):
- indent = " " * indent_level * 4
- print(indent + "Offset timesofar=", self.timesofar, "endtime=", self.endtime)
- print(indent + " polyarea=", self.polyarea.poly, self.polyarea.holes)
- for o in self.inneroffsets:
- o.PrintNest(indent_level+1)
-
- def InitFaceSpokes(self, face_vertices):
- """Initialize the offset representation of a face from vertex list.
-
- If the face has no area or too small an area, don't bother making it.
-
- Args:
- face_vertices: list of int - point indices for boundary of face
- Side effect:
- A new face (list of spokes) may be added to self.facespokes
- """
- n = len(face_vertices)
- if n <= 2:
- return
- points = self.polyarea.points
- area = abs(geom.SignedArea(face_vertices, points))
- if area < AREATOL:
- return
- findex = len(self.facespokes)
- fspokes = [ Spoke(v, face_vertices[(i-1) % n], \
- face_vertices[(i+1) % n], findex, i, points) \
- for i, v in enumerate(face_vertices) ]
- self.facespokes.append(fspokes)
-
- def NextSpokeEvents(self, spoke):
- """Return the OffsetEvents that will next happen for a given spoke.
-
- It might happen that some events happen essentially simultaneously,
- and also it is convenient to separate Edge and Vertex events, so
- we return two lists.
- But, for vertex events, only look at the event with the next Spoke,
- as the event with the previous spoke will be accounted for when we
- consider that previous spoke.
-
- Args:
- spoke: Spoke - a spoke in one of the faces of this object
- Returns:
- (float, list of OffsetEvent, list of OffsetEvent) - time of next event,
- next Vertex event list and next Edge event list
- """
-
- facespokes = self.facespokes[spoke.face]
- n = len(facespokes)
- bestt = 1e100
- bestv = []
- beste = []
- # First find vertex event (only the one with next spoke)
- next_spoke = facespokes[(spoke.index+1) % n]
- ev = spoke.VertexEvent(next_spoke, self.polyarea.points)
- if ev:
- bestv = [ev]
- bestt = ev.time
- # Now find edge events, if this is a reflex vertex
- if spoke.is_reflex:
- prev_spoke = facespokes[(spoke.index-1) % n]
- for f in self.facespokes:
- for other in f:
- if other == spoke or other == prev_spoke:
- continue
- ev = spoke.EdgeEvent(other, self)
- if ev:
- if ev.time < bestt - TOL:
- beste = []
- bestv = []
- bestt = ev.time
- if abs(ev.time - bestt) < TOL:
- beste.append(ev)
- return (bestt, bestv, beste)
-
- def Build(self, target = 2e100):
- """Build the complete Offset structure or up until target time.
-
- Find the next event(s), makes the appropriate inner Offsets
- that are inside this one, and calls Build on those Offsets to continue the
- process until only a single point is left or time reaches target.
+class Offset(object):
+ """Represents an offset polygonal area, and used to construct one.
+
+ Currently, the polygonal area must lie approximately in the XY plane.
+ As well as growing inwards in that plane, the advancing lines also
+ move in the Z direction at the rate of vspeed.
+
+ Attributes:
+ polyarea: geom.PolyArea - the area we are offsetting from.
+ We share the polyarea.points, and add to it as points in
+ the offset polygonal area are computed.
+ facespokes: list of list of Spoke - each sublist is a closed face
+ (oriented CCW); the faces may mutually interfere.
+ These lists are spokes for polyarea.poly + polyarea.holes.
+ endtime: float - time when this offset hits its first
+ event (relative to beginning of this offset), or the amount
+ that takes this offset to the end of the total Build time
+ timesofar: float - sum of times taken by all containing Offsets
+ vspeed: float - speed that edges move perpendicular to offset plane
+ inneroffsets: list of Offset - the offsets that take over after this
+ (inside it)
"""
- bestt = 1e100
- bestevs = [[], []]
- for f in self.facespokes:
- for s in f:
- (t, ve, ee) = self.NextSpokeEvents(s)
- if t < bestt - TOL:
- bestevs = [[], []]
- bestt = t
- if abs(t-bestt) < TOL:
- bestevs[0].extend(ve)
- bestevs[1].extend(ee)
- if bestt == 1e100:
- # could happen if polygon is oriented wrong
- # or in other special cases
- return
- if abs(bestt) < TOL:
- # seems to be in a loop, so quit
- return
- self.endtime = bestt
- (ve, ee) = bestevs
- newfaces = []
- splitjoin = None
- if target < self.endtime:
- self.endtime = target
- newfaces = self.MakeNewFaces(self.endtime)
- elif ve and not ee:
- # Only vertex events.
- # Merging of successive vertices in inset face will
- # take care of the vertex events
- newfaces = self.MakeNewFaces(self.endtime)
- else:
- # Edge events too
- # First make the new faces (handles all vertex events)
- newfaces = self.MakeNewFaces(self.endtime)
- # Only do one edge event (handle other simultaneous edge
- # events in subsequent recursive Build calls)
- splitjoin = self.SplitJoinFaces(newfaces, ee[0])
- nexttarget = target - self.endtime
- if len(newfaces) > 0:
- pa = geom.PolyArea(points = self.polyarea.points)
- pa.color = self.polyarea.color
- newt = self.timesofar+self.endtime
- pa2 = None # may make another
- if not splitjoin:
- pa.poly = newfaces[0]
- pa.holes = newfaces[1:]
- elif splitjoin[0] == 'split':
- (_, findex, newface0, newface1) = splitjoin
- if findex == 0:
- # Outer poly of polyarea was split.
- # Now there will be two polyareas.
- # If there were holes, need to allocate according to
- # which one contains the holes.
- pa.poly = newface0
- pa2 = geom.PolyArea(points = self.polyarea.points)
- pa2.color = self.polyarea.color
- pa2.poly = newface1
- if len(newfaces) > 1:
- # print("need to allocate holes")
- for hf in newfaces[1:]:
- if pa.ContainsPoly(hf, self.polyarea.points):
- # print("add", hf, "to", pa.poly)
- pa.holes.append(hf)
- elif pa2.ContainsPoly(hf, self.polyarea.points):
- # print("add", hf, "to", pa2.poly)
- pa2.holes.append(hf)
- else:
- print("whoops, hole in neither poly!")
- self.inneroffsets = [ Offset(pa, newt, self.vspeed), Offset(pa2, newt, self.vspeed) ]
+ def __init__(self, polyarea, time, vspeed):
+ """Set up initial state of Offset from a polyarea.
+
+ Args:
+ polyarea: geom.PolyArea
+ time: float - time so far
+ """
+
+ self.polyarea = polyarea
+ self.facespokes = []
+ self.endtime = 0.0
+ self.timesofar = time
+ self.vspeed = vspeed
+ self.inneroffsets = []
+ self.InitFaceSpokes(polyarea.poly)
+ for f in polyarea.holes:
+ self.InitFaceSpokes(f)
+
+ def __repr__(self):
+ ans = ["Offset: endtime=%g" % self.endtime]
+ for i, face in enumerate(self.facespokes):
+ ans.append(("<%d>" % i) + str([str(spoke) for spoke in face]))
+ return '\n'.join(ans)
+
+ def PrintNest(self, indent_level=0):
+ indent = " " * indent_level * 4
+ print(indent + "Offset timesofar=", self.timesofar, "endtime=",
+ self.endtime)
+ print(indent + " polyarea=", self.polyarea.poly, self.polyarea.holes)
+ for o in self.inneroffsets:
+ o.PrintNest(indent_level + 1)
+
+ def InitFaceSpokes(self, face_vertices):
+ """Initialize the offset representation of a face from vertex list.
+
+ If the face has no area or too small an area, don't bother making it.
+
+ Args:
+ face_vertices: list of int - point indices for boundary of face
+ Side effect:
+ A new face (list of spokes) may be added to self.facespokes
+ """
+
+ n = len(face_vertices)
+ if n <= 2:
+ return
+ points = self.polyarea.points
+ area = abs(geom.SignedArea(face_vertices, points))
+ if area < AREATOL:
+ return
+ findex = len(self.facespokes)
+ fspokes = [Spoke(v, face_vertices[(i - 1) % n], \
+ face_vertices[(i + 1) % n], findex, i, points) \
+ for i, v in enumerate(face_vertices)]
+ self.facespokes.append(fspokes)
+
+ def NextSpokeEvents(self, spoke):
+ """Return the OffsetEvents that will next happen for a given spoke.
+
+ It might happen that some events happen essentially simultaneously,
+ and also it is convenient to separate Edge and Vertex events, so
+ we return two lists.
+ But, for vertex events, only look at the event with the next Spoke,
+ as the event with the previous spoke will be accounted for when we
+ consider that previous spoke.
+
+ Args:
+ spoke: Spoke - a spoke in one of the faces of this object
+ Returns:
+ (float, list of OffsetEvent, list of OffsetEvent) -
+ time of next event,
+ next Vertex event list and next Edge event list
+ """
+
+ facespokes = self.facespokes[spoke.face]
+ n = len(facespokes)
+ bestt = 1e100
+ bestv = []
+ beste = []
+ # First find vertex event (only the one with next spoke)
+ next_spoke = facespokes[(spoke.index + 1) % n]
+ ev = spoke.VertexEvent(next_spoke, self.polyarea.points)
+ if ev:
+ bestv = [ev]
+ bestt = ev.time
+ # Now find edge events, if this is a reflex vertex
+ if spoke.is_reflex:
+ prev_spoke = facespokes[(spoke.index - 1) % n]
+ for f in self.facespokes:
+ for other in f:
+ if other == spoke or other == prev_spoke:
+ continue
+ ev = spoke.EdgeEvent(other, self)
+ if ev:
+ if ev.time < bestt - TOL:
+ beste = []
+ bestv = []
+ bestt = ev.time
+ if abs(ev.time - bestt) < TOL:
+ beste.append(ev)
+ return (bestt, bestv, beste)
+
+ def Build(self, target=2e100):
+ """Build the complete Offset structure or up until target time.
+
+ Find the next event(s), makes the appropriate inner Offsets
+ that are inside this one, and calls Build on those Offsets to continue
+ the process until only a single point is left or time reaches target.
+ """
+
+ bestt = 1e100
+ bestevs = [[], []]
+ for f in self.facespokes:
+ for s in f:
+ (t, ve, ee) = self.NextSpokeEvents(s)
+ if t < bestt - TOL:
+ bestevs = [[], []]
+ bestt = t
+ if abs(t - bestt) < TOL:
+ bestevs[0].extend(ve)
+ bestevs[1].extend(ee)
+ if bestt == 1e100:
+ # could happen if polygon is oriented wrong
+ # or in other special cases
+ return
+ if abs(bestt) < TOL:
+ # seems to be in a loop, so quit
+ return
+ self.endtime = bestt
+ (ve, ee) = bestevs
+ newfaces = []
+ splitjoin = None
+ if target < self.endtime:
+ self.endtime = target
+ newfaces = self.MakeNewFaces(self.endtime)
+ elif ve and not ee:
+ # Only vertex events.
+ # Merging of successive vertices in inset face will
+ # take care of the vertex events
+ newfaces = self.MakeNewFaces(self.endtime)
else:
- # A hole was split. New faces just replace the split one.
- pa.poly = newfaces[0]
- pa.holes = newfaces[0:findex] + [ newface0, newface1 ] + \
- newfaces[findex+1:]
- else:
- # A join
- (_, findex, othfindex, newface0) = splitjoin
- if findex == 0 or othfindex == 0:
- # Outer poly was joined to one hole.
- pa.poly = newface0
- pa.holes = [ f for f in newfaces if f is not None ]
+ # Edge events too
+ # First make the new faces (handles all vertex events)
+ newfaces = self.MakeNewFaces(self.endtime)
+ # Only do one edge event (handle other simultaneous edge
+ # events in subsequent recursive Build calls)
+ splitjoin = self.SplitJoinFaces(newfaces, ee[0])
+ nexttarget = target - self.endtime
+ if len(newfaces) > 0:
+ pa = geom.PolyArea(points=self.polyarea.points)
+ pa.data = self.polyarea.data
+ newt = self.timesofar + self.endtime
+ pa2 = None # may make another
+ if not splitjoin:
+ pa.poly = newfaces[0]
+ pa.holes = newfaces[1:]
+ elif splitjoin[0] == 'split':
+ (_, findex, newface0, newface1) = splitjoin
+ if findex == 0:
+ # Outer poly of polyarea was split.
+ # Now there will be two polyareas.
+ # If there were holes, need to allocate according to
+ # which one contains the holes.
+ pa.poly = newface0
+ pa2 = geom.PolyArea(points=self.polyarea.points)
+ pa2.data = self.polyarea.data
+ pa2.poly = newface1
+ if len(newfaces) > 1:
+ # print("need to allocate holes")
+ for hf in newfaces[1:]:
+ if pa.ContainsPoly(hf, self.polyarea.points):
+ # print("add", hf, "to", pa.poly)
+ pa.holes.append(hf)
+ elif pa2.ContainsPoly(hf, self.polyarea.points):
+ # print("add", hf, "to", pa2.poly)
+ pa2.holes.append(hf)
+ else:
+ print("whoops, hole in neither poly!")
+ self.inneroffsets = [Offset(pa, newt, self.vspeed), \
+ Offset(pa2, newt, self.vspeed)]
+ else:
+ # A hole was split. New faces just replace the split one.
+ pa.poly = newfaces[0]
+ pa.holes = newfaces[0:findex] + [newface0, newface1] + \
+ newfaces[findex + 1:]
+ else:
+ # A join
+ (_, findex, othfindex, newface0) = splitjoin
+ if findex == 0 or othfindex == 0:
+ # Outer poly was joined to one hole.
+ pa.poly = newface0
+ pa.holes = [f for f in newfaces if f is not None]
+ else:
+ # Two holes were joined
+ pa.poly = newfaces[0]
+ pa.holes = [f for f in newfaces if f is not None] + \
+ [newface0]
+ self.inneroffsets = [Offset(pa, newt, self.vspeed)]
+ if pa2:
+ self.inneroffsets.append(Offset(pa2, newt, self.vspeed))
+ if nexttarget > TOL:
+ for o in self.inneroffsets:
+ o.Build(nexttarget)
+
+ def FaceAtSpokeEnds(self, f, t):
+ """Return a new face that is at the spoke ends of face f at time t.
+
+ Also merges any adjacent approximately equal vertices into one vertex,
+ so returned list may be smaller than len(f).
+ Also sets the destindex fields of the spokes to the vertex they
+ will now end at.
+
+ Args:
+ f: list of Spoke - one of self.faces
+ t: float - time in this offset
+ Returns:
+ list of int - indices into self.polyarea.points
+ (which has been extended with new ones)
+ """
+
+ newface = []
+ points = self.polyarea.points
+ for i in range(0, len(f)):
+ s = f[i]
+ vcoords = s.EndPoint(t, points, self.vspeed)
+ v = points.AddPoint(vcoords)
+ if newface:
+ if v == newface[-1]:
+ s.destindex = len(newface) - 1
+ elif i == len(f) - 1 and v == newface[0]:
+ s.destindex = 0
+ else:
+ newface.append(v)
+ s.destindex = len(newface) - 1
+ else:
+ newface.append(v)
+ s.destindex = 0
+ s.dest = v
+ return newface
+
+ def MakeNewFaces(self, t):
+ """For each face in this offset, make new face extending spokes
+ to time t.
+
+ Args:
+ t: double - time
+ Returns:
+ list of list of int - list of new faces
+ """
+
+ ans = []
+ for f in self.facespokes:
+ newf = self.FaceAtSpokeEnds(f, t)
+ if len(newf) > 2:
+ ans.append(newf)
+ return ans
+
+ def SplitJoinFaces(self, newfaces, ev):
+ """Use event ev to split or join faces.
+
+ Given ev, an edge event, use the ev spoke to split the
+ other spoke's inner edge.
+ If the ev spoke's face and other's face are the same, this splits the
+ face into two; if the faces are different, it joins them into one.
+ We have just made faces at the end of the spokes.
+ We have to remove the edge going from the other spoke to its
+ next spoke, and replace it with two edges, going to and from
+ the event spoke's destination.
+ General situation:
+ __ s ____
+ c\ b\ | /a /e
+ \ \|/ /
+ f----------------g
+ / d \
+ o/ \h
+
+ where sd is the event spoke and of is the "other spoke",
+ hg is a spoke, and cf, fg. ge, ad, and db are edges in
+ the new inside face.
+ What we are to do is to split fg into two edges, with the
+ joining point attached where b,s,a join.
+ There are a bunch of special cases:
+ - one of split fg edges might have zero length because end points
+ are already coincident or nearly coincident.
+ - maybe c==b or e==a
+
+ Args:
+ newfaces: list of list of int - the new faces
+ ev: OffsetEvent - an edge event
+ Side Effects:
+ faces in newfaces that are involved in split or join are
+ set to None
+ Returns: one of:
+ ('split', int, list of int, list of int) - int is the index in
+ newfaces of the face that was split, two lists are the
+ split faces
+ ('join', int, int, list of int) - two ints are the indices in
+ newfaces of the faces that were joined, and the list is
+ the joined face
+ """
+
+ # print("SplitJoinFaces", newfaces, ev)
+ spoke = ev.spoke
+ other = ev.other
+ findex = spoke.face
+ othfindex = other.face
+ newface = newfaces[findex]
+ othface = newfaces[othfindex]
+ nnf = len(newface)
+ nonf = len(othface)
+ d = spoke.destindex
+ f = other.destindex
+ c = (f - 1) % nonf
+ g = (f + 1) % nonf
+ e = (f + 2) % nonf
+ a = (d - 1) % nnf
+ b = (d + 1) % nnf
+ # print("newface=", newface)
+ # if findex != othfindex: print("othface=", othface)
+ # print("d=", d, "f=", f, "c=", c, "g=", g, "e=", e, "a=", a, "b=", b)
+ newface0 = []
+ newface1 = []
+ # The two new faces put spoke si's dest on edge between
+ # pi's dest and qi (edge after pi)'s dest in original face.
+ # These are indices in the original face; the current dest face
+ # may have fewer elements because of merging successive points
+ if findex == othfindex:
+ # Case where splitting one new face into two.
+ # The new new faces are:
+ # [d, g, e, ..., a] and [d, b, ..., c, f]
+ # (except we actually want the vertex numbers at those positions)
+ newface0 = [newface[d]]
+ i = g
+ while i != d:
+ newface0.append(newface[i])
+ i = (i + 1) % nnf
+ newface1 = [newface[d]]
+ i = b
+ while i != f:
+ newface1.append(newface[i])
+ i = (i + 1) % nnf
+ newface1.append(newface[f])
+ # print("newface0=", newface0, "newface1=", newface1)
+ # now the destindex values for the spokes are messed up
+ # but I don't think we need them again
+ newfaces[findex] = None
+ return ('split', findex, newface0, newface1)
else:
- # Two holes were joined
- pa.poly = newfaces[0]
- pa.holes = [ f for f in newfaces if f is not None ] + [ newface0 ]
- self.inneroffsets = [ Offset(pa, newt, self.vspeed) ]
- if pa2:
- self.inneroffsets.append(Offset(pa2, newt, self.vspeed))
- if nexttarget > TOL:
- for o in self.inneroffsets:
- o.Build(nexttarget)
-
- def FaceAtSpokeEnds(self, f, t):
- """Return a new face that is at the spoke ends of face f at time t.
+ # Case where joining two faces into one.
+ # The new face is splicing d's face between
+ # f and g in other face (or the reverse of all of that).
+ newface0 = [othface[i] for i in range(0, f + 1)]
+ newface0.append(newface[d])
+ i = b
+ while i != d:
+ newface0.append(newface[i])
+ i = (i + 1) % nnf
+ newface0.append(newface[d])
+ if g != 0:
+ newface0.extend([othface[i] for i in range(g, nonf)])
+ # print("newface0=", newface0)
+ newfaces[findex] = None
+ newfaces[othfindex] = None
+ return ('join', findex, othfindex, newface0)
+
+ def InnerPolyAreas(self):
+ """Return the interior of the offset (and contained offsets) as
+ PolyAreas.
+
+ Returns:
+ geom.PolyAreas
+ """
+
+ ans = geom.PolyAreas()
+ ans.points = self.polyarea.points
+ _AddInnerAreas(self, ans)
+ return ans
+
+ def MaxAmount(self):
+ """Returns the maximum offset amount possible.
+ Returns:
+ float - maximum amount
+ """
+
+ # Need to do Build on a copy of points
+ # so don't add points that won't be used when
+ # really do a Build with a smaller amount
+ test_points = geom.Points()
+ test_points.AddPoints(self.polyarea.points)
+ save_points = self.polyarea.points
+ self.polyarea.points = test_points
+ self.Build()
+ max_amount = self._MaxTime()
+ self.polyarea.points = save_points
+ return max_amount
+
+ def _MaxTime(self):
+ if self.inneroffsets:
+ return max([o._MaxTime() for o in self.inneroffsets])
+ else:
+ return self.timesofar + self.endtime
- Also merges any adjacent approximately equal vertices into one vertex,
- so returned list may be smaller than len(f).
- Also sets the destindex fields of the spokes to the vertex they
- will now end at.
- Args:
- f: list of Spoke - one of self.faces
- t: float - time in this offset
- Returns:
- list of int - indices into self.polyarea.points (which has been extended with new ones)
- """
+def _AddInnerAreas(off, polyareas):
+ """Add the innermost areas of offset off to polyareas.
- newface = []
- points = self.polyarea.points
- for i in range(0, len(f)):
- s = f[i]
- vcoords = s.EndPoint(t, points, self.vspeed)
- v = points.AddPoint(vcoords)
- if newface:
- if v == newface[-1]:
- s.destindex = len(newface) - 1
- elif i == len(f)-1 and v == newface[0]:
- s.destindex = 0
- else:
- newface.append(v)
- s.destindex = len(newface) - 1
- else:
- newface.append(v)
- s.destindex = 0
- s.dest = v
- return newface
-
- def MakeNewFaces(self, t):
- """For each face in this offset, make new face extending spokes to time t.
-
- Args:
- t: double - time
- Returns:
- list of list of int - list of new faces
- """
+ Assume that polyareas is already using the proper shared points.
- ans = []
- for f in self.facespokes:
- newf = self.FaceAtSpokeEnds(f, t)
- if len(newf) > 2:
- ans.append(newf)
- return ans
-
- def SplitJoinFaces(self, newfaces, ev):
- """Use event ev to split or join faces.
-
- Given ev, an edge event, use the ev spoke to split the
- other spoke's inner edge.
- If the ev spoke's face and other's face are the same, this splits the
- face into two; if the faces are different, it joins them into one.
- We have just made faces at the end of the spokes.
- We have to remove the edge going from the other spoke to its
- next spoke, and replace it with two edges, going to and from
- the event spoke's destination.
- General situation:
- __ s ____
- c\ b\ | /a /e
- \ \|/ /
- f----------------g
- / d \
- o/ \h
-
- where sd is the event spoke and of is the "other spoke",
- hg is a spoke, and cf, fg. ge, ad, and db are edges in
- the new inside face.
- What we are to do is to split fg into two edges, with the
- joining point attached where b,s,a join.
- There are a bunch of special cases:
- - one of split fg edges might have zero length because end points
- are already coincident or nearly coincident.
- - maybe c==b or e==a
-
- Args:
- newfaces: list of list of int - the new faces
- ev: OffsetEvent - an edge event
+ Arguments:
+ off: Offset
+ polyareas: geom.PolyAreas
Side Effects:
- faces in newfaces that are involved in split or join are
- set to None
- Returns: one of:
- ('split', int, list of int, list of int) - int is the index in
- newfaces of the face that was split, two lists are the
- split faces
- ('join', int, int, list of int) - two ints are the indices in
- newfaces of the faces that were joined, and the list is
- the joined face
+ Any non-zero-area faces in the very inside of off are
+ added to polyareas.
"""
- # print("SplitJoinFaces", newfaces, ev)
- spoke = ev.spoke
- other = ev.other
- findex = spoke.face
- othfindex = other.face
- newface = newfaces[findex]
- othface = newfaces[othfindex]
- nnf = len(newface)
- nonf = len(othface)
- d = spoke.destindex
- f = other.destindex
- c = (f-1) % nonf
- g = (f+1) % nonf
- e = (f+2) % nonf
- a = (d-1) % nnf
- b = (d+1) % nnf
- # print("newface=", newface)
- # if findex != othfindex: print("othface=", othface)
- # print("d=", d, "f=", f, "c=", c, "g=", g, "e=", e, "a=", a, "b=", b)
- newface0 = []
- newface1 = []
- # The two new faces put spoke si's dest on edge between
- # pi's dest and qi (edge after pi)'s dest in original face.
- # These are indices in the original face; the current dest face
- # may have fewer elements because of merging successive points
- if findex == othfindex:
- # Case where splitting one new face into two.
- # The new new faces are:
- # [d, g, e, ..., a] and [d, b, ..., c, f]
- # (except we actually want the vertex numbers at those positions)
- newface0 = [ newface[d] ]
- i = g
- while i != d:
- newface0.append(newface[i])
- i = (i+1) % nnf
- newface1 = [ newface[d] ]
- i = b
- while i != f:
- newface1.append(newface[i])
- i = (i+1) % nnf
- newface1.append(newface[f])
- # print("newface0=", newface0, "newface1=", newface1)
- # now the destindex values for the spokes are messed up
- # but I don't think we need them again
- newfaces[findex] = None
- return ('split', findex, newface0, newface1)
+ if off.inneroffsets:
+ for o in off.inneroffsets:
+ _AddInnerAreas(o, polyareas)
else:
- # Case where joining two faces into one.
- # The new face is splicing d's face between
- # f and g in other face (or the reverse of all of that).
- newface0 = [ othface[i] for i in range(0, f+1) ]
- newface0.append(newface[d])
- i = b
- while i != d:
- newface0.append(newface[i])
- i = (i+1) % nnf
- newface0.append(newface[d])
- if g != 0:
- newface0.extend([ othface[i] for i in range(g, nonf) ])
- # print("newface0=", newface0)
- newfaces[findex] = None
- newfaces[othfindex] = None
- return ('join', findex, othfindex, newface0)
-
-
- def InnerPolyAreas(self):
- """Return the interior of the offset (and contained offsets) as PolyAreas.
-
- Returns:
- geom.PolyAreas
- """
-
- ans = geom.PolyAreas()
- ans.points = self.polyarea.points
- _AddInnerAreas(self, ans)
- return ans
-
-
-def _AddInnerAreas(off, polyareas):
- """Add the innermost areas of offset off to polyareas.
-
- Assume that polyareas is already using the proper shared points.
-
- Arguments:
- off: Offset
- polyareas: geom.PolyAreas
- Side Effects:
- Any non-zero-area faces in the very inside of off are
- added to polyareas.
- """
-
- if off.inneroffsets:
- for o in off.inneroffsets:
- _AddInnerAreas(o, polyareas)
- else:
- newpa = geom.PolyArea(polyareas.points)
- for i, f in enumerate(off.facespokes):
- newface = off.FaceAtSpokeEnds(f, off.endtime)
- area = abs(geom.SignedArea(newface, polyareas.points))
- if area < AREATOL:
- if i == 0:
- break
- else:
- continue
- if i == 0:
- newpa.poly = newface
- newpa.color = off.polyarea.color
- else:
- newpa.holes.append(newface)
- if newpa.poly:
- polyareas.polyareas.append(newpa)
-
-
+ newpa = geom.PolyArea(polyareas.points)
+ for i, f in enumerate(off.facespokes):
+ newface = off.FaceAtSpokeEnds(f, off.endtime)
+ area = abs(geom.SignedArea(newface, polyareas.points))
+ if area < AREATOL:
+ if i == 0:
+ break
+ else:
+ continue
+ if i == 0:
+ newpa.poly = newface
+ newpa.data = off.polyarea.data
+ else:
+ newpa.holes.append(newface)
+ if newpa.poly:
+ polyareas.polyareas.append(newpa)
diff --git a/mesh_inset/triquad.py b/mesh_inset/triquad.py
index 628fbdceaf3728f733eaf6162a3c5e34062fdfd9..88affa8985e4b04f4f880da071ab81e7886b572b 100644
--- a/mesh_inset/triquad.py
+++ b/mesh_inset/triquad.py
@@ -16,6 +16,8 @@
#
# ##### END GPL LICENSE BLOCK #####
+# <pep8 compliant>
+
from . import geom
import math
@@ -28,10 +30,10 @@ from math import sqrt
# be tuples instead of lists.
# Vmaps are lists taking vertex index -> Point
-TOL = 1e-7 # a tolerance for fuzzy equality
-GTHRESH = 75 # threshold above which use greedy to _Quandrangulate
-ANGFAC = 1.0 # weighting for angles in quad goodness measure
-DEGFAC = 10.0 # weighting for degree in quad goodness measure
+TOL = 1e-7 # a tolerance for fuzzy equality
+GTHRESH = 75 # threshold above which use greedy to _Quandrangulate
+ANGFAC = 1.0 # weighting for angles in quad goodness measure
+DEGFAC = 10.0 # weighting for degree in quad goodness measure
# Angle kind constants
Ang0 = 1
@@ -42,1112 +44,1129 @@ Ang360 = 5
def TriangulateFace(face, points):
- """Triangulate the given face.
+ """Triangulate the given face.
- Uses an easy triangulation first, followed by a constrained delauney
- triangulation to get better shaped triangles.
+ Uses an easy triangulation first, followed by a constrained delauney
+ triangulation to get better shaped triangles.
- Args:
- face: list of int - indices in points, assumed CCW-oriented
- points: geom.Points - holds coordinates for vertices
- Returns:
- list of (int, int, int) - 3-tuples are CCW-oriented vertices of
- triangles making up the triangulation
- """
+ Args:
+ face: list of int - indices in points, assumed CCW-oriented
+ points: geom.Points - holds coordinates for vertices
+ Returns:
+ list of (int, int, int) - 3-tuples are CCW-oriented vertices of
+ triangles making up the triangulation
+ """
- if len(face) <= 3:
- return [ tuple(face) ]
- tris = EarChopTriFace(face, points)
- bord = _BorderEdges([face])
- triscdt = _CDT(tris, bord, points)
- return triscdt
+ if len(face) <= 3:
+ return [tuple(face)]
+ tris = EarChopTriFace(face, points)
+ bord = _BorderEdges([face])
+ triscdt = _CDT(tris, bord, points)
+ return triscdt
def TriangulateFaceWithHoles(face, holes, points):
- """Like TriangulateFace, but with holes inside the face.
-
- Works by making one complex polygon that has segments to
- and from the holes ("islands"), and then using the same method
- as TriangulateFace.
-
- Args:
- face: list of int - indices in points, assumed CCW-oriented
- holes: list of list of int - each sublist is like face
- but CW-oriented and assumed to be inside face
- points: geom.Points - holds coordinates for vertices
- Returns:
- list of (int, int, int) - 3-tuples are CCW-oriented vertices of
- triangles making up the triangulation
- """
-
- if len(holes) == 0:
- return TriangulateFace(face, points)
- allfaces = [face] + holes
- sholes = [ _SortFace(h, points) for h in holes ]
- joinedface = _JoinIslands(face, sholes, points)
- tris = EarChopTriFace(joinedface, points)
- bord = _BorderEdges(allfaces)
- triscdt = _CDT(tris, bord, points)
- return triscdt
+ """Like TriangulateFace, but with holes inside the face.
+
+ Works by making one complex polygon that has segments to
+ and from the holes ("islands"), and then using the same method
+ as TriangulateFace.
+
+ Args:
+ face: list of int - indices in points, assumed CCW-oriented
+ holes: list of list of int - each sublist is like face
+ but CW-oriented and assumed to be inside face
+ points: geom.Points - holds coordinates for vertices
+ Returns:
+ list of (int, int, int) - 3-tuples are CCW-oriented vertices of
+ triangles making up the triangulation
+ """
+
+ if len(holes) == 0:
+ return TriangulateFace(face, points)
+ allfaces = [face] + holes
+ sholes = [_SortFace(h, points) for h in holes]
+ joinedface = _JoinIslands(face, sholes, points)
+ tris = EarChopTriFace(joinedface, points)
+ bord = _BorderEdges(allfaces)
+ triscdt = _CDT(tris, bord, points)
+ return triscdt
def QuadrangulateFace(face, points):
- """Quadrangulate the face (subdivide into convex quads and tris).
+ """Quadrangulate the face (subdivide into convex quads and tris).
- Like TriangulateFace, but after triangulating, join as many pairs
- of triangles as possible into convex quadrilaterals.
+ Like TriangulateFace, but after triangulating, join as many pairs
+ of triangles as possible into convex quadrilaterals.
- Args:
- face: list of int - indices in points, assumed CCW-oriented
- points: geom.Points - holds coordinates for vertices
- Returns:
- list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
- quadrilaterals and triangles making up the quadrangulation.
- """
+ Args:
+ face: list of int - indices in points, assumed CCW-oriented
+ points: geom.Points - holds coordinates for vertices
+ Returns:
+ list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
+ quadrilaterals and triangles making up the quadrangulation.
+ """
- if len(face) <= 3:
- return [ tuple(face) ]
- tris = EarChopTriFace(face, points)
- bord = _BorderEdges([face])
- triscdt = _CDT(tris, bord, points)
- qs = _Quandrangulate(triscdt, bord, points)
- return qs
+ if len(face) <= 3:
+ return [tuple(face)]
+ tris = EarChopTriFace(face, points)
+ bord = _BorderEdges([face])
+ triscdt = _CDT(tris, bord, points)
+ qs = _Quandrangulate(triscdt, bord, points)
+ return qs
def QuadrangulateFaceWithHoles(face, holes, points):
- """Like QuadrangulateFace, but with holes inside the faces.
-
- Args:
- face: list of int - indices in points, assumed CCW-oriented
- holes: list of list of int - each sublist is like face
- but CW-oriented and assumed to be inside face
- points: geom.Points - holds coordinates for vertices
- Returns:
- list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
- quadrilaterals and triangles making up the quadrangulation.
- """
-
- if len(holes) == 0:
- return QuadrangulateFace(face, points)
- allfaces = [face] + holes
- sholes = [ _SortFace(h, points) for h in holes ]
- joinedface = _JoinIslands(face, sholes, points)
- tris = EarChopTriFace(joinedface, points)
- bord = _BorderEdges(allfaces)
- triscdt = _CDT(tris, bord, points)
- qs = _Quandrangulate(triscdt, bord, points)
- return qs
+ """Like QuadrangulateFace, but with holes inside the faces.
+
+ Args:
+ face: list of int - indices in points, assumed CCW-oriented
+ holes: list of list of int - each sublist is like face
+ but CW-oriented and assumed to be inside face
+ points: geom.Points - holds coordinates for vertices
+ Returns:
+ list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
+ quadrilaterals and triangles making up the quadrangulation.
+ """
+
+ if len(holes) == 0:
+ return QuadrangulateFace(face, points)
+ allfaces = [face] + holes
+ sholes = [_SortFace(h, points) for h in holes]
+ joinedface = _JoinIslands(face, sholes, points)
+ tris = EarChopTriFace(joinedface, points)
+ bord = _BorderEdges(allfaces)
+ triscdt = _CDT(tris, bord, points)
+ qs = _Quandrangulate(triscdt, bord, points)
+ return qs
def _SortFace(face, points):
- """Rotate face so leftmost vertex is first, where face is
- list of indices in points."""
+ """Rotate face so leftmost vertex is first, where face is
+ list of indices in points."""
- n = len(face)
- if n <= 1:
- return face
- lefti = 0
- leftv = face[0]
- for i in range(1, n):
- # following comparison is lexicographic on n-tuple
- # so sorts on x first, using lower y as tie breaker.
- if points.pos[face[i]] < points.pos[leftv]:
- lefti = i
- leftv = face[i]
- return face[lefti:] + face[0:lefti]
+ n = len(face)
+ if n <= 1:
+ return face
+ lefti = 0
+ leftv = face[0]
+ for i in range(1, n):
+ # following comparison is lexicographic on n-tuple
+ # so sorts on x first, using lower y as tie breaker.
+ if points.pos[face[i]] < points.pos[leftv]:
+ lefti = i
+ leftv = face[i]
+ return face[lefti:] + face[0:lefti]
def EarChopTriFace(face, points):
- """Triangulate given face, with coords given by indexing into points.
- Return list of faces, each of which will be a triangle.
- Use the ear-chopping method."""
-
- # start with lowest coord in 2d space to try
- # to get a pleasing uniform triangulation if starting with
- # a regular structure (like a grid)
- start = _GetLeastIndex(face, points)
- ans = []
- incr = 1
- n = len(face)
- while n > 3:
- i = _FindEar(face, n, start, incr, points)
- vm1 = face[(i - 1) % n]
- v0 = face[i]
- v1 = face[(i + 1) % n]
- face = _ChopEar(face, i)
+ """Triangulate given face, with coords given by indexing into points.
+ Return list of faces, each of which will be a triangle.
+ Use the ear-chopping method."""
+
+ # start with lowest coord in 2d space to try
+ # to get a pleasing uniform triangulation if starting with
+ # a regular structure (like a grid)
+ start = _GetLeastIndex(face, points)
+ ans = []
+ incr = 1
n = len(face)
- incr = - incr
- if incr == 1:
- start = i % n
- else:
- start = (i - 1) % n
- ans.append((vm1, v0, v1))
- ans.append(tuple(face))
- return ans
+ while n > 3:
+ i = _FindEar(face, n, start, incr, points)
+ vm1 = face[(i - 1) % n]
+ v0 = face[i]
+ v1 = face[(i + 1) % n]
+ face = _ChopEar(face, i)
+ n = len(face)
+ incr = - incr
+ if incr == 1:
+ start = i % n
+ else:
+ start = (i - 1) % n
+ ans.append((vm1, v0, v1))
+ ans.append(tuple(face))
+ return ans
def _GetLeastIndex(face, points):
- """Return index of coordinate that is leftmost, lowest in face."""
+ """Return index of coordinate that is leftmost, lowest in face."""
- bestindex = 0
- bestpos = points.pos[face[0]]
- for i in range(1, len(face)):
- pos = points.pos[face[i]]
- if pos[0] < bestpos[0] or (pos[0] == bestpos[0] and pos[1] < bestpos[1]):
- bestindex = i
- bestpos = pos
- return bestindex
+ bestindex = 0
+ bestpos = points.pos[face[0]]
+ for i in range(1, len(face)):
+ pos = points.pos[face[i]]
+ if pos[0] < bestpos[0] or \
+ (pos[0] == bestpos[0] and pos[1] < bestpos[1]):
+ bestindex = i
+ bestpos = pos
+ return bestindex
def _FindEar(face, n, start, incr, points):
- """An ear of a polygon consists of three consecutive vertices
- v(-1), v0, v1 such that v(-1) can connect to v1 without intersecting
- the polygon.
- Finds an ear, starting at index 'start' and moving
- in direction incr. (We attempt to alternate directions, to find
- 'nice' triangulations for simple convex polygons.)
- Returns index into faces of v0 (will always find one, because
- uses a desperation mode if fails to find one with above rule)."""
-
- angk = _ClassifyAngles(face, n, points)
- for mode in range(0, 5):
- i = start
- while True:
- if _IsEar(face, i, n, angk, points, mode):
- return i
- i = (i + incr) % n
- if i == start:
- break # try next higher desperation mode
+ """An ear of a polygon consists of three consecutive vertices
+ v(-1), v0, v1 such that v(-1) can connect to v1 without intersecting
+ the polygon.
+ Finds an ear, starting at index 'start' and moving
+ in direction incr. (We attempt to alternate directions, to find
+ 'nice' triangulations for simple convex polygons.)
+ Returns index into faces of v0 (will always find one, because
+ uses a desperation mode if fails to find one with above rule)."""
+
+ angk = _ClassifyAngles(face, n, points)
+ for mode in range(0, 5):
+ i = start
+ while True:
+ if _IsEar(face, i, n, angk, points, mode):
+ return i
+ i = (i + incr) % n
+ if i == start:
+ break # try next higher desperation mode
def _IsEar(face, i, n, angk, points, mode):
- """Return true, false depending on ear status of vertices
- with indices i-1, i, i+1.
- mode is amount of desperation: 0 is Normal mode,
- mode 1 allows degenerate triangles (with repeated vertices)
- mode 2 allows local self crossing (folded) ears
- mode 3 allows any convex vertex (should always be one)
- mode 4 allows anything (just to be sure loop terminates!)"""
-
- k = angk[i]
- vm2 = face[(i - 2) % n]
- vm1 = face[(i - 1) % n]
- v0 = face[i]
- v1 = face[(i + 1) % n]
- v2 = face[(i + 2) % n]
- if vm1 == v0 or v0 == v1:
- return (mode > 0)
- b = (k == Angconvex or k == Angtangential or k == Ang0)
- c = _InCone(vm1, v0, v1, v2, angk[(i + 1) % n], points) and \
- _InCone(v1, vm2, vm1, v0, angk[(i - 1) % n], points)
- if b and c:
- return _EarCheck(face, n, angk, vm1, v0, v1, points)
- if mode < 2:
- return False
- if mode == 3:
- return SegsIntersect(vm2, vm1, v0, v1, points)
- if mode == 4:
- return b
- return True
+ """Return true, false depending on ear status of vertices
+ with indices i-1, i, i+1.
+ mode is amount of desperation: 0 is Normal mode,
+ mode 1 allows degenerate triangles (with repeated vertices)
+ mode 2 allows local self crossing (folded) ears
+ mode 3 allows any convex vertex (should always be one)
+ mode 4 allows anything (just to be sure loop terminates!)"""
+
+ k = angk[i]
+ vm2 = face[(i - 2) % n]
+ vm1 = face[(i - 1) % n]
+ v0 = face[i]
+ v1 = face[(i + 1) % n]
+ v2 = face[(i + 2) % n]
+ if vm1 == v0 or v0 == v1:
+ return (mode > 0)
+ b = (k == Angconvex or k == Angtangential or k == Ang0)
+ c = _InCone(vm1, v0, v1, v2, angk[(i + 1) % n], points) and \
+ _InCone(v1, vm2, vm1, v0, angk[(i - 1) % n], points)
+ if b and c:
+ return _EarCheck(face, n, angk, vm1, v0, v1, points)
+ if mode < 2:
+ return False
+ if mode == 3:
+ return SegsIntersect(vm2, vm1, v0, v1, points)
+ if mode == 4:
+ return b
+ return True
def _EarCheck(face, n, angk, vm1, v0, v1, points):
- """Return True if the successive vertices vm1, v0, v1
- forms an ear. We already know that it is not a reflex
- Angle, and that the local cone containment is ok.
- What remains to check is that the edge vm1-v1 doesn't
- intersect any other edge of the face (besides vm1-v0
- and v0-v1). Equivalently, there can't be a reflex Angle
- inside the triangle vm1-v0-v1. (Well, there are
- messy cases when other points of the face coincide with
- v0 or touch various lines involved in the ear.)"""
- for j in range(0, n):
- fv = face[j]
- k = angk[j]
- b = (k == Angreflex or k == Ang360) \
- and not(fv == vm1 or fv == v0 or fv == v1)
- if b:
- # Is fv inside closure of triangle (vm1,v0,v1)?
- c = not(Ccw(v0, vm1, fv, points) \
- or Ccw(vm1, v1, fv, points) \
- or Ccw(v1, v0, fv, points))
- fvm1 = face[(j - 1) % n]
- fv1 = face[(j + 1) % n]
- # To try to deal with some degenerate cases,
- # also check to see if either segment attached to fv
- # intersects either segment of potential ear.
- d = SegsIntersect(fvm1, fv, vm1, v0, points) or \
- SegsIntersect(fvm1, fv, v0, v1, points) or \
- SegsIntersect(fv, fv1, vm1, v0, points) or \
- SegsIntersect(fv, fv1, v0, v1, points)
- if c or d:
- return False
- return True
+ """Return True if the successive vertices vm1, v0, v1
+ forms an ear. We already know that it is not a reflex
+ Angle, and that the local cone containment is ok.
+ What remains to check is that the edge vm1-v1 doesn't
+ intersect any other edge of the face (besides vm1-v0
+ and v0-v1). Equivalently, there can't be a reflex Angle
+ inside the triangle vm1-v0-v1. (Well, there are
+ messy cases when other points of the face coincide with
+ v0 or touch various lines involved in the ear.)"""
+ for j in range(0, n):
+ fv = face[j]
+ k = angk[j]
+ b = (k == Angreflex or k == Ang360) \
+ and not(fv == vm1 or fv == v0 or fv == v1)
+ if b:
+ # Is fv inside closure of triangle (vm1,v0,v1)?
+ c = not(Ccw(v0, vm1, fv, points) \
+ or Ccw(vm1, v1, fv, points) \
+ or Ccw(v1, v0, fv, points))
+ fvm1 = face[(j - 1) % n]
+ fv1 = face[(j + 1) % n]
+ # To try to deal with some degenerate cases,
+ # also check to see if either segment attached to fv
+ # intersects either segment of potential ear.
+ d = SegsIntersect(fvm1, fv, vm1, v0, points) or \
+ SegsIntersect(fvm1, fv, v0, v1, points) or \
+ SegsIntersect(fv, fv1, vm1, v0, points) or \
+ SegsIntersect(fv, fv1, v0, v1, points)
+ if c or d:
+ return False
+ return True
def _ChopEar(face, i):
- """Return a copy of face (of length n), omitting element i."""
+ """Return a copy of face (of length n), omitting element i."""
- return face[0:i] + face[i + 1:]
+ return face[0:i] + face[i + 1:]
def _InCone(vtest, a, b, c, bkind, points):
- """Return true if point with index vtest is in Cone of points with
- indices a, b, c, where Angle ABC has AngleKind Bkind.
- The Cone is the set of points inside the left face defined by
- segments ab and bc, disregarding all other segments of polygon for
- purposes of inside test."""
-
- if bkind == Angreflex or bkind == Ang360:
- if _InCone(vtest, c, b, a, Angconvex, points):
- return False
- return not((not(Ccw(b, a, vtest, points)) \
- and not(Ccw(b, vtest, a, points)) \
- and Ccw(b, a, vtest, points))
- or
- (not(Ccw(b, c, vtest, points)) \
- and not(Ccw(b, vtest, c, points)) \
- and Ccw(b, a, vtest, points)))
- else:
- return Ccw(a, b, vtest, points) and Ccw(b, c, vtest, points)
+ """Return true if point with index vtest is in Cone of points with
+ indices a, b, c, where Angle ABC has AngleKind Bkind.
+ The Cone is the set of points inside the left face defined by
+ segments ab and bc, disregarding all other segments of polygon for
+ purposes of inside test."""
+
+ if bkind == Angreflex or bkind == Ang360:
+ if _InCone(vtest, c, b, a, Angconvex, points):
+ return False
+ return not((not(Ccw(b, a, vtest, points)) \
+ and not(Ccw(b, vtest, a, points)) \
+ and Ccw(b, a, vtest, points))
+ or
+ (not(Ccw(b, c, vtest, points)) \
+ and not(Ccw(b, vtest, c, points)) \
+ and Ccw(b, a, vtest, points)))
+ else:
+ return Ccw(a, b, vtest, points) and Ccw(b, c, vtest, points)
def _JoinIslands(face, holes, points):
- """face is a CCW face containing the CW faces in the holes list,
- where each hole is sorted so the leftmost-lowest vertex is first.
- faces and holes are given as lists of indices into points.
- The holes should be sorted by softface.
- Add edges to make a new face that includes the holes (a Ccw traversal
- of the new face will have the inside always on the left),
- and return the new face."""
-
- while len(holes) > 0:
- (hole, holeindex) = _LeftMostFace(holes, points)
- holes = holes[0:holeindex] + holes[holeindex + 1:]
- face = _JoinIsland(face, hole, points)
- return face
+ """face is a CCW face containing the CW faces in the holes list,
+ where each hole is sorted so the leftmost-lowest vertex is first.
+ faces and holes are given as lists of indices into points.
+ The holes should be sorted by softface.
+ Add edges to make a new face that includes the holes (a Ccw traversal
+ of the new face will have the inside always on the left),
+ and return the new face."""
+
+ while len(holes) > 0:
+ (hole, holeindex) = _LeftMostFace(holes, points)
+ holes = holes[0:holeindex] + holes[holeindex + 1:]
+ face = _JoinIsland(face, hole, points)
+ return face
def _JoinIsland(face, hole, points):
- """Return a modified version of face that splices in the
- vertices of hole (which should be sorted)."""
+ """Return a modified version of face that splices in the
+ vertices of hole (which should be sorted)."""
- if len(hole) == 0:
- return face
- hv0 = hole[0]
- d = _FindDiag(face, hv0, points)
- newface = face[0:d + 1] + hole + [hv0] + face[d:]
- return newface
+ if len(hole) == 0:
+ return face
+ hv0 = hole[0]
+ d = _FindDiag(face, hv0, points)
+ newface = face[0:d + 1] + hole + [hv0] + face[d:]
+ return newface
def _LeftMostFace(holes, points):
- """Return (hole,index of hole in holes) where hole has
- the leftmost first vertex. To be able to handle empty
- holes gracefully, call an empty hole 'leftmost'.
- Assumes holes are sorted by softface."""
-
- assert(len(holes) > 0)
- lefti = 0
- lefthole = holes[0]
- if len(lefthole) == 0:
+ """Return (hole,index of hole in holes) where hole has
+ the leftmost first vertex. To be able to handle empty
+ holes gracefully, call an empty hole 'leftmost'.
+ Assumes holes are sorted by softface."""
+
+ assert(len(holes) > 0)
+ lefti = 0
+ lefthole = holes[0]
+ if len(lefthole) == 0:
+ return (lefthole, lefti)
+ leftv = lefthole[0]
+ for i in range(1, len(holes)):
+ ihole = holes[i]
+ if len(ihole) == 0:
+ return (ihole, i)
+ iv = ihole[0]
+ if points.pos[iv] < points.pos[leftv]:
+ (lefti, lefthole, leftv) = (i, ihole, iv)
return (lefthole, lefti)
- leftv = lefthole[0]
- for i in range(1, len(holes)):
- ihole = holes[i]
- if len(ihole) == 0:
- return (ihole, i)
- iv = ihole[0]
- if points.pos[iv] < points.pos[leftv]:
- (lefti, lefthole, leftv) = (i, ihole, iv)
- return (lefthole, lefti)
def _FindDiag(face, hv, points):
- """Find a vertex in face that can see vertex hv, if possible,
- and return the index into face of that vertex.
- Should be able to find a diagonal that connects a vertex of face
- left of v to hv without crossing face, but try two
- more desperation passes after that to get SOME diagonal, even if
- it might cross some edge somewhere.
- First desperation pass (mode == 1): allow points right of hv.
- Second desperation pass (mode == 2): allow crossing boundary poly"""
-
- besti = - 1
- bestdist = 1e30
- for mode in range(0, 3):
- for i in range(0, len(face)):
- v = face[i]
- if mode == 0 and points.pos[v] > points.pos[hv]:
- continue # in mode 0, only want points left of hv
- dist = _DistSq(v, hv, points)
- if dist < bestdist:
- if _IsDiag(i, v, hv, face, points) or mode == 2:
- (besti, bestdist) = (i, dist)
- if besti >= 0:
- break # found one, so don't need other modes
- assert(besti >= 0)
- return besti
+ """Find a vertex in face that can see vertex hv, if possible,
+ and return the index into face of that vertex.
+ Should be able to find a diagonal that connects a vertex of face
+ left of v to hv without crossing face, but try two
+ more desperation passes after that to get SOME diagonal, even if
+ it might cross some edge somewhere.
+ First desperation pass (mode == 1): allow points right of hv.
+ Second desperation pass (mode == 2): allow crossing boundary poly"""
+
+ besti = - 1
+ bestdist = 1e30
+ for mode in range(0, 3):
+ for i in range(0, len(face)):
+ v = face[i]
+ if mode == 0 and points.pos[v] > points.pos[hv]:
+ continue # in mode 0, only want points left of hv
+ dist = _DistSq(v, hv, points)
+ if dist < bestdist:
+ if _IsDiag(i, v, hv, face, points) or mode == 2:
+ (besti, bestdist) = (i, dist)
+ if besti >= 0:
+ break # found one, so don't need other modes
+ assert(besti >= 0)
+ return besti
def _IsDiag(i, v, hv, face, points):
- """Return True if vertex v (at index i in face) can see vertex hv.
- v and hv are indices into points.
- (v, hv) is a diagonal if hv is in the cone of the Angle at index i on face
- and no segment in face intersects (h, hv).
- """
-
- n = len(face)
- vm1 = face[(i - 1) % n]
- v1 = face[(i + 1) % n]
- k = _AngleKind(vm1, v, v1, points)
- if not _InCone(hv, vm1, v, v1, k, points):
- return False
- for j in range(0, n):
- vj = face[j]
- vj1 = face[(j + 1) % n]
- if SegsIntersect(v, hv, vj, vj1, points):
- return False
- return True
+ """Return True if vertex v (at index i in face) can see vertex hv.
+ v and hv are indices into points.
+ (v, hv) is a diagonal if hv is in the cone of the Angle at index i on face
+ and no segment in face intersects (h, hv).
+ """
+
+ n = len(face)
+ vm1 = face[(i - 1) % n]
+ v1 = face[(i + 1) % n]
+ k = _AngleKind(vm1, v, v1, points)
+ if not _InCone(hv, vm1, v, v1, k, points):
+ return False
+ for j in range(0, n):
+ vj = face[j]
+ vj1 = face[(j + 1) % n]
+ if SegsIntersect(v, hv, vj, vj1, points):
+ return False
+ return True
def _DistSq(a, b, points):
- """Return distance squared between coords with indices a and b in points."""
+ """Return distance squared between coords with indices a and b in points.
+ """
- diff = Sub2(points.pos[a], points.pos[b])
- return Dot2(diff, diff)
+ diff = Sub2(points.pos[a], points.pos[b])
+ return Dot2(diff, diff)
def _BorderEdges(facelist):
- """Return a set of (u,v) where u and v are successive vertex indices
- in some face in the list in facelist."""
+ """Return a set of (u,v) where u and v are successive vertex indices
+ in some face in the list in facelist."""
- ans = set()
- for i in range(0, len(facelist)):
- f = facelist[i]
- for j in range(1, len(f)):
- ans.add((f[j - 1], f[j]))
- ans.add((f[ - 1], f[0]))
- return ans
+ ans = set()
+ for i in range(0, len(facelist)):
+ f = facelist[i]
+ for j in range(1, len(f)):
+ ans.add((f[j - 1], f[j]))
+ ans.add((f[-1], f[0]))
+ return ans
def _CDT(tris, bord, points):
- """Tris is a list of triangles ((a,b,c), CCW-oriented indices into points)
- Bord is a set of border edges (u,v), oriented so that tris
- is a triangulation of the left face of the border(s).
- Make the triangulation "Constrained Delaunay" by flipping "reversed"
- quadrangulaterals until can flip no more.
- Return list of triangles in new triangulation."""
-
- td = _TriDict(tris)
- re = _ReveresedEdges(tris, td, bord, points)
- ts = set(tris)
- # reverse the reversed edges until done.
- # reversing and edge adds new edges, which may or
- # may not be reversed or border edges, to re for
- # consideration, but the process will stop eventually.
- while len(re) > 0:
- (a, b) = e = re.pop()
- if e in bord or not _IsReversed(e, td, points):
- continue
- # rotate e in quad adbc to get other diagonal
- erev = (b, a)
- tl = td.get(e)
- tr = td.get(erev)
- if not tl or not tr:
- continue # shouldn't happen
- c = _OtherVert(tl, a, b)
- d = _OtherVert(tr, a, b)
- if c is None or d is None:
- continue # shouldn't happen
- newt1 = (c, d, b)
- newt2 = (c, a, d)
- del td[e]
- del td[erev]
- td[(c, d)] = newt1
- td[(d, b)] = newt1
- td[(b, c)] = newt1
- td[(d, c)] = newt2
- td[(c, a)] = newt2
- td[(a, d)] = newt2
- if tl in ts:
- ts.remove(tl)
- if tr in ts:
- ts.remove(tr)
- ts.add(newt1)
- ts.add(newt2)
- re.extend([(d, b), (b, c), (c, a), (a, d)])
- return list(ts)
+ """Tris is a list of triangles ((a,b,c), CCW-oriented indices into points)
+ Bord is a set of border edges (u,v), oriented so that tris
+ is a triangulation of the left face of the border(s).
+ Make the triangulation "Constrained Delaunay" by flipping "reversed"
+ quadrangulaterals until can flip no more.
+ Return list of triangles in new triangulation."""
+
+ td = _TriDict(tris)
+ re = _ReveresedEdges(tris, td, bord, points)
+ ts = set(tris)
+ # reverse the reversed edges until done.
+ # reversing and edge adds new edges, which may or
+ # may not be reversed or border edges, to re for
+ # consideration, but the process will stop eventually.
+ while len(re) > 0:
+ (a, b) = e = re.pop()
+ if e in bord or not _IsReversed(e, td, points):
+ continue
+ # rotate e in quad adbc to get other diagonal
+ erev = (b, a)
+ tl = td.get(e)
+ tr = td.get(erev)
+ if not tl or not tr:
+ continue # shouldn't happen
+ c = _OtherVert(tl, a, b)
+ d = _OtherVert(tr, a, b)
+ if c is None or d is None:
+ continue # shouldn't happen
+ newt1 = (c, d, b)
+ newt2 = (c, a, d)
+ del td[e]
+ del td[erev]
+ td[(c, d)] = newt1
+ td[(d, b)] = newt1
+ td[(b, c)] = newt1
+ td[(d, c)] = newt2
+ td[(c, a)] = newt2
+ td[(a, d)] = newt2
+ if tl in ts:
+ ts.remove(tl)
+ if tr in ts:
+ ts.remove(tr)
+ ts.add(newt1)
+ ts.add(newt2)
+ re.extend([(d, b), (b, c), (c, a), (a, d)])
+ return list(ts)
def _TriDict(tris):
- """tris is a list of triangles (a,b,c), CCW-oriented indices.
- Return dict mapping all edges in the triangles to the containing
- triangle list."""
+ """tris is a list of triangles (a,b,c), CCW-oriented indices.
+ Return dict mapping all edges in the triangles to the containing
+ triangle list."""
- ans = dict()
- for i in range(0, len(tris)):
- (a, b, c) = t = tris[i]
- ans[(a, b)] = t
- ans[(b, c)] = t
- ans[(c, a)] = t
- return ans
+ ans = dict()
+ for i in range(0, len(tris)):
+ (a, b, c) = t = tris[i]
+ ans[(a, b)] = t
+ ans[(b, c)] = t
+ ans[(c, a)] = t
+ return ans
def _ReveresedEdges(tris, td, bord, points):
- """Return list of reversed edges in tris.
- Only want edges not in bord, and only need one representative
- of (u,v)/(v,u), so choose the one with u < v.
- td is dictionary from _TriDict, and is used to find left and right
- triangles of edges."""
-
- ans = []
- for i in range(0, len(tris)):
- (a, b, c) = tris[i]
- for e in [(a, b), (b, c), (c, a)]:
- if e in bord:
- continue
- (u, v) = e
- if u < v:
- if _IsReversed(e, td, points):
- ans.append(e)
- return ans
+ """Return list of reversed edges in tris.
+ Only want edges not in bord, and only need one representative
+ of (u,v)/(v,u), so choose the one with u < v.
+ td is dictionary from _TriDict, and is used to find left and right
+ triangles of edges."""
+
+ ans = []
+ for i in range(0, len(tris)):
+ (a, b, c) = tris[i]
+ for e in [(a, b), (b, c), (c, a)]:
+ if e in bord:
+ continue
+ (u, v) = e
+ if u < v:
+ if _IsReversed(e, td, points):
+ ans.append(e)
+ return ans
def _IsReversed(e, td, points):
- """If e=(a,b) is a non-border edge, with left-face triangle tl and
- right-face triangle tr, then it is 'reversed' if the circle through
- a, b, and (say) the other vertex of tl containts the other vertex of tr.
- td is a _TriDict, for finding triangles containing edges, and points
- gives the coordinates for vertex indices used in edges."""
-
- tl = td.get(e)
- if not tl:
- return False
- (a, b) = e
- tr = td.get((b, a))
- if not tr:
- return False
- c = _OtherVert(tl, a, b)
- d = _OtherVert(tr, a, b)
- if c is None or d is None:
- return False
- return InCircle(a, b, c, d, points)
+ """If e=(a,b) is a non-border edge, with left-face triangle tl and
+ right-face triangle tr, then it is 'reversed' if the circle through
+ a, b, and (say) the other vertex of tl containts the other vertex of tr.
+ td is a _TriDict, for finding triangles containing edges, and points
+ gives the coordinates for vertex indices used in edges."""
+
+ tl = td.get(e)
+ if not tl:
+ return False
+ (a, b) = e
+ tr = td.get((b, a))
+ if not tr:
+ return False
+ c = _OtherVert(tl, a, b)
+ d = _OtherVert(tr, a, b)
+ if c is None or d is None:
+ return False
+ return InCircle(a, b, c, d, points)
def _OtherVert(tri, a, b):
- """tri should be a tuple of 3 vertex indices, two of which are a and b.
- Return the third index, or None if all vertices are a or b"""
+ """tri should be a tuple of 3 vertex indices, two of which are a and b.
+ Return the third index, or None if all vertices are a or b"""
- for v in tri:
- if v != a and v != b:
- return v
- return None
+ for v in tri:
+ if v != a and v != b:
+ return v
+ return None
def _ClassifyAngles(face, n, points):
- """Return vector of anglekinds of the Angle around each point in face."""
+ """Return vector of anglekinds of the Angle around each point in face."""
- return [_AngleKind(face[(i - 1) % n], face[i], face[(i + 1) % n], points) for i in list(range(0, n))]
+ return [_AngleKind(face[(i - 1) % n], face[i], face[(i + 1) % n], points) \
+ for i in list(range(0, n))]
def _AngleKind(a, b, c, points):
- """Return one of the Ang... constants to classify Angle formed by ABC,
- in a counterclockwise traversal of a face,
- where a, b, c are indices into points."""
-
- if Ccw(a, b, c, points):
- return Angconvex
- elif Ccw(a, c, b, points):
- return Angreflex
- else:
- vb = points.pos[b]
- udotv = Dot2(Sub2(vb, points.pos[a]), Sub2(points.pos[c], vb))
- if udotv > 0.0:
- return Angtangential
+ """Return one of the Ang... constants to classify Angle formed by ABC,
+ in a counterclockwise traversal of a face,
+ where a, b, c are indices into points."""
+
+ if Ccw(a, b, c, points):
+ return Angconvex
+ elif Ccw(a, c, b, points):
+ return Angreflex
else:
- return Ang0 # to fix: return Ang360 if "inside" spur
+ vb = points.pos[b]
+ udotv = Dot2(Sub2(vb, points.pos[a]), Sub2(points.pos[c], vb))
+ if udotv > 0.0:
+ return Angtangential
+ else:
+ return Ang0 # to fix: return Ang360 if "inside" spur
def _Quandrangulate(tris, bord, points):
- """Tris is list of triangles, forming a triangulation of region whose
- border edges are in set bord.
- Combine adjacent triangles to make quads, trying for "good" quads where
- possible. Some triangles will probably remain uncombined"""
-
- (er, td) = _ERGraph(tris, bord, points)
- if len(er) == 0:
- return tris
- if len(er) > GTHRESH:
- match = _GreedyMatch(er)
- else:
- match = _MaxMatch(er)
- return _RemoveEdges(tris, match)
+ """Tris is list of triangles, forming a triangulation of region whose
+ border edges are in set bord.
+ Combine adjacent triangles to make quads, trying for "good" quads where
+ possible. Some triangles will probably remain uncombined"""
+
+ (er, td) = _ERGraph(tris, bord, points)
+ if len(er) == 0:
+ return tris
+ if len(er) > GTHRESH:
+ match = _GreedyMatch(er)
+ else:
+ match = _MaxMatch(er)
+ return _RemoveEdges(tris, match)
def _RemoveEdges(tris, match):
- """tris is list of triangles. er is as returned from _MaxMatch or _GreedyMatch.
- Return list of (A,D,B,C) resulting from deleting edge (A,B) causing a merge
- of two triangles; append to that list the remaining unmatched triangles."""
-
- ans = []
- triset = set(tris)
- while len(match) > 0:
- (_, e, tl, tr) = match.pop()
- (a, b) = e
- if tl in triset:
- triset.remove(tl)
- if tr in triset:
- triset.remove(tr)
- c = _OtherVert(tl, a, b)
- d = _OtherVert(tr, a, b)
- if c is None or d is None:
- continue
- ans.append((a, d, b, c))
- return ans + list(triset)
+ """tris is list of triangles.
+ er is as returned from _MaxMatch or _GreedyMatch.
+
+ Return list of (A,D,B,C) resulting from deleting edge (A,B) causing a merge
+ of two triangles; append to that list the remaining unmatched triangles."""
+
+ ans = []
+ triset = set(tris)
+ while len(match) > 0:
+ (_, e, tl, tr) = match.pop()
+ (a, b) = e
+ if tl in triset:
+ triset.remove(tl)
+ if tr in triset:
+ triset.remove(tr)
+ c = _OtherVert(tl, a, b)
+ d = _OtherVert(tr, a, b)
+ if c is None or d is None:
+ continue
+ ans.append((a, d, b, c))
+ return ans + list(triset)
def _ERGraph(tris, bord, points):
- """Make an 'Edge Removal Graph'.
-
- Given a list of triangles, the 'Edge Removal Graph' is a graph whose
- nodes are the triangles (think of a point in the center of them),
- and whose edges go between adjacent triangles (they share a non-border
- edge), such that it would be possible to remove the shared edge
- and form a convex quadrilateral. Forming a quadrilateralization
- is then a matter of finding a matching (set of edges that don't
- share a vertex - remember, these are the 'face' vertices).
- For better quadrilaterlization, we'll make the Edge Removal Graph
- edges have weights, with higher weights going to the edges that
- are more desirable to remove. Then we want a maximum weight matching
- in this graph.
-
- We'll return the graph in a kind of implicit form, using edges of
- the original triangles as a proxy for the edges between the faces
- (i.e., the edge of the triangle is the shared edge). We'll arbitrarily
- pick the triangle graph edge with lower-index start vertex.
- Also, to aid in traversing the implicit graph, we'll keep the left
- and right triangle triples with edge 'ER edge'.
- Finally, since we calculate it anyway, we'll return a dictionary
- mapping edges of the triangles to the triangle triples they're in.
-
- Args:
- tris: list of (int, int, int) giving a triple of vertex indices for
- triangles, assumed CCW oriented
- bord: set of (int, int) giving vertex indices for border edges
- points: geom.Points - for mapping vertex indices to coords
- Returns:
- (list of (weight,e,tl,tr), dict)
- where edge e=(a,b) is non-border edge
- with left face tl and right face tr (each a triple (i,j,k)),
- where removing the edge would form an "OK" quad (no concave angles),
- with weight representing the desirability of removing the edge
- The dict maps int pairs (a,b) to int triples (i,j,k), that is,
- mapping edges to their containing triangles.
- """
-
- td = _TriDict(tris)
- dd = _DegreeDict(tris)
- ans = []
- ctris = tris[:] # copy, so argument not affected
- while len(ctris) > 0:
- (i, j, k) = tl = ctris.pop()
- for e in [(i, j), (j, k), (k, i)]:
- if e in bord:
- continue
- (a, b) = e
- # just consider one of (a,b) and (b,a), to avoid dups
- if a > b:
- continue
- erev = (b, a)
- tr = td.get(erev)
- if not tr:
- continue
- c = _OtherVert(tl, a, b)
- d = _OtherVert(tr, a, b)
- if c is None or d is None:
- continue
- # calculate amax, the max of the new angles that would
- # be formed at a and b if tl and tr were combined
- amax = max(Angle(c, a, b, points) + Angle(d, a, b, points),
- Angle(c, b, a, points) + Angle(d, b, a, points))
- if amax > 180.0:
- continue
- weight = ANGFAC * (180.0 - amax) + DEGFAC * (dd[a] + dd[b])
- ans.append((weight, e, tl, tr))
- return (ans, td)
+ """Make an 'Edge Removal Graph'.
+
+ Given a list of triangles, the 'Edge Removal Graph' is a graph whose
+ nodes are the triangles (think of a point in the center of them),
+ and whose edges go between adjacent triangles (they share a non-border
+ edge), such that it would be possible to remove the shared edge
+ and form a convex quadrilateral. Forming a quadrilateralization
+ is then a matter of finding a matching (set of edges that don't
+ share a vertex - remember, these are the 'face' vertices).
+ For better quadrilaterlization, we'll make the Edge Removal Graph
+ edges have weights, with higher weights going to the edges that
+ are more desirable to remove. Then we want a maximum weight matching
+ in this graph.
+
+ We'll return the graph in a kind of implicit form, using edges of
+ the original triangles as a proxy for the edges between the faces
+ (i.e., the edge of the triangle is the shared edge). We'll arbitrarily
+ pick the triangle graph edge with lower-index start vertex.
+ Also, to aid in traversing the implicit graph, we'll keep the left
+ and right triangle triples with edge 'ER edge'.
+ Finally, since we calculate it anyway, we'll return a dictionary
+ mapping edges of the triangles to the triangle triples they're in.
+
+ Args:
+ tris: list of (int, int, int) giving a triple of vertex indices for
+ triangles, assumed CCW oriented
+ bord: set of (int, int) giving vertex indices for border edges
+ points: geom.Points - for mapping vertex indices to coords
+ Returns:
+ (list of (weight,e,tl,tr), dict)
+ where edge e=(a,b) is non-border edge
+ with left face tl and right face tr (each a triple (i,j,k)),
+ where removing the edge would form an "OK" quad (no concave angles),
+ with weight representing the desirability of removing the edge
+ The dict maps int pairs (a,b) to int triples (i,j,k), that is,
+ mapping edges to their containing triangles.
+ """
+
+ td = _TriDict(tris)
+ dd = _DegreeDict(tris)
+ ans = []
+ ctris = tris[:] # copy, so argument not affected
+ while len(ctris) > 0:
+ (i, j, k) = tl = ctris.pop()
+ for e in [(i, j), (j, k), (k, i)]:
+ if e in bord:
+ continue
+ (a, b) = e
+ # just consider one of (a,b) and (b,a), to avoid dups
+ if a > b:
+ continue
+ erev = (b, a)
+ tr = td.get(erev)
+ if not tr:
+ continue
+ c = _OtherVert(tl, a, b)
+ d = _OtherVert(tr, a, b)
+ if c is None or d is None:
+ continue
+ # calculate amax, the max of the new angles that would
+ # be formed at a and b if tl and tr were combined
+ amax = max(Angle(c, a, b, points) + Angle(d, a, b, points),
+ Angle(c, b, a, points) + Angle(d, b, a, points))
+ if amax > 180.0:
+ continue
+ weight = ANGFAC * (180.0 - amax) + DEGFAC * (dd[a] + dd[b])
+ ans.append((weight, e, tl, tr))
+ return (ans, td)
def _GreedyMatch(er):
- """er is list of (weight,e,tl,tr).
- Find maximal set so that each triangle appears in at most one member of set"""
+ """er is list of (weight,e,tl,tr).
+
+ Find maximal set so that each triangle appears in at most
+ one member of set"""
- er.sort(key = lambda v: v[0], reverse=True) # sort in order of decreasing weight
- match = set()
- ans = []
- while len(er) > 0:
- (_, _, tl, tr) = q = er.pop()
- if tl not in match and tr not in match:
- match.add(tl)
- match.add(tr)
- ans.append(q)
- return ans
+ # sort in order of decreasing weight
+ er.sort(key=lambda v: v[0], reverse=True)
+ match = set()
+ ans = []
+ while len(er) > 0:
+ (_, _, tl, tr) = q = er.pop()
+ if tl not in match and tr not in match:
+ match.add(tl)
+ match.add(tr)
+ ans.append(q)
+ return ans
def _MaxMatch(er):
- """Like _GreedyMatch, but use divide and conquer to find best possible set.
+ """Like _GreedyMatch, but use divide and conquer to find best possible set.
- Args:
- er: list of (weight,e,tl,tr) - see _ERGraph
- Returns:
- list that is a subset of er giving a maximum weight match
- """
+ Args:
+ er: list of (weight,e,tl,tr) - see _ERGraph
+ Returns:
+ list that is a subset of er giving a maximum weight match
+ """
- (ans, _) = _DCMatch(er)
- return ans
+ (ans, _) = _DCMatch(er)
+ return ans
def _DCMatch(er):
- """Recursive helper for _MaxMatch.
-
- Divide and Conquer approach to finding max weight matching.
- If we're lucky, there's an edge in er that separates the edge removal
- graph into (at least) two separate components. Then the max weight
- is either one that includes that edge or excludes it - and we can
- use a recursive call to _DCMatch to handle each component separately
- on what remains of the graph after including/excluding the separating edge.
- If we're not lucky, we fall back on _EMatch (see below).
-
- Args:
- er: list of (weight, e, tl, tr) (see _ERGraph)
- Returns:
- (list of (weight, e, tl, tr), float) - the subset forming a maximum
- matching, and the total weight of the match.
- """
-
- if not er:
- return ([], 0.0)
- if len(er) == 1:
- return (er, er[0][0])
- match = []
- matchw = 0.0
- for i in range(0, len(er)):
- (nc, comp) = _FindComponents(er, i)
- if nc == 1:
- # er[i] doesn't separate er
- continue
- (wi, _, tl, tr) = er[i]
- if comp[tl] != comp[tr]:
- # case 1: er separates graph
- # compare the matches that include er[i] versus those that exclude it
- (a, b) = _PartitionComps(er, comp, i, comp[tl], comp[tr])
- ax = _CopyExcluding(a, tl, tr)
- bx = _CopyExcluding(b, tl, tr)
- (axmatch, wax) = _DCMatch(ax)
- (bxmatch, wbx) = _DCMatch(bx)
- if len(ax) == len(a):
- wa = wax
- amatch = axmatch
- else:
- (amatch, wa) = _DCMatch(a)
- if len(bx) == len(b):
- wb = wbx
- bmatch = bxmatch
- else:
- (bmatch, wb) = _DCMatch(b)
- w = wa + wb
- wx = wax + wbx + wi
- if w > wx:
- match = amatch + bmatch
- matchw = w
- else:
- match = [er[i]] + axmatch + bxmatch
- matchw = wx
- else:
- # case 2: er not needed to separate graph
- (a, b) = _PartitionComps(er, comp, -1, 0, 0)
- (amatch, wa) = _DCMatch(a)
- (bmatch, wb) = _DCMatch(b)
- match = amatch + bmatch
- matchw = wa + wb
- if match:
- break
- if not match:
- return _EMatch(er)
- return (match, matchw)
+ """Recursive helper for _MaxMatch.
+
+ Divide and Conquer approach to finding max weight matching.
+ If we're lucky, there's an edge in er that separates the edge removal
+ graph into (at least) two separate components. Then the max weight
+ is either one that includes that edge or excludes it - and we can
+ use a recursive call to _DCMatch to handle each component separately
+ on what remains of the graph after including/excluding the separating edge.
+ If we're not lucky, we fall back on _EMatch (see below).
+
+ Args:
+ er: list of (weight, e, tl, tr) (see _ERGraph)
+ Returns:
+ (list of (weight, e, tl, tr), float) - the subset forming a maximum
+ matching, and the total weight of the match.
+ """
+
+ if not er:
+ return ([], 0.0)
+ if len(er) == 1:
+ return (er, er[0][0])
+ match = []
+ matchw = 0.0
+ for i in range(0, len(er)):
+ (nc, comp) = _FindComponents(er, i)
+ if nc == 1:
+ # er[i] doesn't separate er
+ continue
+ (wi, _, tl, tr) = er[i]
+ if comp[tl] != comp[tr]:
+ # case 1: er separates graph
+ # compare the matches that include er[i] versus
+ # those that exclude it
+ (a, b) = _PartitionComps(er, comp, i, comp[tl], comp[tr])
+ ax = _CopyExcluding(a, tl, tr)
+ bx = _CopyExcluding(b, tl, tr)
+ (axmatch, wax) = _DCMatch(ax)
+ (bxmatch, wbx) = _DCMatch(bx)
+ if len(ax) == len(a):
+ wa = wax
+ amatch = axmatch
+ else:
+ (amatch, wa) = _DCMatch(a)
+ if len(bx) == len(b):
+ wb = wbx
+ bmatch = bxmatch
+ else:
+ (bmatch, wb) = _DCMatch(b)
+ w = wa + wb
+ wx = wax + wbx + wi
+ if w > wx:
+ match = amatch + bmatch
+ matchw = w
+ else:
+ match = [er[i]] + axmatch + bxmatch
+ matchw = wx
+ else:
+ # case 2: er not needed to separate graph
+ (a, b) = _PartitionComps(er, comp, -1, 0, 0)
+ (amatch, wa) = _DCMatch(a)
+ (bmatch, wb) = _DCMatch(b)
+ match = amatch + bmatch
+ matchw = wa + wb
+ if match:
+ break
+ if not match:
+ return _EMatch(er)
+ return (match, matchw)
def _EMatch(er):
- """Exhaustive match helper for _MaxMatch.
-
- This is the case when we were unable to find a single edge
- separating the edge removal graph into two components.
- So pick a single edge and try _DCMatch on the two cases of
- including or excluding that edge. We may be lucky in these
- subcases (say, if the graph is currently a simple cycle, so
- only needs one more edge after the one we pick here to separate
- it into components). Otherwise, we'll end up back in _EMatch
- again, and the worse case will be exponential.
-
- Pick a random edge rather than say, the first, to hopefully
- avoid some pathological cases.
-
- Args:
- er: list of (weight, el, tl, tr) (see _ERGraph)
- Returns:
- (list of (weight, e, tl, tr), float) - the subset forming a maximum
- matching, and the total weight of the match.
- """
-
- if not er:
- return ([], 0.0)
- if len(er) == 1:
- return (er, er[1][1])
- i = random.randint(0, len(er)-1)
- eri = (wi, _, tl, tr) = er[i]
- # case a: include eri. exlude other edges that touch tl or tr
- a = _CopyExcluding(er, tl, tr)
- a.append(eri)
- (amatch, wa) = _DCMatch(a)
- wa += wi
- if len(a) == len(er)-1:
- # if a excludes only eri, then er didn't touch anything else
- # in the graph, and the best match will always include er
- # and we can skip the call for case b
- wb = -1.0
- bmatch = []
- else:
- b = er[:i] + er[i+1:]
- (bmatch, wb) = _DCMatch(b)
- if wa > wb:
- match = amatch
- match.append(eri)
- matchw = wa
- else:
- match = bmatch
- matchw = wb
- return (match, matchw)
+ """Exhaustive match helper for _MaxMatch.
+
+ This is the case when we were unable to find a single edge
+ separating the edge removal graph into two components.
+ So pick a single edge and try _DCMatch on the two cases of
+ including or excluding that edge. We may be lucky in these
+ subcases (say, if the graph is currently a simple cycle, so
+ only needs one more edge after the one we pick here to separate
+ it into components). Otherwise, we'll end up back in _EMatch
+ again, and the worse case will be exponential.
+
+ Pick a random edge rather than say, the first, to hopefully
+ avoid some pathological cases.
+
+ Args:
+ er: list of (weight, el, tl, tr) (see _ERGraph)
+ Returns:
+ (list of (weight, e, tl, tr), float) - the subset forming a maximum
+ matching, and the total weight of the match.
+ """
+
+ if not er:
+ return ([], 0.0)
+ if len(er) == 1:
+ return (er, er[1][1])
+ i = random.randint(0, len(er) - 1)
+ eri = (wi, _, tl, tr) = er[i]
+ # case a: include eri. exlude other edges that touch tl or tr
+ a = _CopyExcluding(er, tl, tr)
+ a.append(eri)
+ (amatch, wa) = _DCMatch(a)
+ wa += wi
+ if len(a) == len(er) - 1:
+ # if a excludes only eri, then er didn't touch anything else
+ # in the graph, and the best match will always include er
+ # and we can skip the call for case b
+ wb = -1.0
+ bmatch = []
+ else:
+ b = er[:i] + er[i + 1:]
+ (bmatch, wb) = _DCMatch(b)
+ if wa > wb:
+ match = amatch
+ match.append(eri)
+ matchw = wa
+ else:
+ match = bmatch
+ matchw = wb
+ return (match, matchw)
def _FindComponents(er, excepti):
- """Find connected components induced by edges, excluding one edge.
-
- Args:
- er: list of (weight, el, tl, tr) (see _ERGraph)
- excepti: index in er of edge to be excluded
- Returns:
- (int, dict): int is number of connected components found,
- dict maps triangle triple -> connected component index (starting at 1)
- """
-
- ncomps = 0
- comp = dict()
- for i in range(0, len(er)):
- (_, _, tl, tr) = er[i]
- for t in [tl,tr]:
- if t not in comp:
- ncomps += 1
- _FCVisit(er, excepti, comp, t, ncomps)
- return (ncomps, comp)
+ """Find connected components induced by edges, excluding one edge.
+
+ Args:
+ er: list of (weight, el, tl, tr) (see _ERGraph)
+ excepti: index in er of edge to be excluded
+ Returns:
+ (int, dict): int is number of connected components found,
+ dict maps triangle triple ->
+ connected component index (starting at 1)
+ """
+
+ ncomps = 0
+ comp = dict()
+ for i in range(0, len(er)):
+ (_, _, tl, tr) = er[i]
+ for t in [tl, tr]:
+ if t not in comp:
+ ncomps += 1
+ _FCVisit(er, excepti, comp, t, ncomps)
+ return (ncomps, comp)
def _FCVisit(er, excepti, comp, t, compnum):
- """Helper for _FindComponents depth-first-search."""
+ """Helper for _FindComponents depth-first-search."""
- comp[t] = compnum
- for i in range(0, len(er)):
- if i == excepti:
- continue
- (_, _, tl, tr) = er[i]
- if tl == t or tr == t:
- s = tl
- if s == t:
- s = tr
- if s not in comp:
- _FCVisit(er, excepti, comp, s, compnum)
+ comp[t] = compnum
+ for i in range(0, len(er)):
+ if i == excepti:
+ continue
+ (_, _, tl, tr) = er[i]
+ if tl == t or tr == t:
+ s = tl
+ if s == t:
+ s = tr
+ if s not in comp:
+ _FCVisit(er, excepti, comp, s, compnum)
def _PartitionComps(er, comp, excepti, compa, compb):
- """Partition the edges of er by component number, into two lists.
-
- Generally, put odd components into first list and even into second,
- except that component compa goes in the first and compb goes in the second,
- and we ignore edge er[excepti].
-
- Args:
- er: list of (weight, el, tl, tr) (see _ERGraph)
- comp: dict - mapping triangle triple -> connected component index
- excepti: int - index in er to ignore (unless excepti==-1)
- compa: int - component to go in first list of answer (unless 0)
- compb: int - component to go in second list of answer (unless 0)
- Returns:
- (list, list) - a partition of er according to above rules
- """
-
- parta = []
- partb = []
- for i in range(0, len(er)):
-
- if i == excepti:
- continue
- tl = er[i][2]
- c = comp[tl]
- if c == compa or (c != compb and (c&1) == 1):
- parta.append(er[i])
- else:
- partb.append(er[i])
- return (parta, partb)
+ """Partition the edges of er by component number, into two lists.
+
+ Generally, put odd components into first list and even into second,
+ except that component compa goes in the first and compb goes in the second,
+ and we ignore edge er[excepti].
+
+ Args:
+ er: list of (weight, el, tl, tr) (see _ERGraph)
+ comp: dict - mapping triangle triple -> connected component index
+ excepti: int - index in er to ignore (unless excepti==-1)
+ compa: int - component to go in first list of answer (unless 0)
+ compb: int - component to go in second list of answer (unless 0)
+ Returns:
+ (list, list) - a partition of er according to above rules
+ """
+
+ parta = []
+ partb = []
+ for i in range(0, len(er)):
+
+ if i == excepti:
+ continue
+ tl = er[i][2]
+ c = comp[tl]
+ if c == compa or (c != compb and (c & 1) == 1):
+ parta.append(er[i])
+ else:
+ partb.append(er[i])
+ return (parta, partb)
def _CopyExcluding(er, s, t):
- """Return a copy of er, excluding all those involving triangles s and t.
+ """Return a copy of er, excluding all those involving triangles s and t.
- Args:
- er: list of (weight, e, tl, tr) - see _ERGraph
- s: 3-tuple of int - a triangle
- t: 3-tuple of int - a triangle
- Returns:
- Copy of er excluding those with tl or tr == s or t
- """
+ Args:
+ er: list of (weight, e, tl, tr) - see _ERGraph
+ s: 3-tuple of int - a triangle
+ t: 3-tuple of int - a triangle
+ Returns:
+ Copy of er excluding those with tl or tr == s or t
+ """
- ans = []
- for e in er:
- (_, _, tl, tr) = e
- if tl == s or tr == s or tl == t or tr == t:
- continue
- ans.append(e)
- return ans
+ ans = []
+ for e in er:
+ (_, _, tl, tr) = e
+ if tl == s or tr == s or tl == t or tr == t:
+ continue
+ ans.append(e)
+ return ans
def _DegreeDict(tris):
- """Return a dictionary mapping vertices in tris to the number of triangles
- that they are touch."""
+ """Return a dictionary mapping vertices in tris to the number of triangles
+ that they are touch."""
- ans = dict()
- for t in tris:
- for v in t:
- if v in ans:
- ans[v] = ans[v] + 1
- else:
- ans[v] = 1
- return ans
+ ans = dict()
+ for t in tris:
+ for v in t:
+ if v in ans:
+ ans[v] = ans[v] + 1
+ else:
+ ans[v] = 1
+ return ans
def PolygonPlane(face, points):
- """Return a Normal vector for the face with 3d coords given by indexing into points."""
+ """Return a Normal vector for the face with 3d coords given by indexing
+ into points."""
- if len(face) < 3:
- return (0.0, 0.0, 1.0) # arbitrary, we really have no idea
- else:
- coords = [ points.pos[i] for i in face ]
- return Normal(coords)
+ if len(face) < 3:
+ return (0.0, 0.0, 1.0) # arbitrary, we really have no idea
+ else:
+ coords = [points.pos[i] for i in face]
+ return Normal(coords)
# This Normal appears to be on the CCW-traversing side of a polygon
def Normal(coords):
- """Return an average Normal vector for the point list, 3d coords."""
-
- if len(coords) < 3:
- return (0.0, 0.0, 1.0) # arbitrary
-
- (ax, ay, az) = coords[0]
- (bx, by, bz) = coords[1]
- (cx, cy, cz) = coords[2]
-
- if len(coords) == 3:
- sx = (ay - by) * (az + bz) + (by - cy) * (bz + cz) + (cy - ay) * (cz + az)
- sy = (az - bz) * (ax + bx) + (bz - cz) * (bx + cx) + (cz - az) * (cx + ax)
- sz = (ax - bx) * (by + by) + (bx - cx) * (by + cy) + (cx - ax) * (cy + ay)
- return Norm3(sx, sy, sz)
- else:
- sx = (ay - by) * (az + bz) + (by - cy) * (bz + cz)
- sy = (az - bz) * (ax + bx) + (bz - cz) * (bx + cx)
- sz = (ax - bx) * (ay + by) + (bx - cx) * (by + cy)
- return _NormalAux(coords[3:], coords[0], sx, sy, sz)
+ """Return an average Normal vector for the point list, 3d coords."""
+
+ if len(coords) < 3:
+ return (0.0, 0.0, 1.0) # arbitrary
+
+ (ax, ay, az) = coords[0]
+ (bx, by, bz) = coords[1]
+ (cx, cy, cz) = coords[2]
+
+ if len(coords) == 3:
+ sx = (ay - by) * (az + bz) + \
+ (by - cy) * (bz + cz) + \
+ (cy - ay) * (cz + az)
+ sy = (az - bz) * (ax + bx) + \
+ (bz - cz) * (bx + cx) + \
+ (cz - az) * (cx + ax)
+ sz = (ax - bx) * (by + by) + \
+ (bx - cx) * (by + cy) + \
+ (cx - ax) * (cy + ay)
+ return Norm3(sx, sy, sz)
+ else:
+ sx = (ay - by) * (az + bz) + (by - cy) * (bz + cz)
+ sy = (az - bz) * (ax + bx) + (bz - cz) * (bx + cx)
+ sz = (ax - bx) * (ay + by) + (bx - cx) * (by + cy)
+ return _NormalAux(coords[3:], coords[0], sx, sy, sz)
def _NormalAux(rest, first, sx, sy, sz):
- (ax, ay, az) = rest[0]
- if len(rest) == 1:
- (bx, by, bz) = first
- else:
- (bx, by, bz) = rest[1]
- nx = sx + (ay - by) * (az + bz)
- ny = sy + (az - bz) * (ax + bx)
- nz = sz + (ax - bx) * (ay + by)
- if len(rest) == 1:
- return Norm3(nx, ny, nz)
- else:
- return _NormalAux(rest[1:], first, nx, ny, nz)
+ (ax, ay, az) = rest[0]
+ if len(rest) == 1:
+ (bx, by, bz) = first
+ else:
+ (bx, by, bz) = rest[1]
+ nx = sx + (ay - by) * (az + bz)
+ ny = sy + (az - bz) * (ax + bx)
+ nz = sz + (ax - bx) * (ay + by)
+ if len(rest) == 1:
+ return Norm3(nx, ny, nz)
+ else:
+ return _NormalAux(rest[1:], first, nx, ny, nz)
def Norm3(x, y, z):
- """Return vector (x,y,z) normalized by dividing by squared length.
- Return (0.0, 0.0, 1.0) if the result is undefined."""
- sqrlen = x * x + y * y + z * z
- if sqrlen < 1e-100:
- return (0.0, 0.0, 1.0)
- else:
- try:
- d = sqrt(sqrlen)
- return (x / d, y / d, z / d)
- except:
- return (0.0, 0.0, 1.0)
+ """Return vector (x,y,z) normalized by dividing by squared length.
+ Return (0.0, 0.0, 1.0) if the result is undefined."""
+ sqrlen = x * x + y * y + z * z
+ if sqrlen < 1e-100:
+ return (0.0, 0.0, 1.0)
+ else:
+ try:
+ d = sqrt(sqrlen)
+ return (x / d, y / d, z / d)
+ except:
+ return (0.0, 0.0, 1.0)
# We're using right-hand coord system, where
# forefinger=x, middle=y, thumb=z on right hand.
# Then, e.g., (1,0,0) x (0,1,0) = (0,0,1)
def Cross3(a, b):
- """Return the cross product of two vectors, a x b."""
+ """Return the cross product of two vectors, a x b."""
- (ax, ay, az) = a
- (bx, by, bz) = b
- return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
+ (ax, ay, az) = a
+ (bx, by, bz) = b
+ return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
def Dot2(a, b):
- """Return the dot product of two 2d vectors, a . b."""
+ """Return the dot product of two 2d vectors, a . b."""
- return a[0] * b[0] + a[1] * b[1]
+ return a[0] * b[0] + a[1] * b[1]
def Perp2(a, b):
- """Return a sort of 2d cross product."""
+ """Return a sort of 2d cross product."""
- return a[0] * b[1] - a[1] * b[0]
+ return a[0] * b[1] - a[1] * b[0]
def Sub2(a, b):
- """Return difference of 2d vectors, a-b."""
+ """Return difference of 2d vectors, a-b."""
- return (a[0] - b[0], a[1] - b[1])
+ return (a[0] - b[0], a[1] - b[1])
def Add2(a, b):
- """Return the sum of 2d vectors, a+b."""
+ """Return the sum of 2d vectors, a+b."""
- return (a[0] + b[0], a[1] + b[1])
+ return (a[0] + b[0], a[1] + b[1])
def Length2(v):
- """Return length of vector v=(x,y)."""
+ """Return length of vector v=(x,y)."""
- return sqrt(v[0]*v[0] + v[1]*v[1])
+ return sqrt(v[0] * v[0] + v[1] * v[1])
def LinInterp2(a, b, alpha):
- """Return the point alpha of the way from a to b."""
+ """Return the point alpha of the way from a to b."""
- beta = 1-alpha
- return (beta*a[0]+alpha*b[0], beta*a[1]+alpha*b[1])
+ beta = 1 - alpha
+ return (beta * a[0] + alpha * b[0], beta * a[1] + alpha * b[1])
def Normalized2(p):
- """Return vector p normlized by dividing by its squared length.
- Return (0.0, 1.0) if the result is undefined."""
+ """Return vector p normlized by dividing by its squared length.
+ Return (0.0, 1.0) if the result is undefined."""
- (x, y) = p
- sqrlen = x * x + y * y
- if sqrlen < 1e-100:
- return (0.0, 1.0)
- else:
- try:
- d = sqrt(sqrlen)
- return (x / d, y / d)
- except:
- return (0.0, 1.0)
+ (x, y) = p
+ sqrlen = x * x + y * y
+ if sqrlen < 1e-100:
+ return (0.0, 1.0)
+ else:
+ try:
+ d = sqrt(sqrlen)
+ return (x / d, y / d)
+ except:
+ return (0.0, 1.0)
def Angle(a, b, c, points):
- """Return Angle abc in degrees, in range [0,180),
- where a,b,c are indices into points."""
-
- u = Sub2(points.pos[c], points.pos[b])
- v = Sub2(points.pos[a], points.pos[b])
- n1 = Length2(u)
- n2 = Length2(v)
- if n1 == 0.0 or n2 == 0.0:
- return 0.0
- else:
- costheta = Dot2(u, v) / (n1 * n2)
- if costheta > 1.0:
- costheta = 1.0
- if costheta < - 1.0:
- costheta = - 1.0
- return math.acos(costheta) * 180.0 / math.pi
+ """Return Angle abc in degrees, in range [0,180),
+ where a,b,c are indices into points."""
+
+ u = Sub2(points.pos[c], points.pos[b])
+ v = Sub2(points.pos[a], points.pos[b])
+ n1 = Length2(u)
+ n2 = Length2(v)
+ if n1 == 0.0 or n2 == 0.0:
+ return 0.0
+ else:
+ costheta = Dot2(u, v) / (n1 * n2)
+ if costheta > 1.0:
+ costheta = 1.0
+ if costheta < - 1.0:
+ costheta = - 1.0
+ return math.acos(costheta) * 180.0 / math.pi
def SegsIntersect(ixa, ixb, ixc, ixd, points):
- """Return true if segment AB intersects CD,
- false if they just touch. ixa, ixb, ixc, ixd are indices
- into points."""
-
- a = points.pos[ixa]
- b = points.pos[ixb]
- c = points.pos[ixc]
- d = points.pos[ixd]
- u = Sub2(b, a)
- v = Sub2(d, c)
- w = Sub2(a, c)
- pp = Perp2(u, v)
- if abs(pp) > TOL:
- si = Perp2(v, w) / pp
- ti = Perp2(u, w) / pp
- return 0.0 < si < 1.0 and 0.0 < ti < 1.0
- else:
- # parallel or overlapping
- if Dot2(u, u) == 0.0 or Dot2(v, v) == 0.0:
- return False
+ """Return true if segment AB intersects CD,
+ false if they just touch. ixa, ixb, ixc, ixd are indices
+ into points."""
+
+ a = points.pos[ixa]
+ b = points.pos[ixb]
+ c = points.pos[ixc]
+ d = points.pos[ixd]
+ u = Sub2(b, a)
+ v = Sub2(d, c)
+ w = Sub2(a, c)
+ pp = Perp2(u, v)
+ if abs(pp) > TOL:
+ si = Perp2(v, w) / pp
+ ti = Perp2(u, w) / pp
+ return 0.0 < si < 1.0 and 0.0 < ti < 1.0
else:
- pp2 = Perp2(w, v)
- if abs(pp2) > TOL:
- return False # parallel, not collinear
- z = Sub2(b, c)
- (vx, vy) = v
- (wx, wy) = w
- (zx, zy) = z
- if vx == 0.0:
- (t0, t1) = (wy / vy, zy / vy)
- else:
- (t0, t1) = (wx / vx, zx / vx)
- return 0.0 < t0 < 1.0 and 0.0 < t1 < 1.0
+ # parallel or overlapping
+ if Dot2(u, u) == 0.0 or Dot2(v, v) == 0.0:
+ return False
+ else:
+ pp2 = Perp2(w, v)
+ if abs(pp2) > TOL:
+ return False # parallel, not collinear
+ z = Sub2(b, c)
+ (vx, vy) = v
+ (wx, wy) = w
+ (zx, zy) = z
+ if vx == 0.0:
+ (t0, t1) = (wy / vy, zy / vy)
+ else:
+ (t0, t1) = (wx / vx, zx / vx)
+ return 0.0 < t0 < 1.0 and 0.0 < t1 < 1.0
def Ccw(a, b, c, points):
- """Return true if ABC is a counterclockwise-oriented triangle,
- where a, b, and c are indices into points.
- Returns false if not, or if colinear within TOL."""
+ """Return true if ABC is a counterclockwise-oriented triangle,
+ where a, b, and c are indices into points.
+ Returns false if not, or if colinear within TOL."""
- (ax, ay) = (points.pos[a][0], points.pos[a][1])
- (bx, by) = (points.pos[b][0], points.pos[b][1])
- (cx, cy) = (points.pos[c][0], points.pos[c][1])
- d = ax * by - bx * ay - ax * cy + cx * ay + bx * cy - cx * by
- return d > TOL
+ (ax, ay) = (points.pos[a][0], points.pos[a][1])
+ (bx, by) = (points.pos[b][0], points.pos[b][1])
+ (cx, cy) = (points.pos[c][0], points.pos[c][1])
+ d = ax * by - bx * ay - ax * cy + cx * ay + bx * cy - cx * by
+ return d > TOL
def InCircle(a, b, c, d, points):
- """Return true if circle through points with indices a, b, c
- contains point with index d (indices into points).
- Except: if ABC forms a counterclockwise oriented triangle
- then the test is reversed: return true if d is outside the circle.
- Will get false, no matter what orientation, if d is cocircular, with TOL^2.
- | xa ya xa^2+ya^2 1 |
- | xb yb xb^2+yb^2 1 | > 0
- | xc yc xc^2+yc^2 1 |
- | xd yd xd^2+yd^2 1 |
- """
-
- (xa, ya, za) = _Icc(points.pos[a])
- (xb, yb, zb) = _Icc(points.pos[b])
- (xc, yc, zc) = _Icc(points.pos[c])
- (xd, yd, zd) = _Icc(points.pos[d])
- det = xa * (yb * zc - yc * zb - yb * zd + yd * zb + yc * zd - yd * zc) \
- - xb * (ya * zc - yc * za - ya * zd + yd * za + yc * zd - yd * zc) \
- + xc * (ya * zb - yb * za - ya * zd + yd * za + yb * zd - yd * zb) \
- - xd * (ya * zb - yb * za - ya * zc + yc * za + yb * zc - yc * zb)
- return det > TOL * TOL
+ """Return true if circle through points with indices a, b, c
+ contains point with index d (indices into points).
+ Except: if ABC forms a counterclockwise oriented triangle
+ then the test is reversed: return true if d is outside the circle.
+ Will get false, no matter what orientation, if d is cocircular, with TOL^2.
+ | xa ya xa^2+ya^2 1 |
+ | xb yb xb^2+yb^2 1 | > 0
+ | xc yc xc^2+yc^2 1 |
+ | xd yd xd^2+yd^2 1 |
+ """
+
+ (xa, ya, za) = _Icc(points.pos[a])
+ (xb, yb, zb) = _Icc(points.pos[b])
+ (xc, yc, zc) = _Icc(points.pos[c])
+ (xd, yd, zd) = _Icc(points.pos[d])
+ det = xa * (yb * zc - yc * zb - yb * zd + yd * zb + yc * zd - yd * zc) \
+ - xb * (ya * zc - yc * za - ya * zd + yd * za + yc * zd - yd * zc) \
+ + xc * (ya * zb - yb * za - ya * zd + yd * za + yb * zd - yd * zb) \
+ - xd * (ya * zb - yb * za - ya * zc + yc * za + yb * zc - yc * zb)
+ return det > TOL * TOL
def _Icc(p):
- (x, y) = (p[0], p[1])
- return (x, y, x * x + y * y)
+ (x, y) = (p[0], p[1])
+ return (x, y, x * x + y * y)