mesh_bsurfaces.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; version 2
#  of the License.
#
#  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": "Bsurfaces GPL Edition",
    "author": "Eclectiel, Spivak Vladimir(cwolf3d)",
    "version": (1, 5, 2),
    "blender": (2, 80, 0),
    "location": "View3D > EditMode > ToolShelf",
    "description": "Modeling and retopology tool",
    "wiki_url": "https://wiki.blender.org/index.php/Dev:Ref/Release_Notes/2.64/Bsurfaces_1.5",
    "category": "Mesh",
}


import bpy
import bmesh
from bpy_extras import object_utils

import operator
from mathutils import Matrix, Vector
from mathutils.geometry import (
        intersect_line_line,
        intersect_point_line,
        )
from math import (
        degrees,
        pi,
        sqrt,
        )
from bpy.props import (
        BoolProperty,
        FloatProperty,
        IntProperty,
        StringProperty,
        PointerProperty,
        )
from bpy.types import (
        Operator,
        Panel,
        PropertyGroup,
        AddonPreferences,
        )


class VIEW3D_PT_tools_SURFSK_mesh(Panel):
    bl_space_type = 'VIEW_3D'
    bl_region_type = 'UI'
    bl_category = 'Tools'
    #bl_context = "mesh_edit"
    bl_label = "Bsurfaces"

    def draw(self, context):
        layout = self.layout
        scn = context.scene.bsurfaces

        col = layout.column(align=True)
        row = layout.row()
        row.separator()
        col.operator("gpencil.surfsk_init", text="Initialize")
        col.prop(scn, "SURFSK_object_with_retopology")
        col.prop(scn, "SURFSK_object_with_strokes")
        col.separator()
        col.operator("gpencil.surfsk_add_surface", text="Add Surface")
        col.operator("gpencil.surfsk_edit_surface", text="Edit Surface")
        col.operator("gpencil.surfsk_add_strokes", text="Add Strokes")
        col.operator("gpencil.surfsk_edit_strokes", text="Edit Strokes")
        col.prop(scn, "SURFSK_cyclic_cross")
        col.prop(scn, "SURFSK_cyclic_follow")
        col.prop(scn, "SURFSK_loops_on_strokes")
        col.prop(scn, "SURFSK_automatic_join")
        col.prop(scn, "SURFSK_keep_strokes")
        
class VIEW3D_PT_tools_SURFSK_curve(Panel):
    bl_space_type = 'VIEW_3D'
    bl_region_type = 'UI'
    bl_context = "curve_edit"
    bl_category = 'Tools'
    bl_label = "Bsurfaces"

    @classmethod
    def poll(cls, context):
        return context.active_object

    def draw(self, context):
        layout = self.layout

        col = layout.column(align=True)
        row = layout.row()
        row.separator()
        col.operator("curve.surfsk_first_points", text="Set First Points")
        col.operator("curve.switch_direction", text="Switch Direction")
        col.operator("curve.surfsk_reorder_splines", text="Reorder Splines")


# Returns the type of strokes used
def get_strokes_type():
    strokes_type = ""
    strokes_num = 0

    # Check if they are grease pencil
    try:
       gpencil = bpy.context.scene.bsurfaces.SURFSK_object_with_strokes
       layer = gpencil.data.layers[0]
       frame = layer.frames[0]
        
       strokes_num = len(frame.strokes)

       if strokes_num > 0:
           strokes_type = "GP_STROKES"
    except:
        pass
        
    # Check if they are mesh
    try:
       main_object = bpy.context.scene.bsurfaces.SURFSK_object_with_retopology
    except:
        pass

    # Check if they are curves, if there aren't grease pencil strokes
    if strokes_type == "":
        if len(bpy.context.selected_objects) == 2:
            for ob in bpy.context.selected_objects:
                if ob != bpy.context.view_layer.objects.active and ob.type == "CURVE":
                    strokes_type = "EXTERNAL_CURVE"
                    strokes_num = len(ob.data.splines)

                    # Check if there is any non-bezier spline
                    for i in range(len(ob.data.splines)):
                        if ob.data.splines[i].type != "BEZIER":
                            strokes_type = "CURVE_WITH_NON_BEZIER_SPLINES"
                            break

                elif ob != bpy.context.view_layer.objects.active and ob.type != "CURVE":
                    strokes_type = "EXTERNAL_NO_CURVE"
        elif len(bpy.context.selected_objects) > 2:
            strokes_type = "MORE_THAN_ONE_EXTERNAL"

    # Check if there is a single stroke without any selection in the object
    if strokes_num == 1 and main_object.data.total_vert_sel == 0:
        if strokes_type == "EXTERNAL_CURVE":
            strokes_type = "SINGLE_CURVE_STROKE_NO_SELECTION"
        elif strokes_type == "GP_STROKES":
            strokes_type = "SINGLE_GP_STROKE_NO_SELECTION"

    if strokes_num == 0 and main_object.data.total_vert_sel > 0:
        strokes_type = "SELECTION_ALONE"

    if strokes_type == "":
        strokes_type = "NO_STROKES"

    return strokes_type


# Surface generator operator
class GPENCIL_OT_SURFSK_add_surface(Operator):
    bl_idname = "gpencil.surfsk_add_surface"
    bl_label = "Bsurfaces add surface"
    bl_description = "Generates surfaces from grease pencil strokes, bezier curves or loose edges"
    bl_options = {'REGISTER', 'UNDO'}

    is_fill_faces: BoolProperty(
                    default=True
                    )
    selection_U_exists: BoolProperty(
                    default=False
                    )
    selection_V_exists: BoolProperty(
                    default=False
                    )
    selection_U2_exists: BoolProperty(
                    default=False
                    )
    selection_V2_exists: BoolProperty(
                    default=False
                    )
    selection_V_is_closed: BoolProperty(
                    default=False
                    )
    selection_U_is_closed: BoolProperty(
                    default=False
                    )
    selection_V2_is_closed: BoolProperty(
                    default=False
                    )
    selection_U2_is_closed: BoolProperty(
                    default=False
                    )
    
    edges_U: IntProperty(
                    name="Cross",
                    description="Number of face-loops crossing the strokes",
                    default=1,
                    min=1,
                    max=200
                    )
    edges_V: IntProperty(
                    name="Follow",
                    description="Number of face-loops following the strokes",
                    default=1,
                    min=1,
                    max=200
                    )
    cyclic_cross: BoolProperty(
                    name="Cyclic Cross",
                    description="Make cyclic the face-loops crossing the strokes",
                    default=False
                    )
    cyclic_follow: BoolProperty(
                    name="Cyclic Follow",
                    description="Make cyclic the face-loops following the strokes",
                    default=False
                    )
    loops_on_strokes: BoolProperty(
                    name="Loops on strokes",
                    description="Make the loops match the paths of the strokes",
                    default=False
                    )
    automatic_join: BoolProperty(
                    name="Automatic join",
                    description="Join automatically vertices of either surfaces generated "
                                "by crosshatching, or from the borders of closed shapes",
                    default=False
                    )
    join_stretch_factor: FloatProperty(
                    name="Stretch",
                    description="Amount of stretching or shrinking allowed for "
                                "edges when joining vertices automatically",
                    default=1,
                    min=0,
                    max=3,
                    subtype='FACTOR'
                    )
    keep_strokes: BoolProperty(
                name="Keep strokes",
                description="Keeps the sketched strokes or curves after adding the surface",
                default=False
                )
    strokes_type: StringProperty()
    initial_global_undo_state: BoolProperty()
    

    def draw(self, context):
        layout = self.layout
        col = layout.column(align=True)
        row = layout.row()

        if not self.is_fill_faces:
            row.separator()
            if not self.is_crosshatch:
                if not self.selection_U_exists:
                    col.prop(self, "edges_U")
                    row.separator()

                if not self.selection_V_exists:
                    col.prop(self, "edges_V")
                    row.separator()

                row.separator()

                if not self.selection_U_exists:
                    if not (
                          (self.selection_V_exists and not self.selection_V_is_closed) or
                          (self.selection_V2_exists and not self.selection_V2_is_closed)
                          ):
                        col.prop(self, "cyclic_cross")

                if not self.selection_V_exists:
                    if not (
                          (self.selection_U_exists and not self.selection_U_is_closed) or
                          (self.selection_U2_exists and not self.selection_U2_is_closed)
                          ):
                        col.prop(self, "cyclic_follow")

                col.prop(self, "loops_on_strokes")

            col.prop(self, "automatic_join")

            if self.automatic_join:
                row.separator()
                col.separator()
                row.separator()
                col.prop(self, "join_stretch_factor")
            
            col.prop(self, "keep_strokes")

    # Get an ordered list of a chain of vertices
    def get_ordered_verts(self, ob, all_selected_edges_idx, all_selected_verts_idx,
                          first_vert_idx, middle_vertex_idx, closing_vert_idx):
        # Order selected vertices.
        verts_ordered = []
        if closing_vert_idx is not None:
            verts_ordered.append(ob.data.vertices[closing_vert_idx])

        verts_ordered.append(ob.data.vertices[first_vert_idx])
        prev_v = first_vert_idx
        prev_ed = None
        finish_while = False
        while True:
            edges_non_matched = 0
            for i in all_selected_edges_idx:
                if ob.data.edges[i] != prev_ed and ob.data.edges[i].vertices[0] == prev_v and \
                   ob.data.edges[i].vertices[1] in all_selected_verts_idx:

                    verts_ordered.append(ob.data.vertices[ob.data.edges[i].vertices[1]])
                    prev_v = ob.data.edges[i].vertices[1]
                    prev_ed = ob.data.edges[i]
                elif ob.data.edges[i] != prev_ed and ob.data.edges[i].vertices[1] == prev_v and \
                     ob.data.edges[i].vertices[0] in all_selected_verts_idx:

                    verts_ordered.append(ob.data.vertices[ob.data.edges[i].vertices[0]])
                    prev_v = ob.data.edges[i].vertices[0]
                    prev_ed = ob.data.edges[i]
                else:
                    edges_non_matched += 1

                    if edges_non_matched == len(all_selected_edges_idx):
                        finish_while = True

            if finish_while:
                break

        if closing_vert_idx is not None:
            verts_ordered.append(ob.data.vertices[closing_vert_idx])

        if middle_vertex_idx is not None:
            verts_ordered.append(ob.data.vertices[middle_vertex_idx])
            verts_ordered.reverse()

        return tuple(verts_ordered)

    # Calculates length of a chain of points.
    def get_chain_length(self, object, verts_ordered):
        matrix = object.matrix_world

        edges_lengths = []
        edges_lengths_sum = 0
        for i in range(0, len(verts_ordered)):
            if i == 0:
                prev_v_co = matrix @ verts_ordered[i].co
            else:
                v_co = matrix @ verts_ordered[i].co

                v_difs = [prev_v_co[0] - v_co[0], prev_v_co[1] - v_co[1], prev_v_co[2] - v_co[2]]
                edge_length = abs(sqrt(v_difs[0] * v_difs[0] + v_difs[1] * v_difs[1] + v_difs[2] * v_difs[2]))

                edges_lengths.append(edge_length)
                edges_lengths_sum += edge_length

                prev_v_co = v_co

        return edges_lengths, edges_lengths_sum

    # Calculates the proportion of the edges of a chain of edges, relative to the full chain length.
    def get_edges_proportions(self, edges_lengths, edges_lengths_sum, use_boundaries, fixed_edges_num):
        edges_proportions = []
        if use_boundaries:
            verts_count = 1
            for l in edges_lengths:
                edges_proportions.append(l / edges_lengths_sum)
                verts_count += 1
        else:
            verts_count = 1
            for n in range(0, fixed_edges_num):
                edges_proportions.append(1 / fixed_edges_num)
                verts_count += 1

        return edges_proportions

    # Calculates the angle between two pairs of points in space
    def orientation_difference(self, points_A_co, points_B_co):
        # each parameter should be a list with two elements,
        # and each element should be a x,y,z coordinate
        vec_A = points_A_co[0] - points_A_co[1]
        vec_B = points_B_co[0] - points_B_co[1]

        angle = vec_A.angle(vec_B)

        if angle > 0.5 * pi:
            angle = abs(angle - pi)

        return angle

    # Calculate the which vert of verts_idx list is the nearest one
    # to the point_co coordinates, and the distance
    def shortest_distance(self, object, point_co, verts_idx):
        matrix = object.matrix_world

        for i in range(0, len(verts_idx)):
            dist = (point_co - matrix @ object.data.vertices[verts_idx[i]].co).length
            if i == 0:
                prev_dist = dist
                nearest_vert_idx = verts_idx[i]
                shortest_dist = dist

            if dist < prev_dist:
                prev_dist = dist
                nearest_vert_idx = verts_idx[i]
                shortest_dist = dist

        return nearest_vert_idx, shortest_dist

    # Returns the index of the opposite vert tip in a chain, given a vert tip index
    # as parameter, and a multidimentional list with all pairs of tips
    def opposite_tip(self, vert_tip_idx, all_chains_tips_idx):
        opposite_vert_tip_idx = None
        for i in range(0, len(all_chains_tips_idx)):
            if vert_tip_idx == all_chains_tips_idx[i][0]:
                opposite_vert_tip_idx = all_chains_tips_idx[i][1]
            if vert_tip_idx == all_chains_tips_idx[i][1]:
                opposite_vert_tip_idx = all_chains_tips_idx[i][0]

        return opposite_vert_tip_idx

    # Simplifies a spline and returns the new points coordinates
    def simplify_spline(self, spline_coords, segments_num):
        simplified_spline = []
        points_between_segments = round(len(spline_coords) / segments_num)

        simplified_spline.append(spline_coords[0])
        for i in range(1, segments_num):
            simplified_spline.append(spline_coords[i * points_between_segments])

        simplified_spline.append(spline_coords[len(spline_coords) - 1])

        return simplified_spline

    # Cleans up the scene and gets it the same it was at the beginning,
    # in case the script is interrupted in the middle of the execution
    def cleanup_on_interruption(self):
        # If the original strokes curve comes from conversion
        # from grease pencil and wasn't made by hand, delete it
        if not self.using_external_curves:
            try:
                bpy.ops.object.delete({"selected_objects": [self.original_curve]})
            except:
                pass

            bpy.ops.object.delete({"selected_objects": [self.main_object]})
        else:
            bpy.ops.object.select_all('INVOKE_REGION_WIN', action='DESELECT')
            self.original_curve.select_set(True)
            self.main_object.select_set(True)
            bpy.context.view_layer.objects.active = self.main_object

        bpy.ops.object.editmode_toggle('INVOKE_REGION_WIN')

    # Returns a list with the coords of the points distributed over the splines
    # passed to this method according to the proportions parameter
    def distribute_pts(self, surface_splines, proportions):

        # Calculate the length of each final surface spline
        surface_splines_lengths = []
        surface_splines_parsed = []

        for sp_idx in range(0, len(surface_splines)):
            # Calculate spline length
            surface_splines_lengths.append(0)

            for i in range(0, len(surface_splines[sp_idx].bezier_points)):
                if i == 0:
                    prev_p = surface_splines[sp_idx].bezier_points[i]
                else:
                    p = surface_splines[sp_idx].bezier_points[i]
                    edge_length = (prev_p.co - p.co).length
                    surface_splines_lengths[sp_idx] += edge_length

                    prev_p = p

        # Calculate vertex positions with appropriate edge proportions, and ordered, for each spline
        for sp_idx in range(0, len(surface_splines)):
            surface_splines_parsed.append([])
            surface_splines_parsed[sp_idx].append(surface_splines[sp_idx].bezier_points[0].co)

            prev_p_co = surface_splines[sp_idx].bezier_points[0].co
            p_idx = 0

            for prop_idx in range(len(proportions) - 1):
                target_length = surface_splines_lengths[sp_idx] * proportions[prop_idx]
                partial_segment_length = 0
                finish_while = False

                while True:
                    # if not it'll pass the p_idx as an index below and crash
                    if p_idx < len(surface_splines[sp_idx].bezier_points):
                        p_co = surface_splines[sp_idx].bezier_points[p_idx].co
                        new_dist = (prev_p_co - p_co).length

                    # The new distance that could have the partial segment if
                    # it is still shorter than the target length
                    potential_segment_length = partial_segment_length + new_dist

                    # If the potential is still shorter, keep adding
                    if potential_segment_length < target_length:
                        partial_segment_length = potential_segment_length

                        p_idx += 1
                        prev_p_co = p_co

                    # If the potential is longer than the target, calculate the target
                    # (a point between the last two points), and assign
                    elif potential_segment_length > target_length:
                        remaining_dist = target_length - partial_segment_length
                        vec = p_co - prev_p_co
                        vec.normalize()
                        intermediate_co = prev_p_co + (vec * remaining_dist)

                        surface_splines_parsed[sp_idx].append(intermediate_co)

                        partial_segment_length += remaining_dist
                        prev_p_co = intermediate_co

                        finish_while = True

                    # If the potential is equal to the target, assign
                    elif potential_segment_length == target_length:
                        surface_splines_parsed[sp_idx].append(p_co)
                        prev_p_co = p_co

                        finish_while = True

                    if finish_while:
                        break

            # last point of the spline
            surface_splines_parsed[sp_idx].append(
                    surface_splines[sp_idx].bezier_points[len(surface_splines[sp_idx].bezier_points) - 1].co
                    )

        return surface_splines_parsed

    # Counts the number of faces that belong to each edge
    def edge_face_count(self, ob):
        ed_keys_count_dict = {}

        for face in ob.data.polygons:
            for ed_keys in face.edge_keys:
                if ed_keys not in ed_keys_count_dict:
                    ed_keys_count_dict[ed_keys] = 1
                else:
                    ed_keys_count_dict[ed_keys] += 1

        edge_face_count = []
        for i in range(len(ob.data.edges)):
            edge_face_count.append(0)

        for i in range(len(ob.data.edges)):
            ed = ob.data.edges[i]

            v1 = ed.vertices[0]
            v2 = ed.vertices[1]

            if (v1, v2) in ed_keys_count_dict:
                edge_face_count[i] = ed_keys_count_dict[(v1, v2)]
            elif (v2, v1) in ed_keys_count_dict:
                edge_face_count[i] = ed_keys_count_dict[(v2, v1)]

        return edge_face_count

    # Fills with faces all the selected vertices which form empty triangles or quads
    def fill_with_faces(self, object):
        all_selected_verts_count = self.main_object_selected_verts_count

        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='OBJECT')

        # Calculate average length of selected edges
        all_selected_verts = []
        original_sel_edges_count = 0
        for ed in object.data.edges:
            if object.data.vertices[ed.vertices[0]].select and object.data.vertices[ed.vertices[1]].select:
                coords = []
                coords.append(object.data.vertices[ed.vertices[0]].co)
                coords.append(object.data.vertices[ed.vertices[1]].co)

                original_sel_edges_count += 1

                if not ed.vertices[0] in all_selected_verts:
                    all_selected_verts.append(ed.vertices[0])

                if not ed.vertices[1] in all_selected_verts:
                    all_selected_verts.append(ed.vertices[1])

        tuple(all_selected_verts)

        # Check if there is any edge selected. If not, interrupt the script
        if original_sel_edges_count == 0 and all_selected_verts_count > 0:
            return 0

        # Get all edges connected to selected verts
        all_edges_around_sel_verts = []
        edges_connected_to_sel_verts = {}
        verts_connected_to_every_vert = {}
        for ed_idx in range(len(object.data.edges)):
            ed = object.data.edges[ed_idx]
            include_edge = False

            if ed.vertices[0] in all_selected_verts:
                if not ed.vertices[0] in edges_connected_to_sel_verts:
                    edges_connected_to_sel_verts[ed.vertices[0]] = []

                edges_connected_to_sel_verts[ed.vertices[0]].append(ed_idx)
                include_edge = True

            if ed.vertices[1] in all_selected_verts:
                if not ed.vertices[1] in edges_connected_to_sel_verts:
                    edges_connected_to_sel_verts[ed.vertices[1]] = []

                edges_connected_to_sel_verts[ed.vertices[1]].append(ed_idx)
                include_edge = True

            if include_edge is True:
                all_edges_around_sel_verts.append(ed_idx)

            # Get all connected verts to each vert
            if not ed.vertices[0] in verts_connected_to_every_vert:
                verts_connected_to_every_vert[ed.vertices[0]] = []

            if not ed.vertices[1] in verts_connected_to_every_vert:
                verts_connected_to_every_vert[ed.vertices[1]] = []

            verts_connected_to_every_vert[ed.vertices[0]].append(ed.vertices[1])
            verts_connected_to_every_vert[ed.vertices[1]].append(ed.vertices[0])

        # Get all verts connected to faces
        all_verts_part_of_faces = []
        all_edges_faces_count = []
        all_edges_faces_count += self.edge_face_count(object)

        # Get only the selected edges that have faces attached.
        count_faces_of_edges_around_sel_verts = {}
        selected_verts_with_faces = []
        for ed_idx in all_edges_around_sel_verts:
            count_faces_of_edges_around_sel_verts[ed_idx] = all_edges_faces_count[ed_idx]

            if all_edges_faces_count[ed_idx] > 0:
                ed = object.data.edges[ed_idx]

                if not ed.vertices[0] in selected_verts_with_faces:
                    selected_verts_with_faces.append(ed.vertices[0])

                if not ed.vertices[1] in selected_verts_with_faces:
                    selected_verts_with_faces.append(ed.vertices[1])

                all_verts_part_of_faces.append(ed.vertices[0])
                all_verts_part_of_faces.append(ed.vertices[1])

        tuple(selected_verts_with_faces)

        # Discard unneeded verts from calculations
        participating_verts = []
        movable_verts = []
        for v_idx in all_selected_verts:
            vert_has_edges_with_one_face = False

            # Check if the actual vert has at least one edge connected to only one face
            for ed_idx in edges_connected_to_sel_verts[v_idx]:
                if count_faces_of_edges_around_sel_verts[ed_idx] == 1:
                    vert_has_edges_with_one_face = True

            # If the vert has two or less edges connected and the vert is not part of any face.
            # Or the vert is part of any face and at least one of
            # the connected edges has only one face attached to it.
            if (len(edges_connected_to_sel_verts[v_idx]) == 2 and
               v_idx not in all_verts_part_of_faces) or \
               len(edges_connected_to_sel_verts[v_idx]) == 1 or \
               (v_idx in all_verts_part_of_faces and
               vert_has_edges_with_one_face):

                participating_verts.append(v_idx)

                if v_idx not in all_verts_part_of_faces:
                    movable_verts.append(v_idx)

        # Remove from movable verts list those that are part of closed geometry (ie: triangles, quads)
        for mv_idx in movable_verts:
            freeze_vert = False
            mv_connected_verts = verts_connected_to_every_vert[mv_idx]

            for actual_v_idx in all_selected_verts:
                count_shared_neighbors = 0
                checked_verts = []

                for mv_conn_v_idx in mv_connected_verts:
                    if mv_idx != actual_v_idx:
                        if mv_conn_v_idx in verts_connected_to_every_vert[actual_v_idx] and \
                           mv_conn_v_idx not in checked_verts:
                            count_shared_neighbors += 1
                            checked_verts.append(mv_conn_v_idx)

                            if actual_v_idx in mv_connected_verts:
                                freeze_vert = True
                                break

                        if count_shared_neighbors == 2:
                            freeze_vert = True
                            break

                if freeze_vert:
                    break

            if freeze_vert:
                movable_verts.remove(mv_idx)

        # Calculate merge distance for participating verts
        shortest_edge_length = None
        for ed in object.data.edges:
            if ed.vertices[0] in movable_verts and ed.vertices[1] in movable_verts:
                v1 = object.data.vertices[ed.vertices[0]]
                v2 = object.data.vertices[ed.vertices[1]]

                length = (v1.co - v2.co).length

                if shortest_edge_length is None:
                    shortest_edge_length = length
                else:
                    if length < shortest_edge_length:
                        shortest_edge_length = length

        if shortest_edge_length is not None:
            edges_merge_distance = shortest_edge_length * 0.5
        else:
            edges_merge_distance = 0

        # Get together the verts near enough. They will be merged later
        remaining_verts = []
        remaining_verts += participating_verts
        for v1_idx in participating_verts:
            if v1_idx in remaining_verts and v1_idx in movable_verts:
                verts_to_merge = []
                coords_verts_to_merge = {}

                verts_to_merge.append(v1_idx)

                v1_co = object.data.vertices[v1_idx].co
                coords_verts_to_merge[v1_idx] = (v1_co[0], v1_co[1], v1_co[2])

                for v2_idx in remaining_verts:
                    if v1_idx != v2_idx:
                        v2_co = object.data.vertices[v2_idx].co

                        dist = (v1_co - v2_co).length

                        if dist <= edges_merge_distance:  # Add the verts which are near enough
                            verts_to_merge.append(v2_idx)

                            coords_verts_to_merge[v2_idx] = (v2_co[0], v2_co[1], v2_co[2])

                for vm_idx in verts_to_merge:
                    remaining_verts.remove(vm_idx)

                if len(verts_to_merge) > 1:
                    # Calculate middle point of the verts to merge.
                    sum_x_co = 0
                    sum_y_co = 0
                    sum_z_co = 0
                    movable_verts_to_merge_count = 0
                    for i in range(len(verts_to_merge)):
                        if verts_to_merge[i] in movable_verts:
                            v_co = object.data.vertices[verts_to_merge[i]].co

                            sum_x_co += v_co[0]
                            sum_y_co += v_co[1]
                            sum_z_co += v_co[2]

                            movable_verts_to_merge_count += 1

                    middle_point_co = [
                                sum_x_co / movable_verts_to_merge_count,
                                sum_y_co / movable_verts_to_merge_count,
                                sum_z_co / movable_verts_to_merge_count
                                ]

                    # Check if any vert to be merged is not movable
                    shortest_dist = None
                    are_verts_not_movable = False
                    verts_not_movable = []
                    for v_merge_idx in verts_to_merge:
                        if v_merge_idx in participating_verts and v_merge_idx not in movable_verts:
                            are_verts_not_movable = True
                            verts_not_movable.append(v_merge_idx)

                    if are_verts_not_movable:
                        # Get the vert connected to faces, that is nearest to
                        # the middle point of the movable verts
                        shortest_dist = None
                        for vcf_idx in verts_not_movable:
                                dist = abs((object.data.vertices[vcf_idx].co -
                                            Vector(middle_point_co)).length)

                                if shortest_dist is None:
                                    shortest_dist = dist
                                    nearest_vert_idx = vcf_idx
                                else:
                                    if dist < shortest_dist:
                                        shortest_dist = dist
                                        nearest_vert_idx = vcf_idx

                        coords = object.data.vertices[nearest_vert_idx].co
                        target_point_co = [coords[0], coords[1], coords[2]]
                    else:
                        target_point_co = middle_point_co

                    # Move verts to merge to the middle position
                    for v_merge_idx in verts_to_merge:
                        if v_merge_idx in movable_verts:  # Only move the verts that are not part of faces
                            object.data.vertices[v_merge_idx].co[0] = target_point_co[0]
                            object.data.vertices[v_merge_idx].co[1] = target_point_co[1]
                            object.data.vertices[v_merge_idx].co[2] = target_point_co[2]

        # Perform "Remove Doubles" to weld all the disconnected verts
        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='EDIT')
        bpy.ops.mesh.remove_doubles(threshold=0.0001)

        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='OBJECT')

        # Get all the definitive selected edges, after weldding
        selected_edges = []
        edges_per_vert = {}  # Number of faces of each selected edge
        for ed in object.data.edges:
            if object.data.vertices[ed.vertices[0]].select and object.data.vertices[ed.vertices[1]].select:
                selected_edges.append(ed.index)

                # Save all the edges that belong to each vertex.
                if not ed.vertices[0] in edges_per_vert:
                    edges_per_vert[ed.vertices[0]] = []

                if not ed.vertices[1] in edges_per_vert:
                    edges_per_vert[ed.vertices[1]] = []

                edges_per_vert[ed.vertices[0]].append(ed.index)
                edges_per_vert[ed.vertices[1]].append(ed.index)

        # Check if all the edges connected to each vert have two faces attached to them.
        # To discard them later and make calculations faster
        a = []
        a += self.edge_face_count(object)
        tuple(a)
        verts_surrounded_by_faces = {}
        for v_idx in edges_per_vert:
            edges = edges_per_vert[v_idx]
            edges_with_two_faces_count = 0

            for ed_idx in edges_per_vert[v_idx]:
                if a[ed_idx] == 2:
                    edges_with_two_faces_count += 1

            if edges_with_two_faces_count == len(edges_per_vert[v_idx]):
                verts_surrounded_by_faces[v_idx] = True
            else:
                verts_surrounded_by_faces[v_idx] = False

        # Get all the selected vertices
        selected_verts_idx = []
        for v in object.data.vertices:
            if v.select:
                selected_verts_idx.append(v.index)

        # Get all the faces of the object
        all_object_faces_verts_idx = []
        for face in object.data.polygons:
            face_verts = []
            face_verts.append(face.vertices[0])
            face_verts.append(face.vertices[1])
            face_verts.append(face.vertices[2])

            if len(face.vertices) == 4:
                face_verts.append(face.vertices[3])

            all_object_faces_verts_idx.append(face_verts)

        # Deselect all vertices
        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='EDIT')
        bpy.ops.mesh.select_all('INVOKE_REGION_WIN', action='DESELECT')
        bpy.ops.object.mode_set('INVOKE_REGION_WIN', mode='OBJECT')

        # Make a dictionary with the verts related to each vert
        related_key_verts = {}
        for ed_idx in selected_edges:
            ed = object.data.edges[ed_idx]

            if not verts_surrounded_by_faces[ed.vertices[0]]:
                if not ed.vertices[0] in related_key_verts:
                    related_key_verts[ed.vertices[0]] = []

                if not ed.vertices[1] in related_key_verts[ed.vertices[0]]:
                    related_key_verts[ed.vertices[0]].append(ed.vertices[1])

            if not verts_surrounded_by_faces[ed.vertices[1]]:
                if not ed.vertices[1] in related_key_verts:
                    related_key_verts[ed.vertices[1]] = []

                if not ed.vertices[0] in related_key_verts[ed.vertices[1]]:
                    related_key_verts[ed.vertices[1]].append(ed.vertices[0])

        # Get groups of verts forming each face
        faces_verts_idx = []
        for v1 in related_key_verts:      # verts-1 ....
            for v2 in related_key_verts:  # verts-2
                if v1 != v2:
                    related_verts_in_common = []
                    v2_in_rel_v1 = False
                    v1_in_rel_v2 = False
                    for rel_v1 in related_key_verts[v1]:
                        # Check if related verts of verts-1 are related verts of verts-2
                        if rel_v1 in related_key_verts[v2]:
                            related_verts_in_common.append(rel_v1)

                    if v2 in related_key_verts[v1]:
                        v2_in_rel_v1 = True

                    if v1 in related_key_verts[v2]:
                        v1_in_rel_v2 = True

                    repeated_face = False
                    # If two verts have two related verts in common, they form a quad
                    if len(related_verts_in_common) == 2:
                        # Check if the face is already saved
                        all_faces_to_check_idx = faces_verts_idx + all_object_faces_verts_idx

                        for f_verts in all_faces_to_check_idx:
                            repeated_verts = 0

                            if len(f_verts) == 4:
                                if v1 in f_verts:
                                    repeated_verts += 1
                                if v2 in f_verts:
                                    repeated_verts += 1
                                if related_verts_in_common[0] in f_verts:
                                    repeated_verts += 1
                                if related_verts_in_common[1] in f_verts:
                                    repeated_verts += 1

                                if repeated_verts == len(f_verts):
                                    repeated_face = True
                                    break

                        if not repeated_face:
                            faces_verts_idx.append(
                                    [v1, related_verts_in_common[0], v2, related_verts_in_common[1]]
                                    )

                    # If Two verts have one related vert in common and
                    # they are related to each other, they form a triangle
                    elif v2_in_rel_v1 and v1_in_rel_v2 and len(related_verts_in_common) == 1:
                        # Check if the face is already saved.
                        all_faces_to_check_idx = faces_verts_idx + all_object_faces_verts_idx

                        for f_verts in all_faces_to_check_idx:
                            repeated_verts = 0

                            if len(f_verts) == 3:
                                if v1 in f_verts:
                                    repeated_verts += 1
                                if v2 in f_verts:
                                    repeated_verts += 1
                                if related_verts_in_common[0] in f_verts:
                                    repeated_verts += 1

                                if repeated_verts == len(f_verts):
                                    repeated_face = True
                                    break

                        if not repeated_face:
                            faces_verts_idx.append([v1, related_verts_in_common[0], v2])

        # Keep only the faces that don't overlap by ignoring quads
        # that overlap with two adjacent triangles
        faces_to_not_include_idx = []  # Indices of faces_verts_idx to eliminate
        all_faces_to_check_idx = faces_verts_idx + all_object_faces_verts_idx
        for i in range(len(faces_verts_idx)):
            for t in range(len(all_faces_to_check_idx)):
                if i != t:
                    verts_in_common = 0

                    if len(faces_verts_idx[i]) == 4 and len(all_faces_to_check_idx[t]) == 3:
                        for v_idx in all_faces_to_check_idx[t]:
                            if v_idx in faces_verts_idx[i]:
                                verts_in_common += 1
                        # If it doesn't have all it's vertices repeated in the other face
                        if verts_in_common == 3:
                            if i not in faces_to_not_include_idx:
                                faces_to_not_include_idx.append(i)

        # Build faces discarding the ones in faces_to_not_include