mesh_looptools.py

                if vec2.length != 0:
                    vec2 /= vec2.length
            if vec2.length == 0:
                vec2 = mathutils.Vector((1.0, 1.0, 1.0))
            normal = vec2
        normals.append(normal)
    # have plane normals face in the same direction (maximum angle: 90 degrees)
    if ((center1 + normals[0]) - center2).length < \
    ((center1 - normals[0]) - center2).length:
        normals[0].negate()
    if ((center2 + normals[1]) - center1).length > \
    ((center2 - normals[1]) - center1).length:
        normals[1].negate()

    # rotation matrix, representing the difference between the plane normals
    axis = normals[0].cross(normals[1])
    axis = mathutils.Vector([loc if abs(loc) > 1e-8 else 0 for loc in axis])
    if axis.angle(mathutils.Vector((0, 0, 1)), 0) > 1.5707964:
        axis.negate()
    angle = normals[0].dot(normals[1])
    rotation_matrix = mathutils.Matrix.Rotation(angle, 4, axis)

    # if circular, rotate loops so they are aligned
    if circular:
        # make sure loop1 is the circular one (or both are circular)
        if loop2_circular and not loop1_circular:
            loop1_circular, loop2_circular = True, False
            loop1, loop2 = loop2, loop1

        # match start vertex of loop1 with loop2
        target_vector = bm.verts[loop2[0]].co - center2
        dif_angles = [[(rotation_matrix * (bm.verts[vertex].co - center1)
                       ).angle(target_vector, 0), False, i] for
                       i, vertex in enumerate(loop1)]
        dif_angles.sort()
        if len(loop1) != len(loop2):
            angle_limit = dif_angles[0][0] * 1.2 # 20% margin
            dif_angles = [[(bm.verts[loop2[0]].co - \
                bm.verts[loop1[index]].co).length, angle, index] for \
                angle, distance, index in dif_angles if angle <= angle_limit]
            dif_angles.sort()
        loop1 = loop1[dif_angles[0][2]:] + loop1[:dif_angles[0][2]]

    # have both loops face the same way
    if normal_plurity and not circular:
        second_to_first, second_to_second, second_to_last = \
            [(bm.verts[loop1[1]].co - center1).\
            angle(bm.verts[loop2[i]].co - center2) for i in [0, 1, -1]]
        last_to_first, last_to_second = [(bm.verts[loop1[-1]].co - \
            center1).angle(bm.verts[loop2[i]].co - center2) for \
            i in [0, 1]]
        if (min(last_to_first, last_to_second)*1.1 < min(second_to_first, \
        second_to_second)) or (loop2_circular and second_to_last*1.1 < \
        min(second_to_first, second_to_second)):
            loop1.reverse()
            if circular:
                loop1 = [loop1[-1]] + loop1[:-1]
    else:
        angle = (bm.verts[loop1[0]].co - center1).\
            cross(bm.verts[loop1[1]].co - center1).angle(normals[0], 0)
        target_angle = (bm.verts[loop2[0]].co - center2).\
            cross(bm.verts[loop2[1]].co - center2).angle(normals[1], 0)
        limit = 1.5707964 # 0.5*pi, 90 degrees
        if not ((angle > limit and target_angle > limit) or \
        (angle < limit and target_angle < limit)):
            loop1.reverse()
            if circular:
                loop1 = [loop1[-1]] + loop1[:-1]
        elif normals[0].angle(normals[1]) > limit:
            loop1.reverse()
            if circular:
                loop1 = [loop1[-1]] + loop1[:-1]

    # both loops have the same length
    if len(loop1) == len(loop2):
        # manual override
        if twist:
            if abs(twist) < len(loop1):
                loop1 = loop1[twist:]+loop1[:twist]
        if reverse:
            loop1.reverse()

        lines.append([loop1[0], loop2[0]])
        for i in range(1, len(loop1)):
            lines.append([loop1[i], loop2[i]])

    # loops of different lengths
    else:
        # make loop1 longest loop
        if len(loop2) > len(loop1):
            loop1, loop2 = loop2, loop1
            loop1_circular, loop2_circular = loop2_circular, loop1_circular

        # manual override
        if twist:
            if abs(twist) < len(loop1):
                loop1 = loop1[twist:]+loop1[:twist]
        if reverse:
            loop1.reverse()

        # shortest angle difference doesn't always give correct start vertex
        if loop1_circular and not loop2_circular:
            shifting = 1
            while shifting:
                if len(loop1) - shifting < len(loop2):
                    shifting = False
                    break
                to_last, to_first = [(rotation_matrix *
                    (bm.verts[loop1[-1]].co - center1)).angle((bm.\
                    verts[loop2[i]].co - center2), 0) for i in [-1, 0]]
                if to_first < to_last:
                    loop1 = [loop1[-1]] + loop1[:-1]
                    shifting += 1
                else:
                    shifting = False
                    break

        # basic shortest side first
        if mode == 'basic':
            lines.append([loop1[0], loop2[0]])
            for i in range(1, len(loop1)):
                if i >= len(loop2) - 1:
                    # triangles
                    lines.append([loop1[i], loop2[-1]])
                else:
                    # quads
                    lines.append([loop1[i], loop2[i]])

        # shortest edge algorithm
        else: # mode == 'shortest'
            lines.append([loop1[0], loop2[0]])
            prev_vert2 = 0
            for i in range(len(loop1) -1):
                if prev_vert2 == len(loop2) - 1 and not loop2_circular:
                    # force triangles, reached end of loop2
                    tri, quad = 0, 1
                elif prev_vert2 == len(loop2) - 1 and loop2_circular:
                    # at end of loop2, but circular, so check with first vert
                    tri, quad = [(bm.verts[loop1[i+1]].co -
                                  bm.verts[loop2[j]].co).length
                                 for j in [prev_vert2, 0]]

                    circle_full = 2
                elif len(loop1) - 1 - i == len(loop2) - 1 - prev_vert2 and \
                not circle_full:
                    # force quads, otherwise won't make it to end of loop2
                    tri, quad = 1, 0
                else:
                    # calculate if tri or quad gives shortest edge
                    tri, quad = [(bm.verts[loop1[i+1]].co -
                                  bm.verts[loop2[j]].co).length
                                 for j in range(prev_vert2, prev_vert2+2)]

                # triangle
                if tri < quad:
                    lines.append([loop1[i+1], loop2[prev_vert2]])
                    if circle_full == 2:
                        circle_full = False
                # quad
                elif not circle_full:
                    lines.append([loop1[i+1], loop2[prev_vert2+1]])
                    prev_vert2 += 1
                # quad to first vertex of loop2
                else:
                    lines.append([loop1[i+1], loop2[0]])
                    prev_vert2 = 0
                    circle_full = True

    # final face for circular loops
    if loop1_circular and loop2_circular:
        lines.append([loop1[0], loop2[0]])

    return(lines)


# calculate number of segments needed
def bridge_calculate_segments(bm, lines, loops, segments):
    # return if amount of segments is set by user
    if segments != 0:
        return segments

    # edge lengths
    average_edge_length = [(bm.verts[vertex].co - \
        bm.verts[loop[0][i+1]].co).length for loop in loops for \
        i, vertex in enumerate(loop[0][:-1])]
    # closing edges of circular loops
    average_edge_length += [(bm.verts[loop[0][-1]].co - \
        bm.verts[loop[0][0]].co).length for loop in loops if loop[1]]

    # average lengths
    average_edge_length = sum(average_edge_length) / len(average_edge_length)
    average_bridge_length = sum([(bm.verts[v1].co - \
        bm.verts[v2].co).length for v1, v2 in lines]) / len(lines)

    segments = max(1, round(average_bridge_length / average_edge_length))

    return(segments)


# return dictionary with vertex index as key, and the normal vector as value
def bridge_calculate_virtual_vertex_normals(bm, lines, loops, edge_faces,
edgekey_to_edge):
    if not edge_faces: # interpolation isn't set to cubic
        return False

    # pity reduce() isn't one of the basic functions in python anymore
    def average_vector_dictionary(dic):
        for key, vectors in dic.items():
            #if type(vectors) == type([]) and len(vectors) > 1:
            if len(vectors) > 1:
                average = mathutils.Vector()
                for vector in vectors:
                    average += vector
                average /= len(vectors)
                dic[key] = [average]
        return dic

    # get all edges of the loop
    edges = [[edgekey_to_edge[tuple(sorted([loops[j][0][i],
        loops[j][0][i+1]]))] for i in range(len(loops[j][0])-1)] for \
        j in [0,1]]
    edges = edges[0] + edges[1]
    for j in [0, 1]:
        if loops[j][1]: # circular
            edges.append(edgekey_to_edge[tuple(sorted([loops[j][0][0],
                loops[j][0][-1]]))])

    """
    calculation based on face topology (assign edge-normals to vertices)

    edge_normal = face_normal x edge_vector
    vertex_normal = average(edge_normals)
    """
    vertex_normals = dict([(vertex, []) for vertex in loops[0][0]+loops[1][0]])
    for edge in edges:
        faces = edge_faces[edgekey(edge)] # valid faces connected to edge

        if faces:
            # get edge coordinates
            v1, v2 = [bm.verts[edgekey(edge)[i]].co for i in [0,1]]
            edge_vector = v1 - v2
            if edge_vector.length < 1e-4:
                # zero-length edge, vertices at same location
                continue
            edge_center = (v1 + v2) / 2

            # average face coordinates, if connected to more than 1 valid face
            if len(faces) > 1:
                face_normal = mathutils.Vector()
                face_center = mathutils.Vector()
                for face in faces:
                    face_normal += face.normal
                    face_center += face.calc_center_median()
                face_normal /= len(faces)
                face_center /= len(faces)
            else:
                face_normal = faces[0].normal
                face_center = faces[0].calc_center_median()
            if face_normal.length < 1e-4:
                # faces with a surface of 0 have no face normal
                continue

            # calculate virtual edge normal
            edge_normal = edge_vector.cross(face_normal)
            edge_normal.length = 0.01
            if (face_center - (edge_center + edge_normal)).length > \
            (face_center - (edge_center - edge_normal)).length:
                # make normal face the correct way
                edge_normal.negate()
            edge_normal.normalize()
            # add virtual edge normal as entry for both vertices it connects
            for vertex in edgekey(edge):
                vertex_normals[vertex].append(edge_normal)

    """
    calculation based on connection with other loop (vertex focused method)
    - used for vertices that aren't connected to any valid faces

    plane_normal = edge_vector x connection_vector
    vertex_normal = plane_normal x edge_vector
    """
    vertices = [vertex for vertex, normal in vertex_normals.items() if not \
        normal]

    if vertices:
        # edge vectors connected to vertices
        edge_vectors = dict([[vertex, []] for vertex in vertices])
        for edge in edges:
            for v in edgekey(edge):
                if v in edge_vectors:
                    edge_vector = bm.verts[edgekey(edge)[0]].co - \
                        bm.verts[edgekey(edge)[1]].co
                    if edge_vector.length < 1e-4:
                        # zero-length edge, vertices at same location
                        continue
                    edge_vectors[v].append(edge_vector)

        # connection vectors between vertices of both loops
        connection_vectors = dict([[vertex, []] for vertex in vertices])
        connections = dict([[vertex, []] for vertex in vertices])
        for v1, v2 in lines:
            if v1 in connection_vectors or v2 in connection_vectors:
                new_vector = bm.verts[v1].co - bm.verts[v2].co
                if new_vector.length < 1e-4:
                    # zero-length connection vector,
                    # vertices in different loops at same location
                    continue
                if v1 in connection_vectors:
                    connection_vectors[v1].append(new_vector)
                    connections[v1].append(v2)
                if v2 in connection_vectors:
                    connection_vectors[v2].append(new_vector)
                    connections[v2].append(v1)
        connection_vectors = average_vector_dictionary(connection_vectors)
        connection_vectors = dict([[vertex, vector[0]] if vector else \
            [vertex, []] for vertex, vector in connection_vectors.items()])

        for vertex, values in edge_vectors.items():
            # vertex normal doesn't matter, just assign a random vector to it
            if not connection_vectors[vertex]:
                vertex_normals[vertex] = [mathutils.Vector((1, 0, 0))]
                continue

            # calculate to what location the vertex is connected,
            # used to determine what way to flip the normal
            connected_center = mathutils.Vector()
            for v in connections[vertex]:
                connected_center += bm.verts[v].co
            if len(connections[vertex]) > 1:
                connected_center /= len(connections[vertex])
            if len(connections[vertex]) == 0:
                # shouldn't be possible, but better safe than sorry
                vertex_normals[vertex] = [mathutils.Vector((1, 0, 0))]
                continue

            # can't do proper calculations, because of zero-length vector
            if not values:
                if (connected_center - (bm.verts[vertex].co + \
                connection_vectors[vertex])).length < (connected_center - \
                (bm.verts[vertex].co - connection_vectors[vertex])).\
                length:
                    connection_vectors[vertex].negate()
                vertex_normals[vertex] = [connection_vectors[vertex].\
                    normalized()]
                continue

            # calculate vertex normals using edge-vectors,
            # connection-vectors and the derived plane normal
            for edge_vector in values:
                plane_normal = edge_vector.cross(connection_vectors[vertex])
                vertex_normal = edge_vector.cross(plane_normal)
                vertex_normal.length = 0.1
                if (connected_center - (bm.verts[vertex].co + \
                vertex_normal)).length < (connected_center - \
                (bm.verts[vertex].co - vertex_normal)).length:
                # make normal face the correct way
                    vertex_normal.negate()
                vertex_normal.normalize()
                vertex_normals[vertex].append(vertex_normal)

    # average virtual vertex normals, based on all edges it's connected to
    vertex_normals = average_vector_dictionary(vertex_normals)
    vertex_normals = dict([[vertex, vector[0]] for vertex, vector in \
        vertex_normals.items()])

    return(vertex_normals)


# add vertices to mesh
def bridge_create_vertices(bm, vertices):
    for i in range(len(vertices)):
        bm.verts.new(vertices[i])


# add faces to mesh
def bridge_create_faces(object, bm, faces, twist):
    # have the normal point the correct way
    if twist < 0:
        [face.reverse() for face in faces]
        faces = [face[2:]+face[:2] if face[0]==face[1] else face for \
            face in faces]

    # eekadoodle prevention
    for i in range(len(faces)):
        if not faces[i][-1]:
            if faces[i][0] == faces[i][-1]:
                faces[i] = [faces[i][1], faces[i][2], faces[i][3], faces[i][1]]
            else:
                faces[i] = [faces[i][-1]] + faces[i][:-1]
        # result of converting from pre-bmesh period
        if faces[i][-1] == faces[i][-2]:
            faces[i] = faces[i][:-1]

    new_faces = []
    for i in range(len(faces)):
        new_faces.append(bm.faces.new([bm.verts[v] for v in faces[i]]))
    bm.normal_update()
    object.data.update(calc_edges=True) # calc_edges prevents memory-corruption

    return(new_faces)


# calculate input loops
def bridge_get_input(bm):
    # create list of internal edges, which should be skipped
    eks_of_selected_faces = [item for sublist in [face_edgekeys(face) for \
        face in bm.faces if face.select and not face.hide] for item in sublist]
    edge_count = {}
    for ek in eks_of_selected_faces:
        if ek in edge_count:
            edge_count[ek] += 1
        else:
            edge_count[ek] = 1
    internal_edges = [ek for ek in edge_count if edge_count[ek] > 1]

    # sort correct edges into loops
    selected_edges = [edgekey(edge) for edge in bm.edges if edge.select \
        and not edge.hide and edgekey(edge) not in internal_edges]
    loops = get_connected_selections(selected_edges)

    return(loops)


# return values needed by the bridge operator
def bridge_initialise(bm, interpolation):
    if interpolation == 'cubic':
        # dict with edge-key as key and list of connected valid faces as value
        face_blacklist = [face.index for face in bm.faces if face.select or \
            face.hide]
        edge_faces = dict([[edgekey(edge), []] for edge in bm.edges if not \
            edge.hide])
        for face in bm.faces:
            if face.index in face_blacklist:
                continue
            for key in face_edgekeys(face):
                edge_faces[key].append(face)
        # dictionary with the edge-key as key and edge as value
        edgekey_to_edge = dict([[edgekey(edge), edge] for edge in \
            bm.edges if edge.select and not edge.hide])
    else:
        edge_faces = False
        edgekey_to_edge = False

    # selected faces input
    old_selected_faces = [face.index for face in bm.faces if face.select \
        and not face.hide]

    # find out if faces created by bridging should be smoothed
    smooth = False
    if bm.faces:
        if sum([face.smooth for face in bm.faces])/len(bm.faces) \
        >= 0.5:
            smooth = True

    return(edge_faces, edgekey_to_edge, old_selected_faces, smooth)


# return a string with the input method
def bridge_input_method(loft, loft_loop):
    method = ""
    if loft:
        if loft_loop:
            method = "Loft loop"
        else:
            method = "Loft no-loop"
    else:
        method = "Bridge"

    return(method)


# match up loops in pairs, used for multi-input bridging
def bridge_match_loops(bm, loops):
    # calculate average loop normals and centers
    normals = []
    centers = []
    for vertices, circular in loops:
        normal = mathutils.Vector()
        center = mathutils.Vector()
        for vertex in vertices:
            normal += bm.verts[vertex].normal
            center += bm.verts[vertex].co
        normals.append(normal / len(vertices) / 10)
        centers.append(center / len(vertices))

    # possible matches if loop normals are faced towards the center
    # of the other loop
    matches = dict([[i, []] for i in range(len(loops))])
    matches_amount = 0
    for i in range(len(loops) + 1):
        for j in range(i+1, len(loops)):
            if (centers[i] - centers[j]).length > (centers[i] - (centers[j] \
            + normals[j])).length and (centers[j] - centers[i]).length > \
            (centers[j] - (centers[i] + normals[i])).length:
                matches_amount += 1
                matches[i].append([(centers[i] - centers[j]).length, i, j])
                matches[j].append([(centers[i] - centers[j]).length, j, i])
    # if no loops face each other, just make matches between all the loops
    if matches_amount == 0:
        for i in range(len(loops) + 1):
            for j in range(i+1, len(loops)):
                matches[i].append([(centers[i] - centers[j]).length, i, j])
                matches[j].append([(centers[i] - centers[j]).length, j, i])
    for key, value in matches.items():
        value.sort()

    # matches based on distance between centers and number of vertices in loops
    new_order = []
    for loop_index in range(len(loops)):
        if loop_index in new_order:
            continue
        loop_matches = matches[loop_index]
        if not loop_matches:
            continue
        shortest_distance = loop_matches[0][0]
        shortest_distance *= 1.1
        loop_matches = [[abs(len(loops[loop_index][0]) - \
            len(loops[loop[2]][0])), loop[0], loop[1], loop[2]] for loop in \
            loop_matches if loop[0] < shortest_distance]
        loop_matches.sort()
        for match in loop_matches:
            if match[3] not in new_order:
                new_order += [loop_index, match[3]]
                break

    # reorder loops based on matches
    if len(new_order) >= 2:
        loops = [loops[i] for i in new_order]

    return(loops)


# remove old_selected_faces
def bridge_remove_internal_faces(bm, old_selected_faces):
    # collect bmesh faces and internal bmesh edges
    remove_faces = [bm.faces[face] for face in old_selected_faces]
    edges = collections.Counter([edge.index for face in remove_faces for \
        edge in face.edges])
    remove_edges = [bm.edges[edge] for edge in edges if edges[edge] > 1]

    # remove internal faces and edges
    for face in remove_faces:
        bm.faces.remove(face)
    for edge in remove_edges:
        bm.edges.remove(edge)


# update list of internal faces that are flagged for removal
def bridge_save_unused_faces(bm, old_selected_faces, loops):
    # key: vertex index, value: lists of selected faces using it
    vertex_to_face = dict([[i, []] for i in range(len(bm.verts))])
    [[vertex_to_face[vertex.index].append(face) for vertex in \
        bm.faces[face].verts] for face in old_selected_faces]

    # group selected faces that are connected
    groups = []
    grouped_faces = []
    for face in old_selected_faces:
        if face in grouped_faces:
            continue
        grouped_faces.append(face)
        group = [face]
        new_faces = [face]
        while new_faces:
            grow_face = new_faces[0]
            for vertex in bm.faces[grow_face].verts:
                vertex_face_group = [face for face in vertex_to_face[\
                    vertex.index] if face not in grouped_faces]
                new_faces += vertex_face_group
                grouped_faces += vertex_face_group
                group += vertex_face_group
            new_faces.pop(0)
        groups.append(group)

    # key: vertex index, value: True/False (is it in a loop that is used)
    used_vertices = dict([[i, 0] for i in range(len(bm.verts))])
    for loop in loops:
        for vertex in loop[0]:
            used_vertices[vertex] = True

    # check if group is bridged, if not remove faces from internal faces list
    for group in groups:
        used = False
        for face in group:
            if used:
                break
            for vertex in bm.faces[face].verts:
                if used_vertices[vertex.index]:
                    used = True
                    break
        if not used:
            for face in group:
                old_selected_faces.remove(face)


# add the newly created faces to the selection
def bridge_select_new_faces(new_faces, smooth):
    for face in new_faces:
        face.select_set(True)
        face.smooth = smooth


# sort loops, so they are connected in the correct order when lofting
def bridge_sort_loops(bm, loops, loft_loop):
    # simplify loops to single points, and prepare for pathfinding
    x, y, z = [[sum([bm.verts[i].co[j] for i in loop[0]]) / \
        len(loop[0]) for loop in loops] for j in range(3)]
    nodes = [mathutils.Vector((x[i], y[i], z[i])) for i in range(len(loops))]

    active_node = 0
    open = [i for i in range(1, len(loops))]
    path = [[0,0]]
    # connect node to path, that is shortest to active_node
    while len(open) > 0:
        distances = [(nodes[active_node] - nodes[i]).length for i in open]
        active_node = open[distances.index(min(distances))]
        open.remove(active_node)
        path.append([active_node, min(distances)])
    # check if we didn't start in the middle of the path
    for i in range(2, len(path)):
        if (nodes[path[i][0]]-nodes[0]).length < path[i][1]:
            temp = path[:i]
            path.reverse()
            path = path[:-i] + temp
            break

    # reorder loops
    loops = [loops[i[0]] for i in path]
    # if requested, duplicate first loop at last position, so loft can loop
    if loft_loop:
        loops = loops + [loops[0]]

    return(loops)


# remapping old indices to new position in list
def bridge_update_old_selection(bm, old_selected_faces):
    #old_indices = old_selected_faces[:]
    #old_selected_faces = []
    #for i, face in enumerate(bm.faces):
    #    if face.index in old_indices:
    #        old_selected_faces.append(i)

    old_selected_faces = [i for i, face in enumerate(bm.faces) if face.index \
        in old_selected_faces]

    return(old_selected_faces)


##########################################
####### Circle functions #################
##########################################

# convert 3d coordinates to 2d coordinates on plane
def circle_3d_to_2d(bm_mod, loop, com, normal):
    # project vertices onto the plane
    verts = [bm_mod.verts[v] for v in loop[0]]
    verts_projected = [[v.co - (v.co - com).dot(normal) * normal, v.index]
                       for v in verts]

    # calculate two vectors (p and q) along the plane
    m = mathutils.Vector((normal[0] + 1.0, normal[1], normal[2]))
    p = m - (m.dot(normal) * normal)
    if p.dot(p) == 0.0:
        m = mathutils.Vector((normal[0], normal[1] + 1.0, normal[2]))
        p = m - (m.dot(normal) * normal)
    q = p.cross(normal)

    # change to 2d coordinates using perpendicular projection
    locs_2d = []
    for loc, vert in verts_projected:
        vloc = loc - com
        x = p.dot(vloc) / p.dot(p)
        y = q.dot(vloc) / q.dot(q)
        locs_2d.append([x, y, vert])

    return(locs_2d, p, q)


# calculate a best-fit circle to the 2d locations on the plane
def circle_calculate_best_fit(locs_2d):
    # initial guess
    x0 = 0.0
    y0 = 0.0
    r = 1.0

    # calculate center and radius (non-linear least squares solution)
    for iter in range(500):
        jmat = []
        k = []
        for v in locs_2d:
            d = (v[0]**2-2.0*x0*v[0]+v[1]**2-2.0*y0*v[1]+x0**2+y0**2)**0.5
            jmat.append([(x0-v[0])/d, (y0-v[1])/d, -1.0])
            k.append(-(((v[0]-x0)**2+(v[1]-y0)**2)**0.5-r))
        jmat2 = mathutils.Matrix(((0.0, 0.0, 0.0),
                                  (0.0, 0.0, 0.0),
                                  (0.0, 0.0, 0.0),
                                  ))
        k2 = mathutils.Vector((0.0, 0.0, 0.0))
        for i in range(len(jmat)):
            k2 += mathutils.Vector(jmat[i])*k[i]
            jmat2[0][0] += jmat[i][0]**2
            jmat2[1][0] += jmat[i][0]*jmat[i][1]
            jmat2[2][0] += jmat[i][0]*jmat[i][2]
            jmat2[1][1] += jmat[i][1]**2
            jmat2[2][1] += jmat[i][1]*jmat[i][2]
            jmat2[2][2] += jmat[i][2]**2
        jmat2[0][1] = jmat2[1][0]
        jmat2[0][2] = jmat2[2][0]
        jmat2[1][2] = jmat2[2][1]
        try:
            jmat2.invert()
        except:
            pass
        dx0, dy0, dr = jmat2 * k2
        x0 += dx0
        y0 += dy0
        r += dr
        # stop iterating if we're close enough to optimal solution
        if abs(dx0)<1e-6 and abs(dy0)<1e-6 and abs(dr)<1e-6:
            break

    # return center of circle and radius
    return(x0, y0, r)


# calculate circle so no vertices have to be moved away from the center
def circle_calculate_min_fit(locs_2d):
    # center of circle
    x0 = (min([i[0] for i in locs_2d])+max([i[0] for i in locs_2d]))/2.0
    y0 = (min([i[1] for i in locs_2d])+max([i[1] for i in locs_2d]))/2.0
    center = mathutils.Vector([x0, y0])
    # radius of circle
    r = min([(mathutils.Vector([i[0], i[1]])-center).length for i in locs_2d])

    # return center of circle and radius
    return(x0, y0, r)


# calculate the new locations of the vertices that need to be moved
def circle_calculate_verts(flatten, bm_mod, locs_2d, com, p, q, normal):
    # changing 2d coordinates back to 3d coordinates
    locs_3d = []
    for loc in locs_2d:
        locs_3d.append([loc[2], loc[0]*p + loc[1]*q + com])

    if flatten: # flat circle
        return(locs_3d)

    else: # project the locations on the existing mesh
        vert_edges = dict_vert_edges(bm_mod)
        vert_faces = dict_vert_faces(bm_mod)
        faces = [f for f in bm_mod.faces if not f.hide]
        rays = [normal, -normal]
        new_locs = []
        for loc in locs_3d:
            projection = False
            if bm_mod.verts[loc[0]].co == loc[1]: # vertex hasn't moved
                projection = loc[1]
            else:
                dif = normal.angle(loc[1]-bm_mod.verts[loc[0]].co)
                if -1e-6 < dif < 1e-6 or math.pi-1e-6 < dif < math.pi+1e-6:
                    # original location is already along projection normal
                    projection = bm_mod.verts[loc[0]].co
                else:
                    # quick search through adjacent faces
                    for face in vert_faces[loc[0]]:
                        verts = [v.co for v in bm_mod.faces[face].verts]
                        if len(verts) == 3: # triangle
                            v1, v2, v3 = verts
                            v4 = False
                        else: # assume quad
                            v1, v2, v3, v4 = verts[:4]
                        for ray in rays:
                            intersect = mathutils.geometry.\
                            intersect_ray_tri(v1, v2, v3, ray, loc[1])
                            if intersect:
                                projection = intersect
                                break
                            elif v4:
                                intersect = mathutils.geometry.\
                                intersect_ray_tri(v1, v3, v4, ray, loc[1])
                                if intersect:
                                    projection = intersect
                                    break
                        if projection:
                            break
            if not projection:
                # check if projection is on adjacent edges
                for edgekey in vert_edges[loc[0]]:
                    line1 = bm_mod.verts[edgekey[0]].co
                    line2 = bm_mod.verts[edgekey[1]].co
                    intersect, dist = mathutils.geometry.intersect_point_line(\
                        loc[1], line1, line2)
                    if 1e-6 < dist < 1 - 1e-6:
                        projection = intersect
                        break
            if not projection:
                # full search through the entire mesh
                hits = []
                for face in faces:
                    verts = [v.co for v in face.verts]
                    if len(verts) == 3: # triangle
                        v1, v2, v3 = verts
                        v4 = False
                    else: # assume quad
                        v1, v2, v3, v4 = verts[:4]
                    for ray in rays:
                        intersect = mathutils.geometry.intersect_ray_tri(\
                            v1, v2, v3, ray, loc[1])
                        if intersect:
                            hits.append([(loc[1] - intersect).length,
                                intersect])
                            break
                        elif v4:
                            intersect = mathutils.geometry.intersect_ray_tri(\
                                v1, v3, v4, ray, loc[1])
                            if intersect:
                                hits.append([(loc[1] - intersect).length,
                                    intersect])
                                break
                if len(hits) >= 1:
                    # if more than 1 hit with mesh, closest hit is new loc
                    hits.sort()
                    projection = hits[0][1]
            if not projection:
                # nothing to project on, remain at flat location
                projection = loc[1]
            new_locs.append([loc[0], projection])

        # return new positions of projected circle
        return(new_locs)


# check loops and only return valid ones
def circle_check_loops(single_loops, loops, mapping, bm_mod):
    valid_single_loops = {}
    valid_loops = []
    for i, [loop, circular] in enumerate(loops):
        # loop needs to have at least 3 vertices
        if len(loop) < 3:
            continue
        # loop needs at least 1 vertex in the original, non-mirrored mesh
        if mapping:
            all_virtual = True
            for vert in loop:
                if mapping[vert] > -1:
                    all_virtual = False
                    break
            if all_virtual:
                continue
        # loop has to be non-collinear
        collinear = True
        loc0 = mathutils.Vector(bm_mod.verts[loop[0]].co[:])
        loc1 = mathutils.Vector(bm_mod.verts[loop[1]].co[:])
        for v in loop[2:]:
            locn = mathutils.Vector(bm_mod.verts[v].co[:])
            if loc0 == loc1 or loc1 == locn:
                loc0 = loc1
                loc1 = locn
                continue
            d1 = loc1-loc0
            d2 = locn-loc1
            if -1e-6 < d1.angle(d2, 0) < 1e-6:
                loc0 = loc1
                loc1 = locn
                continue
            collinear = False
            break
        if collinear:
            continue
        # passed all tests, loop is valid
        valid_loops.append([loop, circular])
        valid_single_loops[len(valid_loops)-1] = single_loops[i]

    return(valid_single_loops, valid_loops)


# calculate the location of single input vertices that need to be flattened
def circle_flatten_singles(bm_mod, com, p, q, normal, single_loop):
    new_locs = []
    for vert in single_loop:
        loc = mathutils.Vector(bm_mod.verts[vert].co[:])
        new_locs.append([vert,  loc - (loc-com).dot(normal)*normal])

    return(new_locs)


# calculate input loops
def circle_get_input(object, bm, scene):
    # get mesh with modifiers applied
    derived, bm_mod = get_derived_bmesh(object, bm, scene)

    # create list of edge-keys based on selection state
    faces = False
    for face in bm.faces:
        if face.select and not face.hide:
            faces = True
            break
    if faces:
        # get selected, non-hidden , non-internal edge-keys
        eks_selected = [key for keys in [face_edgekeys(face) for face in \
            bm_mod.faces if face.select and not face.hide] for key in keys]
        edge_count = {}
        for ek in eks_selected:
            if ek in edge_count:
                edge_count[ek] += 1
            else:
                edge_count[ek] = 1
        edge_keys = [edgekey(edge) for edge in bm_mod.edges if edge.select \
            and not edge.hide and edge_count.get(edgekey(edge), 1)==1]
    else:
        # no faces, so no internal edges either
        edge_keys = [edgekey(edge) for edge in bm_mod.edges if edge.select \
            and not edge.hide]

    # add edge-keys around single vertices
    verts_connected = dict([[vert, 1] for edge in [edge for edge in \
        bm_mod.edges if edge.select and not edge.hide] for vert in \
        edgekey(edge)])
    single_vertices = [vert.index for vert in bm_mod.verts if \
        vert.select and not vert.hide and not \
        verts_connected.get(vert.index, False)]

    if single_vertices and len(bm.faces)>0:
        vert_to_single = dict([[v.index, []] for v in bm_mod.verts \
            if not v.hide])
        for face in [face for face in bm_mod.faces if not face.select \
        and not face.hide]:
            for vert in face.verts:
                vert = vert.index
                if vert in single_vertices:
                    for ek in face_edgekeys(face):
                        if not vert in ek:
                            edge_keys.append(ek)
                            if vert not in vert_to_single[ek[0]]:
                                vert_to_single[ek[0]].append(vert)
                            if vert not in vert_to_single[ek[1]]:
                                vert_to_single[ek[1]].append(vert)
                    break

    # sort edge-keys into loops
    loops = get_connected_selections(edge_keys)

    # find out to which loops the single vertices belong
    single_loops = dict([[i, []] for i in range(len(loops))])
    if single_vertices and len(bm.faces)>0:
        for i, [loop, circular] in enumerate(loops):
            for vert in loop:
                if vert_to_single[vert]:
                    for single in vert_to_single[vert]:
                        if single not in single_loops[i]:
                            single_loops[i].append(single)

    return(derived, bm_mod, single_vertices, single_loops, loops)


# recalculate positions based on the influence of the circle shape
def circle_influence_locs(locs_2d, new_locs_2d, influence):
    for i in range(len(locs_2d)):
        oldx, oldy, j = locs_2d[i]
        newx, newy, k = new_locs_2d[i]
        altx = newx*(influence/100)+ oldx*((100-influence)/100)
        alty = newy*(influence/100)+ oldy*((100-influence)/100)
        locs_2d[i] = [altx, alty, j]

    return(locs_2d)


# project 2d locations on circle, respecting distance relations between verts
def circle_project_non_regular(locs_2d, x0, y0, r):
    for i in range(len(locs_2d)):
        x, y, j = locs_2d[i]
        loc = mathutils.Vector([x-x0, y-y0])
        loc.length = r
        locs_2d[i] = [loc[0], loc[1], j]

    return(locs_2d)


# project 2d locations on circle, with equal distance between all vertices
def circle_project_regular(locs_2d, x0, y0, r):
    # find offset angle and circling direction
    x, y, i = locs_2d[0]
    loc = mathutils.Vector([x-x0, y-y0])
    loc.length = r
    offset_angle = loc.angle(mathutils.Vector([1.0, 0.0]), 0.0)
    loca = mathutils.Vector([x-x0, y-y0, 0.0])
    if loc[1] < -1e-6:
        offset_angle *= -1
    x, y, j = locs_2d[1]
    locb = mathutils.Vector([x-x0, y-y0, 0.0])
    if loca.cross(locb)[2] >= 0:
        ccw = 1
    else:
        ccw = -1
    # distribute vertices along the circle
    for i in range(len(locs_2d)):
        t = offset_angle + ccw * (i / len(locs_2d) * 2 * math.pi)
        x = math.cos(t) * r