Automated Edge Flip Operator

Last week I ran in to a long lasting Houdini frustration, the lack of proper Edge support.
We get the control over various primitive types, points and verts, but not Edges.

This becomes frustrating when dealing with complex primitive sorting and optimization techniques. For that purpose I decided to create my own Edge object, a rather simple wrapper that let’s me create Edges from points, verts and primitives inside Houdini. An edge can be constructed using two points, a primitive or 2 verts. Functionality includes comparing angles, getting the length, comparing other existing edges on point presence, not order etc.
Just for convenience One of the first things I needed it for was an automatic edge flip operator. The SOP samples primitive data and flips edges based on the optimal angle of neighbouring prims. Similar to the Delauney edge orientation principle. The flipping of an edge is only possible when the union of the two incident faces forms a convex quadrilateral. Flip Constraints are based on the curvature of a reference mesh, optimizing the triangle edge orientation on harsh edges or cutoff points.

To determine if the edge needs to be flipped, both edge lengths are computed, together with the maximum shared angle. If the angle (thus the edge length) is higher for the incident face, the edge can be flipped. Otherwise the angle is compared based on the maximum difference between the point normals and prim to point normals. If the difference comes close to 0, the edge follows curvature. If it is higher then the specified threshold, the edge should be flipped. (conforming the proposed mesh orientation and curvature).

One of the implementation issues: how do you keep track of the mesh when you flip an edge? Flipping an edge makes the next edge (sharing one of the previous edge primitives) most likely invalid, or a possible candidate. You could solve this by recursively update the mesh as you flip an edge, keeping track of the deforming state and marking areas that have been computed. This would be rather slow and tedious to incorporate. I ended up marking the edge that should be flipped, any new edge being sampled is being checked against any shared edge points. If one of the shared edge points match a flipped edge, it isn’t considered. Running the operator a couple of times would asure most edges will align. Most likely ater two computational rounds. I also added a “dirty” mode, not taking into account previously flipped edges. (works surprisingly well)

The actual operator doesn’t flip the edges. Houdini’s Edge Flip Sop does that for me. I simply count the total number of edges to be flipped and save the Edge Points as a Detail String Attribute.
The count is used to specify the number of iterations the loop should make and inside I simply sample each saved edge. It is a limitation I had to work with in Houdini.

A handy tip on speed and flexibility might be to pre-sample your data-set in to clusters. Clusters are formed based on any Vector3 attribute and are weighted over X recursions. Finding the most optimal distribution as it is being computed. Creating clusters (based on Position for example) speeds up the computing of the edges to be flipped, because the scope is limited to only a small island of primitives. In comparison to the entire mesh. Similar to using a Kd-Tree, but simpler. Read up on the algorithm here and here

Now for some code, the edge class should be easy to recreate based on what was told and can be seen in the following snippet:

from Triangulate.hedge import Edge

'''
Globals---------------------------------------------------------------------------------------------------
'''

# Helps me debug my code
debug = hou.pwd().evalParm("debug_mode")
export_name = hou.pwd().evalParm("edge_name")
flip_count_name = hou.pwd().evalParm("count_name")
export_groups = hou.pwd().evalParm("export_groups")
compare_angle = hou.pwd().evalParm("compare_angle")
angle_attrib_name = hou.pwd().evalParm("angle_attribute")
constrain_attrib_name = hou.pwd().evalParm("constrain_attribute")
constrain_angle = hou.pwd().evalParm("constrain_angle")
min_angle = hou.pwd().evalParm("min_angle")
max_angle = hou.pwd().evalParm("max_angle")
constrain_cut = hou.pwd().evalParm("constrain_cutoff")
safe_mode = hou.pwd().evalParm("safe_mode")

# This code is called when instances of this SOP cook.
node = hou.pwd()
geo = node.geometry()
prims = geo.prims()
points = geo.points()
edge_list = []
shared_edges = []
flipped_edges = []
flipped_points = []

'''
Necessary Attributes---------------------------------------------------------------------------------------------------
'''

# Find the count attribute or create it otherwise
count_attrib = geo.findGlobalAttrib(flip_count_name)
if count_attrib is None:
count_attrib = geo.addAttrib(hou.attribType.Global, flip_count_name, 0)
else:
geo.setGlobalAttribValue(count_attrib, 0)

angle_attrib = geo.findPointAttrib(angle_attrib_name)
if angle_attrib is None and compare_angle:
raise hou.NodeError("Can't find angle comparison point attribute: %s" % angle_attrib_name)

constrain_attrib = geo.findPointAttrib(constrain_attrib_name)
if constrain_attrib is None and constrain_angle and compare_angle:
raise hou.NodeError("Can't find constrain attribute: %s" % constrain_attrib_name)

'''
Methods---------------------------------------------------------------------------------------------------
'''

def Console(inString):
"""
Helps me to debug my code
"""

if debug:
print inString

def GetPrimPoints(inPrim):
"""
Get's the points out of a primitive
"""

point_list = []
for vert in inPrim.vertices():
point_list.append(vert.point())
return point_list

def GetHighestAngle(inNormals):
"""
Compares vectors to get the highest angle
"""

max_angle = 0
for n in xrange(len(inNormals)):
n1 = inNormals[n]
for ni in range(n+1, len(inNormals)):
n2 = inNormals[ni]
angle = n1.angleTo(n2)
if angle > max_angle:
max_angle = angle

return max_angle

def GetHighestPrimAngle(inNormals, inPrimNormal):
"""
Compares the normals against the primitive
"""

max_angle = 0
for normal in inNormals:
angle = normal.angleTo(inPrimNormal)
if angle > max_angle:
max_angle = angle

return max_angle

def FindUnsharedPoints(inPrimOne, inPrimTwo, inEdge):
"""
Returns a list of points that isn't shared between the primitives
"""

unshared_points = []
prim_one_points = GetPrimPoints(inPrimOne)
prim_two_points = GetPrimPoints(inPrimTwo)

# compare first primitve
for point in prim_one_points:
if point in inEdge.points:
continue
else:
unshared_points.append(point)

# compare second primitive
for point in prim_two_points:
if point in inEdge.points:
continue
else:
unshared_points.append(point)

return unshared_points

def SampleEdges(inPrim):
"""
Computes the primitive edges
New edges are appended to the edge_list
"""

Console("sampling prim: %d" % inPrim.number())

#Fetch the vertices
point_list = GetPrimPoints(inPrim)

# Now calculate the edge
amp = len(point_list)
for i in range(amp):
if i < (amp-1):
point_one = point_list[i]
point_two = point_list[i+1]
else:
point_one = point_list[-1:]
point_two = point_list

edge = Edge(point_one, point_two)

# if the edge is part of already existing edges, it shared a primitive
# i know this is dirty but i'd rather have this then do 2 lookups :(
try:
Console("shared edge: %s" % edge)
existing_edge = edge_list[edge_list.index(edge)]
shared_edges.append(existing_edge)
except ValueError:
Console("new edge: %s" % edge)
edge_list.append(edge)

def CompareAngle(inPrimOne, inPrimTwo, inEdge, inUnsharedPoints):
"""
Compares the angle of the point normal angle against primitive normal
Returns true if the edge needs to be inverted...
The comparison is based on the maximum difference between the prim point normals,
and the prim point normals compared with the prim normal. If the difference is rather
large on both sides, consider the edge to be flipped

An extra constrain attribute can be specified. Works best with curvature.
If the curvature of the unshared points is higher then the combined edge points,
It is allowed to flip, together with the angle constrain
"""

prim_one_points = GetPrimPoints(inPrimOne)
prim_two_points = GetPrimPoints(inPrimTwo)

# get the point normals
prim_one_normals = []
prim_two_normals = []
for point in prim_one_points:
prim_one_normals.append(hou.Vector3(point.attribValue(angle_attrib)))
for point in prim_two_points:
prim_two_normals.append(hou.Vector3(point.attribValue(angle_attrib)))

# figure out the angles
max_angle_one = GetHighestAngle(prim_one_normals)                                  # the highest measured angle between the point normals
max_angle_two = GetHighestAngle(prim_two_normals)
max_prim_angle_one = GetHighestPrimAngle(prim_one_normals, inPrimOne.normal())     # the highest measued angle between points and prim normal
max_prim_angle_two = GetHighestPrimAngle(prim_two_normals, inPrimTwo.normal())

angle_diff_one = abs(max_angle_one - max_prim_angle_one)                           # difference between prim 1 normals
angle_diff_two = abs(max_angle_two - max_prim_angle_two)                           # difference between prim 2 normals

max_value = max(angle_diff_one, angle_diff_two)                                    # get the higest of the two
min_value = min(angle_diff_one, angle_diff_two)                                    # get the lowest of the two

Console("normal angle one: %s" % max_angle_one)                                    # some debug stuff
Console("normal angle two: %s" % max_angle_two)
Console("prim   angle one: %s" % max_prim_angle_one)
Console("prim   angle two: %s" % max_prim_angle_two)
Console("angle diff one: %s" % angle_diff_one)
Console("angle diff two: %s" % angle_diff_two)

# figure out the additional attribute
if constrain_angle:
com_original = inEdge.pointOne.attribValue(constrain_attrib) + inEdge.pointTwo.attribValue(constrain_attrib)            # original edge point attribute values
com_diff = inUnsharedPoints.attribValue(constrain_attrib) + inUnsharedPoints.attribValue(constrain_attrib)        # possible new edge attribute values
constrain = com_diff / 2.0                                                                                              # average to be compared against attribute cutoff

if com_diff > com_original and constrain > constrain_cut and max_value >= max_angle and min_value >= min_angle:
return True
else:
if max_value >= max_angle and min_value >= min_angle:
return True

return False

def FlipEdge(inEdge):
"""
Checks if the edge needs to be flipped or not
"""

if (safe_mode) and (inEdge.pointOne in flipped_points or inEdge.pointTwo in flipped_points):
Console("ALREADY IN FLIPPED EDGES: %s" % inEdge)
return -1

prim_one = prims[inEdge.edgePrims]
prim_two = prims[inEdge.edgePrims]

# for both primitives find the point that isn't shared and figure out if it should be flipped
unshared_points = FindUnsharedPoints(prim_one, prim_two, inEdge)
if len(unshared_points) != 2:
("%s has more more or less then 2 unshared points" % inEdge)
return -1

Console("unshared points: %s, %s" % (unshared_points.number(), unshared_points.number()))

# check edge length
edge_two = Edge(unshared_points, unshared_points)
edgelength_one = inEdge.edgeLength
edgelength_two = edge_two.edgeLength

if edgelength_two < edgelength_one:
Console("FLIPPING EDGE!")
flipped_edges.append(inEdge)
if safe_mode:
flipped_points.append(inEdge.pointOne)
flipped_points.append(inEdge.pointTwo)
return 1

# check angle
elif (compare_angle) and (CompareAngle(prim_one, prim_two, inEdge, unshared_points)):
Console("FLIPPING EDGE!")
flipped_edges.append(inEdge)
if safe_mode:
flipped_points.append(inEdge.pointOne)
flipped_points.append(inEdge.pointTwo)
return 1

def Finalize():
"""
Creates the detail attributes and groups (if chosen)
The counter specifies the amount of edges to be flipped
"""

count = 0
for edge in flipped_edges:

# set or create the detail attribute
edge_attrib_name = "%s_%d" % (export_name, count)
edge_attrib = geo.findGlobalAttrib(edge_attrib_name)
if edge_attrib is None:
edge_attrib = geo.addAttrib(hou.attribType.Global, edge_attrib_name, "none")
geo.setGlobalAttribValue(edge_attrib,"p%d-%d" % (edge.pointNumberOne, edge.pointNumberTwo))

# see if we export the group
if export_groups:
group = geo.createPointGroup("%s%d_%d" % (export_name, edge.pointNumberOne, edge.pointNumberTwo))
for point in edge.points:

count += 1

# set the global attribute value
geo.setGlobalAttribValue(count_attrib, count)

'''
Compute---------------------------------------------------------------------------------------------------
'''

# For every prim construct 3 edges
# See if the edges have multiple primitives shared
# If multiple primitives are shared, do the comparison

# Calculate shared edges
for prim in prims:
SampleEdges(prim)
Console("\n")

# See if they need to be flipped
for edge in shared_edges:
Console("sampling edge: %s" % edge)
FlipEdge(edge)
Console("\n")

# Set up the edges to be flipped
Finalize()