Removing double edges (reverse edge split) in a script

I have a 2D model of a road that consists of separate rectangles whose ends are approximately lined up. Below an example of a situation where the road turns slightly.

There are thousands of rectangles, and that's a problem for later processing. I need to join ends of most of the roads such that they form a connected mesh, like so:

I can't simply perform a Remove Doubles on all the geometry, because this is done from a script, and the risk of joining vertices in a way that doesn't result in sensible geometry is too high. For example below only the bottom vertices of the edge have been joined. In intersections where roads overlap, completely unrelated roads might become connected in broken ways.

Thus what I need is an operation that joins each edge pair whose ends are sufficiently close to each other. Is there any way easier than writing that entirely by myself? If I were to write that myself, what would be the best strategy?

The end goal is to create a model that can be 3D printed nicely. I extrude the 2D roads, resulting in a tactile map that you get from Touch Mapper. Thousands of objects is a bit much for slicers, which basically need to first separate by loose parts, and then perform an union between every resulting object.

Below is an example of roads in a real-word city centre. Pedestrian roads (being edited in the screenshot) are extruded higher than other roads. I'd like to join eg. the edge whose vertices are hilighted.

Edit: final code

Below is my final code. Changes to Chebhou's answer are:

• Speed up edge proximity searching with mathutils.kdtree.KDTree. Important with over a few thousand edges.
• Merge all verts in one operation using bmesh.ops.weld_verts(). With 1000 verts to join, this shrinks execution time from several seconds to almost nothing.
• Filter out weird edges, and precompute some values in filter_edges().
• Don't merge edges in places where there more than 2 edges meet. This is a good heuristic that avoids all roads becoming a single object, which would not be desirable when 3D printing, because objects shouldn't intersect themselves.

Mesh is assumed to be in object mode initially.

def join_matching_edges(ob):
lt = 0.1  # length difference + -
dt = 0.1  # max distance
at = 0.5  # max sin(angle)  (30°)

bpy.context.scene.objects.active = ob
bpy.ops.object.mode_set(mode = 'EDIT')
bpy.ops.mesh.select_all(action='DESELECT')
from math import sin
bm = bmesh.from_edit_mesh( bpy.context.object.data )
bm.edges.ensure_lookup_table()

center   = lambda e : ( e.verts[0].co + e.verts[1].co ) / 2
length   = lambda e : ( e.verts[0].co - e.verts[1].co ).length
dist     = lambda v1, v2: (  v2 -  v1 ).length
sinAngle = lambda e1, e2: abs(sin((e1.verts[1].co - e1.verts[0].co).angle(e2.verts[1].co - e2.verts[0].co)))

def point_between_edge_neighbor_verts(e):
# Return middle of the verts adjacent to the edge
verts = []
for v in e.verts:
verts.extend((vv for vv in linked_e.verts if vv != e.verts[0] and vv != e.verts[1]))
if len(verts) != 2:
#print("edge has non-2 adjacent verts: " + str(len(verts)))
return None
return ((verts[0].co[0] + verts[1].co[0]) / 2, \
(verts[0].co[1] + verts[1].co[1]) / 2, \
(verts[0].co[2] + verts[1].co[2]) / 2)

class CEdge:
def __init__(self, e, center, into_edge):
self.e = e
self.center = center(e)
self.length = length(e)
self.into_edge = into_edge
self.welded = False

def filter_edges(edges):
out = []
for e in edges:
continue
point_between_edges = point_between_edge_neighbor_verts(e)
if not point_between_edges:
continue
vector_into_edge_face = center(e) - mathutils.Vector(point_between_edges)
out.append(CEdge(e, center, vector_into_edge_face / vector_into_edge_face.length))
return out
candidate_edges = filter_edges(bm.edges)

# Index edges into search tree
kd = mathutils.kdtree.KDTree(len(candidate_edges))
for i, ce in enumerate(candidate_edges):
kd.insert(ce.center, i)
kd.balance()

to_weld = {}
for i, ce in enumerate(candidate_edges[:-1]):
if ce.welded:
continue
ce.welded = True
lmin = ce.length - lt
lmax = ce.length + lt
matches = []
for (_co, oe_index, _dist) in kd.find_range(ce.center, dt):
oe = candidate_edges[oe_index]
if not oe.welded and lmin < oe.length < lmax and sinAngle(ce.e, oe.e) < at:
turn_angle = ce.into_edge.angle(-oe.into_edge)
if turn_angle > math.pi * 0.6: # pi * 0.5 is 90%
#print("not merging edges (%s, %s) pointing to opposite directions, angle is %f" % (ce.e, oe.e, turn_angle))
continue
matches.append(oe)
oe.welded = True

if len(matches) == 1:
# Join nothing where >2 ways meet, else all roads in the scene may become joined and cross itself
ev1, ev2 = ce.e.verts[:]
oev1, oev2 = matches[0].e.verts[:]
if dist(ev1.co, oev1.co) < dist(ev1.co, oev2.co) :
if ev1 != oev1 : to_weld[ev1] = oev1
if ev2 != oev2 : to_weld[ev2] = oev2
else :
if ev1 != oev2 : to_weld[ev1] = oev2
if ev2 != oev1 : to_weld[ev2] = oev1

print("%s: melding %d out of %d edges" % (ob.name, len(to_weld) / 2, len(bm.edges)))
bmesh.ops.weld_verts(bm, targetmap = to_weld)
bmesh.update_edit_mesh(bpy.context.object.data ,True)
bpy.ops.object.mode_set(mode = 'OBJECT')

• You could try select one of these edges and select others by direction from Select Similar operator (Shift+G). Adjusting Threshold in the operator settings, you could select only desired edges. Then remove doubles. This way won't really work if there're important edges similarly directed; probably you'll need to control visibility of the objects. Commented Apr 8, 2016 at 17:26
• There will be all sorts of edges in there, but Select Similar could still be useful as part of the solution. Thanks. Commented Apr 8, 2016 at 17:52

This is modified version of TLousky answer in this script we use three factors to decide merging edges :

• distance between their centers.
• difference in length.
• the angle between these edges.

you can drop any of them if it's not necessary in your scene, and this is the script :

import bpy, bmesh
from math import sin

bm = bmesh.from_edit_mesh( bpy.context.object.data )
bm.edges.ensure_lookup_table()

lt = 0.2 #length difference + -
dt = 0.2  #max distance
at = 0.16  #max sin(angle)  (30°)

center  = lambda e : ( e.verts[0].co + e.verts[1].co ) / 2
length  = lambda e : ( e.verts[0].co - e.verts[1].co ).length
dist    = lambda v1, v2: (  v2 -  v1 ).length
sinAngle   = lambda e1, e2: abs(sin((e1.verts[1].co - e1.verts[0].co).angle(e2.verts[1].co - e2.verts[0].co)))

i =-1
while i < len(bm.edges)-1:
i +=1
e = bm.edges[i]
i = e.index
l = length(e)
lmin = l-lt
lmax = l+lt
c = center(e)
for oe in  bm.edges[i+1:]:
if  lmin < length(oe) < lmax  and dist( center(oe), c) < dt and sinAngle(e, oe)<at :
ev1, ev2 = e.verts[:]
oev1, oev2 = oe.verts[:]
print( ev2, oev2 )
if dist(ev1.co, oev1.co ) < dist(ev1.co, oev2.co ) :
if ev1 != oev1 : bmesh.ops.pointmerge( bm, verts = [ ev1, oev1 ], merge_co = ev1.co )
if ev2 != oev2 : bmesh.ops.pointmerge( bm, verts = [ ev2, oev2 ], merge_co = ev2.co )
else :
if ev1 != oev2 : bmesh.ops.pointmerge( bm, verts = [ ev1, oev2 ], merge_co = ev1.co )
if ev2 != oev1 : bmesh.ops.pointmerge( bm, verts = [ ev2, oev1 ], merge_co = ev2.co )
bm.edges.ensure_lookup_table()
i=-1
break

bmesh.update_edit_mesh(bpy.context.object.data ,True)

• Works perfectly, after I added bm.edges.ensure_lookup_table() after bm's initialization. I added that into the answer's code. Thank you! Commented Apr 26, 2016 at 5:36

This script will delete edges that are closer to each other than the threshold specified edgeThreshold if both their vertices are closer than the vertThreshold.

I use two thresholds in order to avoid iterating over all the edges for each edge, and filtering faraway edges (will make a big difference if the mesh is large and has many edges).

It assumes your mesh is in edit mode.

import bpy, bmesh

bm = bmesh.from_edit_mesh( bpy.context.object.data )

mat           = bpy.context.object.matrix_world
vertThreshold = 0.3
edgeThreshold = 0.4

calcEdgeCenter = lambda e, mat: ( mat * e.verts[0].co + mat * e.verts[1].co ) / 2
dist           = lambda v1, v2, mat: ( mat * v2 - mat * v1 ).length

removedEdges = []
for i, e in enumerate( bm.edges ):
if i in removedEdges: continue

# Find current edge's center point
curEdgeCenter = calcEdgeCenter( e, mat )

# Find close edges (whose center is closer than the edgeThreshold)
otherEdges = [
oe for oe in bm.edges if oe.index != e.index and \
oe.index not in removedEdges and \
dist( curEdgeCenter, calcEdgeCenter( oe, mat ), mat ) < edgeThreshold
]

eV1, eV2 = e.verts[:]

for oe in otherEdges:
oeIdx = oe.index
oeV1, oeV2 = oe.verts[:]

# An edge will be considered redundant if both its verts are closer to
# the current edge's verts. In that case merge the edges' vertices
if dist( eV1.co, oeV1.co, mat ) < vertThreshold and dist( eV2.co, oeV2.co, mat ) < vertThreshold:
if eV1.index != oeV1.index:
bmesh.ops.pointmerge( bm, verts = [ eV1, oeV1 ], merge_co = mat * eV1.co )
if eV2.index != oeV2.index:
bmesh.ops.pointmerge( bm, verts = [ eV2, oeV2 ], merge_co = mat * eV2.co )
removedEdges += [ oeIdx, i ]
elif dist( eV1.co, oeV2.co, mat ) < vertThreshold and dist( eV2.co, oeV1.co, mat ) < vertThreshold:
if eV1.index != oeV2.index:
bmesh.ops.pointmerge( bm, verts = [ eV1, oeV2 ], merge_co = mat * eV1.co )
if eV2.index != oeV1.index:
bmesh.ops.pointmerge( bm, verts = [ eV2, oeV1 ], merge_co = mat * eV2.co )
removedEdges += [ oeIdx, i ]

• If you read my question carefully, I mention that simply removing double verts isn't reliable enough. The geometry is unpredictable, and completely unrelated verts may happen to be arbitrary close to each other. Commented Apr 8, 2016 at 17:42
• Missed that @SamuliKärkkäinen. Will update my answer to devise a better trick. Commented Apr 8, 2016 at 17:45
• @SamuliKärkkäinen, updated my answer with a new script that should do the trick. Commented Apr 8, 2016 at 18:38
• Outstanding answer! For me O(n*n) is too slow, but I'll see about speeding that up using KDTree. Commented Apr 9, 2016 at 7:01
• To make the script work on my scene, I needed to remove "mat" from the "merge_co = mat * eVx.co"'s. After that things work with low thresholds, but when I increase the thresholds, I get every time a crash, or "BMesh data of type BMVert has been removed" on "if dist( eV1.co, oeV1.co, mat ) < vertThreshold and dist( eV2.co, oeV2.co, mat ) < vertThreshold:". In the beginning of the inner loop oe.is_valid is sometimes false, but skipping those oe's doesn't end crashing. I can't figure this out. Commented Apr 9, 2016 at 19:24