Removing double edges (reverse edge split) in a script

Question

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:
          for linked_e in v.link_edges:
              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:
          if len(e.link_faces) != 1:
              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')

Answer 1

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)

Answer 2

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 ]