mesh_looptools.py

# ##### BEGIN GPL LICENSE BLOCK #####
#
#  This program is free software; you can redistribute it and/or
#  modify it under the terms of the GNU General Public License
#  as published by the Free Software Foundation; either version 2
#  of the License, or (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software Foundation,
#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# ##### END GPL LICENSE BLOCK #####

bl_info = {
    "name": "LoopTools",
    "author": "Bart Crouch",
    "version": (4, 5, 2),
    "blender": (2, 69, 3),
    "location": "View3D > Toolbar and View3D > Specials (W-key)",
    "warning": "",
    "description": "Mesh modelling toolkit. Several tools to aid modelling",
    "wiki_url": "http://wiki.blender.org/index.php/Extensions:2.6/Py/"
        "Scripts/Modeling/LoopTools",
    "tracker_url": "https://developer.blender.org/T26189",
    "category": "Mesh"}


import bmesh
import bpy
import collections
import mathutils
import math
from bpy_extras import view3d_utils


##########################################
####### General functions ################
##########################################


# used by all tools to improve speed on reruns
looptools_cache = {}


# force a full recalculation next time
def cache_delete(tool):
    if tool in looptools_cache:
        del looptools_cache[tool]


# check cache for stored information
def cache_read(tool, object, bm, input_method, boundaries):
    # current tool not cached yet
    if tool not in looptools_cache:
        return(False, False, False, False, False)
    # check if selected object didn't change
    if object.name != looptools_cache[tool]["object"]:
        return(False, False, False, False, False)
    # check if input didn't change
    if input_method != looptools_cache[tool]["input_method"]:
        return(False, False, False, False, False)
    if boundaries != looptools_cache[tool]["boundaries"]:
        return(False, False, False, False, False)
    modifiers = [mod.name for mod in object.modifiers if mod.show_viewport \
        and mod.type == 'MIRROR']
    if modifiers != looptools_cache[tool]["modifiers"]:
        return(False, False, False, False, False)
    input = [v.index for v in bm.verts if v.select and not v.hide]
    if input != looptools_cache[tool]["input"]:
        return(False, False, False, False, False)
    # reading values
    single_loops = looptools_cache[tool]["single_loops"]
    loops = looptools_cache[tool]["loops"]
    derived = looptools_cache[tool]["derived"]
    mapping = looptools_cache[tool]["mapping"]

    return(True, single_loops, loops, derived, mapping)


# store information in the cache
def cache_write(tool, object, bm, input_method, boundaries, single_loops,
loops, derived, mapping):
    # clear cache of current tool
    if tool in looptools_cache:
        del looptools_cache[tool]
    # prepare values to be saved to cache
    input = [v.index for v in bm.verts if v.select and not v.hide]
    modifiers = [mod.name for mod in object.modifiers if mod.show_viewport \
        and mod.type == 'MIRROR']
    # update cache
    looptools_cache[tool] = {"input": input, "object": object.name,
        "input_method": input_method, "boundaries": boundaries,
        "single_loops": single_loops, "loops": loops,
        "derived": derived, "mapping": mapping, "modifiers": modifiers}


# calculates natural cubic splines through all given knots
def calculate_cubic_splines(bm_mod, tknots, knots):
    # hack for circular loops
    if knots[0] == knots[-1] and len(knots) > 1:
        circular = True
        k_new1 = []
        for k in range(-1, -5, -1):
            if k - 1 < -len(knots):
                k += len(knots)
            k_new1.append(knots[k-1])
        k_new2 = []
        for k in range(4):
            if k + 1 > len(knots) - 1:
                k -= len(knots)
            k_new2.append(knots[k+1])
        for k in k_new1:
            knots.insert(0, k)
        for k in k_new2:
            knots.append(k)
        t_new1 = []
        total1 = 0
        for t in range(-1, -5, -1):
            if t - 1 < -len(tknots):
                t += len(tknots)
            total1 += tknots[t] - tknots[t-1]
            t_new1.append(tknots[0] - total1)
        t_new2 = []
        total2 = 0
        for t in range(4):
            if t + 1 > len(tknots) - 1:
                t -= len(tknots)
            total2 += tknots[t+1] - tknots[t]
            t_new2.append(tknots[-1] + total2)
        for t in t_new1:
            tknots.insert(0, t)
        for t in t_new2:
            tknots.append(t)
    else:
        circular = False
    # end of hack

    n = len(knots)
    if n < 2:
        return False
    x = tknots[:]
    locs = [bm_mod.verts[k].co[:] for k in knots]
    result = []
    for j in range(3):
        a = []
        for i in locs:
            a.append(i[j])
        h = []
        for i in range(n-1):
            if x[i+1] - x[i] == 0:
                h.append(1e-8)
            else:
                h.append(x[i+1] - x[i])
        q = [False]
        for i in range(1, n-1):
            q.append(3/h[i]*(a[i+1]-a[i]) - 3/h[i-1]*(a[i]-a[i-1]))
        l = [1.0]
        u = [0.0]
        z = [0.0]
        for i in range(1, n-1):
            l.append(2*(x[i+1]-x[i-1]) - h[i-1]*u[i-1])
            if l[i] == 0:
                l[i] = 1e-8
            u.append(h[i] / l[i])
            z.append((q[i] - h[i-1] * z[i-1]) / l[i])
        l.append(1.0)
        z.append(0.0)
        b = [False for i in range(n-1)]
        c = [False for i in range(n)]
        d = [False for i in range(n-1)]
        c[n-1] = 0.0
        for i in range(n-2, -1, -1):
            c[i] = z[i] - u[i]*c[i+1]
            b[i] = (a[i+1]-a[i])/h[i] - h[i]*(c[i+1]+2*c[i])/3
            d[i] = (c[i+1]-c[i]) / (3*h[i])
        for i in range(n-1):
            result.append([a[i], b[i], c[i], d[i], x[i]])
    splines = []
    for i in range(len(knots)-1):
        splines.append([result[i], result[i+n-1], result[i+(n-1)*2]])
    if circular: # cleaning up after hack
        knots = knots[4:-4]
        tknots = tknots[4:-4]

    return(splines)


# calculates linear splines through all given knots
def calculate_linear_splines(bm_mod, tknots, knots):
    splines = []
    for i in range(len(knots)-1):
        a = bm_mod.verts[knots[i]].co
        b = bm_mod.verts[knots[i+1]].co
        d = b-a
        t = tknots[i]
        u = tknots[i+1]-t
        splines.append([a, d, t, u]) # [locStart, locDif, tStart, tDif]

    return(splines)


# calculate a best-fit plane to the given vertices
def calculate_plane(bm_mod, loop, method="best_fit", object=False):
    # getting the vertex locations
    locs = [bm_mod.verts[v].co.copy() for v in loop[0]]

    # calculating the center of masss
    com = mathutils.Vector()
    for loc in locs:
        com += loc
    com /= len(locs)
    x, y, z = com

    if method == 'best_fit':
        # creating the covariance matrix
        mat = mathutils.Matrix(((0.0, 0.0, 0.0),
                                (0.0, 0.0, 0.0),
                                (0.0, 0.0, 0.0),
                                ))
        for loc in locs:
            mat[0][0] += (loc[0]-x)**2
            mat[1][0] += (loc[0]-x)*(loc[1]-y)
            mat[2][0] += (loc[0]-x)*(loc[2]-z)
            mat[0][1] += (loc[1]-y)*(loc[0]-x)
            mat[1][1] += (loc[1]-y)**2
            mat[2][1] += (loc[1]-y)*(loc[2]-z)
            mat[0][2] += (loc[2]-z)*(loc[0]-x)
            mat[1][2] += (loc[2]-z)*(loc[1]-y)
            mat[2][2] += (loc[2]-z)**2

        # calculating the normal to the plane
        normal = False
        try:
            mat = matrix_invert(mat)
        except:
            ax = 2
            if math.fabs(sum(mat[0])) < math.fabs(sum(mat[1])):
                if math.fabs(sum(mat[0])) < math.fabs(sum(mat[2])):
                    ax = 0
            elif math.fabs(sum(mat[1])) < math.fabs(sum(mat[2])):
                ax = 1
            if ax == 0:
                normal = mathutils.Vector((1.0, 0.0, 0.0))
            elif ax == 1:
                normal = mathutils.Vector((0.0, 1.0, 0.0))
            else:
                normal = mathutils.Vector((0.0, 0.0, 1.0))
        if not normal:
            # warning! this is different from .normalize()
            itermax = 500
            iter = 0
            vec = mathutils.Vector((1.0, 1.0, 1.0))
            vec2 = (mat * vec)/(mat * vec).length
            while vec != vec2 and iter<itermax:
                iter+=1
                vec = vec2
                vec2 = mat * vec
                if vec2.length != 0:
                    vec2 /= vec2.length
            if vec2.length == 0:
                vec2 = mathutils.Vector((1.0, 1.0, 1.0))
            normal = vec2

    elif method == 'normal':
        # averaging the vertex normals
        v_normals = [bm_mod.verts[v].normal for v in loop[0]]
        normal = mathutils.Vector()
        for v_normal in v_normals:
            normal += v_normal
        normal /= len(v_normals)
        normal.normalize()

    elif method == 'view':
        # calculate view normal
        rotation = bpy.context.space_data.region_3d.view_matrix.to_3x3().\
            inverted()
        normal = rotation * mathutils.Vector((0.0, 0.0, 1.0))
        if object:
            normal = object.matrix_world.inverted().to_euler().to_matrix() * \
                     normal

    return(com, normal)


# calculate splines based on given interpolation method (controller function)
def calculate_splines(interpolation, bm_mod, tknots, knots):
    if interpolation == 'cubic':
        splines = calculate_cubic_splines(bm_mod, tknots, knots[:])
    else: # interpolations == 'linear'
        splines = calculate_linear_splines(bm_mod, tknots, knots[:])

    return(splines)


# check loops and only return valid ones
def check_loops(loops, mapping, bm_mod):
    valid_loops = []
    for loop, circular in 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
        # vertices can not all be at the same location
        stacked = True
        for i in range(len(loop) - 1):
            if (bm_mod.verts[loop[i]].co - \
            bm_mod.verts[loop[i+1]].co).length > 1e-6:
                stacked = False
                break
        if stacked:
            continue
        # passed all tests, loop is valid
        valid_loops.append([loop, circular])

    return(valid_loops)


# input: bmesh, output: dict with the edge-key as key and face-index as value
def dict_edge_faces(bm):
    edge_faces = dict([[edgekey(edge), []] for edge in bm.edges if \
        not edge.hide])
    for face in bm.faces:
        if face.hide:
            continue
        for key in face_edgekeys(face):
            edge_faces[key].append(face.index)

    return(edge_faces)


# input: bmesh (edge-faces optional), output: dict with face-face connections
def dict_face_faces(bm, edge_faces=False):
    if not edge_faces:
        edge_faces = dict_edge_faces(bm)

    connected_faces = dict([[face.index, []] for face in bm.faces if \
        not face.hide])
    for face in bm.faces:
        if face.hide:
            continue
        for edge_key in face_edgekeys(face):
            for connected_face in edge_faces[edge_key]:
                if connected_face == face.index:
                    continue
                connected_faces[face.index].append(connected_face)

    return(connected_faces)


# input: bmesh, output: dict with the vert index as key and edge-keys as value
def dict_vert_edges(bm):
    vert_edges = dict([[v.index, []] for v in bm.verts if not v.hide])
    for edge in bm.edges:
        if edge.hide:
            continue
        ek = edgekey(edge)
        for vert in ek:
            vert_edges[vert].append(ek)

    return(vert_edges)


# input: bmesh, output: dict with the vert index as key and face index as value
def dict_vert_faces(bm):
    vert_faces = dict([[v.index, []] for v in bm.verts if not v.hide])
    for face in bm.faces:
        if not face.hide:
            for vert in face.verts:
                vert_faces[vert.index].append(face.index)

    return(vert_faces)


# input: list of edge-keys, output: dictionary with vertex-vertex connections
def dict_vert_verts(edge_keys):
    # create connection data
    vert_verts = {}
    for ek in edge_keys:
        for i in range(2):
            if ek[i] in vert_verts:
                vert_verts[ek[i]].append(ek[1-i])
            else:
                vert_verts[ek[i]] = [ek[1-i]]

    return(vert_verts)


# return the edgekey ([v1.index, v2.index]) of a bmesh edge
def edgekey(edge):
    return(tuple(sorted([edge.verts[0].index, edge.verts[1].index])))


# returns the edgekeys of a bmesh face
def face_edgekeys(face):
    return([tuple(sorted([edge.verts[0].index, edge.verts[1].index])) for \
        edge in face.edges])


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

    # calculate selected loops
    edge_keys = [edgekey(edge) for edge in bm_mod.edges if \
        edge.select and not edge.hide]
    loops = get_connected_selections(edge_keys)

    # if only selected loops are needed, we're done
    if input == 'selected':
        return(derived, bm_mod, loops)
    # elif input == 'all':
    loops = get_parallel_loops(bm_mod, loops)

    return(derived, bm_mod, loops)


# sorts all edge-keys into a list of loops
def get_connected_selections(edge_keys):
    # create connection data
    vert_verts = dict_vert_verts(edge_keys)

    # find loops consisting of connected selected edges
    loops = []
    while len(vert_verts) > 0:
        loop = [iter(vert_verts.keys()).__next__()]
        growing = True
        flipped = False

        # extend loop
        while growing:
            # no more connection data for current vertex
            if loop[-1] not in vert_verts:
                if not flipped:
                    loop.reverse()
                    flipped = True
                else:
                    growing = False
            else:
                extended = False
                for i, next_vert in enumerate(vert_verts[loop[-1]]):
                    if next_vert not in loop:
                        vert_verts[loop[-1]].pop(i)
                        if len(vert_verts[loop[-1]]) == 0:
                            del vert_verts[loop[-1]]
                        # remove connection both ways
                        if next_vert in vert_verts:
                            if len(vert_verts[next_vert]) == 1:
                                del vert_verts[next_vert]
                            else:
                                vert_verts[next_vert].remove(loop[-1])
                        loop.append(next_vert)
                        extended = True
                        break
                if not extended:
                    # found one end of the loop, continue with next
                    if not flipped:
                        loop.reverse()
                        flipped = True
                    # found both ends of the loop, stop growing
                    else:
                        growing = False

        # check if loop is circular
        if loop[0] in vert_verts:
            if loop[-1] in vert_verts[loop[0]]:
                # is circular
                if len(vert_verts[loop[0]]) == 1:
                    del vert_verts[loop[0]]
                else:
                    vert_verts[loop[0]].remove(loop[-1])
                if len(vert_verts[loop[-1]]) == 1:
                    del vert_verts[loop[-1]]
                else:
                    vert_verts[loop[-1]].remove(loop[0])
                loop = [loop, True]
            else:
                # not circular
                loop = [loop, False]
        else:
            # not circular
            loop = [loop, False]

        loops.append(loop)

    return(loops)


# get the derived mesh data, if there is a mirror modifier
def get_derived_bmesh(object, bm, scene):
    # check for mirror modifiers
    if 'MIRROR' in [mod.type for mod in object.modifiers if mod.show_viewport]:
        derived = True
        # disable other modifiers
        show_viewport = [mod.name for mod in object.modifiers if \
            mod.show_viewport]
        for mod in object.modifiers:
            if mod.type != 'MIRROR':
                mod.show_viewport = False
        # get derived mesh
        bm_mod = bmesh.new()
        mesh_mod = object.to_mesh(scene, True, 'PREVIEW')
        bm_mod.from_mesh(mesh_mod)
        bpy.context.blend_data.meshes.remove(mesh_mod)
        # re-enable other modifiers
        for mod_name in show_viewport:
            object.modifiers[mod_name].show_viewport = True
    # no mirror modifiers, so no derived mesh necessary
    else:
        derived = False
        bm_mod = bm

    return(derived, bm_mod)


# return a mapping of derived indices to indices
def get_mapping(derived, bm, bm_mod, single_vertices, full_search, loops):
    if not derived:
        return(False)

    if full_search:
        verts = [v for v in bm.verts if not v.hide]
    else:
        verts = [v for v in bm.verts if v.select and not v.hide]

    # non-selected vertices around single vertices also need to be mapped
    if single_vertices:
        mapping = dict([[vert, -1] for vert in single_vertices])
        verts_mod = [bm_mod.verts[vert] for vert in single_vertices]
        for v in verts:
            for v_mod in verts_mod:
                if (v.co - v_mod.co).length < 1e-6:
                    mapping[v_mod.index] = v.index
                    break
        real_singles = [v_real for v_real in mapping.values() if v_real>-1]

        verts_indices = [vert.index for vert in verts]
        for face in [face for face in bm.faces if not face.select \
        and not face.hide]:
            for vert in face.verts:
                if vert.index in real_singles:
                    for v in face.verts:
                        if not v.index in verts_indices:
                            if v not in verts:
                                verts.append(v)
                    break

    # create mapping of derived indices to indices
    mapping = dict([[vert, -1] for loop in loops for vert in loop[0]])
    if single_vertices:
        for single in single_vertices:
            mapping[single] = -1
    verts_mod = [bm_mod.verts[i] for i in mapping.keys()]
    for v in verts:
        for v_mod in verts_mod:
            if (v.co - v_mod.co).length < 1e-6:
                mapping[v_mod.index] = v.index
                verts_mod.remove(v_mod)
                break

    return(mapping)


# calculate the determinant of a matrix
def matrix_determinant(m):
    determinant = m[0][0] * m[1][1] * m[2][2] + m[0][1] * m[1][2] * m[2][0] \
        + m[0][2] * m[1][0] * m[2][1] - m[0][2] * m[1][1] * m[2][0] \
        - m[0][1] * m[1][0] * m[2][2] - m[0][0] * m[1][2] * m[2][1]

    return(determinant)


# custom matrix inversion, to provide higher precision than the built-in one
def matrix_invert(m):
    r = mathutils.Matrix((
        (m[1][1]*m[2][2] - m[1][2]*m[2][1], m[0][2]*m[2][1] - m[0][1]*m[2][2],
        m[0][1]*m[1][2] - m[0][2]*m[1][1]),
        (m[1][2]*m[2][0] - m[1][0]*m[2][2], m[0][0]*m[2][2] - m[0][2]*m[2][0],
        m[0][2]*m[1][0] - m[0][0]*m[1][2]),
        (m[1][0]*m[2][1] - m[1][1]*m[2][0], m[0][1]*m[2][0] - m[0][0]*m[2][1],
        m[0][0]*m[1][1] - m[0][1]*m[1][0])))

    return (r * (1 / matrix_determinant(m)))


# returns a list of all loops parallel to the input, input included
def get_parallel_loops(bm_mod, loops):
    # get required dictionaries
    edge_faces = dict_edge_faces(bm_mod)
    connected_faces = dict_face_faces(bm_mod, edge_faces)
    # turn vertex loops into edge loops
    edgeloops = []
    for loop in loops:
        edgeloop = [[sorted([loop[0][i], loop[0][i+1]]) for i in \
            range(len(loop[0])-1)], loop[1]]
        if loop[1]: # circular
            edgeloop[0].append(sorted([loop[0][-1], loop[0][0]]))
        edgeloops.append(edgeloop[:])
    # variables to keep track while iterating
    all_edgeloops = []
    has_branches = False

    for loop in edgeloops:
        # initialise with original loop
        all_edgeloops.append(loop[0])
        newloops = [loop[0]]
        verts_used = []
        for edge in loop[0]:
            if edge[0] not in verts_used:
                verts_used.append(edge[0])
            if edge[1] not in verts_used:
                verts_used.append(edge[1])

        # find parallel loops
        while len(newloops) > 0:
            side_a = []
            side_b = []
            for i in newloops[-1]:
                i = tuple(i)
                forbidden_side = False
                if not i in edge_faces:
                    # weird input with branches
                    has_branches = True
                    break
                for face in edge_faces[i]:
                    if len(side_a) == 0 and forbidden_side != "a":
                        side_a.append(face)
                        if forbidden_side:
                            break
                        forbidden_side = "a"
                        continue
                    elif side_a[-1] in connected_faces[face] and \
                    forbidden_side != "a":
                        side_a.append(face)
                        if forbidden_side:
                            break
                        forbidden_side = "a"
                        continue
                    if len(side_b) == 0 and forbidden_side != "b":
                        side_b.append(face)
                        if forbidden_side:
                            break
                        forbidden_side = "b"
                        continue
                    elif side_b[-1] in connected_faces[face] and \
                    forbidden_side != "b":
                        side_b.append(face)
                        if forbidden_side:
                            break
                        forbidden_side = "b"
                        continue

            if has_branches:
                # weird input with branches
                break

            newloops.pop(-1)
            sides = []
            if side_a:
                sides.append(side_a)
            if side_b:
                sides.append(side_b)

            for side in sides:
                extraloop = []
                for fi in side:
                    for key in face_edgekeys(bm_mod.faces[fi]):
                        if key[0] not in verts_used and key[1] not in \
                        verts_used:
                            extraloop.append(key)
                            break
                if extraloop:
                    for key in extraloop:
                        for new_vert in key:
                            if new_vert not in verts_used:
                                verts_used.append(new_vert)
                    newloops.append(extraloop)
                    all_edgeloops.append(extraloop)

    # input contains branches, only return selected loop
    if has_branches:
        return(loops)

    # change edgeloops into normal loops
    loops = []
    for edgeloop in all_edgeloops:
        loop = []
        # grow loop by comparing vertices between consecutive edge-keys
        for i in range(len(edgeloop)-1):
            for vert in range(2):
                if edgeloop[i][vert] in edgeloop[i+1]:
                    loop.append(edgeloop[i][vert])
                    break
        if loop:
            # add starting vertex
            for vert in range(2):
                if edgeloop[0][vert] != loop[0]:
                    loop = [edgeloop[0][vert]] + loop
                    break
            # add ending vertex
            for vert in range(2):
                if edgeloop[-1][vert] != loop[-1]:
                    loop.append(edgeloop[-1][vert])
                    break
            # check if loop is circular
            if loop[0] == loop[-1]:
                circular = True
                loop = loop[:-1]
            else:
                circular = False
        loops.append([loop, circular])

    return(loops)


# gather initial data
def initialise():
    global_undo = bpy.context.user_preferences.edit.use_global_undo
    bpy.context.user_preferences.edit.use_global_undo = False
    object = bpy.context.active_object
    if 'MIRROR' in [mod.type for mod in object.modifiers if mod.show_viewport]:
        # ensure that selection is synced for the derived mesh
        bpy.ops.object.mode_set(mode='OBJECT')
        bpy.ops.object.mode_set(mode='EDIT')
    bm = bmesh.from_edit_mesh(object.data)

    return(global_undo, object, bm)


# move the vertices to their new locations
def move_verts(object, bm, mapping, move, influence):
    for loop in move:
        for index, loc in loop:
            if mapping:
                if mapping[index] == -1:
                    continue
                else:
                    index = mapping[index]
            if influence >= 0:
                bm.verts[index].co = loc*(influence/100) + \
                    bm.verts[index].co*((100-influence)/100)
            else:
                bm.verts[index].co = loc
    bm.normal_update()
    object.data.update()


# load custom tool settings
def settings_load(self):
    lt = bpy.context.window_manager.looptools
    tool = self.name.split()[0].lower()
    keys = self.as_keywords().keys()
    for key in keys:
        setattr(self, key, getattr(lt, tool + "_" + key))


# store custom tool settings
def settings_write(self):
    lt = bpy.context.window_manager.looptools
    tool = self.name.split()[0].lower()
    keys = self.as_keywords().keys()
    for key in keys:
        setattr(lt, tool + "_" + key, getattr(self, key))


# clean up and set settings back to original state
def terminate(global_undo):
    # update editmesh cached data
    obj = bpy.context.active_object
    if obj.mode == 'EDIT':
        bmesh.update_edit_mesh(obj.data, tessface=True, destructive=True)

    bpy.context.user_preferences.edit.use_global_undo = global_undo


##########################################
####### Bridge functions #################
##########################################

# calculate a cubic spline through the middle section of 4 given coordinates
def bridge_calculate_cubic_spline(bm, coordinates):
    result = []
    x = [0, 1, 2, 3]

    for j in range(3):
        a = []
        for i in coordinates:
            a.append(float(i[j]))
        h = []
        for i in range(3):
            h.append(x[i+1]-x[i])
        q = [False]
        for i in range(1,3):
            q.append(3.0/h[i]*(a[i+1]-a[i])-3.0/h[i-1]*(a[i]-a[i-1]))
        l = [1.0]
        u = [0.0]
        z = [0.0]
        for i in range(1,3):
            l.append(2.0*(x[i+1]-x[i-1])-h[i-1]*u[i-1])
            u.append(h[i]/l[i])
            z.append((q[i]-h[i-1]*z[i-1])/l[i])
        l.append(1.0)
        z.append(0.0)
        b = [False for i in range(3)]
        c = [False for i in range(4)]
        d = [False for i in range(3)]
        c[3] = 0.0
        for i in range(2,-1,-1):
            c[i] = z[i]-u[i]*c[i+1]
            b[i] = (a[i+1]-a[i])/h[i]-h[i]*(c[i+1]+2.0*c[i])/3.0
            d[i] = (c[i+1]-c[i])/(3.0*h[i])
        for i in range(3):
            result.append([a[i], b[i], c[i], d[i], x[i]])
    spline = [result[1], result[4], result[7]]

    return(spline)


# return a list with new vertex location vectors, a list with face vertex
# integers, and the highest vertex integer in the virtual mesh
def bridge_calculate_geometry(bm, lines, vertex_normals, segments,
interpolation, cubic_strength, min_width, max_vert_index):
    new_verts = []
    faces = []

    # calculate location based on interpolation method
    def get_location(line, segment, splines):
        v1 = bm.verts[lines[line][0]].co
        v2 = bm.verts[lines[line][1]].co
        if interpolation == 'linear':
            return v1 + (segment/segments) * (v2-v1)
        else: # interpolation == 'cubic'
            m = (segment/segments)
            ax,bx,cx,dx,tx = splines[line][0]
            x = ax+bx*m+cx*m**2+dx*m**3
            ay,by,cy,dy,ty = splines[line][1]
            y = ay+by*m+cy*m**2+dy*m**3
            az,bz,cz,dz,tz = splines[line][2]
            z = az+bz*m+cz*m**2+dz*m**3
            return mathutils.Vector((x, y, z))

    # no interpolation needed
    if segments == 1:
        for i, line in enumerate(lines):
            if i < len(lines)-1:
                faces.append([line[0], lines[i+1][0], lines[i+1][1], line[1]])
    # more than 1 segment, interpolate
    else:
        # calculate splines (if necessary) once, so no recalculations needed
        if interpolation == 'cubic':
            splines = []
            for line in lines:
                v1 = bm.verts[line[0]].co
                v2 = bm.verts[line[1]].co
                size = (v2-v1).length * cubic_strength
                splines.append(bridge_calculate_cubic_spline(bm,
                    [v1+size*vertex_normals[line[0]], v1, v2,
                    v2+size*vertex_normals[line[1]]]))
        else:
            splines = False

        # create starting situation
        virtual_width = [(bm.verts[lines[i][0]].co -
                          bm.verts[lines[i+1][0]].co).length for i
                          in range(len(lines)-1)]
        new_verts = [get_location(0, seg, splines) for seg in range(1,
            segments)]
        first_line_indices = [i for i in range(max_vert_index+1,
            max_vert_index+segments)]

        prev_verts = new_verts[:] # vertex locations of verts on previous line
        prev_vert_indices = first_line_indices[:]
        max_vert_index += segments - 1 # highest vertex index in virtual mesh
        next_verts = [] # vertex locations of verts on current line
        next_vert_indices = []

        for i, line in enumerate(lines):
            if i < len(lines)-1:
                v1 = line[0]
                v2 = lines[i+1][0]
                end_face = True
                for seg in range(1, segments):
                    loc1 = prev_verts[seg-1]
                    loc2 = get_location(i+1, seg, splines)
                    if (loc1-loc2).length < (min_width/100)*virtual_width[i] \
                    and line[1]==lines[i+1][1]:
                        # triangle, no new vertex
                        faces.append([v1, v2, prev_vert_indices[seg-1],
                            prev_vert_indices[seg-1]])
                        next_verts += prev_verts[seg-1:]
                        next_vert_indices += prev_vert_indices[seg-1:]
                        end_face = False
                        break
                    else:
                        if i == len(lines)-2 and lines[0] == lines[-1]:
                            # quad with first line, no new vertex
                            faces.append([v1, v2, first_line_indices[seg-1],
                                prev_vert_indices[seg-1]])
                            v2 = first_line_indices[seg-1]
                            v1 = prev_vert_indices[seg-1]
                        else:
                            # quad, add new vertex
                            max_vert_index += 1
                            faces.append([v1, v2, max_vert_index,
                                prev_vert_indices[seg-1]])
                            v2 = max_vert_index
                            v1 = prev_vert_indices[seg-1]
                            new_verts.append(loc2)
                            next_verts.append(loc2)
                            next_vert_indices.append(max_vert_index)
                if end_face:
                    faces.append([v1, v2, lines[i+1][1], line[1]])

                prev_verts = next_verts[:]
                prev_vert_indices = next_vert_indices[:]
                next_verts = []
                next_vert_indices = []

    return(new_verts, faces, max_vert_index)


# calculate lines (list of lists, vertex indices) that are used for bridging
def bridge_calculate_lines(bm, loops, mode, twist, reverse):
    lines = []
    loop1, loop2 = [i[0] for i in loops]
    loop1_circular, loop2_circular = [i[1] for i in loops]
    circular = loop1_circular or loop2_circular
    circle_full = False

    # calculate loop centers
    centers = []
    for loop in [loop1, loop2]:
        center = mathutils.Vector()
        for vertex in loop:
            center += bm.verts[vertex].co
        center /= len(loop)
        centers.append(center)
    for i, loop in enumerate([loop1, loop2]):
        for vertex in loop:
            if bm.verts[vertex].co == centers[i]:
                # prevent zero-length vectors in angle comparisons
                centers[i] += mathutils.Vector((0.01, 0, 0))
                break
    center1, center2 = centers

    # calculate the normals of the virtual planes that the loops are on
    normals = []
    normal_plurity = False
    for i, loop in enumerate([loop1, loop2]):
        # covariance matrix
        mat = mathutils.Matrix(((0.0, 0.0, 0.0),
                                (0.0, 0.0, 0.0),
                                (0.0, 0.0, 0.0)))
        x, y, z = centers[i]
        for loc in [bm.verts[vertex].co for vertex in loop]:
            mat[0][0] += (loc[0]-x)**2
            mat[1][0] += (loc[0]-x)*(loc[1]-y)
            mat[2][0] += (loc[0]-x)*(loc[2]-z)
            mat[0][1] += (loc[1]-y)*(loc[0]-x)
            mat[1][1] += (loc[1]-y)**2
            mat[2][1] += (loc[1]-y)*(loc[2]-z)
            mat[0][2] += (loc[2]-z)*(loc[0]-x)
            mat[1][2] += (loc[2]-z)*(loc[1]-y)
            mat[2][2] += (loc[2]-z)**2
        # plane normal
        normal = False
        if sum(mat[0]) < 1e-6 or sum(mat[1]) < 1e-6 or sum(mat[2]) < 1e-6:
            normal_plurity = True
        try:
            mat.invert()
        except:
            if sum(mat[0]) == 0:
                normal = mathutils.Vector((1.0, 0.0, 0.0))
            elif sum(mat[1]) == 0:
                normal = mathutils.Vector((0.0, 1.0, 0.0))
            elif sum(mat[2]) == 0:
                normal = mathutils.Vector((0.0, 0.0, 1.0))
        if not normal:
            # warning! this is different from .normalize()
            itermax = 500
            iter = 0
            vec = mathutils.Vector((1.0, 1.0, 1.0))
            vec2 = (mat * vec)/(mat * vec).length
            while vec != vec2 and iter<itermax:
                iter+=1
                vec = vec2
                vec2 = mat * vec